1 Introduction
The origin of magnetic activity in stellar interiors is a fundamental problem of magnetohydrodynamics. The global solar magnetic field oscillates with a mean period of 22 years (leading to an 11-year activity cycle) and is believed to be generated via a dynamo acting (at least in part) deep within the Sun. The Sun’s magnetic field is largely dipolar; i.e. the mean azimuthal field that leads to the formation of active regions is generally antisymmetric about the equator. However, when this field is weak at the end of a cycle, it takes on a more mixed character, with a quadrupole component that becomes significant (Sokoloff & Nesme-Ribes Reference Sokoloff and Nesme-Ribes1994). Furthermore, direct observations and proxy data demonstrate that the amplitude of the solar cycle is modulated on longer time scales. There is indeed a period of reduced activity between 1645 and 1715 – the Maunder minimum – when the occurrence of sunspots was much reduced (Eddy Reference Eddy1976; Usoskin et al.
Reference Usoskin, Arlt, Asvestari, Hawkins, Käpylä, Kovaltsov, Krivova, Lockwood, Mursula and O’Reilly2015). Analysis of the abundances of the cosmogenic isotopes
$^{10}\text{Be}$
in polar ice and
$^{14}\text{C}$
in tree rings reveals 27 grand minima in the past 11 000 years, separated by aperiodic intervals of approximately 200 years (McCracken et al.
Reference McCracken, Beer, Steinhilber and Abreu2013; Usoskin Reference Usoskin2013). A key observation for our understanding of the processes leading to modulation is that as the Sun emerged from the Maunder minimum, sunspots were largely restricted to the southern hemisphere, showing that the magnetic field emerged with a mixed character with both dipole and quadrupole components (Sokoloff & Nesme-Ribes Reference Sokoloff and Nesme-Ribes1994). Moreover, between 1750 and 1775, the solar magnetic field took on a more quadrupolar character, with sunspots appearing at the equator (Arlt Reference Arlt2009). There is now evidence from the cosmogenic isotope records that the Sun switches on a long time scale between strong modulation with clusters of deep grand minima and weaker modulation, which can be associated with symmetry breaking. This ‘supermodulation’ is an example of chaotic (though deterministic) modulational effects (Weiss & Tobias Reference Weiss and Tobias2016). Evidence for modulation in other stars arises from the long-term monitoring of the CaII H
$+$
K flux of solar-type stars started by Wilson in 1968. The so-called Mount Wilson Observatory survey provides a panel of different stellar activities, in which 60 % of stars exhibit periodic variations, and 25 % show irregular or aperiodic variability (Baliunas et al.
Reference Baliunas, Donahue, Soon and Henry1998; Oláh et al.
Reference Oláh, Kolláth, Granzer, Strassmeier, Lanza, Järvinen, Korhonen, Baliunas, Soon and Messina2009). Evidence for changes of symmetry in young rapidly rotating stars is also now beginning to emerge (Hackman et al.
Reference Hackman, Lehtinen, Rosén, Kochukhov and Käpylä2016).
Stellar magnetic fields are thought to be maintained against ohmic dissipation by dynamo action through the flow of an electrically conducting fluid (Moffatt Reference Moffatt1978). Although it is known that systematic activity can be generated through the interaction of turbulent flows with rotation, shear and magnetic fields, no satisfactory self-consistent nonlinear model of dynamo action is currently available (Jones, Thompson & Tobias Reference Jones, Thompson and Tobias2010; Charbonneau Reference Charbonneau2014). Direct numerical simulations aimed at understanding these interactions are restricted to parameters well away from those pertaining to stellar interiors (with Reynolds numbers and magnetic Reynolds numbers (
$Rm$
) orders of magnitudes smaller than would be realistic). For this reason, much attention has been focused on mean-field models of dynamo action (Krause & Rädler Reference Krause and Rädler1980). In this paradigm, only the large-scale flows and magnetic fields are modelled, with small-scale interactions being parameterized via transport coefficients such as the
${\it\alpha}$
-tensor and the turbulent diffusivity. Although there are many issues with the mean-field formalism – primary among these is whether mean fields can ever be seen at high
$Rm$
or whether the solution is dominated by the fluctuations – these models are of use in describing the dynamics of mean fields once they have been generated. In particular, the mean-field equations naturally respect the symmetries of the underlying rotating spherical system (Knobloch Reference Knobloch, Proctor and Gilbert1994), and capture the nonlinear interactions between magnetic modes of different symmetries and the underlying large-scale velocity field that is driving the dynamo.
Mean-field dynamo models have demonstrated that modulation of the basic cycle may occur through stochastic fluctuations in the underlying transport coefficients (Schmitt, Schuessler & Ferriz-Mas Reference Schmitt, Schuessler and Ferriz-Mas1996; Choudhuri & Karak Reference Choudhuri and Karak2012; Hazra, Passos & Nandy Reference Hazra, Passos and Nandy2014) or more naturally via nonlinear interactions inherent in the dynamo equations leading to chaotic (though deterministic) modulation (Pipin Reference Pipin1999; Bushby & Mason Reference Bushby and Mason2004). The type of modulation can be classified according to the key nonlinear interactions that are primarily responsible (Tobias Reference Tobias2002). In the first (type 1 modulation), magnetic modes of different symmetry (e.g. dipole and quadrupole modes) interact to produce modulation of the basic cycle, with significant changes in the symmetry (parity) of solutions. This behaviour is similar to that seen in the sunspot record over the past 300 years. In the second (imaginatively termed type 2 modulation), a magnetic mode with a given symmetry undergoes modulation via interaction with a large-scale velocity field; here, changes in the amplitude of the basic cycle occur with no significant changes in the symmetry of solutions. Recently, Weiss & Tobias (Reference Weiss and Tobias2016) have argued from analysis of cosmogenic isotope records that both of these modulational mechanisms have been at play in the solar dynamo, leading to ‘supermodulation’ on long time scales. The precise modulational effects important in the system are sometimes model-dependent. For this reason, progress can also be made by considering low-order systems based on symmetry considerations (Knobloch & Landsberg Reference Knobloch and Landsberg1996; Knobloch, Tobias & Weiss Reference Knobloch, Tobias and Weiss1998; Weiss Reference Weiss2011). These models demonstrate that the dynamics found in the ad hoc mean-field models is robust and may be expected in simulations of the full three-dimensional dynamo system. Symmetry arguments have also proved useful in explaining the dynamics of dynamo experiments and the geodynamo where the first bifurcation is stationary (Pétrélis et al. Reference Pétrélis, Fauve, Dormy and Valet2009).
In this paper, we present the results of three-dimensional numerical solutions of dynamos driven by anelastic convection in a spherical shell. We do not attempt to model solar convection directly, as this is well beyond the scope of modern-day computations. Rather, we focus on the symmetries and nonlinear interactions that lead to modulation in dynamos, and provide examples of the basic types of modulation and supermodulation. These results are important for our understanding of magnetic field generation via dynamo action, not only in late-type stars, but also in other astrophysical objects such as planets.
2 Governing equations
We consider electrically conducting fluid in a spherical shell rotating at angular velocity
${\it\Omega}\,\boldsymbol{e}_{z}$
. The shell is bounded by two concentric spheres of radius
$r_{i}$
and
$r_{o}$
and we define the shell width as
$d=r_{o}-r_{i}$
and aspect ratio as
${\it\chi}=r_{i}/r_{o}$
. We rely on the LBR anelastic approximation (Braginsky & Roberts Reference Braginsky and Roberts1995; Lantz & Fan Reference Lantz and Fan1999) to model a perfect gas with kinematic viscosity
${\it\nu}$
, turbulent entropy diffusivity
${\it\kappa}$
, specific heat
$c_{p}$
and magnetic diffusivity
${\it\eta}$
(all assumed to be constant). The gravity is given by
$\boldsymbol{g}=-GM\boldsymbol{e}_{r}/r^{2}$
, where
$G$
is the gravitational constant and
$M$
is the central mass. The equilibrium polytropic solution of the anelastic system defines the reference-state pressure
$\overline{P}=P_{c}{\it\zeta}^{n+1}$
, density
$\overline{{\it\varrho}}={\it\varrho}_{c}{\it\zeta}^{n}$
and temperature
$\overline{T}=T_{c}{\it\zeta}$
, with
${\it\zeta}=c_{0}+c_{1}d/r$
,
$c_{0}=(2{\it\zeta}_{0}-{\it\chi}-1)/(1-{\it\chi})$
,
$c_{1}=(1+{\it\chi})(1-{\it\zeta}_{o})/(1-{\it\chi})^{2}$
and
${\it\zeta}_{0}=({\it\chi}+1)/({\it\chi}\exp (N_{{\it\varrho}}/n)+1)$
. The constants
$P_{c}$
,
${\it\varrho}_{c}$
and
$T_{c}$
are the reference-state pressure, density and temperature mid-way between the inner and outer boundaries. These reference values serve as units for these variables, while length is scaled by
$d$
, time by
$d^{2}/{\it\eta}$
, entropy by
${\rm\Delta}s$
(the entropy drop across the layer) and magnetic field by
$\sqrt{{\it\Omega}{\it\varrho}_{c}{\it\mu}{\it\eta}}$
, where
${\it\mu}$
is the magnetic permeability. Then, the governing equations are (Jones et al.
Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011)
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}\boldsymbol{v}}{\text{D}t}=Pm\left[-\frac{1}{E}\boldsymbol{{\rm\nabla}}\frac{P^{\prime }}{{\it\zeta}^{n}}+\frac{Pm}{Pr}Ra\frac{s}{r^{2}}\boldsymbol{e}_{r}-\frac{2}{E}\,\boldsymbol{e}_{z}\times \boldsymbol{v}+\boldsymbol{F}_{{\it\nu}}+\frac{1}{E\,{\it\zeta}^{n}}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{B}}{\partial t}=\boldsymbol{{\rm\nabla}}\times (\boldsymbol{v}\times \boldsymbol{B})+{\rm\nabla}^{2}\boldsymbol{B}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{D}s}{\text{D}t}={\it\zeta}^{-n-1}\frac{Pm}{Pr}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\zeta}^{n+1}\,\boldsymbol{{\rm\nabla}}s)+\frac{Di}{{\it\zeta}}[E^{-1}{\it\zeta}^{-n}(\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})^{2}+Q_{{\it\nu}}], & \displaystyle\end{eqnarray}$$
with the constraints
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\zeta}^{n}\boldsymbol{v})=0$
and
$\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$
. In the Navier–Stokes equation (2.1),
$P^{\prime }$
denotes the pressure perturbation, and the viscous force
$\boldsymbol{F}_{{\it\nu}}$
is given by
$\boldsymbol{F}_{{\it\nu}}={\it\zeta}^{-n}\boldsymbol{{\rm\nabla}}\unicode[STIX]{x1D64E}$
, with
$S_{ij}=2{\it\zeta}^{n}(e_{ij}-(1/3){\it\delta}_{ij}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})$
and
$2e_{ij}=\partial _{j}v_{i}+\partial _{i}v_{j}$
. The expressions for the dissipation parameter
$Di$
and the viscous heating
$Q_{{\it\nu}}$
in (2.3) are
$Di=c_{1}Pr/(PmRa)$
and
$Q_{{\it\nu}}=2[e_{ij}e_{ij}-(1/3)(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})^{2}]$
. Following Jones et al. (Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011), we impose stress-free boundary conditions for the velocity field, and the magnetic field matches a potential field inside and outside the fluid shell. The convection is driven by an imposed entropy difference
${\rm\Delta}s$
between the inner and outer boundaries. The above system involves seven control parameters: the Rayleigh number
$Ra=GMd{\rm\Delta}s/({\it\nu}{\it\kappa}c_{p})$
, the Ekman number
$E={\it\nu}/({\it\Omega}d^{2})$
, the Prandtl number
$Pr={\it\nu}/{\it\kappa}$
, the magnetic Prandtl number
$Pm={\it\nu}/{\it\eta}$
, together with the aspect ratio
${\it\chi}$
, the polytropic index
$n$
and the number of density scale heights
$N_{{\it\varrho}}\equiv \ln [\overline{{\it\varrho}}(r_{i})/\overline{{\it\varrho}}(r_{o})]$
. We set
$E=10^{-4}$
,
$Pr=1$
,
$Pm=1$
,
${\it\chi}=0.35$
,
$n=2$
and choose a relatively weak density stratification
$N_{{\it\varrho}}=0.5$
, to limit the computational time. The critical Rayleigh number for the linear onset of convection is then
$Ra_{c}=3.34\times 10^{5}$
(after Schrinner et al.
Reference Schrinner, Petitdemange, Raynaud and Dormy2014).
The anelastic equations are integrated for between five and 60 magnetic diffusion times, which is certainly long enough to establish dynamo action, utilizing the pseudo-spectral code parody (Dormy, Cardin & Jault Reference Dormy, Cardin and Jault1998), whose anelastic version (Schrinner et al.
Reference Schrinner, Petitdemange, Raynaud and Dormy2014) reproduces the anelastic dynamo benchmark proposed by Jones et al. (Reference Jones, Boronski, Brun, Glatzmaier, Gastine, Miesch and Wicht2011). Typical resolutions use 288 points in the radial direction and a spherical harmonic decomposition truncated at degree
$l_{max}\sim 80$
and order
$m_{max}\sim 60$
. As an empirical validation of convergence, we ensure for both spectra a decrease of more than three orders of magnitude over the range of
$l$
and
$m$
. We define the kinetic energy
$E_{k}=(1/2)\int {\it\zeta}^{n}\boldsymbol{v}^{2}\,\text{d}V$
and the magnetic energy
$E_{b}=Pm/(2E)\int \boldsymbol{B}^{2}\,\text{d}V$
. With our choice of units, a non-dimensional measure of the velocity amplitude is naturally given by the magnetic Reynolds number
$Rm=\sqrt{2E_{k}/V}$
,
$V$
being the volume of the fluid shell. Crucially for this investigation, which is concerned with the symmetries of the solutions about the equatorial plane, we also decompose both the kinetic and magnetic energies according to their symmetry about the equator (
$E_{k}^{S}$
,
$E_{k}^{A}$
,
$E_{b}^{S}$
and
$E_{b}^{A}$
respectively). For clarity, we prefer to avoid the terms dipole and quadrupole families which are also in use to denote the different parities of the magnetic field; we further adopt the same definition as Knobloch et al. (Reference Knobloch, Tobias and Weiss1998), according to which symmetric refers to an overall field with dipole symmetry (and vice versa). This choice is consistent with the properties of pseudo-vectors: the dipole family is invariant under reflection with respect to the equatorial plane, but the quadrupole family is not.
Figure 1. (a) Three-dimensional phase portrait showing the projection of the system trajectory onto the space
$(E_{b}^{S},E_{b}^{A},E_{Z})$
. Solid circles denote antisymmetric solutions while crosses denote symmetric solutions. The limit cycle at
$Ra=1.40\times 10^{6}$
(solid blue line) has been obtained from a similar mixed-mode solution after decreasing the value of
$Ra$
. (b) Energy time series for the limit cycle at
$Ra=1.47\times 10^{6}$
.
3 Results
In this paper, we aim to study the symmetry interactions and low-frequency modulations of the dynamo waves that are characteristics of the so-called multipolar dynamo branch (Gastine, Duarte & Wicht Reference Gastine, Duarte and Wicht2012; Schrinner et al.
Reference Schrinner, Petitdemange, Raynaud and Dormy2014). This branch is the only one that can be sustained at low magnetic Reynolds number
$Rm\sim 40$
(Raynaud, Petitdemange & Dormy Reference Raynaud, Petitdemange and Dormy2015). We stress at the outset that the dynamo magnetic fields we consider here, independent of their symmetry about the equator, are dominated by their
$m=1$
component, and note that this differs from the Sun – although the Sun does show a tendency for active longitudes. By considering almost Boussinesq models with
$N_{{\it\varrho}}=0.1$
, Raynaud, Petitdemange & Dormy (Reference Raynaud, Petitdemange and Dormy2014) showed that the non-axisymmetry is related to the choice of a gravity profile corresponding to a central mass distribution. It should be noted that at the low values of
$Rm$
considered here the advective time is comparable with the ohmic diffusive time (in contrast to stars). These solutions are usually interpreted in terms of Parker (Reference Parker1955) waves, in both the Boussinesq (Busse & Simitev Reference Busse and Simitev2006; Schrinner, Petitdemange & Dormy Reference Schrinner, Petitdemange and Dormy2011; Dietrich, Schmitt & Wicht Reference Dietrich, Schmitt and Wicht2013) and anelastic frameworks (Gastine et al.
Reference Gastine, Duarte and Wicht2012), although this interpretation relies on crude estimates of the
${\it\alpha}$
-effect via the flow helicity. Following the methodology of Schrinner, Petitdemange & Dormy (Reference Schrinner, Petitdemange and Dormy2012), we confirm the key role played by differential rotation in the generation of the toroidal magnetic field in our sample of models. It is well known that the
${\it\alpha}{\it\Omega}$
dynamo instability generically sets in as a Hopf bifurcation leading to oscillatory solutions. Our aim here is to identify the changes in the symmetry of the solutions as
$Ra$
is increased with other parameters held fixed. In practice, we easily distinguish different branches of solution by restarting from the closest simulations performed with other parameters; in a few cases, we also tested their stability by restarting the simulation after killing one or the other parity of the magnetic field. At
$Ra=1.39\times 10^{6}$
, the flow does not break the equatorial symmetry and magnetic modes of different parity are linearly decoupled. Depending on the choice of the initial conditions, we effectively observe a bistability between symmetric and antisymmetric solutions, illustrated by the red cross and the red dot in figure 1(a). In this figure, the trajectory of the system is projected for different Rayleigh numbers onto the space spanned by the symmetric and antisymmetric magnetic energies
$E_{b}^{S}$
and
$E_{b}^{A}$
, and the zonal wind energy measured by the axisymmetric toroidal kinetic energy
$E_{Z}$
– a projection introduced by Knobloch et al. (Reference Knobloch, Tobias and Weiss1998). In spite of a misleading effect of perspective, it can be noted that the contribution of the antisymmetric magnetic field does reduce to a negligible fraction for the symmetric solutions (crosses), for which
$E_{b}^{A}/E_{b}^{S}\leqslant 10^{-2}$
. We further stress that the aforementioned bistability must not be confused with the hysteretic transition between the dipolar and multipolar branches resulting from the use of stress-free boundary conditions (Schrinner et al.
Reference Schrinner, Petitdemange and Dormy2012). When the magnetic field is predominantly antisymmetric, the flow is characterized by an
$m=8$
convection mode; on the other hand, when the magnetic field is predominantly symmetric, the flow is then characterized by an
$m=9$
convection mode and larger fluctuations of the kinetic energy. We also note – although it is difficult to see from figure 1(a) – that at these parameters the symmetric mode is quasiperiodic, having undergone a bifurcation from the periodic state, while the antisymmetric mode is strictly periodic (taking the form of a dynamo wave).
Figure 2. Snapshots of
$B_{r}$
at the surface of the model for
$Ra=1.47\times 10^{6}$
taken at the times highlighted by the dashed (a) and solid (b) vertical lines in figure 1(b). It should be noted that when
$E_{b}^{S}>E_{b}^{A}$
, the magnetic field displays an overall dipole symmetry (and conversely).
Increase of the Rayleigh number from
$1.40\times 10^{6}$
to
$1.45\times 10^{6}$
leads to the destabilization of the antisymmetric solution (blue dot in figure 1
a) and the discovery of an asymmetric solution that takes the form of a limit cycle in this phase space (green solid line); the basic dynamo wave is modulated by change in the underlying parity of the solution. This solution coexists with the symmetric solution (green cross). An example of the limit cycle at
$Ra=1.47\times 10^{6}$
is given in figure 1(b) which shows the time series of the antisymmetric kinetic energy
$E_{k}^{A}$
(green solid line), together with those for the symmetric and antisymmetric magnetic energies
$E_{b}^{S}$
and
$E_{b}^{A}$
(represented by the solid red and dashed black lines respectively). This solution is characterized by weak symmetry breaking of the flow coupling magnetic modes of different parity. Indeed, we clearly see a periodic exchange of energy between modes of opposite parity, which could be described as a type 1 modulation, in reference to the terminology introduced by Knobloch et al. (Reference Knobloch, Tobias and Weiss1998). Figure 2 shows projections of the radial magnetic field at times when the solution is mixed and antisymmetric. In figure 2(a), we note that when the solution is a mixed mode the magnetic field tends to be localized in one hemisphere (Grote & Busse Reference Grote and Busse2000; Gastine et al.
Reference Gastine, Duarte and Wicht2012). We believe that this mixed-mode solution is born in a subcritical secondary Hopf bifurcation from the antisymmetric state. Evidence for this arises from the hysteresis that can be identified. As we can see in figure 1(a), this state indeed coexists with both the symmetric and antisymmetric states down to
$Ra=1.40\times 10^{6}$
, below which it disappears (presumably in a saddle-node bifurcation). We stress that what sets the dependence of the period of both the basic cycle and the modulation of the dynamos for strongly nonlinear solutions is an open problem and one that is important for understanding stellar activity (Tobias Reference Tobias1998; Dubé & Charbonneau Reference Dubé and Charbonneau2013). In our sample of models, the modulation period
$T_{mod}$
is sensitive to the value of the Rayleigh number but tends towards a constant when it approaches the critical bifurcation value: for
$Ra\in [1.40\times 10^{6},1.43\times 10^{6}]$
, we have
$T_{mod}\simeq 3.0\pm 0.2$
; in contrast, for
$Ra\in [1.45\times 10^{6},1.55\times 10^{6}]$
, it seems that
$T_{mod}\propto 1/\sqrt{Ra}$
, although accurate measurements are compromised by the fact that we just have six data points, four of which are only metastable for
$Ra\geqslant 1.49\times 10^{6}$
. Figure 1(a) indeed shows that this mixed-mode limit cycle eventually loses its stability when the Rayleigh number is increased to
$1.49\times 10^{6}$
(dashed black line), and the solution develops more sign of spatio-temporal complexity, as described below.
When the Rayleigh number is further increased, we find that the dynamics of the magnetic field progressively switches from parity to amplitude modulations, i.e. from type 1 to type 2 modulation. This transition is particularly clear when comparing the three-dimensional phase portraits represented for increasing values of the Rayleigh number in figure 3, in which the trajectory of the system has been smoothed by applying a moving average that removes the basic dynamo cycle and short-period oscillations. For
$Ra=1.55\times 10^{6}$
(see figure 3
a), the dynamics is mainly governed by the energy exchange between
$E_{b}^{S}$
and
$E_{b}^{A}$
(i.e. type 1 modulation), and we only distinguish the first signs of the type 2 modulation through intermittent decays of the magnetic energy, always followed by an increase of the zonal wind. In stark contrast, we see in figure 3(c) that the system trajectory in the space
$(E_{b}^{S},E_{b}^{A},E_{Z})$
is actually confined near the antisymmetric subspace (i.e.
$E_{b}^{S}\ll E_{b}^{A}$
) and characterized by the strong amplitude modulation of the antisymmetric energy by the zonal wind for
$Ra=1.85\times 10^{6}$
. This is clear type 2 modulation. Most interesting, however, is the attractor for
$Ra=1.65\times 10^{6}$
in figure 3(b). This clearly shows the solution exhibiting both types of modulation; type 1 modulation where there are no minima in activity but energy transfer between the modes of different symmetries and type 2 modulation where the antisymmetric solution regularly visits grand minima in activity through interactions with the zonal wind. The transition between these two types of modulation has been termed supermodulation and is believed to be prominent in solar activity records (Weiss & Tobias Reference Weiss and Tobias2016).
Figure 3. Three-dimensional phase portraits showing the projection of the system trajectory onto the space
$(E_{b}^{S},E_{b}^{A},E_{Z})$
.
Figure 4. Time series corresponding to the phase portraits in figure 3. (a–c) Time series of zonal wind energy (blue) and total magnetic energy (red). (d–f) Time series of symmetric magnetic energy (red) and antisymmetric magnetic energy (black). The shaded regions in (b) and (e) highlight the occurrences of type 1 modulation.
The time series corresponding to the trajectories in figure 3 are shown in figure 4. For
$Ra=1.55\times 10^{6}$
, figures 4(a) and 4(d) reveal that the magnetic energy does not exhibit any deep minima although four dips in total energy are visible, and that there are periods when the symmetric energy is greater than the antisymmetric energy and periods when they are comparable. In contrast, figures 4(c) and 4(f) show the strong amplitude modulation of both the zonal wind and the magnetic energy, leading to grand minima observed at
$Ra=1.85\times 10^{6}$
. This temporal evolution is reminiscent of the relaxation oscillations that can be observed in turbulent hydrodynamic convection (Grote, Busse & Tilgner Reference Grote, Busse and Tilgner2000; Christensen Reference Christensen2002). In that model, the relaxation phenomenon originates from the fact that the columnar convection feeds the differential rotation through the action of Reynolds stresses but also tends to be disrupted by the shear due to differential rotation. This competition between the convection and the zonal wind is present in our models, since we note, for instance, that the Nusselt number is always minimum when the zonal wind reaches its maximum. However, a hydrodynamic simulation performed at
$Ra=1.75\times 10^{6}$
– where the system tends to switch from supermodulation to pure type 2 modulation – demonstrates that the flow does not continue to break the equatorial symmetry when the magnetic field is turned off, and also that there is no amplitude modulation without the backreaction of the Lorentz force; therefore, as expected, the magnetic field is playing a key role here. The sudden growth of the zonal wind results thus from the decrease of the magnetic field. In general, we observe that the Nusselt number is higher when the magnetic field is present, which confirms that the magnetic field promotes the columnar convection and thus the heat transport by reduction of the zonal wind (Grote et al.
Reference Grote, Busse and Tilgner2000; Yadav et al.
Reference Yadav, Gastine, Christensen, Duarte and Reiners2016). At
$Ra=1.85\times 10^{6}$
, the antisymmetric magnetic energy is always larger than that for the symmetric field. A closer examination suggests that minima are caused by interactions with the zonal wind, and we report that the typical time scale between two minima is affected by only 10 % variations of the Prandtl numbers. It tends to increase at lower
$Pm$
(or higher
$Pr$
), and the modulation even disappears if one of these parameter values is lowered from 1 to of the order of 0.7 (not shown). More generally, these observations are compatible with the mean-field dynamo results of Tobias (Reference Tobias1996), who explained the dependence on the magnetic Prandtl number in terms of the time taken for the zonal velocity to decay once its energy source has been diminished by the interaction of magnetic fields with convection.
The supermodulation is shown in figures 4(b) and 4(e), which highlight how the nonlinear solution naturally transitions between the different types of modulational processes. For example, in the shaded regions, the solution undergoes changes in symmetry with no deep minima (type 1 modulation), while between
$t\approx 6$
and
$t\approx 9$
clusters of grand minima are found. From a mathematical perspective, it is no surprise that a chaotic nonlinear dynamo solution exhibits such behaviour, which has indeed been predicted (Weiss & Tobias Reference Weiss and Tobias2016).
Figure 5. Butterfly diagrams representing the axisymmetric azimuthal component of the magnetic field just below the outer surface (
$r=0.998\,r_{o}$
) at
$Ra=1.65\times 10^{6}$
. The time intervals in (b) and (c) correspond to the period between the vertical lines in (a) – dotted for (b) and solid for (c).
Finally, if we examine more closely the evolution of the axisymmetric magnetic field as a function of colatitude and time, we see in figure 5 that both type 1 and type 2 modulations affect the so-called butterfly diagrams in the form of interesting patterns. These butterfly diagrams correspond to the model with
$Ra=1.65\times 10^{6}$
, whose phase portrait is shown in figure 3(b), and for which supermodulation is present. In addition to underlining the oscillatory nature of these dynamos, they demonstrate that both modulational processes occur on time scales that are not comparable to the period of the dynamo wave, which is in general of the order of 0.1 magnetic diffusion times in our sample of models. Figure 5(a) illustrates type 2 modulation for
$t\in [6,8]$
, and figure 5(b) emphasizes the change in amplitude of the magnetic field when the system emerges from a grand minimum. In contrast, the characteristic features of type 1 modulation are shown in figure 5(c), with the hemispherical magnetic field undergoing a change of parity at constant amplitude. Furthermore, the comparison between figures 5(b) and 5(c) (which both cover the same time span) indicates that the period of the basic cycle is likely to be affected by the superimposed modulation process, which could be reminiscent of the
$\pm$
30 % variability in the duration of the sunspot cycle (McCracken et al.
Reference McCracken, Beer, Steinhilber and Abreu2013). Although this point deserves further study, we mention that it may not be in contradiction with the simplified Parker wave dispersion relation, which predicts, for instance, the scaling
${\it\omega}\propto {E_{Z}}^{1/4}$
for the frequency of an
${\it\alpha}{\it\Omega}$
dynamo (Busse & Simitev Reference Busse and Simitev2006; Schrinner et al.
Reference Schrinner, Petitdemange and Dormy2011; Gastine et al.
Reference Gastine, Duarte and Wicht2012). To conclude, we underline that the magnetic activity appears to be concentrated at high latitudes, which is probably related to the much smaller aspect ratio of our models (we recall that we set
${\it\chi}=0.35$
, whereas the solar convective zone has an aspect ratio closer to 0.7); this is also consistent with the fact that the surface magnetic field is predominantly non-axisymmetric at low latitudes, which results in the low values displayed by the butterfly diagrams close to the equator. It should be noted that in all cases the axisymmetric toroidal field migrates towards the poles, which is reminiscent of the poleward branch of solar magnetic activity but contrasts with the equatorward migration of the active latitudes displayed by the solar butterfly diagrams.
4 Conclusion
In this paper, we have examined the hydromagnetic interactions between dynamo modes generated by rotating anelastic convection in a spherical shell. Motivated by direct and indirect observations of solar magnetic activity, our primary aim was to investigate the interactions between modes with different equatorial symmetries. Mathematically, these dynamos display a dynamical behaviour reminiscent of the results obtained with (axisymmetric) mean-field models or low-order systems, with the caveat for the comparison being that the dynamo solutions presented here are dominated by a non-axisymmetric (
$m=1$
) mode. Hemispheric dynamos of the type reported by Grote & Busse (Reference Grote and Busse2000), and studied in more detail by Gallet & Pétrélis (Reference Gallet and Pétrélis2009), have also been found. The present study demonstrates that this hemispheric configuration is also pertinent to the understanding of the dynamics of oscillatory dynamos, and thus could be relevant to explaining the hemispheric magnetic configuration that has been observed on the Sun at the end of the Maunder minimum (Sokoloff & Nesme-Ribes Reference Sokoloff and Nesme-Ribes1994; Beer, Tobias & Weiss Reference Beer, Tobias and Weiss1998; Knobloch et al.
Reference Knobloch, Tobias and Weiss1998).
We stress again that all current direct numerical simulations of convective dynamos – including those here – are far away from what one can imagine as a ‘realistic’ parameter regime. There is, therefore, the question of the robustness of these results. Of course, increasing
$Ra$
for fixed Ekman number should lead to more disordered states, gradually breaking all symmetries. What happens after this is a matter of conjecture/debate. It is possible to argue that for very high
$Ra$
, symmetry is re-established on average in the turbulent state and then similar symmetry-breaking interactions may occur in the averaged equations. Support for this comes from the finding of such interactions in mean-field models, which (despite all their drawbacks) retain the symmetry properties of the underlying system. We note that symmetry arguments are therefore very powerful, and we expect similar behaviour to be observed in Boussinesq and indeed fully compressible models.
Our primary result is that we have demonstrated that the interactions between such modes can lead naturally to a pattern of supermodulation (Arlt & Weiss Reference Arlt and Weiss2014; Weiss & Tobias Reference Weiss and Tobias2016) where the system alternates between modulation with little change of symmetry (with clusters of deep minima) and modulation that involves significant changes in symmetry. We believe that this is the first demonstration of such an interaction between the two types of modulation, leading to supermodulation in the full partial differential equations for convective dynamos.
Acknowledgements
This study was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. Numerical simulations were also carried out at the TGCC computing center (GENCI project x2013046698). R.R. thanks E. Dormy, C. Gissinger, L. Petitdemange and F. Pétrélis for various discussions. The authors thank N. O. Weiss for helpful comments.
