{'fraction_non_alphanumeric': 0.03337554269175109, 'fraction_numerical': 0.014969247467438496, 'mean_word_length': 4.6654710052096675, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the sensitivity of magnetic cycles in global simulations of solar-like stars

The periods of magnetic activity cycles in the Sun and solar-type stars do not exhibit a simple or even single trend with respect to rotation rate or luminosity. Dynamo models can be used to interpret this diversity, and can ultimately help us understand why some solar-like stars do not exhibit a magnetic cycle, whereas some do, and for the latter what physical mechanisms set their magnetic cycle period. Three-dimensional non-linear magnetohydrodynamical simulations present the advantage of having only a small number of tunable parameters, and produce in a dynamically self-consistent manner the flows and the dynamo magnetic fields pervading stellar interiors. We conducted a series of such simulations within the EULAG-MHD framework, varying the rotation rate and luminosity of the modeled solar-like convective envelopes. We find decadal magnetic cycles when the Rossby number near the base of the convection zone is moderate (typically between 0.25 and 1). Secondary, shorter cycles located at the top of the convective envelope close to the equator are also observed in our numerical experiments, when the local Rossby number is lower than 1. The deep-seated dynamo sustained in these numerical experiments is fundamentally non-linear, in that it is the feedback of the large-scale magnetic field on the large-scale differential rotation that sets the magnetic cycle period. The cycle period is found to decrease with the Rossby number, which offers an alternative theoretical explanation to the variety of activity cycles observed in solar-like stars.


INTRODUCTION
Over fifty years ago, Olin Wilson initiated the Mt Wilson stellar activity monitoring program relying on a measure of emission in the core of the Calcium H+K lines. This remarkable observational effort was further extended by other similar observational programs, some still ongoing (see Hall 2008 and references therein). Spatially-resolved solar observations had already long shown that emission in the Ca H+K lines is closely associated with surface magnetism, in particular active regions and surrounding plages. A large fraction of the solar-type stars in the Mt Wilson sample were found to exhibit temporal variability in H+K emission, often in the form of more or less regular defined cycles, with periods ranging from 5 to some 20 yr (Wilson 1978;Baliunas et al. 1995; see also Hall et al. 2007). Such cyclic activity is mostly observed in stars of masses, rotation rates and evolutionary states similar to the Sun. The natural interpretation is to consider these as stellar analogs of the 11yr solar magnetic activity cycle. However, even amongst antoine.strugarek@cea.fr solar-type stars in the sample, some show highly irregular, non-cyclic variability or very little variability over decades; magnetism is clearly ubiquitous among solartype stars, but regular activity cycles are not.
Efforts at understanding these stellar activity cycles have been hampered by the fact that most stars in which cycles have been measured are field stars, for which the precise evolutionary status (age, luminosity, etc.) was difficult to pin down accurately. This situation has changed in the past decade, following accurate parallax determination by space missions such as Hipparcos and GAIA, and in some cases by asteroseismic analyses made possible by photometric observations from the Corot and Kepler missions (e.g. Mathur et al. 2012;Chaplin et al. 2014;do Nascimento et al. 2014).
Again quite naturally, physical interpretation of observed stellar cycles has been sought within solar dynamo theory (e.g. Noyes et al. 1984;Baliunas et al. 1996;Saar & Brandenburg 1999). Working on 13 Mt-Wilson stars with well-defined cycle periods, Noyes et al. (1984) used simple mean-field dynamo models to establish empirical relationship linking rotation period (P rot ), cycle period (P cyc ), and spectral type. The Rossby number (Ro = P rot /τ c ), a dimensionless quantity defined as the ratio of rotation period to the convective turnover time (τ c ), emerged from their analysis as the key parameter. From the dynamo point of view this is satisfying, as in this context the Rossby number measures the influence of rotation on convection. Large-scale differential rotation and cyclonic turbulence, in turn, are two key inductive processes in solar dynamo theory, and both require this influence (see, e.g., Miesch & Toomre 2009; Charbonneau 2014 and references therein). Noyes et al. (1984) showed that cycle periods in their reduced stellar sample could be tolerably well-represented by the relationship P cyc ∝ Ro 5/4 ; yet their modeling work also showed that any such relationship is sensitive to modeling choices made in setting up the dynamo model.
In a similar vein, but working with an expanded sample and longer datasets, Saar & Brandenburg (1999) showed that the significant scatter about the Noyes et al. (1984) empirical relationship could be markedly reduced by first dividing the sample into "active" and "inactive" stars (based on the overall emission levels), leading to two "branches" in the P cyc vs Ro diagram. In their analyses both branches showed an increase of P cyc with Ro, but with different slopes. Bohm Vitense (2007) further suggested an evolutionary scenario whereby each branch is associated with a distinct dynamo mode within the star, a proposal buttressed by the fact that some stars showed dual-period cyclic variability, with the two periods each lying on one of the two activity branches. As they age and spin down, young rapidly rotating stars initially populating the active branch eventually switch to a less efficient dynamo mode, and transit to the inactive branch. This scenario, plausible and compatible with observed stellar activity cycles, remains unsatisfactory in one important aspect: the Sun, arguably the best data point of the lot, is left hanging far in between the two stellar activity branches, forcing one to assume that our star just happens to be in its transition phase between the two branches (for a new-and-improved version of this type of scenario, see Metcalfe & van Saders 2017). Larger samples of stars with a well-defined cyclic activity are today emerging, suggesting that this gap between the branches is actually a selection bias and questioning the rationale behind this interpretation (e.g., see Boro Saikia et al. 2018).
The use of dynamo models to interpret stellar activity cycles, as pioneered by Noyes et al. (1984), requires a number of input "ingredients". Broadly defined and at minimum: (1) Choosing a specific class of solar/stellar dynamo model; (2) Choosing a nonlinear amplitudelimiting mechanism; (3) Specifying how relevant inductive flows (differential rotation, meridional circulation, cyclonic turbulence) vary as a function of rotation rate, structural evolution, etc. Even in the solar case there is currently no consensus on the first two items (see, e.g., discussions in Charbonneau 2010; Brun & Browning 2017). For the present day Sun the differential rotation (point 3) is well-constrained by helioseismology. However, in the case of stars, where comparable asteroseismic constraints are difficult to obtain, even under physically reasonable assumptions there remains sufficient leeway to produce almost any P cyc vs P rot relationship (see, e.g., Charbonneau & Saar 2001;Jouve et al. 2010).
Global magnetohydrodynamical (MHD) simulations of solar/stellar convection and associated dynamo action in principle bypass these problems. Available computational resources limit the range of spatial and temporal scales that can be captured by such simulations, with the consequence that some potentially important physical mechanisms cannot be included, for example the surface decay of active regions. Nevertheless, the generation of large-scale flows by rotationally-influenced convection as well as the non-linear backreaction of magnetism on internal flows are both captured in a dynamically correct manner at all resolved scales.
As a case in point, the transition from "solar-like" differential rotation (equatorial regions rotating faster than the poles) to "anti-solar" differential rotation (slower equator) is now well-understood. Not surprisingly, the Rossby number emerges again as a key parameter, with high-Ro ( 1) convection leading to anti-solar differential rotation, and intermediate-Ro (∼ 0.1) convection producing solar-like differential rotation (see Guerrero et al. 2013;Gastine et al. 2014;Featherstone & Miesch 2015;Mabuchi et al. 2015;. Interestingly, these simulation results also indicate that the Sun may be near the tipping point between these two rotational regimes. It must be noted however that the limited range of spatial scales in numerical simulations may lead one to over-estimate the amplitude of convective flows (see Hanasoge et al. 2015;Greer et al. 2015 and references therein), which would in turn change the critical stellar Rossby number at which this transition occur. One could nevertheless infer from these simulation results that a simple, monotonic relationship between cycle and rotation periods is not to be expected, considering that differential rotation is a key inductive process in the vast majority of extant solar/stellar dynamo scenarios.
Finding counterparts of stellar magnetic cycles in global MHD simulations has proven a challenging undertaking. While thermally-driven turbulent convection has no difficulty producing magnetic fields (e.g. Brun et al. 2004), their solar-like organization into a well-structured and cyclically varying large-scale component has only been achieved relatively recently (Ghizaru et al. 2010;Brown et al. 2011;Racine et al. 2011;Nelson et al. 2011Nelson et al. , 2013Passos & Charbonneau 2014;Fan & Fang 2014;Masada et al. 2013;Käpylä et al. 2012Käpylä et al. , 2016Augustson et al. 2015). There are actually few adjustable physical parameters in such simulations once a background solar/stellar structure is set: typically luminosity and rotation are the two primary physical knobs, and both are very well-constrained in the solar case. While low-Ro simulations tend to produce more cyclic solution, the overall variety of cyclic (and non-cyclic) behaviors uncovered in the afore-cited simulations likely reflects, at least in part, the key influence of small-scale dissipative processes, which are handled quite differently across computational frameworks.
In a recent paper , we have presented a set of global MHD simulations of solar convection produced using the EULAG-MHD framework (see Smolarkiewicz & Charbonneau 2013;Strugarek et al. 2016 and references therein), with the rotation rates varying between 0.55 and 1.1 times solar, and the convective luminosity between 0.2 and 0.61 solar luminosity.
Most simulations in the set exhibit fairly regular largescale magnetic cycles with decadal-type periods, the latter showing systematic variations with rotation rate and luminosity. We revisit the dependency of cycle periods on stellar characteristics with this simulation set. A striking result from our analyses is the decrease of the magnetic cycle period with increasing rotation period. Indeed, cycles in the simulation set are well represented by the relationship P cyc ∝ Ro −1.8 ; here the negative exponent implies a trend opposite to that suggested by Noyes et al. (1984). A key aspect of this new trend is the dependency of cycle period on luminosity; when accounted for in estimating Rossby numbers for the Mt Wilson stellar sample, all solar twins with well-defined periods are back on a single branch in the P cyc -Ro diagram. More importantly perhaps, the Sun is now on the same branch, finally back amongst other solar-type stars.
In this paper we extend the Rossby number range of the Strugarek et al. (2017) simulation set, and show the existence of three well-defined dynamo regimes as a function of the Rossby number: low-Ro with short magnetic cycles appearing close to the surface at low latitude; intermediate-Ro where deep-seated solar-like cycles are observed; and high-Ro where stable wreaths of largescale magnetic field develop (see Fig. 2 for a summary). We further carry out a detailed analysis of the nonlinear driving of large-scale magnetic polarity reversals by magnetically-mediated alterations of the internal differential rotation.
The remainder of the paper is structured as follows. We present the numerical EULAG-MHD model on which this study is based in §2. The dynamo states achieved in the series of 17 simulations analyzed in this paper are summarized in §3. One prototype cyclic solution is then detailed in §4, and the non-linear dynamo acting in these simulations is presented in §5. The trends of the magnetic cycles periods are finally discussed in §6, and we conclude in §7.

TURBULENT CONVECTION ZONE MODEL
Building on the preliminary work of Strugarek et al. (2016), we consider again here spherical fully convective shells with a solar-like aspect ratio (the inner and outer radii are defined by R bot = 0.7 R and R top = R ). We detail in what follows the set of anelastic equations we use to describe the turbulent state of the convective layers of solar-like stars ( §2.1), and the background and ambient state around which the aforementioned perturbed equations are derived ( §2.2). We also explain the modelling technique used to drive turbulent convective motions as well as boundary conditions in §2.3.

Formulation of the anelastic equations
We consider the Lantz-Braginsky-Roberts (LBR, 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 the subscript a) that theoretically may differ from the background state (hereafter denoted with overbars) 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 Ω are where the perturbed quantities are denoted without prime for the sake of simplicity, and D t is the material derivative. We recall that we use standard notation for the basic quantities, i.e. u is the fluid velocity, ρ its density, p its pressure, Θ its potential temperature, and B is the magnetic field. In addition, we use a standard perfect gas equation of state which 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 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). It ensures that on long time-scales, 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 models directly the divergence of the diffusive flux across the domain, and can be shown to be mathematically and physically equivalent to boundary forcing, provided the ambiant 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 towards the prescribed ambient state. The convective luminosity ends up being set by different control parameters under each of these two modelling approaches, hence any direct comparison must be carried out with care.

Background and ambient states
Our numerical setup closely follows the anelastic benchmark of Jones et al. (2011). We consider a spherical shell of aspect ratio β = R top /R bot , and we set d = R top − R bot = R top (1 − β). 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) where n is the polytropic index, ρ i , P i , T i are the density, pressure and temperature at the bottom of the domain, and the constants c 0 and c 1 are given by with α = β + 1 β exp (N ρ /n) + 1 , and where N ρ = ln(ρ i /ρ 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 3 in all the models. The background potential temperature is given bȳ where the standard adiabatic exponent for a perfect gas is γ = c p /c v = 5/3. We chose a polytropic exponent n = 3/2 to naturally ensure an isentropic background state withΘ = P 1/γ i /ρ i , and the equivalent background entropy profile isS = c p ln Θ (we consider here c p = 3.4 × 10 4 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 which we recall is related to the ambient potential temperature profile by Θ a = exp (S a /c p ).

Numerical method
The Eulerian-Lagrangian (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, Eqs. (1-4) are written as a set of Eulerian conservation laws and projected on a spherical coordinate system (Prusa & Smolarkiewicz 2003). EU-LAG solves the evolution equations using MPDATA (multidimensional positive definite advection transport algorithm), 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 stressfree, impermeable boundaries at top and bottom of the domain such that 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 perfect conductor at the bottom and purely radial field at the top, such that The simulations presented in this work are initialized with random, small-amplitude perturbations around the background profiles defined in §2.2. The originality of this numerical method lies in the ILES approach used in the EULAG-MHD code. While this was shown to minimize dissipation at large-scales, it also renders comparison with other simulation results less straightforward. Following Strugarek et al. (2016), we characterized our simulations with effective dissipation coefficients deduced from a spectral analysis of the simulations results. The details of this technical procedure are given in Appendix B. The numerical dissipation of the various MHD fields fit 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 8 m 2 /s, which is what is typically expected for the relatively coarse numerical resolution adopted in this work (Strugarek et al. 2016). The Prandtl number P r = ν eff /κ eff and magnetic Prandtl number P m = ν eff /η eff are generally slightly larger than one, 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) characterize well 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 than that of a Laplacian operator (for more on this see Strugarek et al. 2016).

OVERVIEW OF THE SIMULATION SET
3.1. Dynamo regimes in the simulation set 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 anti-solar differential rotation when the Rossby number reaches values higher than 1 (see panel h in Fig. 1 and . Several dynamo states are achieved in this series with distinctive properties as shown in Fig. 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 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 of model s01, shown in the bottom panels. Notice also that this particular case possesses an anti-solar differential rotation profile. Finally, some models (not shown here, see e.g. Fig. 10) also present stochastic reversals of the deep-seated magnetic field.
The various dynamo states realized in our set of simulations occur in particular ranges of the adimensional  -Three representative cases of the dynamo states achieve in this work. From left to right, we represent the azimuthal magnetic field averaged over longitude Bϕ ϕ at r = 0.75Rtop as a function of latitude and time; the differential rotation profile averaged over hundreds of years; and an instantaneous radial magnetic field in an orthographic projection close to the top boundary of the domain. Panels a to c represent a typical cyclic solution (model O2). Panels d1, e and f represent model 06, and in panel d2 we zoom for this model on the time window labeled by the black box on panel d1. In this panel, we represent Bϕ ϕ at r = 0.98Rtop for which we filtered signals with periodicities larger than 5 years to make the short cycle appear more clearly. Finally, panels g to i represent model sO1, which is non-cyclic and close to the anti-solar differential rotation transition.
fluid Rossby number. It can be defined as where θ,ϕ,t stands for the average over time, latitude and longitude. The Rossby number quantifies the importance of rotation in the dynamics of the system. We define here two Rossby numbers averaged over the where the short cycle is observed in the simulations developing it). Both Rossby numbers are reported in Table 1, and note that the averaging intervals were 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 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πr 2ρ c P u r T θ,ϕ,t at these two locations as well (with T representing the temperature fluctuations with respect to the ambient state).
We summarize in Fig. 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 the function of the Rossby number Ro b . The color of the labels indicates the dynamo state each model achieves, which are 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 06 and S1) possess both a deep-seated cycle and a short cycle at the same time, which was also reported in the millenium 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 Fig. 1), which is consistent here with the Parker-Yoshimura (Parker 1955;Yoshimura 1975) rule of a dynamo wave propagating along the near-cylindrical iso-contours of differential rotation. The deep-seated cycle does not follow such a rule, and originates from the non-linear feedback of the Lorentz force on the differential rotation itself, as will be made clear in §4 and §5.
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 (around 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 Fig. 1) to anti-solar (slow equator, fast higher latitudes, e.g. panel h in Fig. 1). We discuss in more details this transition in §6.4. We now briefly characterize the energetics of our simulation set before detailing the magnetic cycles themselves in subsequent sections.

