Powering Stellar Magnetism: Energy Transfers in Cyclic Dynamos of Sun-like Stars

We use the ASH code to model the convective dynamo of solar-type stars. Based on a series of 15 3-D MHD simulations spanning 4 bins in rotation and mass, we show what mechanisms are at work in these stellar dynamos with and without magnetic cycles and how global stellar parameters affect the outcome. We also derive scaling laws for the differential rotation and magnetic field based on these simulations. We find a weaker trend between differential rotation and stellar rotation rate, ($\Delta \Omega \sim (|\Omega|/\Omega_{\odot})^{0.46}$) in the MHD solutions than in their HD counterpart $(|\Omega|/\Omega_{\odot})^{0.66})$, yielding a better agreement with the observational trends based on power laws. We find that for a fluid Rossby number between $0.15 \lesssim Ro_f \lesssim 0.65$ the solutions possess long magnetic cycle, if $Ro_f \lesssim 0.42$ a short cycle and if $Ro_f \gtrsim 1$ (anti-solar-like differential rotation) a statistically steady state. We show that short-cycle dynamos follow the classical Parker-Yoshimura rule whereas the long-cycle period ones do not. We further demonstrate that the Rossby number dependency of the large-scale surface magnetic field in the simulation ($B_{L,surf} \sim Ro_{f}^{-1.26}$) agrees better with observations ($B_{V} \sim Ro_{s}^{-1.4 \pm 0.1}$) and differs from dynamo scaling based on the global magnetic energy ($B_{bulk} \sim Ro_{f}^{-0.5}$). We also show that up to few percents of the stellar luminosity can be channelled into the star's magnetism, hence providing a large energy reservoir for possible surface eruptive events.


