On the Sensitivity of Magnetic Cycles in Global Simulations of Solar-like Stars

We consider the Lantz–Braginsky–Roberts (see Lantz & Fan 1999; Braginsky & Roberts 1995)—or, equivalently, the Lipps–Hemler (Lipps & Hemler 1982, 1985)—set of anelastic equations (see Vasil et al. 2013 for a review on various approximations of anelastic equations). The perturbed equations are written with respect to an ambient state (hereafter denoted by a subscript a) that theoretically may differ from the background state (hereafter denoted by an overbar), around which the generic anelastic equations are derived. The background state is chosen to be isentropic, and both the ambient and background states are assumed to satisfy hydrostatic equilibrium. The anelastic equations written in the stellar rotating frame ${{\boldsymbol{\Omega }}}_{\star }$ are

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\bar{\rho }{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{D}_{t}{\boldsymbol{u}}=-{\rm{\nabla }}\left(\displaystyle \frac{p}{\bar{\rho }}\right)-\displaystyle \frac{{\rm{\Theta }}}{\bar{{\rm{\Theta }}}}{\boldsymbol{g}}-2{{\boldsymbol{\Omega }}}_{\star }\times {\boldsymbol{u}}+\displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{{\mu }_{0}\bar{\rho }},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{D}_{t}{\rm{\Theta }}=-({\boldsymbol{u}}\cdot {\rm{\nabla }}){{\rm{\Theta }}}_{{\rm{a}}}-\displaystyle \frac{{\rm{\Theta }}}{\tau },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{D}_{t}{\boldsymbol{B}}=({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{u}}-({\rm{\nabla }}\cdot {\boldsymbol{u}}){\boldsymbol{B}},\end{eqnarray} \tag{ 4 }$$

where the perturbed quantities are denoted without a prime for the sake of simplicity and D_t is the material derivative. We use standard notation for the basic quantities, i.e., ${\boldsymbol{u}}$ is the fluid velocity, ρ the fluid density, p the fluid pressure, Θ the potential temperature, and ${\boldsymbol{B}}$ the magnetic field. In addition, we use a standard perfect gas equation of state that is linearized around the background state.

In the preceding equations, convection is forced by the combined advection of the unstable ambient entropy profile S_a and a Newtonian cooling term of a characteristic timescale τ (for details, see Prusa et al. 2008; Smolarkiewicz & Charbonneau 2013). The Newtonian cooling damps entropy perturbations over the timescale τ, which is always chosen to exceed the convective overturning time (here, τ = 600 days). This ensures that on long timescales, the model mimics a stellar convection zone remaining in thermal equilibrium (e.g., Cossette et al. 2017). The volumetric forcing of convection is an alternative to forcing via an imposed heat flux across the boundaries of the computational domain. Volumetric forcing directly models the divergence of the diffusive flux across the domain and can be shown to be mathematically and physically equivalent to boundary forcing, provided the ambient state is adequately specified (see Appendix A.1 in Cossette et al. 2017). A relaxed state is then achieved, which results from the equilibrium between the convective heat transport and the volumetric forcing toward the prescribed ambient state. The convective luminosity ends up being set by different control parameters under each of these two modeling approaches; hence any direct comparison must be carried out with care.

Our numerical setup closely follows the anelastic benchmark of Jones et al. (2011). We consider a spherical shell of aspect ratio $\beta ={R}_{\mathrm{top}}/{R}_{\mathrm{bot}}$, and we set $d={R}_{\mathrm{top}}-{R}_{\mathrm{bot}}={R}_{\mathrm{top}}(1-\beta )$. In all simulations presented here, the aspect ratio is solar-like, i.e., β = 0.7. We assume a gravity profile $g={GM}/{r}^{2}$, for which the anelastic equations admit an equilibrium (denoted by overbars) polytropic solution (see, e.g., Jones et al. 2011):

$$\begin{eqnarray}&&\bar{\rho }={\rho }_{i}{\xi }^{n},\bar{P}={P}_{i}{\xi }^{n+1},\bar{T}={T}_{i}\xi ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\xi ={c}_{0}+\displaystyle \frac{{c}_{1}d}{r},\end{eqnarray} \tag{ 6 }$$

where n is the polytropic index, ${\rho }_{i},{P}_{i},{T}_{i}$ are the density, pressure, and temperature at the bottom of the domain, and the constants c₀ and c₁ are given by

$$\begin{eqnarray}&&{c}_{0}=\displaystyle \frac{2\alpha -\beta -1}{1-\beta },{c}_{1}=\displaystyle \frac{(1+\beta )(1-\alpha )}{{(1-\beta )}^{2}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\mathrm{with}\ \alpha =\displaystyle \frac{\beta +1}{\beta \exp ({N}_{\rho }/n)+1},\end{eqnarray} \tag{ 8 }$$

where ${N}_{\rho }=\mathrm{ln}({\rho }_{i}/{\rho }_{o})$ is the number of density scale heights in the layer. The density at the base of the domain is fixed to ρ_i = 200 kg m⁻³ in all the models. The background potential temperature is given by

$$\begin{eqnarray}&&\bar{{\rm{\Theta }}}=\displaystyle \frac{{\bar{P}}^{1/\gamma }}{\bar{\rho }}=\displaystyle \frac{{P}_{i}^{1/\gamma }}{{\rho }_{i}}{\xi }^{(n+1-n\gamma )/\gamma },\end{eqnarray} \tag{ 9 }$$

where the standard adiabatic exponent for a perfect gas is $\gamma ={c}_{p}/{c}_{v}=5/3$. We choose a polytropic exponent n = 3/2 to naturally ensure an isentropic background state with $\bar{{\rm{\Theta }}}={P}_{i}^{1/\gamma }/{\rho }_{i}$, and the equivalent background entropy profile is $\bar{S}={c}_{p}\mathrm{ln}(\bar{{\rm{\Theta }}})$ (we here consider c_p = 3.4 × 10⁴ J kg⁻¹ K⁻¹).

The ambient state needs to be specified only in terms of potential temperature. The entropy jump throughout the domain, ΔS, is used to define the ambient entropy profile by

$$\begin{eqnarray}&&{S}_{{\rm{a}}}(r)=\bar{S}+{\rm{\Delta }}S\displaystyle \frac{{\xi }^{-n}({R}_{\mathrm{top}})-{\xi }^{-n}(r)}{{\xi }^{-n}({R}_{\mathrm{top}})-{\xi }^{-n}({R}_{\mathrm{bot}})},\end{eqnarray} \tag{ 10 }$$

which is related to the ambient potential temperature profile by ${{\rm{\Theta }}}_{{\rm{a}}}=\exp ({S}_{{\rm{a}}}/{c}_{p})$.

The EULAG code is designed to use either Eulerian (flux form) or semi-Lagrangian (advective form) integration schemes (see Prusa et al. 2008; Smolarkiewicz & Charbonneau 2013). In the case presented here, Equations (1)–(4) are written as a set of Eulerian conservation laws and projected on a spherical coordinate system (Prusa & Smolarkiewicz 2003). EULAG solves the evolution equations using the multidimensional positive definite advection transport algorithm (MPDATA), which belongs to the class of nonoscillatory Lax–Wendroff schemes (Smolarkiewicz 2006) and is more specifically a second-order-accurate nonoscillatory forward-in-time template. Since all dissipation is delegated to MPDATA, this provides an implicit turbulence model (Domaradzki et al. 2003). In EULAG, all linear forcing terms are integrated in time using a second-order Crank–Nicholson scheme.

In terms of boundary conditions, we consider stress-free, impermeable boundaries at the top and bottom of the domain such that

$$\begin{eqnarray}&&{u}_{r}=0;{\partial }_{r}({u}_{\theta }/r)={\partial }_{r}({u}_{\varphi }/r)=0\end{eqnarray} \tag{ 11 }$$

at both the top and bottom boundaries. The potential temperature gradient is set to zero on the upper and lower boundaries. The magnetic boundary conditions are chosen to be a perfect conductor at the bottom and a purely radial field at the top, such that

$$\begin{eqnarray}&&{B}_{r}=0;{\partial }_{r}({{rB}}_{\theta })={\partial }_{r}({{rB}}_{\varphi })=0\ \mathrm{at}\ r={R}_{\mathrm{bot}},\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{B}_{\theta }={B}_{\varphi }=0\ \mathrm{at}\ r={R}_{\mathrm{top}}.\end{eqnarray} \tag{ 13 }$$

The simulations presented in this work are initialized with random, small-amplitude perturbations around the background profiles defined in Section 2.2.

The originality of this numerical method lies in the implicit large-eddy simulation (ILES) approach used in the EULAG-MHD code. While this is shown to minimize dissipation at large scales, it also renders comparison with other simulation results less straightforward. Following Strugarek et al. (2016), we characterize our simulations with effective dissipation coefficients deduced from a spectral analysis of the simulation results. The details of this technical procedure are given in Appendix B. The numerical dissipation of the various MHD fields fits a classical Laplacian operator very well for the smallest scales of the simulation (spherical harmonics l > 25, typically). We find that the effective dissipation coefficients are of the order of 10⁸ m² s⁻¹, which is what is typically expected for the relatively coarse numerical resolution adopted in this work (Strugarek et al. 2016). The Prandtl number ${\rm{\Pr }}={\nu }_{\mathrm{eff}}/{\kappa }_{\mathrm{eff}}$ and magnetic Prandtl number ${\rm{Pm}}={\nu }_{\mathrm{eff}}/{\eta }_{\mathrm{eff}}$ are generally slightly larger than 1, which is what can be expected from such a numerical approach when using the same subgrid-scale model for all quantities (Domaradzki et al. 2003). Comparison with other simulation results must be carried out with care, though, as the dissipation coefficients reported in Table 4 (see Appendix B) well characterize the small scales in the simulations but are by no means representative of the dissipation of the largest spatial scales associated with the magnetic cycles. Moreover, locally, the dissipation can be significantly different from that of a Laplacian operator (for more on this see Strugarek et al. 2016).

We present here a series of 17 simulations using the same numerical grid and varying the rotation rate, entropy contrast, and number of density scale heights of the convection zone. The parameters of the models are listed in Table 1. All the simulations presented here exhibit either a solar-like differential rotation self-sustained by convective turbulence under the influence of rotation (e.g., Brun & Toomre 2002) or an antisolar differential rotation when the Rossby number reaches values higher than 1 (see panel (h) in Figure 1 and Brun et al. 2017). Several dynamo states are achieved in this series with distinctive properties, as shown in Figure 1. Some of the simulations display a magnetic cycle deep-seated in the convective envelope, such as model O2 (panel a). A short cycle located close to the equator and in subsurface layers is also observed in some simulations, such as model O6 (panel d2). Finally, in some simulations a stable wreath of magnetic field builds up in the bottom part of the convection zone, showing only little time variability due to the convective motions perturbing it (see also Brown et al. 2010; Nelson et al. 2013). For instance, this is the case for model sO1, shown in the bottom panels. Notice also that this particular case possesses an antisolar differential rotation profile. Finally, some models (not shown here; see, e.g., Figure 10) also present stochastic reversals of the deep-seated magnetic field.

Zoom In Zoom Out Reset image size

Table 1.  Input Parameters and Global Properties of the Models

Case	${{\rm{\Omega }}}_{\star }/{{\rm{\Omega }}}_{\odot }$	ΔS (J kg–1 K–1)	Nρ	${\mathrm{Ro}}_{{\rm{b}}}$	${L}_{{\rm{b}}}/{L}_{\odot }$	Cycle	${\mathrm{Ro}}_{{\rm{t}}}$	${L}_{{\rm{t}}}/{L}_{\odot }$	Short Cycle
sO1	0.33	1.00	3.22	1.39 ± 0.12	0.29 ± 0.01	No	4.86 ± 0.69	1.17	No
O1	0.55	1.00	3.22	0.81 ± 0.09	0.31 ± 0.01	Yes	2.73 ± 0.38	1.21	No
O2	0.69	1.00	3.22	0.63 ± 0.06	0.33 ± 0.01	Yes	2.12 ± 0.30	1.26	No
O3	0.83	1.00	3.22	0.50 ± 0.05	0.34 ± 0.01	Yes	1.69 ± 0.23	1.28	No
O4	1.00	1.00	3.22	0.39 ± 0.05	0.31 ± 0.01	Yes	1.26 ± 0.16	1.09	No
O5	1.10	1.00	3.22	0.34 ± 0.05	0.30 ± 0.01	Yes	1.11 ± 0.14	1.07	No
O6	1.25	1.00	3.22	0.29 ± 0.05	0.27 ± 0.01	Yes	0.93 ± 0.10	0.92	Yes
O7	1.66	1.00	3.22	0.21 ± 0.04	0.27 ± 0.01	No	0.64 ± 0.07	0.86	Yes
O8	3.00	1.00	3.22	0.10 ± 0.02	0.28 ± 0.02	No	0.32 ± 0.03	0.78	Yes
O9	5.52	1.00	3.22	0.04 ± 0.01	0.20 ± 0.02	No	0.12 ± 0.02	0.39	Yes
S1	1.10	0.80	3.22	0.31 ± 0.05	0.19 ± 0.01	Yes	0.92 ± 0.10	0.60	Yes
S2a	0.55	1.50	3.22	1.01 ± 0.07	0.55 ± 0.01	Yes	3.88 ± 0.52	2.60	No
S2b	1.10	1.50	3.22	0.44 ± 0.04	0.56 ± 0.02	Yes	1.54 ± 0.20	2.37	No
R0	1.10	1.00	1.00	0.87 ± 0.13	11.34 ± 2.82	No	0.82 ± 0.04	17.46	No
R1	1.10	1.00	2.00	0.59 ± 0.05	1.89 ± 0.22	Yes	1.08 ± 0.09	4.58	No
R2	1.10	1.00	3.50	0.31 ± 0.05	0.19 ± 0.01	Yes	1.05 ± 0.13	0.66	No
R3	1.10	1.00	4.00	0.25 ± 0.05	0.10 ± 0.01	No	0.94 ± 0.13	0.32	Yes

Note. The subscripts b and t indicate that the quantities are derived from averages in longitude and latitude and over radial bands $r\in [0.75{R}_{\mathrm{top}},0.8{R}_{\mathrm{top}}]$ and $r\in [0.95{R}_{\mathrm{top}},0.97{R}_{\mathrm{top}}]$.

Download table as:  ASCIITypeset image

The various dynamo states realized in our set of simulations occur in particular ranges of the adimensional fluid Rossby number. It can be defined as

$$\begin{eqnarray}&&\mathrm{Ro}=\displaystyle \frac{{\langle | {\boldsymbol{\nabla }}\times {\boldsymbol{u}}| \rangle }_{\theta ,\varphi ,t}}{2{{\rm{\Omega }}}_{\star }},\end{eqnarray} \tag{ 14 }$$

where ${\langle \rangle }_{\theta ,\varphi ,t}$ stands for the average over latitude, longitude, and time. The Rossby number quantifies the importance of rotation in the dynamics of the system. We define here two Rossby numbers averaged over the $[0.75{R}_{\mathrm{top}}$,$0.8{R}_{\mathrm{top}}]$ domain (${\mathrm{Ro}}_{{\rm{b}}}$, where the deep-seated magnetic field lies) and over the $[0.95{R}_{\mathrm{top}}$,$0.97{R}_{\mathrm{top}}]$ domain (${\mathrm{Ro}}_{{\rm{t}}}$, where a short cycle is observed in the simulations developing it). Both Rossby numbers are reported in Table 1, and note that the averaging intervals are chosen to safely exclude the radial boundaries. Our series of simulations covers a range of Rossby numbers Ro_b between 0.04 and 1.39 and Ro_t between 0.12 and 4.86 (see Table 1). The error bars in the Rossby numbers are taken from the standard deviation of Ro over the radial averaging interval. For the sake of completeness, we also report in Table 1 the convective luminosities, defined as $L(r)=4\pi {r}^{2}\bar{\rho }{c}_{P}{\langle {u}_{r}T\rangle }_{\theta ,\varphi ,t}$, at these two locations (with T representing the temperature fluctuations with respect to the ambient state).

We summarize in Figure 2 our ensemble of dynamo states, where the ratio between the time-averaged kinetic energy contained in the differential rotation (DRKE; see Appendix A) and the total kinetic energy (KE) is shown as a function of the Rossby number Ro_b. The color of the symbols indicates the dynamo state each model achieves, which is also reported at the top of the figure. We immediately see that three regimes of Rossby numbers exist in our set of simulations. When Ro_b is below 0.3, a short magnetic cycle is observed near the surface and close to the equator (blue symbols). For 0.25 ≲ Ro_b ≲ 1, a deep-seated magnetic cycle is realized (red symbols). Finally, as soon as Ro_b > 1, the magnetic cycles disappear, and stable wreaths of magnetic field are observed at the base of the convective envelope (black symbols). Note that two of the simulations in our set (namely, O6 and S1) possess both a deep-seated cycle and a short cycle at the same time, a feature also reported in the millennium EULAG simulation (see, e.g., Beaudoin et al. 2016 and references therein). The short cycle close to the surface presents a poleward migration (panel (d2) in Figure 1), which is consistent here with the Parker–Yoshimura (Parker 1955; Yoshimura 1975) rule of a dynamo wave propagating along the near-cylindrical isocontours of differential rotation. The deep-seated cycle does not follow such a rule and originates from the nonlinear feedback of the Lorentz force on the differential rotation itself, as will be made clear in Sections 4 and 5.

Zoom In Zoom Out Reset image size

It has been noted before that cyclic dynamos are observed over a limited range of Rossby numbers, starting with the pioneering work of Gilman (1983; see, e.g., his Figure 31). We observe here a similar trend, where the low Rossby number transition to cyclic solutions corresponds to a threshold value of DRKE/KE of 0.6 (Ro_b ∼ 0.25). The high Rossby number transition is found at a lower threshold (DRKE/KE ∼ 0.1 to 0.2) and is found to occur when the differential rotation switches from being solar-like (fast equator, slow poles; see panel (e) in Figure 1) to being antisolar (slow equator, fast higher latitudes; e.g., panel (h) in Figure 1). We discuss this transition in more detail in Section 6.4. We now briefly characterize the energetics of our simulation set before detailing the magnetic cycles themselves in subsequent sections.

The simulations in our set strongly vary in their kinetic energy and magnetic energy (ME) properties, which we list in Table 2 (see Appendix A for the definitions of the various energies). We find that the differential rotation (DRKE, third column) accounts for a significant part of the kinetic energy in most cases—between 10% and 60% of the total energy (KE, second column). The convective kinetic energy (CKE, fourth column) accounts for more than 40% of the KE in all cases.

The majority of our simulations develop a dynamo in subequipartition, with ME (fifth column) of the order of 10% of the KE. The two simulations (O8 and O9) at the highest rotation rates, for which the specifics of the dynamo solution are discussed in Section 6.3, achieve a superequipartition regime up to ME ∼ 5 KE (magnetostrophic state; see Brun & Browning 2017). In all other cases, 30% to 50% of the total ME is distributed roughly equally over the large-scale axisymmetric toroidal (TME, sixth column) and poloidal (PME, seventh column) fields. The fluctuating ME spectrum peaks at mid-to-small scale (typically around the spherical harmonic degree l = 30). As a result, the large-scale field is dominated by the axisymmetric dipole and quadrupole, while the small-scale field is dominated by nonaxisymmetric modes.

In all simulations presenting a deep-seated primary cycle, we observe temporal modulations of the mean differential rotation correlated with the magnetic cycle (see Figure 4). We propose here that these modulations play a major role in allowing the magnetic cycle to operate and in setting its period. Let us first characterize in more detail this modulation, shown in the top panels of Figure 5 for model O5 (the left panels represent time–latitude diagrams at r = 0.75R_top; the right panels display radius–latitude diagrams at latitude 44°). The modulation of the differential rotation $\delta {\rm{\Omega }}/{\langle {\rm{\Omega }}\rangle }_{t}$ amounts to ∼1% of the total differential rotation in case O5, which corresponds to ∼50% when the rotating frame Ω_⋆ is subtracted. Systematic accelerations (in red, contoured for 0.25% in black, and highlighted by white arrows to guide the eye) are correlated with the magnetic cycle, propagating toward the equator (left panel) and toward the surface (right panel).

Zoom In Zoom Out Reset image size

We first demonstrate that these modulations originate from the feedback of the magnetic field on the differential rotation. In the second row of Figure 5, we show the energy transfer due to the Lorentz force toward the toroidal kinetic energy (last term in Equation (2); see also Equations (B.12) and (B.13) in Nelson et al. 2013). A positive (red) transfer means that energy is effectively transferred from magnetic to kinetic, and negative (blue) indicates the opposite. We observe that locally, at midlatitudes and in the lower part of the convection zone, energy is transferred from ME to kinetic energy. The energy transfer furthermore peaks just before the cyclic acceleration in the upper panel (the contours of $\delta {\rm{\Omega }}$ are overlaid here for clarity). The positive transfer that is sustained over roughly 100 solar days is indeed able to induce a flow of magnitude of ∼10 m s⁻¹, which is on par with the acceleration displayed in the top panels. Thus, the acceleration phases observed in our simulations are indeed of magnetic origin. Note that these modulations differ from the ones observed in the millennium EULAG simulation (Beaudoin et al. 2016), which were shown to originate indirectly from magnetically induced meridional circulation modulations (Passos et al. 2012). Analogous modulations of Ω were also reported in the simulations of Augustson et al. (2015) and were shown to be at the origin of the equatorward migration of the magnetic structures. Here equatorial migration is very weak, but the direct feedback of the magnetic field on the differential rotation is stronger (see energies in Figure 3).

However, do these modulations affect dynamo action? We display the equivalent of the dynamo omega effect in our three-dimensional (3D) nonlinear simulations in the third row of Figure 5. The mean shearing of the magnetic field (or Ω-effect), ${({\langle {\boldsymbol{B}}\rangle }_{p}\cdot {\boldsymbol{\nabla }}){\langle {\boldsymbol{u}}\rangle }_{\varphi }| }_{\varphi }$, is represented in time–latitude and time–radius diagrams, with the same contours of δΩ overlaid in white. We observe that during the acceleration phases of the differential rotation, the mean shear severely drops in amplitude and ultimately locally changes signs (going from yellow to blue or blue to yellow) due to a change of sign in the latitudinal gradient of the differential rotation (see, e.g., Figure S8 in Strugarek et al. 2017). The azimuthal field consequently ceases to be sustained (last row in Figure 5), and as a result the poloidal field (not shown here; see Strugarek et al. 2017) also loses its source and rapidly decreases.

As the magnetic field wanes, the acceleration of the differential rotation fades away as the energy transfer mediated by the Lorentz force severely weakens (see second row of panels). During this phase, the mean shear (third panel) is still of opposite polarity, and thus a weak azimuthal field of opposite polarity is generated. This new seed azimuthal field is strong enough that a poloidal field (also of opposite polarity) is slowly generated as well. The differential rotation finally recovers its ground state (blue phases in the first panel), and the dynamo process is able to amplify again both components of the field, albeit with opposite polarity to the previous configuration. Eventually, the field becomes strong enough such that an efficient energy transfer again locally perturbs the differential rotation, which triggers a new reversal. This mechanism repeats over and over, and does so in all simulations of our set that generate a deep-seated magnetic cycle.

Most of the dynamo states achieved in this series of simulations present a quadrupolar symmetry—e.g., an azimuthal field symmetric with respect to the equator (see Figure 4). This appears to be contrary to the solar dynamo, which is dominated by a dipolar symmetry. We characterize the dynamo solutions by calculating the energy repartition between dipolar ("antisymmetric") and quadrupolar ("symmetric") families (Roberts & Stix 1972; McFadden et al. 1991; Gubbins & Zhang 1993). Projecting a vector on the vectorial spherical harmonic basis (Rieutord 1987; Mathis & Zahn 2005; Strugarek et al. 2013)

$$\begin{eqnarray}&&{\boldsymbol{X}}=\displaystyle \sum _{l,m}{A}_{m}^{l}{{\boldsymbol{ \mathcal R }}}_{l}^{m}+{B}_{m}^{l}{{\boldsymbol{ \mathcal S }}}_{l}^{m}+{C}_{m}^{l}{{\boldsymbol{ \mathcal T }}}_{l}^{m},\end{eqnarray} \tag{ 15 }$$

we can separate its quadrupolar and dipolar components as

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{{\rm{Q}}}\,=\,\cdots \,+\,{A}_{m}^{m}{{\boldsymbol{ \mathcal R }}}_{m}^{m}+{B}_{m}^{m}{{\boldsymbol{ \mathcal S }}}_{m}^{m}+{C}_{m}^{m+1}{{\boldsymbol{ \mathcal T }}}_{m+1}^{m}\,+\,\cdots ,\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{{\boldsymbol{X}}}_{{\rm{D}}}\,=\,\cdots \,+\,{A}_{m}^{m+1}{{\boldsymbol{ \mathcal R }}}_{m+1}^{m}+{B}_{m}^{m+1}{{\boldsymbol{ \mathcal S }}}_{m+1}^{m}+{C}_{m}^{m}{{\boldsymbol{ \mathcal T }}}_{m}^{m}+\cdots .\end{eqnarray} \tag{ 17 }$$

Based on this decomposition, we further define the parity of the magnetic field (Tobias 1997) with the ratio of magnetic energies E_Q and E_D:

$$\begin{eqnarray}&&P=\displaystyle \frac{{E}_{{\rm{Q}}}-{E}_{{\rm{D}}}}{{E}_{{\rm{Q}}}+{E}_{{\rm{D}}}}.\end{eqnarray} \tag{ 18 }$$

The parity P is equal to +1 when the field is mostly quadrupolar and to −1 when it is mostly dipolar. The parity of two representative simulations (O5 and O2) at the top of the domain is shown in Figure 6, along with the total ME (Equation 23).

Zoom In Zoom Out Reset image size

Case O5 (Figure 4) exhibits a dynamo solution dominated by quadrupolar symmetry at all times. When the total ME peaks, the parity falls to zero in this case (dipoles and quadrupoles have the same energy). This is somehow the opposite of what is observed in the Sun, with a dipolar symmetry dominating at cycle minimum and a quadrupolar symmetry dominating at cycle maximum (DeRosa et al. 2012). Case O2 presents an interesting beating between the two families over a timescale longer than the magnetic cycle itself. Between t = 590 yr and t = 650 yr, quadrupolar symmetry dominates the field, as in model O5. However, model O2 behaves more like the Sun during this period, with dipolar symmetry most prominent during cycle minimum (see, e.g., DeRosa et al. 2012). After t = 650 yr, the polarity trend is reversed, and the field is dominated by dipolar symmetry.

Coexisting modulation of magnetic parity and large-scale flows is a well-known property of nonlinear dynamos, which has been explored at length using dynamical system approaches and geometrically simplified mean-field and mean-field-like dynamo models (e.g., Beer et al. 1998; Knobloch et al. 1998; Weiss & Tobias 2000, and references therein) and to some extent in nonlinear 3D simulations (Augustson et al. 2015; Raynaud & Tobias 2016). The development of joint modulations of parity and large-scale flow in our simulations is consistent with the magnetic polarity reversal scenario proposed in the following section, under which energy exchange between ME and kinetic energy reservoirs is a key ingredient in the evolution of the magnetic cycle.

In order to characterize the cycle period, we need a robust methodology to extract the most significant period out of a potentially multiperiodic signal. One such method was recently successfully applied to similar simulations by Käpylä et al. (2016) based on the so-called empirical-mode decomposition method (we refer the reader to Appendix C and Käpylä et al. 2016 for a more in-depth discussion of this method). We base our analysis on the libeemd library (Luukko et al. 2015) to make use of the complete variant of ensemble empirical-mode decomposition (CEEMDAN). This method automatically decomposes a signal on a series of intrinsic mode functions, leveraging the addition of noise to better separate different frequencies in the spectra.

Based on our analysis of the cases presented in the present study, we perform CEEMDAN analysis at two representative radii and latitudes, $[{R}_{{\rm{b}}}=0.72{R}_{\mathrm{top}},\theta =45^\circ ]$ and $[{R}_{{\rm{t}}}=0.98{R}_{\mathrm{top}},\theta =25^\circ ]$, on the longitudinally averaged azimuthal component of the magnetic field. At the top position, the CEEMDAN analysis is applied to the detrended ${\langle {B}_{\varphi }\rangle }_{\varphi }$ (see panel (d2) in Figure 1). For each position, we identify the most energetic mode as the representative cyclic mode (see Appendix C). The period (P), energy (E), and robustness (R) of the modes are given in Table 3. Their energy is calculated relative to the total energy of all the decomposed modes, and the robustness is estimated by comparing the results with the results of CEEMDAN applied to a white-noise time series with the same length, mean, and standard deviation as the analyzed simulation time series (robustness increases as ${\text{}}R\to$ 1). The signal of the shorter cycle at R_t is generally weaker than the signal of the primary cycle at R_b. As a result, any cycle for which all modes have either (1) a relative energy less than 0.5 for the primary cycle or 0.3 for the short cycle or (2) a robustness smaller than 0.2 is considered undetected (symbol "..." in Table 3). The error in the cycle period is estimated based on the standard deviation of the corresponding empirical mode. We have reproduced this analysis with the latitude slightly varied at the chosen depths $({R}_{{\rm{b}}},{R}_{{\rm{t}}})$ without finding any significant changes in our results, all detected cycle periods being well located within their error bars.

We find that, in our sample of models, deep-seated magnetic cycles materialize when 0.25 ≲ Ro_b ≲ 1, and short cycles occur when Ro_b ≲ 0.3 (or, equivalently, Ro_t ≲ 1), as shown in Figure 2. Only two models in our set manage to simultaneously sustain the two types of cycles (models O6 and S1). When the Rossby number is significantly larger than 1, no cyclic magnetism could be realized in our simulations (we will come back to this point in Section 6.4).

The period of the deep-seated magnetic cycle decreases with an increasing Rossby number (Tables 1 and 3). Furthermore, the cycle period decreases like a power law with an exponent of −1.6 ± 0.14 when normalized to the rotation period of the star, as shown in Figure 7. This exponent is slightly smaller than the one derived in Strugarek et al. (2017), where a smaller sample of cyclic models was used. A slope smaller than −1 nevertheless remains a robust feature of our nonlinear 3D simulations. All models fit relatively well a single power-law dependency on the Rossby number Ro_b.

Zoom In Zoom Out Reset image size

Model R1 (green dot in the middle of Figure 7) is the most distant point to the power-law fit, suggesting that the density contrast in the convective envelope slightly affects the power-law dependency of the cycle period on the Rossby number found in this work. The strong effect of stratification on the dynamo state was already reported by Käpylä et al. (2014), who showed that the transition to a cyclic dynamo occurs at higher Rossby numbers as density stratification is increased. In the present case, additional cyclic solutions with larger density contrasts would be required to precisely assess the effect of this parameter on the power-law fit. However, the density contrast in stellar convective envelopes changes according to their aspect ratio. As a result, both effects ultimately need to be studied jointly to fully characterize how magnetic cycles depend on density contrast, which we leave here for future work.

The nonlinear character of the dynamo action detailed in Section 5 was proposed to be at the origin of the negative slope observed in Figure 7 by Strugarek et al. (2017). Using this interpretation, the magnetic cycle length is set by the Lorentz force feedback on the large-scale differential rotation: the stronger the feedback is, the shorter the cycle will be. This interpretation can be further tested by calculating the timescale τ_M associated with the Maxwell torque in our simulations. By equating the time derivative in Equation (2) with the right-hand-side large-scale Lorentz force, one may estimate τ_M by

$$\begin{eqnarray}&&{\tau }_{{\rm{M}}}\approx \sqrt{\displaystyle \frac{\mathrm{DRKE}}{\mathrm{TME}\cdot \mathrm{PME}}}\cdot d\cdot \sqrt{\displaystyle \frac{{{ \mathcal V }}_{\mathrm{CZ}}\hat{\rho }}{2}},\end{eqnarray} \tag{ 19 }$$

where ${{ \mathcal V }}_{\mathrm{CZ}}$ is the volume of the convection zone, $d\,={R}_{\mathrm{top}}-{R}_{\mathrm{bot}}$ is the depth of the convection zone, and $\hat{\rho }$ is the background density averaged over the convection zone. We display in Figure 8 the cycle period of the cyclic models as a function of τ_M estimated with Equation (19). We find a positive correlation between the two timescales: if the Maxwell torque timescale is short, the magnetic cycle is short as well. In addition, τ_M is systematically shorter than P_cyc in our simulations, showing that the magnetic feedback is indeed fast enough to act over the magnetic cycle timescale. This further strengthens our interpretation of the nonlinear dynamo action occurring in our set of simulations.

Zoom In Zoom Out Reset image size

Two additional remarkable features can be noted from Figure 8. First, there seems to be a saturation of the cycle period as τ_M increases past a threshold. This is not necessarily a surprise: as τ_M increases, the feedback of the large-scale magnetic field on the differential rotation weakens to ultimately become negligible (or, equivalently, the feedback occurs on a timescale too long compared to the convective turnover timescale). As a result, we expect our model to possess a maximal cycle length, leading to such a thresholding effect. Second, one simulation clearly stands out of this P_cyc–τ_M relationship (red circle at the bottom of Figure 8). This particular simulation is model S2a, which is the cyclic solution with the largest Rossby number in our sample (Ro_b ∼ 1; see Table 1). This particular model is at the edge of the high Rossby number transition and is starting to exhibit an antisolar differential rotation profile. It is well known that magnetic torques and stresses can affect the delicate balance establishing the differential rotation profile (e.g., Brun et al. 2004; Fan & Fang 2014; Gastine et al. 2014; Karak et al. 2015). As a result, models lying very close to the Rossby transition (such as model S2a) are not necessarily expected to fall on the same trend as the other cyclic models in our sample. Nonetheless, model S2a follows the robust Rossby number trend identified in Figure 7.

Several of our models do not produce a deep-seated magnetic cycle (see Figure 2), which is found to materialize only in a relatively narrow Rossby number range. A sharp transition is observed at Ro_b ∼ 0.25 and Ro_b ∼ 1.0, where the Maxwell torque timescale diverges and becomes too long to participate in setting the period of the deep-seated cycle. We now detail what occurs at both of these transitions.

As we increase the rotation rate of the convective envelope or its density contrast, the Rossby number decreases (Table 1). We observe that when Ro_b ≲ 0.25, the dynamo loses its deep-seated cycle (see Table 3 and Figure 2). We display in the left panels of Figure 9 three cases (R3, O7, and O8) tracing the transition of our dynamo solutions from a cyclic solution (Ro > 0.25) to a randomly reversing dynamo (cases R3 and O7) to solutions developing stable wreaths of toroidal field at the base of the convective envelope (cases O8 and O9). The latter cases are morphologically similar to the simulations of Brown et al. (2010), who originally showed that strong azimuthal magnetic structures could be sustained within turbulent convective envelopes.

Zoom In Zoom Out Reset image size

As the Rossby number decreases, the total ME increases considerably in proportion to the DRKE (see right panel in Figure 9). Meanwhile, the azimuthal magnetic field energy markedly decreases relative to the total ME (see Table 2). As a result, the fluctuating (i.e., nonaxisymmetric) ME dominates the total ME for these models. The large-scale feedback from the Maxwell torque is thus much less efficient in these cases, which is why the dynamo does not operate in the same way as for Ro_b ≳ 0.25.

Simulations operating near the low Rossby number transition develop a short-period magnetic cycle (blue circles in the right panel of Figure 9) near the top of the domain. We display the short cycle period normalized to the rotation rate of the model as a function of the Rossby number Ro_t in Figure 10. Even though we only have six models exhibiting such a short cycle, a clear decreasing trend with the local Rossby number is also observed (compared to the deep-seated magnetic cycle in Figure 7). The power law is shallower for the short cycle (keeping in mind that the exact slope value is strongly influenced by a few simulations in the sample) but is still characterized by a slope steeper than −1, which suggests a trend similar to the that of the deep-seated cycle. Like Käpylä et al. (2014), we find that stronger stratification favors such short cycles near the top of the domain (compare, e.g., models R2 and R3), which here simply reflect a decrease of the local Rossby number Ro_t when the number of density scale heights embedded in the model is increased. Interestingly, the short cycles are observed as soon as Ro_t is smaller than 1, and we see no hint of short-cycle disappearance at low Ro_t.

Zoom In Zoom Out Reset image size

In our set of simulations, the deep-seated magnetic cycles also disappear when the Rossby number exceeds unity. A structural change in the differential rotation occurs when Ro ∼ 1 (e.g., Brun et al. 2017), with the differential rotation switching from solar-like to antisolar (fast to slow equator, slow to fast high latitudes). This transition is observed, for instance, in panel (h) of Figure 1 for model sO1, for which a strong antisolar differential rotation starts to appear. Model sO1 interestingly develops a stable azimuthal wreath at the base of the convection zone, akin to the low Rossby number cases. In this case, though, the differential rotation and convective energies are balanced, and the ME is less than 10% of the KE (see Table 2). The ME is dominated by the nonaxisymmetric components of the field, which points to a decrease of the large-scale field as the Rossby number of our models is increased. This is further confirmed when computing the dipole strength f_dip, defined here as the energy of the axisymmetric dipole divided by the total ME, and shown in Figure 11. We observe that the dipole strength decreases as a function of the Rossby number Ro_b, and report a hint of saturation at a high Rossby number around f_dip ∼ 0.1–0.2. Additional models at higher Rossby numbers are nevertheless required to fully investigate this trend in the context of the effect of a dynamo transition for old, slowly rotating stars (e.g., Metcalfe & van Saders 2017), which we leave here for future work.

Zoom In Zoom Out Reset image size