Energetics in the simulation set
The simulations in our set strongly vary in their kinetic and magnetic energy 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 total kinetic energy in all cases.
The majority of our simulations develop a dynamo in sub-equipartition, with magnetic energy (ME, fifth column) of the order of 10% of the total kinetic energy. The two simulations (O8 and O9) at the higher rotation rates, for which the specifics of the dynamo solution are discussed in §6.3, achieve a super-equipartition regime up to ME∼ 5 KE (magnetostrophic state, see Brun & Browning 2017). In all other cases, 30% to 50% of the total magnetic energy is distributed roughly equally over the large-scale axisymmetric toroidal (TME, sixth column) and poloidal (PME, seventh column) fields. The -Ratio of the differential rotation kinetic energy to the total kinetic energy (third column in Table 2) as a function of the Rossby number Ro b for all the models in our set of simulations The errorbars are given in Table 1. The color bars at the top label the Rossby number regimes where we observe a short cycle (blue), deep-seated cycle (red) and no cycle at all (black). fluctuating magnetic energy spectrum peaks at mid-tosmall scale (typically around the spherical harmonics 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 non-axisymmetric modes.

A REPRESENTATIVE CYCLIC SOLUTION
We now focus on a reference solution displaying features found, to some extent, in all the cyclic simulations presented above. We opted to show simulation O5, which is rotating slightly faster than the Sun and has a solarlike convective luminosity at its top (as shown in Table  1). The magnetic energy amounts to a non-negligible fraction of the kinetic energy (∼10%), the later being dominated by the fluctuating components ('CKE'). Fig. 3 shows time series of the energy contributions  Fig. 3.-Temporal series of energies taken over a magnetic cycle and integrated over the whole domain in model O5. Blue and red are respectively the total kinetic and magnetic energies. Orange is the kinetic energy related to convection, and green is the differential rotation energy. Purple is the magnetic energy in the mean toroidal field, while brown and pink are respectively energies in the mean poloidal field and in the non-axisymmetric components of the magnetic field.
in simulation O5 over a limited time interval spanning two magnetic reversals. Almost all energies are temporally modulated on both long (of the order of 20 years) and short (less than a year) timescales. The total magnetic energy (red) shows a cyclic modulation anti-correlated with the kinetic energy of the differential rotation (green), indicating a coupling between the large-scale flows and fields through the Lorentz force (as further investigated in §5). Note that the energy variation in both of them is of comparable amplitude, which suggests that the energy transfer between the two plays a significant role in the dynamo action. The energy of the turbulent, convective motions (orange) appears to be insensitive to these decadal variations, similarly to previously published EULAG and ASH global MHD simulations presenting magnetic cycles (e.g. Brun et al. 2005;Racine et al. 2011;Augustson et al. 2015). The toroidal (purple) and fluctuating (pink) magnetic energies are found to vary in phase with the total energy. The poloidal magnetic energy (brown) is also essentially varying in phase with the total magnetic energy, albeit small time-lags are observed for some reversal (around t ∼ 334 yr in Fig. 3) which can be traced to hemispheric decorrelations.
An example of a magnetic cycle is shown in Fig. 4 for model O5, along with the mean differential rotation profile (panel c) and its variations during one magnetic reversal (panels d to g). The two top panels show the longitudinally averaged toroidal magnetic fields in latitudetime representations, taken in the lower half of the convection zone (r = 0.75 R top ). Panel a covers the full length of the simulation, with a rectangular box delimiting the smaller time interval displayed in panel b, covering one magnetic cycle. Some solar-like features in the large-scale magnetic field are apparent here, such as regular inversions of polarities, a slight propagation towards the equator during a half-cycle, and magnetic intensities comparable to the ones expected at the base of the solar convection zone. However, the field is quadrupolar, not entirely well synchronized between the hemispheres and located at mid-to high latitudes (some cases exhibit a beating between dipolar and quadrupolar symmetry, see §5.2). This magnetic cycle is the source of the signatures seen in Fig. 3 in both magnetic and differential rotation energies. Panel c of Fig. 4 shows the time-averaged, solarlike differential rotation profile Ω t of model O5 with a fast equator and slow poles, the rotation being defined as Ω = Ω + u ϕ ϕ / (r sin θ). At the bottom right, the departure from the mean profile δΩ = Ω − Ω t is shown at four time steps corresponding to the dashed lines in panel b. These time steps span different phases of the magnetic cycle: from a maximum (t 1 , panel d) to a minimum (t 2 and t 3 ; panels e and f), then to a maximum again (t 4 , panel g). Fast (slow) angular velocities compared to the mean are represented in red (blue). Isocontours of positive (negative) toroidal magnetic field are overlaid in black solid (dashed) lines. During the phases of maximum (t 1 , t 4 ; panels d and g), strong accelerations (red) of the differential rotation are observed (see also §5) on both sides of the magnetic structures (black lines). During the miminum of the cycle (t 2 , panel e), the overall differential rotation is reinforced (acceleration at the equator, deceleration at mid-latitudes) as the simulation is closer to a purely hydrodynamic state (Beaudoin et al. 2018). This phenomenon is at the heart of the dynamo action taking place in these simulations, which we now turn to.

DYNAMO ACTION
The dynamo action sustaining the cycling magnetic field described in §4 requires further investigation. We first examine the polarity inversion mechanism ( §5.1) and then characterize the symmetry properties of the dynamo solution ( §5.2).

Polarity inversion mechanism
In all simulations presenting a deep-seated primary cycle, we observe temporal modulations of the mean differential rotation correlated with the magnetic cycle (see Fig. 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 details this modulation, shown in the top panels of Fig.  5 for model O5 (the left panels represent time-latitude diagrams at r = 0.75R top , right panel display radiuslatitude diagrams at latitude 44 • ). The modulation of the differential rotation δΩ/ Ω 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 towards the equator (left panel) and towards the surface (right panel).
We first demonstrate that these modulations originate from the feedback of the magnetic field on the differential rotation. In the second row of Fig. 5, we show the energy transfer due to the Lorentz force towards the toroidal kinetic energy (last term in Eq. 2, see also Eqs 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 mid-latitude and in the lower part of the convection zone, energy is transferred from magnetic to kinetic energy. The energy transfer furthermore peaks just before the cyclic acceleration in the upper panel (the contours of δΩ have been 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. As a result, the acceleration phases observed in our simulations are indeed of magnetic origin. Note that these modulations differ from the one observed in the millenium EULAG simulation , which were shown to originate indirectly from the 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 their magnetic structures. Here the equatorial migration is very weak, but the direct feedback of the magnetic field on the differential rotation is stronger (see energies in Fig. 3).
However, are these modulations affecting dynamo action? We display the equivalent of the dynamo Omegaeffect in our 3D non-linear simulations in the third row of Fig. 5. The mean shearing of the magnetic field (or Ω-effect), B p · ∇ u ϕ ϕ , is represented in timelatitude and time-radius diagrams, with the same contours of δΩ overlaid in white. We observe that immediately during the acceleration phases of the differential rotation, the mean shear severely drops in amplitude, and ultimately locally changes sign (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., Fig. S8 in Strugarek et al. 2017). The azimuthal field consequently ceases to be sustained (last row in Fig. 5), and as a result the poloidal field (not shown here, see Strugarek et al. 2017) also loses its source and rapidly decreases following the azimuthal field decrease. 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 the opposite polarity and thus a weak azimuthal field of opposite polarity is generated. This new seed azimuthal field is strong enough so 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 compared to the previous configuration. Eventually, the field becomes large enough such that an efficient energy transfer again perturbs locally the differential rotation, which triggers a new reversal. This mechanism repeats over and over, and similarly so in all the simulations of our set which generate a deep-seated magnetic cycle.

Dynamo families
Most of the dynamo states achieved in this series of simulations present a quadrupolar symmetry, with e.g. an azimuthal field symmetric with respect to the equator (see Fig. 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 ('anti-symmetric') and quadrupolar ('symmetric') families (Roberts & Stix 1972;Gubbins & Zhang 1993;McFadden et al. 1991). Projecting a vector on the vectorial spherical harmonics basis (Rieutord 1987;Mathis & Zahn 2005;Strugarek et al. 2013 we can separate its quadrupolar and dipolar components as Based on this decomposition, we further define the parity of the magnetic field (Tobias 1997)

magnetic energies E Q and E
The parity P is equal to +1 when the field is mostly quadrupolar, and −1 when it is mostly dipolar. The parity of two representative simulations (O5 and O2) is shown in Fig. 6 at the top of the domain, along with the total magnetic energy ME (Eq. A4). Case O5 (Fig. 4) exhibits a dynamo solution dominated by the quadrupolar symmetry at all times. When the total magnetic energy peaks, the parity falls down to zero in this case (dipoles and quadrupoles have the same energy). This is somehow opposite to what is observed in the Sun, with a dipolar symmetry dominating during 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 time-scale larger than the magnetic cycle itself. Between t = 590 yr and t = 650 yr, quadrupolar symmetry dominates the field, as in model O5. Nevertheless, model O2 behaves more like the Sun during this period, with the dipolar symmetry being 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 the dipolar symmetry.
Co-existing modulation of magnetic parity and largescale 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 also recently in non-linear three dimensional 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 what follows, under which energy exchange between magnetic and kinetic energy reservoirs is a key ingredient in the evolution of the magnetic cycle.

TRENDS IN THE CYCLE PERIODS AND
TRANSITIONS IN DYNAMO STATES 6.1. Characterization of the cycle periods In order to characterize the cycle period, we need a robust methodology to extract the most significant period out of a potentially multi-periodic signal. One such method was recently successfully applied to similar sim- shown. The first row shows the differential rotation perturbation δΩ (see text) from the time-averaged differential rotation, with white arrow labeling the acceleration phases. The black contours trace the 0.25% variation level and are overlayed in black in the second and fourth rows, and in white in the third row. The second row shows the energy transfer from magnetic energy to azimuthal kinetic energy, red denoting a positive transfer from magnetic to kinetic. The third row shows the mean shear of the poloidal field by the differential rotation (see text), yellow being positive and blue negative. Finally, the fourth row show the mean azimuthal field with blue denoting negative Bϕ and red positive Bϕ.
ulations by Käpylä et al. (2016) based on the so-called empirical-mode decomposition (EMD) method (we refer the reader to Appendix C and Käpylä et al. 2016 for a more in depth discussion about this method). We based our analysis on the libeemd library (Luukko et al. 2015) to make use of the complete variant ensemble empirical mode decomposition (CEEMDAN) method. 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 the CEEMDAN analysis at two representative radii and latitudes, [R b = 0.72R top , θ = 45 • ] and [R t = 0.98R top , θ = 25 • ] on the longitudinally-averaged azimuthal component of the magnetic field. At the top position, the CEEMDAN analysis is applied to the detrended B ϕ ϕ (see panel d2 in Fig. 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 which has the same length, mean and standard deviation as the analyzed simulation time series (robustness increases as R→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, are 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 varying slightly the latitude at the chosen depths (R b , R 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 when Ro b 0.3 (or equivalently Ro t 1), as shown in Fig. 2. Only two models in our set manage to sustain simultaneously 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 §6.4).

Trends of cyclic solutions
The period of the deep-seated magnetic cycle decreases with increasing Rossby number (viz . Tables 1 and 3). Furthermore, the cycle period decreases like a powerlaw with an exponent of −1.6 ± 0.14 when normalized to the rotation period of the star, as shown in Fig. 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 non-linear 3D simulations. All models fit relatively well one single power-law dependency with the Rossby number Ro b . Deep-seated magnetic cycles O models S models R models Fig. 7.-Cycle period normalized to the stellar rotation period as a function of the Rossby number Ro b in logarithmic scale. The different models are color-coded as models where the rotation rate was changed ('O' models, in blue), the entropy contrast in the convective envelope was changed ('S' models, in red), and the number of density scale heights was changed ('R' models, in green). Orthogonal distance regression was performed to find the powerlaw fitted indicated by the gray dashed line, using the error-bars for both the estimated cycle period (Table 3) and Rossby number (Table 1).
Model R1 (green dot in the middle of Fig. 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 to 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 occurred at higher Rossby numbers as the density stratification was increased. In the present case, additional cyclic solutions with larger density contrasts would be required to assess precisely the effect of this parameter on the power-law fit. However, the density contrast in stellar convective envelopes changes accordingly with 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 non-linear character of the dynamo action detailed in §5 was proposed to be at the origin of the negative slope observed in Fig. 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 Eq. 2 with the right-hand side large-scale Lorentz force, one may estimate τ M by where V CZ is the volume of the convection zone, d = R top − R bot is the depth of the convection zone, andρ is the background density averaged over the convection zone. We display in Fig. 8 the cycle period of the cyclic models as a function of τ M estimated with Eq. 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 smaller 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 non-linear dynamo action occurring in our set of simulations. O models S models R models Two additional remarkable features can be noted from Fig. 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 Fig. 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 starts to exhibit an anti-solar 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;Gastine et al. 2014;Fan & Fang 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 still follows the robust Rossby number trend identified in Fig. 7.
Several of our models do not produce a deep-seated magnetic cycle (see Fig. 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 be-comes too long to participate in setting the period of the deep-seated cycle. We now detail what occurs at both of these transitions.

Low Rossby number regime
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 Fig. 2). We display in the left panels of Fig. 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.
As the Rossby number decreases, the total magnetic energy strongly increases in proportion of the differential rotation kinetic energy (see right panel in Fig. 9). In the meanwhile, the azimuthal magnetic field energy strongly decreases relatively to the total magnetic energy (see Table 2). As a result, the fluctuating (i.e. nonaxisymmetric) magnetic energy dominates the total magnetic energy 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 Fig. 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 Fig. 10. Even though we only have 6 models exhibiting such a short cycle, a clear decreasing trend with the local Rossby number is also observed (compare to the deep-seated magnetic cycle in Fig. 7). The powerlaw 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 deep-seated cycle. Similarly to Käpylä et al. (2014), we find that a stronger stratification favors such short cycles near the top of the domain (compare, e.g., models R2 and R3), which here simply reflects 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 one, and we see no hint of a short-cycle disappearance at low Ro t .

Large Rossby number regime
In our set of simulations, the deep-seated magnetic cycles are also lost when the Rossby number exceeds unity. A structural change in the differential rotation occurs when Ro ∼ 1 (e.g. , 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 Fig. 1 for model sO1 for which a strong antisolar differential rotation starts to appear. Model sO1 interestingly devel-  (Table 3) and Rossby number (Table 1). ops 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 magnetic energy is less than 10% of the total kinetic energy (see Table 2). The magnetic energy is dominated by the non-axisymmetric 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 magnetic energy, and shown in Fig. 11. We observe that the dipo-lar strength decreases as a function of the Rossby number Ro b , and report a hint of saturation at 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. 7. CONCLUSIONS In this paper, we offered a detailed presentation of an extended set of global MHD simulations of solar convection and dynamo action. Generated using the EU-LAG code, this simulation set extends the subset originally published in Strugarek et al. (2017). Considering a spherical convective shell with a solar-like aspect ratio, we varied the rotation rate of the shell, its convective luminosity and the number of density scale-heights spanning it. We found that, for adequate parameters well characterized by a range in Rossby number (Eq. 14), the dynamo sustained in the turbulent spherical shell generates a spatially well-organized large-scale magnetic field exhibiting a fairly regular cyclic behavior. When Ro b 0.25, the dynamo loses its regular cycles and irregularly reverses. Stable large-scale magnetic structures resembling the magnetic wreaths found by Brown et al. (2010) are observed when the Rossby number is even further decreased. When Ro b 1, the differential rotation switches from solar-like (fast equator, slow poles) to antisolar-like (slow equator, fast high latitudes), as in . The cyclic dynamo also disappears, and the magnetic field becomes more and more nonaxisymmetric while retaining a predominantly azimuthal component at the base of the convection zone. We recall here that the EULAG code uses an implicit large-eddy simulation (ILES) approach, where all the dissipation is handled by the numerical scheme itself. This has the important consequence that the large spatial scales in this set of simulations are essentially inviscid, and the dissipation at mid-to small scales is approximated by standard Laplacian-like dissipative operators, with dissipation coefficients reported in Table 4 (see also Strugarek et al. 2016).
The cyclic dynamo states found in our 3D, non-linear turbulent simulations are in a fundamentally non-linear regime where the feedback of the Lorentz force on the differential rotation plays an essential role in setting the cycle period. The feedback induces variations of only ∼ 1% of the total differential rotation, but is strong enough to locally affect the latitudinal gradient of the differential rotation and trigger the reversal of the global magnetic field. This interpretation, supported by the detailed analysis of dynamo action in §5, is further confirmed by the correlation between the cycle period and the characteristic timescale associated with the Maxwell torques ( §6.2).
If the dynamo action observed in this series of simulations does occur in the convective envelope of solar-like stars, several important conclusions for stellar cycles can be drawn from this study.
First, the cycle period of such stars is set by the nonlinear feedback of Maxwell torque deep inside the convection zone. As a result, the surface effects such as the Babcock-Leighton mechanism, while certainly possible and in fact observed at the solar surface, are not essential to the dynamo nor to its cyclic character.
Second, close to the low-Rossby transition, two types of cycles can easily be mixed when comparing observational results. These cycles, which are non-linearly coupled in turbulent convective envelopes, are not necessarily expected to follow the same trends and dependencies, nor even to originate from the exact same dynamo mechanism. Indeed, in our set of simulations we find two different trends for the magnetic cycle periods, and we also find different propagation properties of the cyclic fields (see panels d1 and d2 in Fig. 1). Particular care is thus in order when determining trends (or branches) from observed cycles on a large sample of stars.
Third, at high Rossby numbers, our results suggest a change in the dynamo state promoting smaller scale, less axisymmetric fields. This transition at Ro ∼ 1 is interesting in several aspects. First it has been shown to likely be the transition of the state of differential rotation from solar to anti-solar  and references therein). Second, it has been argued by (van Saders et al. 2016;Metcalfe & van Saders 2017) that such a change could impact the type of dynamo occurring (by modifying, e.g. the omega-effect) and hence explain a possible break of stellar rotation-age relationship (so-called gyrochronology, see Barnes 2003). In our set of simulations, at the transition Ro ∼ 1 the large-scale field does not disappear and the dipolar fraction remains within 10% of the total magnetic energy. Hence the change of the properties of the dynamo field we observe at Ro ∼ 1 does not seem to be large enough to explain by itself a strong decrease of the torque exerted by the magnetized wind of the star. We intend to explore further dynamo action near the Ro ∼ 1 transition to quantitatively assess how the dynamo and magnetic activity is modified and how this affect the stellar wind properties.
In the regime where a cyclic dynamo is found, Strugarek et al. (2017) showed that the solar cycle period fits very well the Rossby trend we found. However, some solar-type stars of the analyzed sample diverged significantly from our single-branch theoretical trend. These differences could be ascribed to varying aspect ratios and absolute thickness of the convection zone, thus spanning a different range of density scale heights. The differential rotation and Rossby number achieved at the base of such convection zones are thus expected to have different dependency to the stellar parameters (rotation, luminosity) than in our study where we kept the aspect ratio constant. Other stellar internal structures (as in , eventually including the coupling with an underlying stable layer (see Beaudoin et al. 2018 for a first step on this aspect), are now needed to assess the robustness of our scaling-law for solar-like stars of different spectral types and metallicities.
Finally, we want to stress again that the simulations presented in this study operate in an overall parameter regime very far from that expected to characterize physical conditions in solar and stellar interiors. Furthermore, a systematic comparison of the different ways of forcing convection in stellar envelopes and how surface layers and stratification may influence the description of stellar turbulent transport of heat, angular momentum and magnetic field is in order in the community, and some attempts in these directions have been initiated (see e.g. Lord et al. 2014;Strugarek et al. 2016;Cossette & Rast 2016;Cossette et al. 2017;Nelson et al. 2018). The dynamo mechanism identified in this series of simulations nevertheless gives a novel robust mechanism explaining the origin of magnetic cycles in turbulent convective envelopes. These simulations are able to reproduce the solar cycle period, in spite of our relatively coarse spatial discretization grid. We hope in a near future that better resolved simulations, and ultimately comparison with results obtained with other numerical techniques, will help enhance the realism of the dynamo mechanism unveiled in this study, and further validate it applicability to stellar dynamo and cycles.
We thank S. Tobias and J.-D. do Nascimento Jr.
for useful discussions on magnetic cycles and dynamos, and an anonymous referee for suggestions that improved the presentation of our results. The authors acknowledge support from the Natural Sciences and Engineering Research Council of Canada. A. Strugarek ac-knowledges support from CNES in preparation of the Solar Orbiter mission. This work was also supported by the INSU/PNST. We acknowledge access to supercomputers through GENCI (project 1623) and Calcul Québec, a member of the Compute Canada consortium.

A. DEFINITIONS OF THE GLOBAL ENERGIES
We quantify the energies of each simulation by defining the total kinetic energy (KE), along with its differential rotation (DRKE) and convective (CKE) components, as follows where t stands for the temporal average, ϕ for the azimuthal average,X = X − X ϕ and dV = r 2 sin θdrdθdϕ.
We furthermore define the total magnetic energy (ME) and its mean toroidal (TME), mean poloidal (PME) and fluctuating (FME) components as All these energies are listed in Table 2 for the whole simulation set.
B. DISSIPATIVE PROPERTIES OF THE SIMULATIONS We characterize the dissipative properties of the models following the methodology detailed in Strugarek et al. (2016). The basic principle is as follows: the set of anelastic equations (1-4) is projected on the spherical harmonics basis. By projecting each spectral equation (for example, we project Equation 2 on the vectorial spherical harmonics and take the dot-product with u) and integrating them over concentric spheres, we obtain evolution equations for the kinetic energy, potential temperature, and magnetic energy spectra (see Strugarek et al. 2016). The evolution equation for the magnetic energy spectrum can be written, using the same notation as in Strugarek et al. 2016, where T L is the energy transfer from the kinetic energy reservoir (see Strugarek et al. 2013). Formally, EULAG-MHD solves the induction equation with no ohmic dissipation. As a result, the residualĖ M L − T L is an estimate of the level of numerical dissipation added by the MPDATA algorithm, which is generically scale and time-dependent. We match here this residual (averaged over time) to an explicit ohmic dissipation with an effective coefficient η eff , namely The procedure was illustrated at length in Strugarek et al. (2016), we do not repeat it here and simply present in Fig.  12 the resulting η eff as a function of depth and spherical harmonic degree L for model O5. We see that the effective ohmic dissipation is of the order of 1-2×10 8 m 2 /s in the bulk of the convective envelope for modes L > 25. The grey area corresponds to a region where the effective dissipation cannot be matched to a classical laplacian operator and is smaller than the dissipation at larger scales. The spatial distribution with slightly higher η eff in the upper part of the convection zone is a direct consequence of the scale distribution of the magnetic structures as a function of depth. We note one more time that the largest scales in the domain (typically L < 10) are almost not dissipated compared to the smallest scales, which are in turn the scales well characterized by the dissipation coefficients derived through this methodology. Finally, the overall procedure is repeated for all models and also for the effective viscosity ν eff and heat dissipation coefficient κ eff , which are given in Table 4.  C. EMPIRICAL MODE DECOMPOSITION ALGORITHM In order to robustly extract the periodicity of the dynamo cycle, we make use in this work the Complete variant of the Ensemble Empirical Mode Decomposition with Adaptative Noise (CEEMDAN algorithm). Empirical Mode Decomposition (EMD) is a method to decompose an input signal into a series of Intrinsic Mode Functions (IMFs), which are oscillatory modes that are not required to be simple sinusoidal functions but still retain meaningful local frequencies (for an in-depth discussion of the EMD method, see Luukko et al. 2015;Käpylä et al. 2016 and references therein). The CEEMDAN variant was further developed by Torres et al. (2011), which achieves completeness of the decomposition and improves the robustness of the method for noisy signals. The CEEMDAN method has been made available in the open-source libeemd library (Luukko et al. 2015), which was used in this work. We show in Fig. 13 one example of a CEEMDAN run applied to B ϕ ϕ at r = 0.72 R and latitude 57 • for model O5. The original signal is shown at the top panel in blue, in physical units. Given the length, mean and standard deviation of the signal, the CEEMDAN method decomposes it automatically into 11 IMFs that are shown in the other panels. The residual of the decomposition is shown in green in the bottom panel. For each IMF, one can compute its relative energy E to assess how much its contributes to the total signal. In this particular example, IMF 10 (in red) completely dominates the decomposition with a relative energy of 76%. It is also important to assess the significance of the decomposition, characterized here by the robustness parameter R indicated in each panel. To have an estimate of the robustness of each IMF, we build a white noise signal of the same length as the original signal with the same mean and standard deviation. We perform the CEEMDAN on this white noise signal, and calculate the robustness of each IMF of the original signal by its relative energy divided by the relative energy of its white noise IMF counterpart. The parameter R is then the normalized so that the sum of the robustnesses is 1, and we find that IMF 10 is also very robust (87%) giving us confidence that it captures accurately the main periodicity of the original signal. Finally, a mean period (P) with its error-bars can be easily computed for each IMF, giving us a robust estimate of the main periodicity of the dynamo at that depth and latitude for model O5.