INTRODUCTION
Sun-like stars go through various magnetic activity phases in their lives. From young very active TTauri stars rotating much faster than our Sun to old stars that are less active, it is key to understand how convection, rotation, turbulence, magnetism and surface activity evolve and feedback on one another over secular time. Of particular interest is the generation of magnetic field via dynamo action, because it is both as the source of key temporal variabilities like the Schwabe 11yr or Gleissberg 90yr magnetic cycles in the Sun and at the heart of a complex feedback loop between stel-lar magnetism and rotation via wind braking and the loss of mass and angular momentum by the star Brun & Browning 2017;Vidotto 2021). It is also key in providing the free energy reservoir needed to power eruptive events such as flares or CME's (Shibata et al. 2013;Aschwanden et al. 2015;Maehara et al. 2017). In this work, we seek to assess how solar-like stars with different masses and rotation rates can power their magnetism by means of dynamo action in their convective envelopes.
Various activity indicators have been derived observationally over the last 50 years using for instance photometric and spectroscopic variability (Baliunas et al. 1995;Oláh et al. 2009;Egeland 2017;Boro Saikia et al. 2018), and more recently through asteroseismology (García & Ballot 2019) to connect the spectral class 2 A.S. Brun et al. and age of a star to its dynamical properties and activity level. Turning specifically to solar type stars, spectropolarimetric studies have revealed several interesting properties (Marsden et al. 2014). In Vidotto et al. (2014), it was shown that the large-scale magnetic field is following a scaling law with the stellar Rossby number B V ∝ Ro −1.38±0.14 s for stars with Ro s > 0.1 (here the stellar Rossby number is defined as the ratio between the rotation period and the convective turnover time). More recently, See et al. (2019a) have revisited this trend and found a similar result with B V ∝ Ro −1.40±0.10 s . It was also proposed in Petit et al. (2008) and later by See et al. (2015) that toroidal magnetic field dominates over poloidal field for fast rotators. It was further shown that no significant collapse of the large-scale field with respect to higher multipole moments was observed as the star evolved and are found less active (Vidotto et al. 2016). Recent work by Lehtinen et al. (2021) including more evolved stars seems to help constrain better the rotation-activity relationship, confirming that considering the Rossby number is better than the rotation period alone. In Karoff et al. (2018) the possibility that larger metallicity increases the activity level of solar analogs was also proposed.
Moreover, long observational studies based on Ca II H&K chromospheric observations have shown that magnetic activity of solar-like stars (Wilson 1978;Saar 1990;Plachinda & Tarasova 1999;Hall 2008;Hall et al. 2007) can be found to be either irregular with no obvious cyclic activity or to possess activity cycles with short magnetic cycle periods (Metcalfe et al. 2010;Jeffers et al. 2018) or long (decadal) ones (Noyes et al. 1984;Baliunas et al. 1995) as in the Sun. Such studies have further indicated the existence of a relation between Rossby number and magnetic cycle periods, its exact nature being still debated given the relatively small numbers of truly confirmed cyclic magnetic stars (do Nascimento et al. 2014;Egeland 2017).
A puzzling property regarding stellar magnetic cycles has been the existence (or not) of active and inactive branches of stellar activity, as proposed by Saar & Brandenburg (1999) and Böhm-Vitense (2007). Recently it has been argued that stars may be transiting from one state to the other as they evolve and that such distinct activity branches do not exist. Instead, activity would be decreasing while rotation would be almost unchanged beyond a certain age or stellar internal dynamo state. A key quantity to characterize this activity state transition is again the Rossby number. Metcalfe & van Saders (2017) proposed that once their Rossby number becomes large, stars stop braking through their stellar wind, hence departing from the classical Skumanich law Ω(t) ∝ t −0.5 and gyrochronology trend (Skumanich 1972; Barnes 2003Barnes , 2007. This is still highly debated in the community as some observers find stars older than the Sun still following Skumanich's law (Meibom et al. 2015; Barnes et al. 2016;Lorenzo-Oliveira et al. 2018do Nascimento et al. 2020) while others do not Metcalfe & Egeland 2019;Hall et al. 2021). The disagreement could also be due to the observation techniques (photometric versus chromospheric studies for instance) and observational data set (Kepler data vs long-term monitoring of individual stars) used, since each have rotation rates and ages determination that differ sometime significantly (Lorenzo-Oliveira et al. 2016;do Nascimento et al. 2020). Another alternative would be that stars temporarily stop spinning down before starting again (Spada & Lanzafame 2020;Curtis et al. 2020) or that the coronal temperature drops, yielding a smaller mass loss for older stars (Ó Fionnagáin & Vidotto 2018). Thus, understanding what happens from a theoretical point of view to stellar dynamo and magnetic field geometry for large Rossby numbers is crucial in helping to interpret the most recent observations. This is one of the goals of this study. Given the close link between surface activity and stellar magnetism, a key aspect to characterize is the amount of magnetic energy made available in a given solar-like star by dynamo action. We know that flare intensity is linked to the magnetic energy made available to the magnetic structures. It is thus crucial to better characterize energy transfers in solar-type star dynamos for a wide range of Rossby numbers.
Characterizing the differential rotation realised at the base and in the convective envelope of solar-type stars is central to the understanding of their magnetic field generation, activity level and rotation, as it is directly linked to the Ω-effect (e.g. stretching of the poloidal magnetic field lines by large scale shear). Hence, the role of the differential rotation (DR) in driving the star's magnetic activity level and field properties should be clarified (Donahue et al. 1996). Doppler imaging (Donati & Collier Cameron 1997; Barnes et al. 2005), asteroseismology (Gizon & Solanki 2004;Reinhold et al. 2013;García et al. 2014), classical spot models (Lanza et al. 2014) and short-term Fourier-transform (Vida et al. 2014) are methods to infer differential rotation. The combination of all these observations on stellar rotation and magnetism helps constrain the trends linking rotation with stellar differential rotation and magnetic activity. Various analyses of stellar differential rotation revealed different dependencies between DR and star's rotation (∆Ω ∝ Ω n ), with n varying between 0.15 and 0.7 (Barnes et al. 2005;Reiners 2006;Reinhold et al. 2013).

3
There is no clear consensus in the community for now, some authors are even advocating that such laws should be derived per spectral stellar classes and that the confusion comes from mixing together F and K stars (Balona & Abedigamba 2016). Saar (2011) ;Brandenburg & Giampapa (2018) also propose that the dependency of the differential rotation with the rotation rate may not be monotonic, with a break near Rossby equals unity. By contrast, a more systematic and stronger dependency is observed with the star's temperature (∆Ω ∝ T 8.92 eff Barnes et al. 2005;Reinhold et al. 2013 and ∆Ω ∝ T 8.6 eff Collier Cameron 2007). Hence, we expect large-scale shear to vary both in amplitude and profile (as a function of latitude and depth), as the global stellar parameters change. Some recent studies have confirmed this is happening in solar-type stars by inverting seismically their profile (Benomar et al. 2018), pointing to a possible anti-solar differential rotation state (e.g. slow equator/fast poles) which was possibly already guessed in F stars (Reiners 2007) and advocated to exist in numerical simulations (Matt et al. 2011;Gastine et al. 2014;Brun et al. 2015, see below).
Considering the large number of global stellar parameters probed by these different observational studies, it is expected that the excitation of various types of convective dynamos may occur (Weiss 1994;Tobias 1998;Brun & Browning 2017;Brandenburg & Giampapa 2018;Charbonneau 2020). In order to quantify the influence of key parameters such as rotation and mass in characterizing the dynamo and magnetic level achieved in solar-like stars and given the intrinsic nonlinear mechanisms at work in stellar dynamos, multi-D numerical simulations have been developed over the years in an attempt to provide more quantitative answers.
Some studies have used the 2.5D mean field dynamo approach to do so, extending solar mean field dynamo models to other stellar spectral types (Chabrier & Küker 2006;Jouve et al. 2010;Küker et al. 2011;Kitchatinov et al. 2018, and references therein). While these studies are very helpful, most of them lack the full nonlinearity and genuine parametric dependence of 3D magnetohydrodynamic (MHD) simulations. Recent developments by Pipin (2021) are starting to overcome these limits and have extended the work of Rempel (2006) on the Sun to solar-type stars with various rotation rates. Nevertheless, with the arrival of more powerful supercomputers, other authors have used instead global 3D MHD simulations to model differential rotation and stellar magnetism in the convection zone of solar-like stars (Glatzmaier & Gilman 1982;Miesch et al. 2000;Brun et al. 2004;Miesch et al. 2006;Brown et al. 2008;Ghizaru et al. 2010;Brown et al. 2010;Brun et al. 2011;Käpylä et al. 2011Käpylä et al. , 2014Gastine et al. 2014;Augustson et al. 2015;Karak et al. 2015). These studies pointed out the large magnetic temporal variability and the critical effect of stellar rotation and mass on magnetic field generation through dynamo mechanism, leading in some parameter regimes to configurations with cyclic activity (Gilman & Miller 1981;Gilman 1983;Glatzmaier 1985a;Racine et al. 2011;Brown et al. 2011;Nelson et al. 2013;Käpylä et al. 2013;Augustson et al. 2013Augustson et al. , 2015Beaudoin et al. 2016;Guerrero et al. 2016;Strugarek et al. 2017Strugarek et al. , 2018Warnecke 2018;Viviani et al. 2018Viviani et al. , 2019Guerrero et al. 2019;Matilsky & Toomre 2020). Several studies pointed out the positive effect of a stable region underneath the convection zone (Parker 1993) on the efficient storage of intense toroidal field and the lengthening of the stellar dynamo cycle period (Glatzmaier 1985b;Browning et al. 2006;Lawson et al. 2015;Beaudoin et al. 2016;Guerrero et al. 2016Guerrero et al. , 2019Käpylä et al. 2019;Bice & Toomre 2020). Over the last decade significant progress has been made in successfully simulating large-scale mean flows and stellar activity cycle using different numerical codes and methods (Jones et al. 2011). This is quite reassuring that a global consensus is growing on the nature of solar-like star dynamos. It is common knowledge that there are still key transitions in Rossby number (at low and high values of this parameter) that need to be understood further, as well as what is the exact type of convective dynamos realized in solar-like stars as their global parameters are varied. This study continues this effort by doing an even broader systematic parametric study of solar-like star dynamos coupled to a stably stratified layer below than what have been published so far. It extends the work published in Varela et al. (2016) and  with the MHD anelastic spherical harmonic code (ASH) . In particular, we wish to better characterize energy transfers and how much of a star's energy (luminosity) is converted into magnetic energy by nonlinear global convective dynamos over a wide range of Rossby numbers, generalizing to solar-like stars the work by Starr & Gilman (1966) and Rempel (2006).
In the following sections, we analyze how differential rotation and magnetism feedback on one another (Brun 2004;Fan & Fang 2014) as well as how kinetic and magnetic energies flow within a stellar magnetized rotating convective envelope, using 15 convective dynamo MHD simulations for model stars with different masses and rotation rates (hence Rossby numbers) in order to achieve this goal. In section §2 we present the equations and model setup. In section §3 we make a quick overview of one of the dynamo solutions emphasizing the main properties of a cyclic solution. In §4 we discuss the 4 A.S. Brun et al. various DR profiles obtained in our parametric studies, expanding Varela et al. (2016) to include 15 models. We discuss angular momentum and various scaling laws of the differential rotation contrast ∆Ω. In §5 we analyze our dynamo solutions for various key properties as a function of the Rossby number, such as their activity level, the amount of magnetic flux generated by the dynamo, the existence or not of an activity cycle and torsional oscillations, how the cycle period for cyclic solutions changes, what is the relative contribution of dipolar and quadrupolar magnetic fields in the overall dynamo generated magnetic field, and interpret our simulations in terms of mean field α-Ω classification. We further expand our data set with the 17 simulations published previously Strugarek et al. (2017Strugarek et al. ( , 2018 with the Eulag-MHD code (Smolarkiewicz & Charbonneau 2013), in order to improve the statistics. In §6 we perform an extensive study of energy transfer between various reservoirs in stellar dynamos, assessing how much magnetic energy is accessible to stars like our Sun to power eruptive events. We compute all MHD transfers between kinetic and energy reservoirs for the large-scale flows and magnetic fields. In §7 we reflect on our findings in an astrophysical context, comparing our results with recent observational results and then conclude.

NUMERICAL SETUP
In this section we present the main features of the ASH code, describing the boundary and initial conditions of the numerical models and our choice of global parameters.

Set of equations solved
We perform 3D MHD simulations of convective dynamo action coupled to a stable radiative interior where the anelastic MHD equations are solved for the motions of a conductive plasma in a rotating sphere (Jones et al. 2011). The anelastic approximation captures the effects of density stratification without having to resolve sound waves, which would severely limit the time step (Brown et al. 2012). In the MHD context, the anelastic approximation filters out fast magneto-acoustic waves but retains Alfvén waves.
The code ASH uses a pseudo-spectral method (Clune et al. 1999). The velocity (v), magnetic (B), and thermodynamic variables (entropy S, pressure P ) are expanded in spherical harmonics Y m (θ, φ) for their horizontal structure and in Chebyshev polynomials T n (r) for their radial structure ). The density (ρ), entropy, pressure and temperature (T ) are linearized about the spherically symmetric background values, denoted by the symbol (ˆ). The equations solved by ASH are : with the velocity field v = v rêr + v θêθ + v ϕêϕ , the magnetic field B = B rêr + B θêθ + B ϕêϕ , the angular velocity in the rotation frame Ω * = Ω * êz ,ê z the unit vector along the rotation axis, g the magnitude of the gravitational acceleration. A volumetric heating term ρ is also taken into account to approximate generation of energy by nuclear reactions in the stellar core. The nuclear reactions are modeled very simply by assuming that = 0T nc . By enforcing that the integrated luminosity of the star matches its known surface value, we can determine 0 and n c as listed in Table 7 of . Note that only the low-mass star series of models (e.g. 0.5 and 0.7 M ) require that heating source term, since their computational domain includes a portion of the nuclear energy generation core.
The diffusion tensor D and the dissipative term Φ d are defined as: with e ij the stress tensor and J = c/4π∇×B the current density. The energy flux q is the sum of a radiation flux and of a turbulent entropy diffusion flux: with ν, κ and η the effective eddy diffusivities of the momentum, heat and magnetic field transport, κ r the atomic radiation diffusion coefficient, κ 0 the effective thermal diffusivity acting only on the spherically symmetric (l = 0) entropy gradient and c p the specific heat at constant pressure. Due to limitations in computing resources, current numerical simulations cannot capture all scales of solar convective motions and magnetic fields from global to atomic dissipation scales. The simulations described in this study resolve nonlinear interactions among a large range of scales but motions and magnetic field still exist Dynamo action in G and K stars 5 in solar-like stars on scales smaller than our grid resolution. Hence, our models should be considered as large-eddy simulations (LES) with parameterization to account for subgrid-scale (SGS) motions. The effective eddy diffusivities ν, κ, and η represent momentum, heat, and magnetic field transport by motions which are not resolved by the simulation. They are allowed to vary with radius but are independent of latitude, longitude, and time for a given simulation. In the simulations reported here, ν, κ, and η have the following profile: and with ν top in cm 2 s −1 and r t and σ t in cm, as provided in Table 7 of , α is -0.5 for all cases. All models assumed a Prandtl number P r = ν/κ of 0.25, so that κ can be directly obtained from the amplitude and profile of ν. The magnetic Prandtl number P m = ν/η is equal to 1 or 2 depending on the case considered (see Table 3), so that η can as well be deduced from ν. These tapered profiles are chosen in order to take into account the much smaller sub-grid scale transport expected in the stably stratified radiative interior. A representative profile is shown in Figure 1. Their amplitudes are adapted for each rotation rates and stellar masses considered in order to achieve the best turbulent convective dynamo states while retaining a reasonable numerical resolution and computing effort (still, each model has used of the order of 8 to 10 million cpu hours spread over several years). The diffusivity κ 0 is set such as to have the unresolved eddy flux carrying the solar flux outward at the top of the domain (see Figure 2). It drops off exponentially with depth in order to avoid a large inward heat flux in the stable zone (see Miesch et al. 2000). Of course there is some arbitrariness in choosing the exact shape and amplitude of our diffusivity profiles, and we do our best to limit their influence on the results reported here.
The mass flux and magnetic vector fields are maintained divergenceless by a streamfunction formalism : A perfect ideal gas equation is used for the mean state and the fluctuations are linearized: Typical radial profile of kinematic viscosity ν used in this study, here for case M11R3m. Profiles of κ and η are the same, but their amplitude depends respectively for each cases on the chosen Prandtl and magnetic Prandtl numbers (see Table 3). ρ/ρ = P/P − T /T = P/γP − S/c p with γ = 5/3 the adiabatic exponent. The anelastic MHD system of equations requires 12 boundary conditions. We use an impenetrable and stress-free boundary conditions at the top and bottom of the domain, i.e.: Magnetic boundary conditions are perfectly conducting at the lower radial boundary and the magnetic field matches to a potential field in the upper boundary: Bϕ r | r bot = 0 and B r | rtop = ∇Ψ ⇒ ∆Ψ = 0, with r top , r bot respectively the radius of the top and bottom of the numerical domain and r bcz that of the base of the convective layer (cf . Table 1). Finally, we maintain the entropy flux at the top and bottom. Keeping the values of dŜ/dr| rtop,r bot fixed at all times in the simulations further implies that the fluctuating dS/dr is set to zero at both boundary conditions.

Model structure and initialization
The simulation is focused on the bulk convection zone, avoiding regions too close to the stellar surface. We include a stably stratified layer below the convective envelope, hence providing a realistic bottom boundary condition for the fields and flow that are allowed to be pumped down and to penetrate into the radiative interior. The code uses a realistic background stratification for the profiles of entropy (Ŝ), density (ρ ), temperature (T ) derived from a one-dimensional solar structure model CESAM (Morel 1997;. Our starting point are the G and K star rotating convective 3D models published in  (see also Matt et al. 2011;Brun et al. 2015 and Table 1).
The MHD models are initialized from their equivalent progenitor hydrodynamical models in which a small magnetic field perturbation is introduced in the convective envelope (many orders smaller than the final magnetic field observed in the simulation). In that hydrodynamic study we published 15 simulations covering 4 mass bins and 4 rotation rates. We have models for stellar masses 0.5, 0.7, 0.9, 1.1 M and rotation rates ranging from 1/8 to 5 times the solar rotation rates. In keeping with the naming nomenclature of , we name our model such as to indicate the mass of the star and its rotation rate. The models are named MAxrm, where 'A' is the mass of the star and 'r' the rotation rate of the star (in solar rotation rate). The index 'x' indicates slow/anti-solar (x = s) and prograde (x = R) differential rotation models (except model M11R1m that is also anti-solar) and m stands for magnetism, to distinguish between the hydrodynamic progenitor published in 2017 and their MHD dynamo counterparts considered in this study.
The models have a numerical resolution of (N r , N θ , N φ ) 769 x 256 x 512 except for few cases in the M09m & M11m series rotating at Ω * = 3 or 5Ω where N θ is 512 and N φ is 1024.
In Fig. 2 we show in (a) an example of the radial dependence of the entropy gradient, and in (b) the temporal and azimuthally averaged radial energy fluxes bal-ance as luminosities for the model with 1.1 solar mass and 3 times the rotation rate of the Sun. We note the sharp increase of the stratification at the base of the convective envelope that is coherent with the stiff radiative interior found in main sequence solar-like stars. Such a realistic interface, as opposed to an impenetrable wall, allows the convective motions to overshoot beyond the radius where the entropy gradient changes sign. By doing so they generate internal gravity waves and pump magnetic field. Since we are in this study mostly interested with the magnetic state of our simulation, we refer the reader to the following multidimensional studies of internal gravity waves generation in solar-like stars (Rogers & Glatzmaier 2006;Brun et al. 2011;Alvan et al. 2014Alvan et al. , 2015. Turning to the radial flux balance, we note that the enthalpy flux (dash-triple-dotted line) dominates energy transport in most of the convective envelope. The diffusive fluxes (radiative in the bottom half of the computational domain (long-dash) and unresolved near the top (dotted line)) carry the stellar luminosity at each end of the domain. We note an inward kinetic energy flux (dash dot line) reaching about 10% of the star's luminosity, as is common to find in stratified convection simulations. The Poynting and viscous fluxes account for less than 1% of the radial energy balance. Finally, we note the negative enthalpy flux near the base of the convective envelope, that is compensated by a local increase of the radiative flux, such as to reach a satisfactory radial energy balance and thermal equilibrium. All the listed values were computed with the CESAM stellar evolution code (Morel 1997). We adopt M = 1.989 × 10 33 g, R = 6.9599 × 10 10 cm, and L = 3.846 × 10 33 erg · s −1 . The density ratios ∆czρ and ∆ fρ are evaluated by forming the ratio between the value of the density respectively at the base of the convection and the top of the domain and at the bottom and the top of the domain.  Table 2. Models dimensional characteristics, averaged over a small interval of 0.01R at the middle of the convective envelopes (unless stated otherwise). Characteristic velocities, differential rotation, magnetic fields, and time scales are listed. The differential rotation is taken between latitude 60 • and the equator at the surface of the models (see §4.1). Likewise, the total magnetic flux is computed at the surface of the models, and averaged over at least one magnetic cycle for the cyclic cases (see §5.4).

G and K stars parametric study
As indicated above, we initiate each of the 15 dynamo simulations from mature, relaxed hydrodynamics convective states and introduce a seed magnetic field in the convective envelope only. These hydrodynamical progenitors have been run long enough to reach a statistically stationary state in the convection zone and a well established rotation profile. They possess a gen-uinely established tachocline, defined as the transition between differential rotation in the outer convective envelope to solid-body rotation in their stable radiative interior, leading to regions with strong shear (Spiegel & Zahn 1992). The tachocline plays an important role in the dynamo process of magnetic field generation in solarlike stars as reported in simulations performed by several authors (Glatzmaier 1985b;Browning et al. 2006 Table 3. Model characteristics non-dimensional numbers, averaged over a small interval of 0.01R at the middle of the convective envelopes. Re =ṽD/ν is the Reynolds number. The Prandtl number P r = ν/κ = 1/4 in all cases. Pm = ν/η is the magnetic Prandtl number, Rm = RePm is the magnetic Reynolds number, and P e = ReP r is the Péclet number. Ra = (−∂ρ /∂S)∆SgD 3 /ρ νκ is the Rayleigh number, and Ra /Rac is the modified Rayleigh number as computed by Takehiro et al. (2020). T a = 4Ω 2 * D 4 /ν 2 is the Taylor number. We also list three Rossby numbers: the fluid Rossby number Ro f =ω/2Ω , the convective Rossby number Roc = Ra/T aP r, and the stellar Rossby number Ros = Prot/τ CS c . The latter is useful for comparison from observationnally-derived Rossby numbers. For Ros we have therefore considered the empirical convective turnover time derived by Cranmer & Saar (2011) + 0.002 days. We note that it correlates well with our fluid Rossby number and find Ros 2.26Ro f . The Ekman number is defined as Ek = ν/(Ω D 2 ), and the Elsasser number as Λ =B 2 /8πρ DṽΩ . Racine et al. 2011;Masada et al. 2013;Lawson et al. 2015;Guerrero et al. 2016).
The main parameters of the models are listed in Tables 2 and 3. The density scale heights between the top and the base of the convection zone and between the top and the bottom of the model are defined as Nρ bcz = ln(ρ out /ρ bcz ) and Nρ tot = ln(ρ out /ρ in ). For the M05 model Nρ bcz = 3.25 and Nρ tot = 4.70, M07 model Nρ bcz = 3.48 and Nρ tot = 5.78, M09 model Nρ bcz = 3.31 and Nρ tot = 5.99, M11 model Nρ bcz = 3.28 and Nρ tot = 5.60. The convective flows at the middle of the convective envelope vary from 5 m/s up to about 300 m/s in our sample, and the convective turnover time from 7 (case M11R1m) to 222 (case M05R3m) days. The surface differential rotation between the equator and latitude 60 • varies from -102 to +278 nHz in our sample, and we will study its maintenance in detail in §4.1. In the middle of the convection zone, the rootmean-squared magnetic field typically varies from 1 G to 70 G in our models. It is found to be maximum close to the bottom of the convective envelope, where the large-scale shear efficiently powers the dynamo, and the magnetic field can be stored in the tachocline close to the convective-radiative boundary.
Our sample of simulations were designed to operate in a relatively homogeneous turbulent Reynolds number regime, as seen in the first column of Table 3. The supercriticality degree can be characterized by the Rayleigh number achieved in our models, compare to a critical Rayleigh number for the onset of convection. Such a modified Rayleigh number was proposed by Takehiro et al. (2020) and is listed in column 6 of Table 3. All our models exhibit a Rayleigh number at least five time larger than the critical Rayleigh number. We have run these models as long as we could while maintaining a reasonable numerical cost to achieve the large parameter study presented here, and while computing the models showing a magnetic cycle over several decades. In this study we use the Rossby number as a measure of the influence of rotation on the flows maintaining the differential rotation as well as to power dynamos. Several Rossby number definitions have been proposed in the community, and we have computed the fluid Rossby number Ro f , the convective Rossby number Ro c , and the stellar Rossby number Ro s in Table 3. We refer the reader to  and their Appendix for a more in-depth discussion of these various definitions of the Rossby number. We will focus here on the fluid Rossby number Ro f and note that the other two are related to Ro f through a nearly linear relationship. The fluid Rossby number decreases with rotation rate and increases with mass, and varies from 0.07 to 2.35 in our sample of models. This range nicely covers the transition from solar-like to anti-solar differential rotation regimes (the transition is nearly at Ro f 1), and our smallest Rossby number (namely model M05R5m) is close to the expected fast-rotators saturation regime. We now briefly present one representative cyclic dynamo solution before entering a more detailed analysis of our dynamo simulations ensemble in sections 4 to 6. . Temporal evolution in case M09R3m of kinetic (KE) and magnetic (ME) energies. We also show their axisymmetric toroidal TKE, TME and poloidal PKE and PME components and their fluctuating components FKE and FME. We note the rise over about 500 days of ME just after having introduced a weak seed field. Then follows a modulation of ME with a 9 year period. Case M09R3m is indeed one of our cyclic cases (see also Figure 4).

OVERVIEW OF ONE CYCLIC DYNAMO CASE
To illustrate the richness of the dynamo solutions discussed in this study, it is key to show how the subtle nonlinear interplay between convection, rotation, and turbulence leads to the generation of time dependent complex magnetic fields. All 15 models discussed in detail in this study successfully generate and maintain a dynamo-generated magnetic field against Ohmic dissipation.
We defer the systematic comparison between all 15 models to the next sections and focus here on the representative case M09R3m. Indeed, M09R3m is in an intermediate Rossby number regime (R of = 0.27) and therefore lies in the middle of our sample of models. The temporal evolution of kinetic and magnetic energies is shown in Fig. 3. The magnetic energies first rise very fast to then saturate after about 1000 days in this case, and exhibit long-term oscillations over a decadal timescale reminiscent of a solar-like magnetic cycle. All components of the magnetic energy (toroidal, poloidal, fluctuating) oscillate in phase in this model. The mean toroidal kinetic energy also presents oscillation of the same amplitude, albeit anti-correlated with the magnetic ones. These energy trends are similar to the ones found in the magnetic cycles obtained with the EULAG code in Strugarek et al. (2018) (see their Fig. 3) and points toward a similar dynamo mechanism involving a strong feedback of the magnetic field on the differential rotation within the convective envelope. We will perform a detailed analysis of this mechanism in §5. Here, we first illustrate the dynamics of the dynamo achieved in model M09R3m in Fig. 4. The top row shows the 3D structure of our model by means of a potential field extrapolation outside our computational domain, at three different instances covering a magnetic reversal. We see that the field at the South Pole changes from blue to black, showing the polarity reversal. The strong toroidal field at the base of the convective envelope can be seen through transparency. In the leftmost panel, this deep wreath is mainly blue (westward oriented). Its polarity is reversed in the rightmost panel (red, eastward oriented), showing that the polarity reversal takes places over the full convective domain. The subsequent rows show spherical slices of B φ at the base of the convective envelope (second row), B r at the top of the domain (third row), entropy fluctuations (S ) at the top of the domain (fourth row), and enstrophy |∇ × v| 2 (last row).
We recover in the second row the magnetic wreath located at the base of the convective envelope and at midlatitude that changes polarity as the cycle progresses. The toroidal field reaches high values up to 1.5 × 10 4 G, with a strong temporal variation, as seen in the middle panel during the reversal. The surface radial field (third row) reaches values of about 100G and exhibits a complex topology, mixing dipolar and quadrupolar symmetries. We see again here that both fields oscillate in phase and reach a minimum in the midst of the magnetic reversal (middle panels). Finally, the two last rows show the thermal (entropy) fluctuations and the vortical motions (enstrophy) in our simulations. The first striking aspect is that these two quantities vary very little along the magnetic cycle. Indeed, the magnetic field modifies the large-scale motions and the average convective state in our models. Yet the magnetic cycles (when present)  The upper row shows a 3D potential extrapolation of the modelled magnetic field, with blue lines denoting field lines oriented outward and black lines oriented inward. Behind a semi-transparent representation, the radial velocity close to the surface, while deeper below the magnetic wreaths are shown by red (oriented eastward) and blue (oriented westward) lines. The second row shows the azimuthal field at the bottom of the convection zone, and the third row the radial field at the top of the domain. The fourth row shows the entropy fluctuations at the top of the domain, and the lowest row the enstrophy at the same depth.
show little imprint on the convective flows themselves and mainly act on the mean flows (see §4.1). The specific entropy fluctuations have two distinctive features. Firstly, a mean pole-to-equator contrast is well established in the model, with higher entropy fluctuations at the poles. Such contrast is expected in models with a solar-like differential rotation (Brun & Toomre 2002) and can be generally related to the pressure field required to drive the observed meridional flow. Secondly, patterns in S are imprinted by the non axisymmetric convective motions themselves, which are also recognizable in the enstrophy in the lower panel. The enstrophy is concentrated at the boundaries of the so-called banana cells at low latitude (Miesch et al. 2000), and is distributed between convective cell centers and boundaries at high latitude. We now turn to the detailed analysis of the large-scale flows ( §4), magnetic properties ( §5), and energetic balances ( §6) achieved in the 15 models. The reader mostly interested in the astrophysical consequences of our study may consider going directly to §7 for a summary.

LARGE-SCALE FLOWS IN THE MODELS
In this section we analyze the differential rotation profiles of the models including both a stable subadiabatic layer and magnetic field self-consistently generated by dynamo action. The aim of the study is to compare the differential rotation profiles of the hydrodynamical and MHD models. We confirm our preliminary results (Varela et al. 2016) and that of others (Karak et al. 2015;Guerrero et al. 2016;Viviani et al. 2018) that the presence of magnetic fields leads to different trends for the differential rotation with stellar rotation rate and mass when compared to their hydrodynamical counterpart. We further discuss how the meridional circulation is impacted by the presence of a magnetic field and discuss the main mechanisms acting to redistribute angular momentum within the convective shell. We also observe torsional oscillations in our set of dynamo simulations but delay their discussion to section 5.3.

Differential rotation profiles as a function of Rossby number
We analyze the differential rotation of the simulations that results from the angular momentum redistribution occurring mostly in the convection zone. The panels of Figure 5 show a meridional cut of the axisymmetric differential rotation averaged over 10 overturning convective times, defined as τ c = rtop r bcz dr/ṽ r (see Table 2). We observe that for the simulations M05Sm, M07Sm, M09Sm and M11R1m there is an anti-solar differential rotation, with the poles rotating faster than the equator, like their hydrodynamical counterparts (see Figure   6 and also ). The cases rotating at an intermediate rotation rate show a solar-like differential rotation. Finally, the cases rotating the fastest (R5 series) show almost no differential rotation (in particular for cases M05Sm and M07Sm). This constitutes a big difference with their hydrodynamical counterpart cases. The magnetic field here had a major impact, with almost solid body rotation imposed throughout the convective envelope. There is little asymmetry in the profiles between the northern and southern hemispheres, as expected when the average is performed over an interval long enough with respect to the convective overturning time (except for M05Sm for which the rotational constraint is the weakest and the longitudinal average less meaningful). Figure 6 also displays radial cuts of the rotation for the MHD cases (blue lines) and hydrodynamic progenitor cases (gray lines). In cases rotating 1, 3 and 5 times the solar rotation rate (bottom three rows), the velocity range in latitude (different styles of line) is generally reduced in the presence of magnetic field. This effect is observed to be stronger as the rotation rate increases. Conversely, the effect of magnetic field is mild for the slowly-rotating cases (upper row), except on the slowly-rotating case M11R1m which still shows some degree of magnetic feedback on its differential rotation. As one may expect, the radial gradient of the differential rotation nearby the tachocline is generally weaker in all the MHD simulations compared with hydrodynamic progenitors. This points to a magnetic feedback of the dynamo field on differential rotation itself, a feedback that is observed to strengthen as the rotation rate increases.
We have calculated the surface latitudinal differential rotation ∆Ω for each model, defined as the difference between the equator and 60 o latitude. A positive value thus denotes a solar-like differential rotation, and a negative value an anti-solar differential rotation. We report these values for the magnetic cases as well as the hydrodynamic progenitors in Table 2 (fourth column).
The differential rotation of our sample spans a range between -102 and +278 nHz, with some fast-rotating models presenting an extremely weak DR, like M05R5m with ∆Ω = 9 nHz. We find that the absolute differential rotation generally weakens in MHD models compared to their hydrodynamic progenitors, as expected from the radial profiles shown in Fig. 6. This is particularly striking for fast rotators such as M09R5m that goes from 338 nHz in hydro to 76 nHz in MHD.
We investigate in Fig. 7 the differential rotation trends with respect to the rotation rate (top panel) and rotation period (bottom panel). The differential rotation of the hydrodynamic progenitors and of the MHD  cases are respectively shown in small semi-transparent and large opaque symbols. The shape of the symbol labels the rotation of the model, and the color the mass of the modelled star, as indicated in the legend. In the bottom panel, we compare the model differential rotation to the differential rotation in the Kepler sample obtained by Reinhold & Gizon (2015) (shown as black dots). The dotted lines correspond to their estimated observational detection limits. We first note that the absolute value of our differential rotations agree well with the observed values. In addition, the differential rotation range in our sample increases as the rotation period decreases, like what is observed in the Kepler satellite sample. Sev-eral of our models nevertheless lie outside the observed values: the three anti-solar differential rotations on the right (triangles), and two of our fast-rotating models. Several reasons can explain this discrepancy. Slowlyrotating stars could produce very few starspots, or even no starspots at all (see for instance van Saders et al. 2019), making their differential rotation impossible to detect with photometry. Another possibility is that they lie outside the presently detectable limit with the Kepler data, due to their long rotational period (up to about 200 days for our most slowly rotating model). Finally, the two fast rotating models (M05R5m and M07R5m) show very weak differential rotations due to magnetic Dynamo action in G and K stars feedback, which are outside the detection limits of Kepler (dotted black lines). The top panel of Fig. 7 shows the differential rotation trend with the rotation rate. Using only the hydrody-namic progenitors, we previously showed that the differential rotation scales as Ω 0.66 . Blindly trying to fit such a power-law to the MHD sample, we find that the exponent reduces to Ω 0.46 . This weaker de-pendency is expected due to the magnetic feedback on the differential rotation through the Lorentz force. It also agrees better with the observational trends, which are still quite uncertain and were found to vary from 0.2 (Balona & Abedigamba 2016 for G stars), 0.3 (Reinhold et al. 2013 for cool stars) to even 0.7 (Donahue et al. 1996 for F-K stars). Looking more closely at our sample on the top panel of Fig. 7, it clearly appears that a power-law fit is a poor representation of the differential rotation in our sample. Rather, we see that ∆Ω increases with Ω for slow rotators, while it dramatically drops for fast rotators due to the magnetic feedback. Following Saar (2011), we recast in Fig. 8 the differential rotation trend in terms of relative differential rotation ∆Ω/Ω with respect to the fluid Rossby number R of (for the different definitions of Rossby number used in this work, see the caption of Table 3 or the appendix A of . We find a trend that is very similar to the observational trend reported by Saar (2011) (shown by the dashed line in Fig. 8): ∆Ω/Ω is roughly constant for inverse Rossby numbers lower than a certain threshold (here Ro −1 f 5), and it drops for fast rotators as ∆Ω/Ω ∝ Ro p f . Saar (2011) proposed that p = 2, but here our sample agree with a somewhat large range p ∈ [2, 6]. Additional models with even higher turbulence level are required to confirm the exact amplitude of the drop in differential rotation contrast found in the fast rotating cases. Finally, our sample also shows some hint of an increase of ∆Ω/Ω at large Rossby numbers, which is outside the observable constraints for now. It would be interesting to search observationally for candidate solar-like stars possibly possessing such anti-solar rotation states.
In , we have proposed that the differential rotation could follow two power-laws with respect to the Rossby number and the stellar mass. Here, we find that the differential rotation is weakened at high Rossby number, and therefore we do not recover a simple power-law trend, as we saw in Fig. 7. We can nevertheless attempt to fit such a combined power-law on a sub-sample of our models, excluding the fast rotating case but retaining the slow rotators. We obtain in this way In the MHD case, we find again that the differential rotation is less sensitive to both the Rossby number and the stellar mass. The power-law fit is nevertheless questionable here, as the range covered by our Rossby numbers and masses is quite small. We have nevertheless included the results of the fit here to compare with the purely hydrodynamic case . We can conclude here that the clear trend in stellar mass and effective temperature found in the hydrodynamic study ) is less significant when magnetism is taken into account, but overall we see a better agreement with observations of the dynamo models compared to their hydrodynamical progenitors.  On the bottom panel, the differential rotation in the Kepler sample obtained by Reinhold & Gizon (2015) is shown as black dots, and the observational detection limit by the two dotted black lines.
The MHD simulations therefore show that the magnetic field changes the angular momentum redistribu- The symbols shape and color are the same as in Fig. 7. The trend found in the MHD sample is highlighted by the grey area, and the observational trend reported by Saar (2011) is shown by the dashed black line.
tion, especially for fast rotating stars. In the next section, we perform a detailed analysis of this balance for 4 representative models.

Angular momentum transfer
We can better understand how the differential rotation profiles are achieved in the dynamo models by identifying the main physical processes responsible for redistributing angular momentum within rotating convective shells. Our choice of stress-free and potential-field boundary conditions at the top and stress-free and perfect conductor boundary conditions at the bottom of the computational domain have the advantage that no net external torque is applied, and thus angular momentum is conserved. We can assess the transport of angular momentum by considering the mean radial (F r ) and latitudinal (F θ ) angular momentum fluxes, applying the procedure used in Brun et al. (2004). Starting from the φ-component of the momentum equation expressed in conservative form and averaged in time and longitude: involving the mean radial angular momentum flux and the mean latitudinal angular momentum flux In these equations, the terms on the right-hand-side represent for both fluxes contributions respectively from viscous diffusion (which we denote as F VD In Figure 9 we show the components of F r and F θ for M07 case series, having integrated over co-latitude and radius as follows: Thus Φ r represents the net angular momentum flux through horizontal shells at different radii and Φ θ represents the net flux through cones at different latitudes. This representation is helpful in assessing the direction and amplitude of angular momentum transport within the computational domain by each component of F r and F θ . For each of the four cases we display Φ r on the left panel and Φ θ on the right panel both normalized by R 2 * . Turning to the radial angular momentum transfer, we first note a very good overall radial balance. We find that the Reynolds stresses (green dash-dotted curves) transport angular momentum outward in all the low Rossby number models. By contrast, M07S the slowly rotating case has the Reynolds stresses transporting angular momentum inward. The viscous diffusion and Maxwell stresses oppose this transport, tending to rigidify the rotation state in the radial direction. The meridional circulation has one large cell per hemisphere for the M07Sm case (see §4.3). It opposes the Reynolds stresses, but as the rotation rate increases and the Maxwell stresses gain in amplitude, it changes in profiles and direction to yield a radial balance of angular momentum, from the angular momentum equation. Note that the mean large-scale magnetic torques (black Figure 9. Angular momentum transport in M07 case series. We display cases M07Sm (left) and M07R1m (right) in the top row and M07R3m (left) and M07R5m (right) in the bottom row. For each model radial (left) and latitudinal (right) angular momentum balance is shown, with Reynolds stresses contribution shown as green dash-dotted lines, viscous stresses as dashtriple-dotted yellow lines, meridional circulation as dashed cyan lines, Maxwell stresses as long dashed magenta lines, the large scale magnetic torques as black dotted line and the sum of all contributions as solid black line. We note that a very good angular momentum balance is achieved in most models in both directions with the sum being close to zero. dotted line) have little influence in the overall radial angular momentum balance.
Considering now Φ θ , we can assess the balance of latitudinal angular momentum transport. We first notice that the Reynolds stresses (green curves) are systematically equatorward in both hemisphere (positive in the northern hemisphere and negative in the southern one). Since most cases exhibit a very good latitudinal balance, as demonstrated by the solid black curve, these Reynolds stresses must be nicely counterbalanced. A quick survey of the right panels for all four models indicates that many contributors act depending on the Rossby number of the simulations. For the slowly rotating case (upper left corner), we see that it is mostly the meridional circulation (cyan dashed curve) that does most of the work (we defer the reader to section 4.3 for a discussion of the meridional circulation patterns in the various dynamo cases). By contrast, magnetic terms do not play much role in the case of M07Sm. For M07R1m (right top corner) it is now the viscous diffusion that plays that role of opposing the Reynolds stresses. For that case the meridional circulation is not doing much, but we do see a 20% contribution of the large scale magnetic torques, the Maxwell stresses being still weak. As the Rossby number decreases and the dynamo action becomes more intense, we see that the magnetic terms start influencing the latitudinal angular momentum transport more and more, tending to oppose the Reynolds stresses. It is particularly noticeable for the Maxwell stresses. They are the dominant player for M07R3m case (bottom left corner), helped by the large scale magnetic torques. In that case the meridional circulation is somewhat helping the Reynolds stresses, notably at low latitudes near the equator.
For the M07R5m case, the story becomes less clear, except for the Reynolds stresses all terms fluctuate and sometimes oppose or reinforce the turbulent stresses. Maxwell stresses still play an important role as does the meridional circulation. In that model, the differential rotation has been so significantly quenched by dynamo action, that it is not surprising the trends are less clear and systematic. In summary, in most cases the transport of angular momentum by Reynolds stresses are opposed by a combination of meridional circulation, viscous stresses and Maxwell stresses.

Meridional circulation profiles
The meridional flow patterns are also affected by the presence of magnetism in our set of models, especially for the fast rotating cases. We immediately note that the meridional circulation is indirectly modified by magnetism (as will be made clear in §6.2). Indeed, magnetic stresses play a negligible role in setting the meridional flows in our models, and the differences we observe compared to the hydrodynamical counterparts originate from changes in the differential rotation (see e.g. Passos & Charbonneau 2014).
We illustrate the meridional flow pattern achieved in the M07m set of simulations in Fig. 10. The slow rotating case (first panel) is very similar to its hydrodynamic progenitor, with a well-defined circulation cell in each hemisphere. Both cells circulate from the equator to the pole at the surface, and from the pole to the equator at the base of the convective envelope. The second model rotating at the solar rate (second panel) is also similar to its hydrodynamical progenitor and shows a more complex circulation profile. These are consistent with previous numerical experiments by e.g. Karak et al. (2015). It consists of stacked cells elongated along the rotation axis outside the tangent cylinder, and two counter-rotating cells in each hemisphere at high latitude. Finally, the fast-rotating models (third and fourth panels) exhibit a peculiar meridional circulation pattern concentrated at the equator, with two stacked transequatorial cells (see e.g. Simitev & Busse 2009). These profiles can be understood as follows. In these models, the differential rotation is strongly quenched by magnetic feedback as seen in the previous section. In particular, the radial shear of differential rotation vanishes at the equator as seen in Fig. 6. As a result, gyroscoping pumping McIntyre 2007;Featherstone & Miesch 2015) dramatically weakens along the equator and the resulting meridional circulation is both very weak (this can be seen in the drop of meridional flow kinetic energy in Table 5) and mainly driven by the remaining latitudinal shear. This leads to two meridional cells crossing the equator, as seen in the last panels of Fig. 10. Having presented the large-scale flows achieved in the simulations, we now turn to discussing their magnetic properties.

MAGNETIC PROPERTIES
In this section, we discuss in more details various aspects of our dynamo simulations, such as their type, their temporal variability, the amount of magnetic flux they generate and the distribution in space and size of their magnetic fields.  Table 4. Magnetic properties of the modelled dynamos. The first column indicates the presence or absence of a long (decadal), deeply-seated magnetic cycle. When the time-series were long-enough to identify unambiguously a cycle period, its value is given with error bars. Otherwise, the existence of such a cycle is indicated by a yes ('y'), and its absence by a no ('n'). The second column shows the same for the short magnetic cycle that we identify in the upper convection zone near the equator. We do not indicate this information for model M11R5m that was not run long enough to determines the existence or absence of magnetic cycles. The third column indicates the amplitude of the torsional oscillations at the surface in nHz (see §5.3). The fourth column shows the active latitudinal band at the bottom of the convection zone, based on the azimuthally-averaged and temporally-varying azimuthal field straddling the base of the convection zone. The fourth column shows the total magnetic flux at the surface, in units of 10 24 Mx, with minimum and maximum as subscript and superscript (see §5.4). The three next columns show the root-mean-squared surface field in Gauss, the surface dipole in Gauss, and the surface large-scale radial field B L,surf (taken for l < 5) in G with the same layout (see §7). Finally, the last two columns show the fractions of the large-scale dipole (f dip ) and quadrupole (f quad ), as defined in §5.5.
magnetic cycles, and stable magnetic wreaths concentrated close to the bottom of the convection zones. These three states are illustrated in Fig. 11 with models M07R5m, M09R3m and M09R1m. Let us first focus on the decadal cycles, as the one found for M09R3m (see middle left panels in Fig. 11). In this model, we find that the global magnetic field of the star reverses with a period of 10 years (see first column in Table 4). The averaged azimuthal field at the bottom of the convection zone presents a solar-like butterfly diagram, with both a polar and an equatorial branches. The magnetic field is generally consistent with dipolar symmetry, with azimuthal field of inverse polarities in each hemisphere. We also see some departures from hemispheric symmetry (for instance around t=42 years). The azimuthal field is found to be concentrated at the base of the convective envelope and in the tachocline, where the radial shear of Ω is maximized, as shown in the time-radius and meridional diagrams. It develops over a relatively large latitudinal extent, as shown by the active latitudinal band reported in the fourth column of Table 4. We find this band to be centered at higher latitudes the slower the model rotates for low and intermediate Rossby numbers. Conversely, this activity band moves to high latitudes for models with anti-solar differential rotation. The averaged radial magnetic field at the surface is also found to reverse with the same timescale. At the surface, the migration branches are nevertheless not as clear as deep within the convection zone in this case. Strugarek et al. (2017Strugarek et al. ( , 2018 also found similar deeply-seated cycles using the EULAG code, as well as Augustson et al. (2015) using the ASH code. The cyclic behavior in their results originate from the nonlinear magnetic feedback of the large-scale Lorentz force onto the differential rotation. This weakens the source of mean toroidal field that decreases and reverses, while the associated poloidal field closely follows due to the sign inversion of the electromotive force. We find the same mechanism in this new sample of simulation with the ASH code. Indeed, the DRKE cyclic variation observed in Fig. 3 Figure 11. Various dynamo states achieved in our sample our models, as illustrated by M07R5m (top panels), M09R3m (middle panels) and M09R1m (lower panels). In the top panel we show time-latitude, time-radius, and instantaneous meridional plane of Br (red denotes positive values and blue negative values), with sampling times indicated by a vertical dashed line. These illustrate the short magnetic cycles achieved by our models. The four middle panels illustrate both the short and long cycles achieved in model M09R3m. The first panels show latitude-time (at the base of the convection zone) and radius-time (at mid-latitude) diagrams of the mean azimuthal magnetic field that reverses on a decadal timescale. The two panels below show the mean radial field at the top, which also shows the same cyclicity. Once the long cycle is filtered, the short cycle appears in the zoomed panel on the right at a particular epoch and around the equator. The lowest panels show the mean azimuthal field for model M09R1m and illustrate a dynamo with no cycles but which sustains strong stable wreaths at the base of the convective envelope.
of the existence of a subtle nonlinear feedback loop between the large-scale shear and the toroidal magnetic field, is therefore confirmed by the present study using a different numerical code than Strugarek et al. (2017). We stress that its existence can be unveiled here only because we consider fully nonlinear convective dynamos, with a self-consistent differential rotation maintenance and magnetic field generation.
Still, we have attempted to interpret our simulations through mean-field dynamo theory by inverting the α tensor and its antisymmetric part γ through the means of singular-value decomposition (SVD) technique (see Augustson et al. 2015;Simard et al. 2016). The details of this procedure are given in Appendix B. One can then use the derivedᾱ profile to compute the Parker-Yoshimura rule (Parker 1955;Yoshimura 1975) and assess the consistency of a mean-field approach with our 3D turbulent model. We therefore compute where λ = r sin θ. The time-latitude variations of S θ are shown at the base of the convection zone of M09R3m in the top panel of Fig. 12, with red/white denoting a southward migration rule and blue/black a northward migration rule. We overlay contours of B φ as black contours (plain/dashed denoting positive/negative contours) in the top panel. We see that the derived Parker-Yoshimura dynamo wave rule does not agree with the observed latitudinal propagation, which strengthen our interpretation in terms of a cycle dominated by the Lorentz-force feedback on the differential rotation itself. We also find another type of cyclical behavior in our sample of models: short cycles, that seem to preferentially be sited close the equator and in the upper part of the convection zone. Such types of cycles have already been reported in previous publications with numerical models (Käpylä et al. 2016;Beaudoin et al. 2016;Strugarek et al. 2018) and could be reminiscent of the possible quasi-biennial oscillations observed in the Sun (Broomhall et al. 2012;Simoniello et al. 2013). They oscillate on a yearly timescale as shown in the second column of Table 4. Short cycles are interestingly found in almost all of our models except the slowly-rotating cases. Two short cycles are illustrated in Fig. 11 for case M07R5m and M09R3m. In the former fast rotating case, no long deeply-seated cycle is observed, and the short cycle clearly appears in both the latitude-time and radius-time diagrams. In the case of M09R3m, both types of cycle are found at the same time, and the short cycle appears clearly once the signal of the long cycle is removed (see zoomed panel). The short cycles are found to always show a poleward propagation branch, and to be concentrated close to the equator. We have performed the same SVD analysis and show the Parker-Yoshimura rule S θ (Eq. 13) for model M07R5m, which is shown in the bottom panel of Figure 12. In this case, the analysis is carried out in the upper part of the convective envelope, and contours of B r are overlaid above the propagation rule. The Parker-Yoshimura rule is found here to be in good qualitative agreement with the poleward branch, suggesting that an α − Ω or an α 2 − Ω dynamo could be at the source of this type of cycles. The short cycles furthermore embed much less magnetic energy than the deeply-seated ones, and we do not find any clear DRKE beating associated with them. As a result, we find that the two types of cyclical behaviors likely originate from two different dynamo processes: the deep-seated cycle from the large-scale feedback loop between the magnetic field and the differential rotation through Maxwell torques, and the short cycles from the standard α−Ω or α 2 −Ω dynamo loop. Finally, short cycles were also reported in the study of Strugarek et al. (2018), with the same type of localization within the convective envelope. In this former study, the short cycles were only found at small Rossby number, i.e. for the fastly rotating cases. Here we find short magnetic cycles much more ubiquitously in our sample models, as they only disappear at large Rossby numbers. It is possible that the coarse resolution used in Strugarek et al. (2018) with the EULAG code prevented models at intermediate Rossby number to develop such magnetic cycles. Additional modelling effort pushing the turbulence level of the simulations are required to properly assess this point, which is left for future work.
Finally, some models in our sample do not present any cyclical behavior. Instead, they sustain a steady dynamo with stable magnetic wreaths within their convective envelope and tachocline. This is the case for instance of model M09R1m shown in the lower panels of Fig. 11. We obtain such solutions only in the high Rossby number regime, close and above the transition toward an anti-solar differential rotation.
To summarize, we find that the different types of cyclical behaviors exist in specific Rossby number ranges in our sample. We illustrate this in Fig. 13 where we follow Gilman (1983) and show DRKE/KE as a function of Ro f in our set of models. Short cycles are found for Ro f 0.42, deeply-seated solar-like cycle for 0.15 Ro f 0.65, and steady magnetic fields for Ro f 1.0. The exact boundaries between these cyclical behaviors regimes are not precisely defined and may depend on a number of factors. First, let us note that the same trend was found in Strugarek et al. (2018) with the EULAG code, as shown by the colored stars also plotted in Fig. 13. This is very important because it again demonstrates that the results discussed in this study are not code or setup dependent, but the results of genuine nonlinear convective dynamo action in a rotating spherical shell. It confirms that the Rossby number is one of the key parameters to characterize the various dynamo states found in the literature, and that cyclic convective dynamo solutions clearly exist in a parameter regime that our study helps to refine. The transitions between the different types of cycles were found at slightly different Rossby numbers, possibly due to different Reynolds, Prandtl and Rayleigh numbers regimes achieved in the two ensemble of simulations. Indeed, Nelson et al. (2013) showed that fast-rotating models exhibiting stable wreaths of magnetism ) could produce reversals when the Reynolds and Rayleigh numbers are increased. Since the Rossby num-ber of the more turbulent models nevertheless change significantly as well, it is therefore unclear whether this can be attributed to a fundamental change in the dynamo action or if it is the consequence of a change in Rossby number. Fundamental exploration aimed at predicting the Rossby number of turbulent numerical experiments such as Anders et al. (2019) are very promising in that respect, and need now to be extended to the full MHD regime. For the time being, we can conclude here that qualitatively the different regimes highlighted by our simulations are robust, yet simulations at much higher turbulent levels are required to assess the exact regime boundaries. Please note that case M11R5m is sometimes omitted in ensemble analysis in §5 and §6 because it is not as well numerically converged as all the other cases and can sometimes be an outlier in some analysis. This does not impact our conclusions in any of the results reported in the paper.

Short cycles Cycles No cycles
This work S17 Figure 13. Summary of the dynamo states found in our study (circles) and in the previous study of Strugarek et al. (2017) (stars). In both studies, we find a clear trend in the type of cyclical behavior that models tend to produce as a function of the Rossby number. They are shown here by the ratio between the differential rotation and total kinetic energies. For small Rossby numbers, only short cycles are found. At intermediate Rossby numbers, decade long cycles resembling the solar cycle start to appear on a relatively narrow parameter space. At high Rossby numbers, magnetic cycles disappear and our models produce energetic stable wreaths of magnetic field in their convective envelopes.

Dependencies of the cycle periods
We have calculated the period of the short and long cycles and reported their values in the second and third columns of Table 4. We use the approach initially followed by Käpylä et al. (2016); Strugarek et al. (2018) and rely on an empirical-mode decomposition method ( Luukko et al. 2015) to identify quasi-periodic signals. Five of our models clearly exhibit a deeply-seated long cycle that can be identified by eye. We were nevertheless able to calculate accurately the associated period for four of them. The cycle period of the fourth model would require at least twice as long integration times to be identified. This would require an even more massive numerical effort and will be explored in future work. Still we can deduce with some confidence what characterizes this long cycle nonlinear dynamo case. Conversely, the short cycles take place higher up in the convective envelope and their short periods allow us to determine the cycle periods for all the models exhibiting them. The error bars on the cycle periods are directly estimated with the empirical-mode decomposition method, as explained in Strugarek et al. (2018). The left panel of Fig. 14 shows the cycle periods (in years) as a function of the rotation period (in days) of our models. We report both short and long cycles here, respectively in blue circles and red circles. We have also added the cycles found with the EULAG code and reported in Strugarek et al. (2017Strugarek et al. ( , 2018 as red and blue stars. Finally, we have overlaid the detected cycles of distant stars reported by Böhm-Vitense (2007) as gray squares, as well as the Sun right in the middle of the figure. Our three identified long cycles are achieved by models with different masses, which makes their direct comparison subject to caution in a (P cyc , P rot ) diagram. Overall, we do not recover the dichotomy between active and inactive branches as initially proposed by Saar & Brandenburg (1999) and Böhm-Vitense (2007). Rather, our sample of models combining ASH and EULAG simulation spans the whole diagram, including the hypothetical gap where the Sun stands.
Using the EULAG sample of simulations only, we have previously shown that the cycle period is controlled by the effective Rossby number achieved by the simulated convection zone . This is shown for the long and short cycles in the middle and right panels of Fig. 14. Here we find that our new ASH simulations are compatible with the trends obtains with the EULAG sample, which strengthens the similarities between the modelled dynamos in our two studies. This is moreover remarkable as the ASH simulations include a tachocline and a deeper radiative layer, whereas the EULAG sample considered only an isolated convective shell.
The fact that the cycle period seems to decrease with the Rossby number has also been reported by other research groups using yet another code (see e.g. Warnecke 2018). So far only one study relying on 3D turbulent simulations (Guerrero et al. 2019) has shown some evidence for a cycle period increasing with rotation period. We believe this is due to how their differential rotation scales with rotation rate. Indeed, their simulations exhibit a differential rotation that strengthens as the rotation rate decreases (i.e. the rotation period increases). This is at odds with all the aforementioned studies (including the present work), where we find it to increase with the rotation rate up to a point where magnetic feedback strongly back-reacts to suppress it. We suspect that the thermal treatment of the radiativeconvective interface may produce this effect in the work of Guerrero et al. (2019), albeit additional analyses are required to confirm this interpretation.
Finally, it is worth noting that more complex dynamo states have also been reported in a similar Rossby number regime with the PENCIL code by Viviani et al. (2019). This warrants again caution in the interpretation of simulations results at moderate Reynolds number, and highlights the need of achieving more turbulent regimes in future work to confirm our trends.

Torsional Oscillations in Cyclic Solutions
We observe clear and strong torsional oscillations δ t Ω (Basu & Antia 2019) in all our models that exhibit a long, deeply-seated cycle. Torsional oscillations take the form of a modulation of the azimuthally-averaged rotation rate Ω(r, θ, t) in both depth, latitude, and time. We illustrate the torsional oscillations at the base of the convection zone of model M09R3m in Fig. 15. The torsional positive/negative oscillations are shown in red/blue in nHz as a function of time and latitude. We have overlaid iso-contours of B φ in black (plain lines correspond to 4000 G, dashed lines to -4000 G). The torsional oscillations are observed to be in phase with long magnetic cycle. At cycle minimum (in between black contours), the poles are rotating slower (blue) and the equator faster (red), meaning that the latitudinal differential rotation is strengthened as the magnetic field weakens and the associated magnetic torque stops inhibiting it. During cycle maximum, the opposite situations occur, and the differential rotation is found to decrease substantially. We observe torsional oscillations very similar to what was found with EULAG simulations by Strugarek et al. (2017Strugarek et al. ( , 2018 and previous ASH simulations by Nelson et al. (2013); Augustson et al. (2015). In all these simulations, the torsional oscillations are found to play a major role in producing the deeply-seated cycle. This is reassuring because such nonlinear interplay between the flow and field seems independent of setup details such as BC's or numerical schemes. Moreover, torsional oscillations in our models are very energetic: they reach more the 20 nHz at the base of the convective envelope in model M09R3m, and their energy corresponds to the energy variations in the total magnetic energy (ME) seen in Fig. 3. As a result, we find they play an active role in allowing deeply-seated cycles by reversing locally ∂Ω/∂θ and hence generating a toroidal field of opposite sign.
We have also searched for torsional oscillations at the locations of short magnetic cycles, i.e. at the surface and close to the equator of fast rotating models. We find a temporal modulation of the local rotation rate at the surface in all our models. We have nevertheless not found any evidence for a correlation between these temporal variations and the short cycles themselves. This confirms that a different dynamo process sustains the short cycles, which is likely related to a more standard α − Ω mechanism as we have seen in Sec. 5.1.
Finally, we have characterized the surface torsional oscillations in all our models and reported in Table 4 the average values of δ t Ω within the activity band identified in Table 4. The surface torsional oscillations range from about 1 to 39 nHz in our sample of simulations, which corresponds to 0.4 to 6% of the model rotation rates. Torsional oscillations associated with short cycles are found to be very weak, and the ones associated with the long cycle to be prominent deep inside the convective envelope. As a result, we do not observe any strong correlation between the amplitude of the surface torsional oscillations and the Rossby number of our models: a linear regression gives δ t Ω/Ω ∝ Ro 1.1±0.15 f Ro f .

Magnetic flux budget
To further assess the magnetic properties of a dynamo solution, we display in Figure 16  ous measures of the magnetic flux available at the top boundary layer, that is Φ + , Φ − , the magnetic fluxes for B r | r=rtop > 0 and B r | r=rtop < 0 respectively, the total flux Φ tot = |Φ + |+|Φ − |, the net flux Φ net = Φ + +Φ − and the southern Φ S and northern Φ N hemispheric fluxes (i.e. integrated only over the northern and southern hemisphere respectively). First, we see the very good conservation of divergenceless nature of the magnetic field, with Φ net being systematically null (so implying that Φ − = −Φ + , as clearly evident). This is the direct consequence of solving the induction equation via a poloidal-toroidal field decomposition (see Eq. 6). Likewise, the two hemispherical measures of Φ have opposite signs, but a much smaller amplitude than Φ + , Φ − by about a factor 10. This is likely due to a highly structured magnetic field, since for an axial dipole they are expected to be equal. When adding up the absolute value of Φ + , Φ − , we can assess the total amount of magnetic flux generated by the dynamo. We find fluxes from 10 24 to 10 25 Mx, which are in good agreement with values observed in the Sun (see for instance Fig. 3 of Schrijver & Harvey 1994). We also note that in M05R1m and M09R3m cases, both of which possess a clear and long magnetic cycle, the temporal modulation of the magnetic fluxes is obvious. In M05R1m case, the modulation is about a factor of 2 from minimum to maximum of activity. In case M09R3m it reaches almost a factor of 8 (compared to 5 for the Sun). Here again the larger mass (luminosity) of M09R3m and its higher rotation rate leads to larger temporal modulation of the magnetic energy and hence the magnetic flux. Finally, for the steady dynamo case M11R1m, possessing an anti-solar differential rotation, a very small magnetic flux variability is observed. However, it is the model with the highest value of total magnetic flux, reaching about 10 times what is observed in the present Sun.
We furthermore see a tendency for Φ tot to increase with both stellar mass and rotation rate, in good agreement with the level of magnetic energy found in the simulations. However, more robust tendencies appear on the rotation when one considers only model with Ro f < 1. They are interestingly compatible with a simple linear dependency, with Φ tot 2.3 Ω 0.84±0.42 for the rotation rate. When considering how the total magnetic flux scales with rotation rate Φ tot ∝ Ω n * , different values from n = 1.2 (Saar 2001) to n = 2.8 (Schrijver et al. 2003) have been proposed (Rempel 2008). In our study we find a tentative scaling with the fluid Rossby number as Φ tot 1.19 Ro −0.88±0.31 as shown in Fig. 17, where the time-averaged total flux of each model is considered (see also Table 2). Our models depart significantly from this trend when their Rossby number exceeds one, indicating a possible change for very slowly rotating stars. In this regime, our sample of models suggests that the total magnetic flux increases with Rossby number, as shown by the dashdotted line. Additional models at large Rossby numbers are required to fully characterize this regime properly, which we leave for future work. To summarize, we find that the total magnetic flux follows a trend compatible with the one from Saar (2001) for intermediate and small Rossby numbers, and that this trend reverses for slow rotators (Ro f > 1).

Dynamo families and f dip
We now turn to considering how the change of differential rotational state as a function of the Rossby number may influence the relative amplitude of the dynamo modes. We have seen in the previous sections that as we vary the Rossby number the type of dynamo solution  Hall et al. (2021) that the Sun and solar-like stars older than the Sun may be undergoing a magnetic activity transition around a Rossby number of 1 (see Lorenzo-Oliveira et al. (2018) for an alternative view). In particular, they argue that the wind braking efficiency may be collapsing around that rotational state transition. This would result in stars rotating more rapidly than what Skumanich law or gyrochronology would have predicted (Skumanich 1972;Barnes 2003). If for instance a collapse of the large-scale dynamo modes (mainly dipole and quadrupole) would occur after transiting to anti-solar differential rotation this would provide a very simple explanation, as it is well known that the most efficient wind braking for sun-like stars is found for the simplest magnetic field geometry (Kawaler 1988;Réville et al. 2015;Finley & Matt 2018). In order to assess if such a change of magnetic geometry occurs at or near the R of ∼ 1 limit, we will use a measure called f dip , that was introduced by Christensen & Aubert (2006), and that permits the assessment of < Ω Ω 3Ω 5Ω Figure 18. f dip and f quad in all 15 models (color symbols) and those published in Strugarek et al. (2017Strugarek et al. ( , 2018) (gray stars). Models with Rossby number greater than 1 possess an anti-solar differential rotation. We see only a weak decreasing trend of f dip and f quad with Rossby number (for the parameter space explored). In addition, there does not seem to be a collapse of the large-scale magnetic field for slowly rotating stars. the energy content of the dipolar field with respect to the first 12 magnetic modes. We also introduce f quad , using the same principle, as a quadrupolar field configuration is still quite efficient at spinning down a star via its associated wind braking. Both are defined as where a l,m are the spherical harmonics coefficient of the radial magnetic field at the upper boundary (surface) of our models.
In Figure 18 we show f dip (top panel) and f quad (bottom panel) from 32 dynamo cases: the 15 cases analyzed in details in this paper, to which we add the 17 published in Strugarek et al. (2018), using the EULAG-MHD code (Smolarkiewicz & Charbonneau 2013). This allows us to extend our database and to compare nonlinear dynamo solutions obtained with two different MHD codes using very different numerical techniques, hence giving us confidence that the trend found in our simulations is not due to a given code. We observe a relatively good agreement between the ASH and EULAG databases for f dip , and surprisingly find that the EULAG set of simulations produces systematically a weaker f quad compare to the ASH database. In both series we find a weak trend for a decrease of f dip and f quad with Rossby number. Nevertheless, we do not find any hint of a collapse of f dip or f quad when the Rossby number exceeds 1 and the differential rotation realized in the simulations becomes anti-solar. The weak decreasing trend is not significant enough to explain the stalling of stellar wind braking advocated by van Saders et al. (2016); Metcalfe et al. (2016). Hence, it seems unlikely that field geometry is the source of the wind braking regime change for old solar-type stars. This is in agreement with the observational study of Vidotto et al. (2016), whom have analyzed spectro-polarimetric inversion for a suite of sunlike stars, and they too did not find a collapse of the dipole strength as they crossed the Ro f = 1 limit. So if such a stalling of stellar spin down occurs, it must come from another mechanism (see §7).
In summary, we have shown in §5 that the dynamo solutions presented in this study possess very interesting magnetic properties that agree very well with observations and other theoretical studies. In particular, we have confirmed the key role of the Rossby number (and magnetic Reynolds number) in determining the type of dynamo realized. Now we wish to characterize better their energy content and how energies flows back and forth from kinetic to magnetic reservoirs.

ENERGY CONTENT AND TRANSFERS IN STELLAR CONVECTIVE DYNAMOS
In the following section we analyze the kinetic and magnetic energy contained in the models and how they are distributed between their various components.

Global measure of kinetic and magnetic energies
We now turn to discussing the global energy content in the convective envelope of the 15 dynamo cases presented in this study. In Table 5 we list the kinetic (KE) and magnetic (ME) energy densities and their axisymmetric and non-axisymmetric components (see their definition in Appendix A and Brun et al. 2004). We first notice that as we increase the stellar mass, the KE is found to slightly decrease. This is due to the lower averaged mean density due to the shallower convective envelope in more massive stars. The averaged density over the simulated convective envelopes varies from 4 to 0.05 g/cm 3 when going from models M05m to M11m, so a drop by a factor of 80. This is in part compensated by the higher luminosity (convective velocity) of the more massive stars, leading to values of KE in the range of 10 6 to 10 7 erg/cm 3 . Note that the total kinetic energy (i.e. the energy density multiplied by the volume) increases with stellar mass due to the much larger volume occupied by larger-mass stars. If we now decompose KE into its axisymmetric poloidal (MCKE) and toroidal (DRKE) and non-axisymmetric (CKE) components, we can further understand how the energy is being distributed in the various models.
First, as it is often the case, MCKE is found to play a minor role in all models independently of their mass or rotation rates. In most cases MCKE is of the order of 10 4 erg/cm 3 so about 1% or less of KE. This results in DRKE and CKE being the dominant components. Analyzing these two components, a clear trend is observed common to all masses. As the rotation rate is increased, going from Rossby number greater than 1 to value less than about 0.1, we note that DRKE first increases to constitute up to 96% of KE. This means that most of the kinetic energy is in the differential rotation with both strong latitudinal and radial shear across the convective envelope and at its base (we refer the reader to §4.1 where the angular velocity profiles of each model is discussed in details). Such a behavior is similar to what was observed in the purely hydrodynamic progenitors published in . Hence, up to a certain rotational influence, the presence of dynamo generated magnetic fields in the simulations does not modify significantly the trends observed before in the hydrodynamic cases. As a direct consequence, CKE is found to contribute less and less to the overall dynamics. CKE is found to be dominant for the slowly rotating cases, their convective motions having little azimuthal mean. As the Rossby number is decreased and the rotational influence on convective motions made stronger, we see that CKE drops to be less than a few percents of the total kinetic energy. However, this is not the case when the rotational influence increases even further. For all the fastest cases with the smallest Rossby numbers, we notice a sudden drop of DRKE both in percentage and absolute value, while CKE contributes relatively more to KE (but KE also undergoes a decrease of its amplitude). This is due to the strong feedback of the Lorentz force on the differential rotation, a phenomenon often called Ω-quenching (Glatzmaier 1985a;Brun 2004;Brun et al. 2005;Karak et al. 2015) and seen only in global spherical rotating models by similitude to α-quenching (Blackman & Field 2001;Subramanian & Brandenburg 2004;Brun et al. 2004) found in most local dynamo simulations (at the origin of the interface dynamo paradigm Parker 1993; Mason et al. 2008) and characterized in our simulations by the absolute concomitant drop of CKE. This significant drop of DRKE or "Ω-quenching", accompanied by a smaller decrease of CKE or "α-quenching", leads to a strong decrease of KE. This confirms that dynamo simulations do not have the same rotational dependence as the purely hydrodynamic cases. Since most solar-like stars are likely to have magnetic fields, such a finding indicates that scaling laws derived in this work will likely be more accurate when compared to observations. Since the influence of magnetic field becomes more and more dominant as we lower Ro f , it is also instructive to analyze how the magnetic energy content evolves as well.
In Table 5, we also provide the value of the magnetic energy densities (total magnetic energy [ME], axisymmetric poloidal [PME] and toroidal [TME] components, and non axisymmetric components [FME]). Here there are some surprises given what we just discussed for their kinetic energy counterparts MCKE, DRKE and FKE. First, the axisymmetric poloidal component PME con-tributes more to total ME than MCKE contributes to KE. It often represents few % of ME and in one case M11R5m it is even found to be dominant. Interestingly, PME is found to reach its lowest values for intermediate rotators close to the Ro f = 1 regime. In the Ro f > 1 we find that PME rises again, confirming the trend we observed on the total magnetic flux in §5.4. TME somewhat follows DRKE, it first increases with rotation rate, more and more energy being pumped by the large-scale shear into toroidal magnetic energy via the dynamo Ωeffect and also via complex convective motions. TME can reach values between 80 to 85% of the total ME. However, the Lorentz force feedback is so strong past a certain point that the large shear is quenched (the feedback destroying its generating source). In most of these highly rotationally constrained cases, the magnetic energy is found in the non-axisymmetric magnetic field. These trends are also illustrated in Fig. 19. It is worth noting that the three magnetic energies show an overall similar trend: the total energy density ME decreases with an increasing Rossby number until Ro f 1. The four models at Ro f > 1 then exhibit a large scatter, and only PME shows an unambiguous increase with Rossby number in this regime. We also see a hint of a saturation and possibly a slight decrease of TME at low Rossby numbers. Additional simulations at even lower Rossby numbers are required to confirm this trend, which is to be expected based on the observed saturation of magnetic activity for fast rotators (see e.g Wright A.S. Brun et al.  Reiners et al. 2014). In all panels, we have indicated the inverse Rossby number trend as a gray dashed line. We remark that the three magnetic energy densities are all compatible with Ro −1 f trend at intermediate Rossby number, as expected from standard dynamo scaling laws in this regime (see Augustson et al. 2019). This translates into a bulk magnetic field B bulk ∝ ME 1/2 ∝ Ro −0.5 f . We note that this scaling does not necessarily translate into the same scaling for the surface large-scale magnetic field, as will be made clear in §7.
The relative energies (shown as percentages in Table  5) also present interesting trends. We note first that in the slowly rotating cases, ME is only a few percent of KE. As we lower Ro f , this value increases to reach equipartition by a subtle combination of both ME increasing while KE first increases and then decreases, as we have just seen. These variations inevitably lead to the fact that for the fastest rotating cases ME is even larger than KE and the simulations are in a so-called global super-equipartition state. This is very interesting, because it means that the kinetic energy in the convective envelope is not the maximal value that the magnetic energy can reach. This is due to a change in the force balance in the Navier-Stokes equation between Lorentz, inertia, buoyancy and Coriolis forces. As the rotation rate is increased and the Coriolis force becomes stronger and stronger, the balance at first shifts from being between mostly inertia and Lorentz force to a magnetostrophic state that implies a balance between Lorentz and Coriolis forces. We refer the reader to these following studies for more detailed discussions of dynamo scaling laws (Christensen 2010;Davidson 2014;Oruba & Dormy 2014;Brun et al. 2015;Augustson et al. 2019).
Overall we see that the dynamo states reached in our 15 cases do not show a strong difference as a function of mass, at least in the range studied here. However, both in terms of amplitude of the magnetic field and in the time variability of the magnetic field (cyclic, unsteady or steady solutions), we confirm that rotation plays a key role in determining the type of dynamo found in our simulations. We also note that the mean axisymmetric magnetic fields are not negligible in most of the models, often reaching values of 5 % of the total energy content for the poloidal field and a large fraction for the toroidal magnetic field. For the latter, this has important consequences for the energy made available for the various magnetic phenomena occurring at the surface of solar-like stars (see §7).
Note that we did not look for hysteresis around the Ro f = 0.1 limit, by running various cases with different value of the seed magnetic field, as was done in some geophysical dynamo studies (Schrinner et al. 2014). We consider that stars acquired their magnetic field through a complex formation process, in which the seed magnetic field is likely very weak (interstellar medium magnetic field amplitude are on averaged about 10-100 microGauss) and that starting the dynamo process with a weak seed field is the most likely scenario (Emeriau-Viard & Brun 2017). However, some studies have shown that weak and strong dynamo branches may exist under certain initial conditions (weak or strong seed magnetic field Charbonneau 2004) or parameters such as the magnetic Prandtl number (Simitev & Busse 2009;Petitdemange 2018). Such weak or strong dynamo branches may explain some observed magnetic and rotational states seen in M dwarfs (Morin et al. 2011). Since this would depend on the local astrophysical context, we have decided to focus on the most common case of a weak seed magnetic field and refer the reader to these other complementary studies.
Having discussed how the kinetic and magnetic energies are distributed in our various models, we wish to go further in understanding exactly how these subtle balances come about. For this purpose we have computed the details of the energy transfers in our models, fo-cusing on the mean axisymmetric components MCKE, DRKE, PME and TME, since large-scale fields and flows are of key astrophysical interest.

Main transfer mechanisms between energy reservoirs
In this section we discuss the various energy transfers occurring in a rotating magnetized convective envelope. We refer the reader to Appendix A for the detailed derivation of the energy transfer equations, in which we have followed Starr & Gilman (1966); Rempel (2006), generalizing their derivation to global 3D spherical geometry. We focus here on the energy budget for the mean (axisymmetric) fields in the convective envelope of our models. We decompose energies into toroidal (along the azimuth) and poloidal (in the meridional plane) components. The budgets can be summarized as where all the different terms are detailed in Appendix A. We have computed individually each of the terms and show them normalized to the stellar luminosity in Fig. 20, as a function of the fluid Rossby number of the models. For each model, we have averaged the balances (17-20) over typically one hundred convective turnover time τ c such that the sum of the terms is close to zero. Cyclic cases show large variations of the energy balance (we return to this point hereafter), in these cases we averaged on a shorter time span chosen at cycle maximum. In addition, we have tabulated the transfers for three representative cases in Table 6 in units of both %L and %L . The differential rotation (upper left panel of Figure  20) is always sustained primarily by Reynolds stresses in the models (as discussed in §4.2), with a dominant contribution of the radial component v r v φ over the latitudinal component v θ v φ . The cases exhibiting anti-solar differential rotation (Rossby number larger than 1) present a reversal of the latter term, showing that the latitudinal component of the Reynolds stress is detrimental to the differential rotation kinetic energy in these cases. The magnetic contributions Q Ω (blue) and Q DR MS (red) start playing a significant role for fast rotating cases (low Rossby numbers, see model M07R5m in Table 6), sometimes even dominating completely viscous dissipation (Q DR ν , purple). In all cases the magnetic contributions tend to oppose differential rotation, as seen in §4.1. The power associated with the maintenance of differential rotation can reach about 30% of the stellar luminosity, and drops at minimum to about 4% in our sample of models. We remark that simulations with fluid Rossby numbers around Ro f ∼ 0.2 achieve the most powerful maintenance of differential rotation that can reach values up to 17% of the solar luminosity. At larger Rossby numbers, the star does not rotate fast enough and the differential rotation is weakly maintained. At lower Rossby numbers, the magnetic feedback from the dynamo field is so efficient that the power associated with the maintenance of differential rotation decreases significantly.
The meridional circulation energy balance (upper right panel of Figure 20) is dominated by a balance between the work of pressure (Q ∇P , cyan), buoyancy (Q b , blue-green) and Coriolis (Q c , green) forces (see also Table 6 where the dominant transfer terms are highlighted in bold font). The latter almost always remains negative, indicating an energy transfer from the meridional flow to the differential rotation when models are in a steady-state. Viscous dissipation (purple) plays a much lesser role for MCKE compared to DRKE, and magnetic contributions can be considered as negligible, except maybe for small Rossby number cases possessing trans-equatorial meridional cells (see Fig. 10). We find that the relative contribution of buoyancy and pressure gradients vary from model to model, and also vary in time for each model. We believe that is due to the anelastic approximation used in this study, and expect that a Lantz-Braginsky formulation (Brown et al. 2012) would lead to more systematic relative contributions of these two important terms for MCKE. Finally, we note that the power associated with the meridional circulation maintenance increases with Rossby number, and does not go above 15% of the stellar luminosity in our sample.
Let us now turn to the power sustaining magnetism in our models. The toroidal (TME) and poloidal (PME) magnetic energy budgets are shown in the left and right lower panels of Figure 20. We immediately note that the power sustaining magnetism corresponds at maximum to 3% of the stellar luminosity in our sample for TME. This corresponds to an absolute maximum of 6% of the solar luminosity. A very large amount of power is therefore indeed channeled to sustain the large toroidal magnetic energy reservoir that the dynamo builds up in the simulations. Hence, it is expected that a significant proportion of this large magnetic energy reservoir will be accessible to trigger various surface magnetic activity events (Shibata et al. 2013). The power associated with PME is a bit weaker, but still reaches up to 0.4% of the stellar luminosity. We find again that the most powerful transfers occur around R of ∼ 0.2. The power involved saturates for lower Rossby numbers, which is reminiscent of the saturation of magnetic activity observed in the X-ray luminosity of fast-rotating stars (e.g. Wright et al. 2011). It slowly drops for large Rossby numbers, but the power maintains a value of at least 0.01% of the stellar luminosity even in our most slowly rotating models. These figures are in good qualitative agreement with the value of 0.1% found for the Sun by Rempel (2006) using 2.5D mean field dynamo models. Let us stress again that with values ranging in our sample between 0.01% and 3% of the star's luminosity, this is a massive reservoir of magnetic energy extracted by dynamo action.
The poloidal magnetic energy balance is relatively straightforward: it is sustained primarily by the turbulent electromotive force originating from the convective motions (Q PM EMF , yellow) and opposed by Ohmic dissipa-  tion (purple). Mixed stresses involving the mean meridional flow (Q MC PM , salmon) are not observed to play any major role here. The toroidal magnetic energy balance is slightly more complex. In most of our models, it is primarily sustained by the Omega-effect (Q Ω , blue), and saturated by Ohmic dissipation (purple). Interestingly, we find that the role of the turbulent electromotive force can change from one model to the other (see Table 4), and it can even change sign with time in our cyclic solutions. This is highlighted in Fig. 21 where we observe how the various transfer terms for TME vary during one long cycle for model M09R3m in the left panel (TME is overplotted in black), and one short cycle for model M07R5m in the right panel. First, we observe that the amplitude of the transfers vary by an order of magnitude along the long cycle (left panel), being maximum when the magnetic energy is maximum as one may expect. We also see that electromotive force (yellow) plays a dominant role when TME increases right after cycle minimum, and then switches sign and draws energy from TME when TME decreases. This striking behavior is at odds with the classical picture of constant-in-time parametrization of mean-field coefficients. It furthermore supports our interpretation that the dynamo processes behind the decadal magnetic cycles observed in some models involve a complex interplay between sources and sinks of magnetic energy that vary at different stages of the cycle. This is important because it reinforces the conclusions drawn in §5 about the special nature of the long cycle period dynamo simulations presented in this study. We also see that the short cycle (right panel) behaves differently than the long cycle on the left. In model M07R5m, the electromotive force sometimes equates or even dominates the Ω−effect while still being balanced by Ohmic dissipation. In this case, the amplitude of the transfer terms vary much less with time, and we recover a behavior expected for α 2 − Ω dynamos. These simulations could therefore be categorized either as α − Ω or α 2 − Ω dynamos depending on the phases of evolution. We observe that the SVD analysis discussed in §5.1 and Appendix B shows coherent results when we take into account these temporal variations of the production terms, as shown in Fig. 21. Given the highly time-dependent nature of these nonlinear convective dynamo simulations, the analysis presented in this section about their dynamical properties is more robust than the SVD decomposition we performed in Appendix B as a companion analysis, since it does not assume any scale-separation approximation.

ASTROPHYSICAL IMPLICATIONS AND CONCLUSION
We have shown in the previous sections how various magnetic properties of solar-type stellar dynamo simulations change as a function of stellar mass and rotation. Often such variations can be understood using the Rossby number as a key control parameter. We here wish to reflect upon these findings and what are their astrophysical implications. There are several properties of solar-like stars such as their convective power and spectra, rotation profile, level of activity and presence of a magnetic cycle to cite only a few, that are of keen interest to be characterized. Our set of simulations can help us discuss some of these properties and provide clues to understand the physical mechanisms acting behind them.
Take for instance their interior rotation profile, we have seen in §4 that various states can be achieved in our set of simulations. We have further confirmed that such states depend on the Rossby numbers of the simulations. In  it was advocated, based on the hydrodynamic counterpart of the dynamo cases studied here, that three states of internal rotation could be found: solar-like (fast equator-slow poles), Jupiterlike (cylindrical profile with alternations of prograde and retrograde zonal jets) and anti-solar like (slow equator, fast poles). How is the presence of a self-sustained dynamo field changing this statement? We find that two states are retained: solar-like and anti-solar, and that the third one found for small Rossby numbers has been replaced by a new state. Indeed, we find that as the Rossby number decreases the feedback of the Lorentz force on the convective motion (via Maxwell stresses opposing Reynolds stresses in the angular momentum transport balance) yields smaller angular velocity contrast. This comes about because the rotation state tends towards uniform rotation (see §4). So for very small Rossby numbers, cases such as M11R5m or M09R5m are mostly showing a solid body rotation in their convective envelope, in sharp contrast with the banded profile of their hydrodynamics counterpart. However, the disappearance of cylindrical banded differential rotation state may be due to the range of Reynolds and magnetic Reynolds numbers considered in our study. The strong Lorentz force feedback may be due to our moderate state of turbulent convection. It is possible that at higher Reynolds numbers a cylindrical state would be retained even for a state near super-equipartition between kinetic and magnetic energy. This is a point that needs to be investigated further with a dedicated low Rossby/high Reynolds numbers study. Said differently: is there a level at which the magnetic energy contained in the convective envelope is so high that quasiuniform internal rotation is inevitable? We believe this is a reasonable assumption given the tendency of magnetic field to quench differential rotation as identified by many authors (Glatzmaier & Gilman 1982;Charbonneau 2004;Brun et al. 2005;Karak et al. 2015;Warnecke & Käpylä 2020, and references therein). So in summary, we find that the likely rotation states of solartype stars depend on their increasing Rossby number: quasi-uniform, banded/cylindrical, solar-like and antisolar. Such variations of the differential rotation states translated into an overall variation of surface angular velocity contrast being less sensitive to the bulk rotation rate, with ∆Ω ∝ Ω 0.46 , down from Ω 0.66 as in . We also find another potential interesting property for the differential rotation of solar-like stars: a scaling law may not be the best fit to our simulations database. As in Saar (2011), we find that there is a clear change of trend for small Rossby number (see Fig. 8). This is interestingly the change of rotation state from solar to almost uniform rotation. Determining for these various rotation states the exact transition in Rossby number will require more numerical study at higher levels of turbulence and continued dedicated observations. We intend to contribute to this effort with dedicated new simulations but also in preparing the scientific exploitation of PLATO (Rauer et al. 2014).
These various transitions of rotation profiles must impact the resulting dynamo and field properties. We have shown in the paper (sections 5 and 6) that this is indeed the case. Going from low to high Rossby number we find that dynamo action yields short cycle, long cycle and statistically steady (yet irregular) magnetic field evolution. This is very interesting because we can guide observations to search for these transitions in rotation state or temporal variability of the magnetic field. This will also help us discriminate between various dynamo scenarios.
Our set of dynamo solutions can help us characterize the mechanisms at work to generate and maintain magnetic fields for different sets of global stellar parameters. The rich range of magnetic phenomena occurring in stars rely on the free energy available in magnetic structures created by dynamo mechanism. In this study we have focused our analysis on a key aspect of the convective dynamo: energy transfers. We have done an extensive study on how the energy flows to and from the kinetic and magnetic energy reservoirs, separating them into their toroidal and poloidal components. The first key result is that a significant amount of the star's luminosity is being transferred into kinetic and magnetic energies. In Table 6 we listed as a function of the star luminosity (also with respect to the solar one) the amount of accessible energy. We demonstrated that typical numbers for the kinetic energy contained in the differential rotation are of order 10%, for the meridional circulation 1% of the star's luminosity. We also showed that for the toroidal magnetic energy, the energy available is also around 1% (with a maximum of 3%) and of the order of 0.1% for the poloidal magnetic energy. Having access to 1% of the star luminosity to power stellar magnetism via collective emergence of toroidal structures is significant. This means that there is large reservoir of magnetic energy accessible for the manifestation of various magnetic phenomena at the star's surface. We find for instance that our modelled stars can power dynamos such that they reach a global magnetic energy content from 10 37 to 10 39 erg. Part of this energy is found to be stored in the mean toroidal magnetic field (up to 6 × 10 38 erg), and the mean poloidal magnetic field is generally found to be much less energetic (reaching at most 4×10 37 erg). The corresponding total (unsigned) magnetic flux Φ tot is found to vary between 10 24 to 10 25 Mx over the range of mass and rotation covered by our study, thus very similar to observations of the Sun and other solar-type stars. In dynamo cases with long cycles such as case M09R3m, Φ tot is found to vary by a factor between 7 and 8 (see Figure 16) which is slightly more than what is found for the Sun (a factor of about 5 has been found for cycle 21 Schrijver & Harvey 1994).
We also found that Φ tot follows a scaling law with the Rossby number Φ tot ∼ R −0.88 of in qualitative agreement with observations (see Figure 17).
Another interesting finding of our study, which confirms results published in Augustson et al. (2015) with the same ASH code and in Strugarek et al. (2017Strugarek et al. ( , 2018 with the Eulag-MHD code, is the existence of a so-called nonlinear cyclic dynamo. Of course, convective dynamos are nonlinear in essence but what is meant here is that through the feedback of the Lorentz force on the flow, a cyclic behavior of the dynamo arises. Standard kinematic α − Ω mean field dynamos follow the Parker-Yoshimura (P-Y) rule (Parker 1955;Yoshimura 1975) and do not take into account nonlinear retroaction or do so in a limited way via the so-called Malkus-Proctor approach (Covas et al. 2005;Bushby 2006;Lopes et al. 2014, and references therein). By contrast, more and more 3D MHD convective dynamo simulations find that in a limited range of the parameter space, the P-Y rule does not apply anymore. This is the case in this study, where we find that for intermediate values of the Rossby number, typically 0.15 Ro f 0.65, the long cycle periods are due to a subtle interplay between the large-scale flow and the field. As the rotation rate is increased and the toroidal component of the dynamo generated magnetic field becomes more and more dominant via an efficient Ω-effect acting on the large scale poloidal field, the associated Lorentz force starts to quench the differential rotation via the action of Maxwell stresses opposing Reynolds stresses. This quenching of the differential rotation in turn implies that the Ω-effect is modified to the point that locally its latitudinal variation ∂Ω/∂θ reverses sign, leading to the generation of a toroidal field of opposite polarity, and through the action of turbulent convection, a reversed poloidal field. This nonlinear cyclic dynamo behavior is in sharp contrast with P-Y mechanism. Note that this is a delicate dynamo state to achieve, as the magnetic energy needs to be neither too weak nor too strong as discussed in Gilman (1983) (see for instance their Figure 31 or in Brun et al. (2005) where such a modulated dynamo state was also found in stellar core dynamos). To demonstrate that further, we have computed in Figure 12 the P-Y rule for one typical long cycle period dynamo case of our study and confirm that it is unable to explain the dynamo wave and cyclic behavior of this subset of dynamo cases (M09R3m, M11R3m for instance). However, we do find that for low Rossby number (Ro f < 0.42), the P-Y rule still works, and for instance in a case such as M07d5m also shown in Figure 12, we clearly have poleward dynamo waves compatible with the radial shear and the α-effect. Hence, we may have been able in this study to identify when P-Y vs nonlinear cyclic dynamos (in the sense defined in this study, e.g. feedback of the magnetic field on the local shear) take place. This is very important as it tells us how to reconcile various recent publications in the community that sometimes were finding that global convective dynamo could be interpreted as classical α-Ω dynamos (Warnecke 2018;Viviani et al. 2018Viviani et al. , 2019, and references therein), whereas others did not (Augustson et al. 2015;Strugarek et al. 2017Strugarek et al. , 2018. We propose that it is linked to different effective values of the Rossby number used in these various dynamo simulations. As we have seen above, it is instructive to make the link between full 3D MHD convective dynamo simulation and mean field dynamo concepts. Mean field dynamo theory usually uses the α-effect to parameterize turbulent magnetic field generation. In this study, we have estimated it through both the kinetic helicity (see §C and Pouquet et al. 1976) and an SVD decomposition (see §5.1 and Racine et al. 2011;Dubé & Charbonneau 2013;Augustson et al. 2015;Emeriau-Viard & Brun 2017). In the former case, we do not find a significant change of sign nor amplitude in the kinetic helicity of models possessing an anti-solar differential rotation. In the range of parameters considered in this study, this means that anti-solar-like stars need to be modeled with an α-effect similar to solar-like stars at least in their radial dependency, if not in amplitude. In the mean field α − Ω dynamo paradigm this implies that anti-solarlike stars will have a dynamo wave with a propagation reversed to that of the Sun, e.g. poleward from the equator to mid-latitudes as imposed by the P-Y rule. In our 3D simulations, we do not find such cyclic poleward dynamos for slowly rotating simulations, instead we find that they are statistically steady (but highly time dependent on short time scales). This is likely due to a less favorable phasing between poloidal and toroidal magnetic field generation in the convective envelope of these slowly rotating case that develops via complex nonlinear interactions between the fields and flows, which are not fine-tuned but instead evolves depending on the global parameters considered. Another interesting aspect is to assess how the dynamo-generated magnetic field is distributed over spatial scales. It is well known that there is a nonlinear feedback loop between rotation, dynamo, stellar wind and magnetic braking over secular time scales (Skumanich 1972;Brown 2014;Matt et al. 2015;Brun & Browning 2017;Metcalfe & van Saders 2017;Brun 2020;Vidotto 2021). It has been demonstrated that the magnetic torque provided by stellar winds are mostly controlled by the dipolar and quadrupolar modes (Réville et al. 2015;Garraffo et al. 2015;Finley & Matt 2018). Hence, one key question is to assess what happens with dipolar and quadrupolar modes when the dynamo changes its properties. To this end we showed in Figure 18 how magnetic geometry changes by computing quantities known as f dip and f quad . This allows us to assess the overall contribution of these two dynamo modes to the overall magnetic energy spectra. We found that they are key contributors to the overall magnetic energy with values ranging from 0.05 to 0.6, with most of the cases studied possessing f dip and f quad around 0.2 -0.3. We do not see any clear trend with Rossby number. Fast rotators and slows rotators both possess large dipolar and quadrupolar components. So from a stellar dynamo point of view it is difficult to invoke a drop in the largescale magnetic field to explain a possible break of stellar spin down for slow rotators as proposed by Metcalfe & van Saders (2017). Similar findings are obtained from observations of magnetic fields in cool stars as shown in Vidotto et al. (2016). The advocated Rossby number transition in magnetic field geometry to explain a collapse of magnetic breaking is thus unlikely. This study suggests that we must find a different explanation, maybe a less efficient heating mechanism inducing a sudden drop of coronal temperature and wind mass loss (Ó Fionnagáin & Vidotto 2018) which directly impacts angular momentum loss. Self-consistent rotating wind models with detailed treatment of the coronal heating mechanism are needed (see for instance Shoda et al. 2020;Hazra et al. 2021) in order to confirm the existence or not of such a transition in mass loss at slow rotation rates.
We have focused our analysis on the global energetics of the dynamo, and showed that the global dynamo field followed roughly a B bulk Ro −0.5 f trend (see §6.1) in agreement with previously published dynamo scaling laws (Augustson et al. 2019). It is also useful to interpret our simulations only considering the top of the dynamo domain, making a more direct link with stellar observations of surface magnetism. In this context,   Figure 22. Large-scale field at the surface of our modelled convective envelopes as a function of the Rossby number. The first panel shows the total dipole field, the second panel the large-scale fields (spherical harmonics < 5, including the nonaxisymmetric components (m = 0)). The third panel shows the ratio between the total root-mean-square (rms) field at the surface, and the equilibrium field based on the gas pressure at the photosphere. It can be considered as a measure of the filling factor f (see Cranmer & Saar 2011;See et al. 2019a). The symbols used in the panels are the same as in Fig. 7 we show in Fig. 22 the trend in Rossby number for the surface dipole field (first panel), surface large-scale field (second panel, see Table 5), and the ratio of the rootmean-square (rms) surface field to the equipartition field (third panel). The error-bars were deduced from the temporal variability of the fields, and the values are reported in Both trends are compatible with the trends deduced from Zeeman-Doppler Imaging surveys, that generally find the large-scale surface magnetic field to follow a Ro −1.3 f trend at intermediate Rossby numbers (See et al. 2019b). Finally, it is also instructive to assess the level of equipartition at the surface through the ratio between the surface rms field B rms and the equipartition field B eq (as defined in Johns-Krull & Valenti 2000) deduced here from the gas pressure at the surface of the stellar models we considered. Indeed, Cranmer & Saar (2011) have proposed that this ratio measures the filling-factor f of the large-scale field that shapes the lower stellar corona and ultimately determines the angular momentum loss rate of stars. See et al. (2019a) have found observationally that this ratio decreases with Rossby number. We find a similar trend here as seen on the third panel of Fig. 22, with f B rms /B eq 0.03 Ro −0.97±0.27 f . Finally, we note that the three magnetic field measures shown in Fig. 22 all exhibit an increase in amplitude at high Rossby number. This again strengthens the case that dynamo action within cool-stars does not exhibit any significant decrease of the large-scale magnetic field for slow rotators.
How are these results informing us about our star, the Sun? First, we note that the study of Strugarek et al. (2017Strugarek et al. ( , 2018 is about 1 solar mass stars and is taken into account in the analysis presented in this study. Given the good agreement seen in many of the plots discussed in §5 between the study done with the Eulag-MHD code and the one presented here with the ASH code (independently of models details), we are confident that the dynamo solutions discussed in this study are useful to understand the physical nature of the cyclic activity of a 1 solar mass star such as the Sun. Second, in this parametric stellar dynamo study we are proposing that in order to get both a solar-like conical differential rotation and a deep slow decadal-long magnetic cycles, the Rossby number must be between 0.15 and 0.65. Hence, we here acknowledge that cases M09R1m and M11R1m rotating at the solar rate do not show behaviors that are sun-like with respect to their magnetic activity (no cycles present) because their Rossby number is not falling in the 0.15-0.65 range. Instead, we believe that M09R3m or M11R3m are better, closer representations of the Sun even though their rotation rate is faster than the Sun, because their Rossby number is in the correct range of values. This means that while the overall trends found in our study are robust, the specific location of any given star must be thought with extreme care due to the so-called convective-conundrum, i.e a mismatch between global convection simulations and solar helioseismic inversion regarding the amplitude of giant convection cells (Hanasoge et al. 2016;Hotta & Kusano 2021). This is likely due to the fact that for any given rotation rate because of the convective-conundrum, the Rossby number achieved in the rotating convection simulation is slightly too large. So in order to be likely closer to the solar state and to aim for the correct value of the solar Rossby number, models rotating faster such as M09R3m or M11R3m cases are somewhat a better match to model the Sun than M09R1m or M11R1m. Thanks to this knowledge, we will next build a new global convective dynamo model of the Sun with an improved set of parameters by keeping the rotation rate to the solar one while controlling the effective Rossby number achieve in the simulation to be in the right range of values. We will report our finding in a future work.
To conclude, our study has confirmed the richness of dynamo solutions in parameter regimes that are likely to be found in solar-like stars and the large amount of magnetic energy and flux made available to the star and its surface activity by dynamo action. We have also identified the Rossby number regimes for different realization of differential rotation profiles and magnetic temporal modulations (cyclic or not), generalizing in an MHD context what we published in . Two key transitions in parameter space seem to be present, one at low Rossby number (Ro f < 0.1), another at high Rossby number (Ro f > 1). We need to study them with even more detail and at higher resolution and turbulence level to confirm the trends and scaling laws we have reported here. We intend to do so in the near future as well as study in more details the influence of a realistic atmosphere and of a wind ) on the dynamo properties.
We acknowledge partial financial supports by: DIM-ACAV+ ANAIS2 project, ERC Whole Sun Synergy grant #810218 and STARS2 starting grant #207430, ANR Toupies, INSU/PNST, Solar Orbiter and PLATO CNES funds. We thank GENCI via project 1623 for having provided the massive computing resources needed to perform this extensive study. J.V. acknowledge Project 2019-T1/AMB-13648 and multiannual agreement with UC3M ("Excelencia para el Profesorado Universitario" -EPUC3M14) -Fifth regional research plan 2016-2020 founded by the Comunidad de Madrid. This work benefited from discussions within the international team "The Solar and Stellar Wind Connection: Heating processes and angular momentum loss", supported by the International Space Science Institute (ISSI). A.S.B is thankful to B.P. Brown In this appendix we list the set of equations describing the energy transfer occurring in a star, focusing on mean energy quantities such as the poloidal and toroidal mean axisymmetric kinetic and magnetic energies. Following Starr & Gilman (1966); Brandenburg et al. (1996);De Rosa et al. (2002); Rempel (2006), we derive the set of equations of full energy transfers in spherical MHD configurations.
Let us denote the azimuthal average by a bar, and the derivation from it by a prime. For example, the radial velocity component will be written as v r = v r + v r . In order to characterize the axisymmetric magnetic (E m ) and kinetic (E k ) energy transfers between the various reservoirs of energy (thermal, potential, kinetic and magnetic) we will split E m and E k into three components: with DRKE and TME the mean axisymmetric toroidal energies, MCKE and PME the mean axisymmetric poloidal energies and FKE and FME the non-axisymmetric energies. To find the energy transfer equation for these various components we project the Navier-Stokes or induction equation onto the direction we wish to write the energy equation for, e.g. φ for TME for instance and inject the decomposition between mean and prime quantities. Then we perform an azimuthal average, thereby eliminating all terms that are linear in prime quantities. For each energy equations, we then multiply by a bar quantity (for instance B φ for TME) and rearrange the terms. For MCKE and PME, we combine the radial and latitudinal equations. Doing so systematically leads to the following set of equations 1 .

A.1. Overall Energy budgets
We follow the approach of Starr & Gilman (1966) and write the energy budgets in the following way (see the schematic in Fig. 23):  Figure 23. Global energy budget schematic. In red, we list all the key energy transport terms (see appendix B). Black arrows correspond to the direction of transport between the various energy reservoirs. Surface terms are indicated as black disks. We omit curvature terms to avoid crowding the figure.
In all that follows, quantities are separated into mean and fluctuating components through and the corresponding terms in the original derivation of Starr & Gilman (1966) are given by the blue 'SG66: [XX]' labels at the end of each equation, where XX is the term or equation number in Starr & Gilman (1966). Note that we have extra curvature terms C x due to our choice of spherical coordinates.

A.4. Axisymmetric toroidal magnetic equation (TME)
Q Ω and Q MC TM were defined previously in A9 and A16. The remaining terms in Eq. A5 are A.5. Axisymmetric poloidal magnetic equation (PME) Q PM MC have already been defined in Eq. A19. The remaining terms in Eq. A6 are

B. MEAN FIELD SVD DECOMPOSITION OF DYNAMO SOLUTION
It is instructive to compare our 3D simulation results with the concepts used in mean field dynamo theory (see §5). For instance, the generation of poloidal magnetic field in the simulation is dominated by the action of the fluctuating EMF: E FI = E = v × B . This process can also be interpreted through the α-effect approximation, which is a first order expansion of E around the mean magnetic field and its gradient: with α ij a rank-two pseudo-vector and β ijk a rank-three tensor. In the following, we will neglect the β term. However, this will increase the systematic error when estimating the α term. Thus, a single-value decomposition (SVD) including the β-effect has been calculated in order to provide a lower-bound on the systematic error as discussed in Augustson et al. (2015). In the following analysis, α has been decomposed into its symmetric and antisymmetric components Thanks to the SVD decomposition we can quantify the relative efficiency of the α-effect in generating the mean magnetic field and characterize the type of dynamo through the relative influence of its regenerating terms. We can start by evaluating how the convective flows regenerate mean magnetic fields. This can be determined by finding the amplitude of an estimated α-effect relative to the rms value of the non-axisymmetric velocity field E α v rms = 3 2(r 3 top − r 3 bcz ) × i,j drdθr 2 sin θ α ij α ij {v · v } (B35) and E γx = γ x E . By calculating this matrix, see Table 7, we notice that for the antisymmetric part γ, the predominant term is γ ϕ that impacts the poloidal component of the magnetic field. Only for M07R5m are the three components of the same order of magnitude. In the three other cases shown, γ r and γ θ have roughly the same order of magnitude and are smaller by a factor 2 to 3 compared to γ ϕ . By looking at the symmetric part α S , we see the same trend. The predominant term is α (rr) with α (rθ) and α (θθ) close second. They all act on the poloidal component of the magnetic field. The smallest term is in most cases α (ϕϕ) which is at least 5 times smaller than the predominant term except once more in case M07R5m where it is of the same order of magnitude. The sum of all α terms varies between 51% in case M07R5m up to 73% in case M09R3m. Hence, the γ terms (the antisymmetric part of the alpha-tensor) account for 49% in case M07R5m down to 27% in case M09R3m.
In order to better quantify this relative influence we can compute the α P /α ϕ ratio: Looking at Table 7 where we report the value of this ratio for all 4 representative models, we note the predominance of the poloidal field regeneration over the toroidal field regeneration for all models as the ratio α P /α ϕ is always above 1. This ratio varies from 1.59 in M09R3m up to 12.4 in case M07R5m.
Turning now to the regeneration of the toroidal field, we know from mean-field dynamo theory that it can be due to either the α-effect, coming from the fluctuating emf E , or from the Ω effect that acts on the poloidal field through differential rotation. In all our models, we note that the regeneration of B ϕ by the α-effect is small, compared to the one of B pol . Therefore, we now want to measure the relative influence of the Ω-effect to that of the α-effect, since the toroidal magnetic field can be regenerated through both effects: Ω α ϕ = 3 2(r 3 top − r 3 bcz ) × drdθr 2 sin θ r sin θ B ϕ B P · ∇ Ω B ϕ φ · ∇ × E .
We note that in all models the Ω-effect is much stronger than the α-effect in generating the toroidal magnetic field (the ratio Ω/α ϕ is greater than 5), except for case M09R1m for which it is closer to 1. This confirms that most of the dynamo models considered in this study can be classified as α-Ω dynamos rather than α 2 -Ω. Statistically steady simulations such as M09R1m on the contrary are closer to be classified as α 2 -Ω. Of course, this mean field dynamo classification is mostly useful for short magnetic cycle period cases (illustrated in the table with case M07R5m) as they also follow Parker-Yoshimura rule (see §5). For long magnetic cycle period cases such as M09R3m and M11R3m this is less significant, as we observe a complex nonlinear feedback that leads to a different type of cyclic dynamo. Further, we have shown in section 6 and Fig 21 that these dynamo mechanisms are highly variable in time, and can sometimes be quenched while at other times they become dominant. Hence, a mean field classification on such solutions could vary depending on the dynamo phase considered.

C. KINETIC HELICITY IN SOLAR AND ANTI SOLAR CASES
In Figure 24 we display several realizations of the horizontally-averaged radial profile of the kinetic helicity H k = v·ω in our set of convective dynamo models. These profiles have been averaged over the northern hemisphere only. Note: Results of the mean field SVD dynamo analysis on four representative models (M07R5m, M09R1m, M09R3m, M11R3m) ordered from top to bottom in increasing Rossby number values. The first column represents the α tensor with its symmetric: αs and antisymmetric: γ (italic) portions (see Eq B33). The middle column gives the relative importance of the Ω-effect to the α-effect for the regeneration of the toroidal field. The last column quantifies the ratio of the α-effect used for the regeneration of the poloidal magnetic field to the one used for the regeneration of the toroidal field. On the left panel we display the kinetic helicity profiles for the M05m series. We first note that the kinetic helicity is negative in most of the domain and changes sign at the bottom of the convective envelope and is close to zero in the deep radiative interior below. This sign reversal of H k is understood by the change of sign of the vorticity field in the downward plumes. As they splash onto the top of the radiative zone (whose realistic stiffness we recall is directly taken