<figure><img src="image_42.png" /></figure>

FIGURE 1.33 – Évolution temporelle du taux de rotation de l'enveloppe convective (traits pleins) et du cœur radiatif (pointillés) pour des étoiles de 1 $M_\odot$. Les différentes couleurs représentent des rotateurs lents (orange), intermédiaires (bleu clair) et rapides (bleu foncé), par rapport à leur rotation initiale $P_{\text{rot,ini}} = 1.4\,5$ et 8 jours. Les points avec les barres d'erreurs correspondent aux $25\%$, $50\%$ et $90\%$ percentiles des distributions de rotation au sein d'amas stellaires observés par Gallet & Bouvier (2015). La loi de Skumanich est représenté par la ligne noire, normalisée par rapport au Soleil (⊙). Adapté de Ahuir et al. (2021).

deux couches au sein des étoiles sont alors ajustés pour correspondre aux observations, et dépendent de la masse des étoiles considérées. On voit alors apparaître différentes phases rotationnelles durant la vie d'une étoile. Lors de sa formation, l'étoile subit une phase d'effondrement tout au long de la PMS. Cependant, sa vitesse de rotation reste d'abord constante pendant quelques millions d'années, ici dû à l'hypothèse d'interaction avec le disque proto-stellaire environnant (voir la revue de Zhao et al. 2020). Une fois le disque dissipé par accrétion, éjection et formation de planètes, l'étoile accélère en se contractant par conservation de son moment cinétique, jusqu'à atteindre la ZAMS. Une fois rentrée sur sa séquence principale, la rotation ralentit en raison de la perte de moment cinétique par l'intermédiaire du vent stellaire magnétisé, comme nous l'avons vu à la section précédente (voir équation 1.44). La variation de cette rotation va alors influencer le mécanisme dynamo à l'origine du champ magnétique de l'étoile, qui va lui-même influencer la structuration du vent (Réville et al., 2015; Finley & Matt, 2018). Le vent est ensuite à l'origine de la perte de masse, qui va donc influencer l'évolution de la rotation. Cette boucle de processus physiques non-linéaires fait alors converger les étoiles de masses similaires vers un comportement commun durant la séquence principale. Pour comprendre cela, on peut considérer par approximation qu'une étoile en rotation plus rapide (lente) aura tendance à créer un champ magnétique plus (moins) intenses, produisant alors un freinage par vent magnétisé plus (moins) fort. Ce comportement commun avait initialement été observé par Skumanich (1972), donnant son nom à la loi empirique du même nom, indiquant que $\Omega_* \propto t^{-1/2}$. Il est ainsi devenu possible d'avoir une idée de l'âge d'une étoile en fonction de sa masse et de sa période de rotation, ce que l'on appelle la gyrochronologie (Barnes, 2003, 2007; Meibom et al., 2015).

On peut noter ici que ces résultats ont été obtenus en ne considérant que la contribution continue du vent magnétisé pour la perte de masse. Cette dernière peut cependant se faire aussi de manière intermittente par le
<figure><img src="image_43.png" /></figure>

FIGURE 1.34 – Gauche : Évolution du champ magnétique moyen à grandes échelles non signé $\langle|B_V|\rangle$ dérivé par ZDI en fonction de l’âge de l’étoile. Crédits : Vidotto et al. (2014). Droite : Plot en densité du traceur chromosphérique R'IST de GAIA (similaire au traceur R_HK) en fonction de la température effective. Les cibles dont la gravité de surface log₁₀g< 4.0 ont été ici tronquée pour mettre en évidence le comportement sur la séquence principale. Crédits : Lanzafame et al. (2022).

biais de phénomènes éruptifs de surface, telles que les éjections de masse coronales (voir par exemple Aarnio et al. 2012; Drake et al. 2013). Il est possible que les événements transitoires jouent un rôle plus important que le vent dans la perte de masse de jeunes étoiles très actives (Osten & Wolk, 2015; Jardine et al., 2020), cependant il apparaît que le vent reste la contribution majoritaire pour le cas des étoiles de type solaire (Howard et al., 1985; Güdel, 2019).

Ainsi, on a pu voir que le champ magnétique joue un rôle important pour la perte de masse des étoiles de type solaire, et est amené à évoluer durant la vie d'une étoile. Nous allons maintenant voir les caractéristiques de cette évolution.

### 1.3.3 Évolution du magnétisme

Par l'intermédiaire des différentes techniques de détection mentionnées en Section 1.2.3.3, il est aussi possible d'étudier l'évolution du magnétisme durant la vie des étoiles de type solaire, en observant des étoiles à différents âges. Vidotto et al. (2014) ont ainsi relevé les champs magnétiques de 73 étoiles de 0.1 à 2 M☉, dont les âges varient entre 1 million et 10 milliards d'années. Leurs mesures ZDI du coefficient Stokes V de polarisation circulaire permettent de remonter aux composantes à grandes échelles du champ magnétique, et montrent que ce dernier est décroissant en fonction de l'âge, tel que $\langle|B_V|\rangle \propto t^{-0.655\pm0.045}$, illustré sur le panneau gauche de la Figure 1.34 (voir aussi Folsom et al. 2018). Cette relation similaire à ce qui est observé pour la rotation permet alors une détermination similaire de l'âge des étoiles, cette-fois à partir de leur champ magnétique, ce qui est appelé *magnétochronologie*. Cependant, la détermination des âges pour les cibles de la calibration initiale souffre pour l'instant d'incertitudes rendant son utilisation difficile (Lebreton et al., 2014). Les récentes données de la mission Gaia (Lanzafame et al., 2022) ainsi que la future mission PLATO (Rauer et al., 2014) devraient en ce sens apporter de nouvelles contraintes pour améliorer cette précision.

L'origine du champ magnétique stellaire provient du mécanisme dynamo créé par l'action conjointe des mouvements convectifs et de la rotation (voir Section 1.1.6). Pour le premier mécanisme, il est possible de déterminer le temps de retournement convectif $\tau_c$ à partir de la température effective (voir par exemple Cranmer & Saar 2011). Aucune relation claire entre le magnétisme et la température effective n'a pour l'instant été relevée, bien que les récentes données de GAIA semble montrer différents comportements pour l'activité chromosphérique
<figure><img src="image_44.png" /></figure>

FIGURE 1.35 – Gauche : Champ magnétique moyen à grandes échelles non signé $\langle B_V \rangle$ en fonction du nombre de Rossby stellaire Ros, où la masse est représentée par la barre de couleurs. L'intervalle magenta représente la variation de la valeur solaire durant le cycle 24. Crédits : See et al. (2019a). Droite : Temps caractéristiques de freinage dans les régimes saturé (ligne noire inférieure) des rotations rapides, et non saturé (ligne noire supérieure) des rotations lentes, en fonction de la masse (Matt et al., 2015). Les lignes pointillées montrent les temps similaires pour les modèles de Barnes (2010) (B10, en vert), van Saders & Pinsonneault (2013) (vSP13, en bleu), Reiners & Mohanty (2012) (RM12, en rouge), ainsi qu'une version modifiée de Kawaler (1988). Crédits : Matt et al. (2015)

maximum qui peut être atteinte pour une température (voir panneau droit de la Figure 1.35). Un comportement similaire est observé pour le maximum d'activité avec la période de rotation, sans avoir non plus une tendance clairement définie (Santos et al., 2021). Cependant, lorsque convection et rotation sont combinées par l'intermédiaire du nombre de Rossby (voir Section 1.2.2.4), des tendances claires apparaissent (Folsom et al., 2016; Lehtinen et al., 2020) et sont illustrées sur le panneau de gauche de la Figure 1.35 (See et al., 2019a). On observe alors deux régimes pour le champ magnétique à grande échelle en fonction du nombre de Rossby stellaire, avec premièrement les jeunes étoiles en rotation rapide possédant une composante magnétique $\langle B_V \rangle \sim 250$ G invariante, dans un régime dit saturé. Au-delà d'un Rossby de transition Ros 0.06, la tendance du champ est décroissante dans un régime dit non-saturé, telle que $\langle B_V \rangle \propto R_o^{-1.40 \pm 0.10}$ (voir aussi Figure 1.28). Ces différents régimes en fonction du nombre de Rossby ont aussi été relevés dans les émissions chromosphériques (Noyes et al., 1984; Lehtinen et al., 2021), les émissions en rayons X (Pizzolato et al., 2003; Wright et al., 2011), ou encore le flux magnétique caractérisant le champ magnétique à petites échelles (Reiners et al., 2009, 2022). En fonction du régime considéré, les modèles d'évolution stellaire montrent que l'impact du champ magnétique sur la perte de moment cinétique va être modifiée (Barnes, 2010; van Saders & Pinsonneault, 2013; Matt et al., 2015), inversant alors la dépendance des temps caractéristiques de freinage par le vent magnétisé en fonction de la masse (voir panneau droit de la Figure 1.35). Les étoiles de faibles masses restent donc des rotateurs rapides plus longtemps que les étoiles plus massives, tandis qu'une fois le régime de rotations lentes atteint, elles vont ralentir plus efficacement au point de passer en rotation plus lente que les étoiles plus massives de même âge.

Ainsi, le champ magnétique tend à diminuer durant la séquence principale, tout comme la rotation de l'étoile, qui participe à son origine. Grâce à la caractérisation de cette évolution, les modèles actuels tendent à relativement bien expliquer les observations des étoiles de type solaire sur la première partie de leur séquence principale. Cependant, l'évolution magnéto-rotationnelle des vieilles étoiles de type solaire est actuellement sujet à controverses et discussions, comme nous allons le voir dans la Section suivante.
### 1.3.4 Changement de comportement en fin de séquence principale

Avec l'augmentation récente du nombre d'étoiles pour lesquelles des périodes de rotation sont mesurées (Santos et al., 2019; Santos et al., 2021; Reinhold & Hekker, 2020; Claytor et al., 2022), il a été observé que la perte de moment cinétique des étoiles sur la séquence principale subit des époques de perte en efficacité (van Saders et al., 2016) ou même d'arrêt (Curtis et al., 2020; David et al., 2022). Les évolutions rotationnelles s'écartent alors de la loi classique de Skumanich, impliquant un besoin de révisions pour les techniques telles que la gyrochronologie.

En particulier, les observations de van Saders et al. (2016) ont montré que les analogues solaires plus vieux que le Soleil possèdent un taux de rotation plus rapide que celui attendu par les prédictions de Skumanich (1972), comme illustré au panneau droit de la Figure 1.36. Les mesures de rotations obtenues principalement par photométrie (et astérosismologie pour certains analogues solaires) pointent alors vers une baisse significative de la perte de moment cinétique, après avoir passé un certain âge. Ce changement de comportement semble apparaître vers 2 à 3 milliards d'années pour les étoiles de types F, 4 à 5 milliards d'années pour les étoiles de type G (soit environ l'âge du Soleil), et 6 à 7 milliards d'années pour les étoiles de type K (Metcalfe & van Saders, 2017). Cette dépendance évolutive à la température effective pointe alors vers une dépendance au nombre de Rossby, pour lequel une déviation à la loi de Skumanich s'enclencherait une fois qu'une valeur critique de $Ro_s \sim 2$ est franchie. Enfin, la possibilité d'un biais d'observation a également été proposée, car les étoiles à rotation lente pourraient perdre leurs taches de surface, rendant la détection des modulations rotationnelles stellaires relativement difficile avec les techniques photométriques. Cette question a longuement été débattue depuis, notamment car l'étude de (Lorenzo-Oliveira et al., 2019) n'observe pas ce signe de freinage magnétique affaibli au sein de leur échantillon de jumeaux solaires. Cependant, l'étude de Hall et al. (2021) semble confirmer l'existence de cet affaiblissement du freinage en utilisant seulement des périodes astérosismiques, qui ne sont donc à priori pas sensible à un tel biais. De plus, des discussions récentes de travaux à paraître prochainement (Metcalfe, private communication) semblerait pointer vers le fait que l'échantillon de Lorenzo-Oliveira et al. (2019) serait finalement cohérent avec un scénario de freinage moins efficace. Ces récents éléments apportés au débat pointent donc vers une probabilité non négligeable que l'écart à la gyrochronologie observé viendrait à être confirmé et d'origine physique.

Comme nous l'avons vu à l'équation 1.44, la perte de moment cinétique $\dot{J}$ est principalement influencée par le taux de rotation, le magnétisme ainsi que la perte de masse. Pour interpréter un tel changement de régime du freinage aux âges tardifs, deux scénarios sont alors envisagés. Le premier est un changement dans le processus dynamo des étoiles en rotation lente, comme le suggèrent des observations récentes (Metcalfe & Egeland, 2019). En effet, les modèles montrent que le rayon d'Alfvén $r_A$, dont dépend le couple magnétique fourni par les vents, dépend principalement de la topologie magnétique à grande échelle, et plus particulièrement de l'amplitude de la composante dipolaire (Réville et al., 2015; Finley & Matt, 2018). Cette hypothèse est actuellement la plus soutenue puisqu'une dépendance de $\dot{J}$ au carré de $r_A$ rend la chute de freinage d'autant plus sensible à la chute de la topologie à grande échelle. Cependant, la caractérisation observationnelle des topologies magnétiques stellaires s'avère difficile pour les régimes de rotation lente (Donati et al., 2006; See et al., 2019b). Nous reviendrons sur ce scénario en Section 6.1.2.

Le deuxième scénario envisagé concernerait cette fois-ci un changement de régime des propriétés de chauffage de la couronne, impliquant alors une chute du taux de perte de masse $\dot{M}$ qui rendrait à son tour le vent moins efficace pour extraire le moment angulaire. En reprenant les observations de température coronale $T_c$ (Johnstone & Güdel (2015) représentées au panneau droit de la Figure 1.36, Ó Fionnagáin & Vidotto (2018) ont calculé la perte de masse résultante, à partir d'un modèle polytropique 1D. En se basant sur un changement de régime suggéré par les observations, passant de $T_c \propto \Omega_*^{0.38}$ à $T_c \propto \Omega_*^{1.14}$ (courbes pointillées rouges), l'évolution de la perte de masse chute alors en passant de $\dot{M} \propto t^{-0.74}$ à $\dot{M} \propto t^{-3.89}$. Il est intéressant de noter que dans ce cas, le Soleil aurait déjà passé cette transition de régime, apparaissant pour $1.4 \Omega_*$. Cependant,
<figure><img src="image_45.png" /></figure>

FIGURE 1.36 – *Gauche* : Évolution temporelle de la période de rotation. La ligne pleine montre un modèle d'évolution rotationnelle standard en accord avec la loi de Skumanich (1972). La ligne en pointillées montre le modèle modifié de van Saders et al. (2016), visant à reproduire l'échantillon présenté. La zone grise représente la nouvelle tendance une fois les barres d'erreurs prises en compte. Les points jaunes représentent des analogues solaires bien caractérisés. Crédits : Metcalfe & van Saders (2017). *Droite* : Évolution de la température coronale en fonction du taux de rotation Ω*. Les points indiquent les observations de Johnstone & Güdel (2015), et les lignes pointillées représentent des régressions linéaires faites sur différents intervalles de ces points. La zone grise représente le régime où le comportement a changé durant l'évolution, ici à partir de 1.4 Ω⊙. Crédits : Ó Fionnagáin & Vidotto (2018)

nous manquons encore de contraintes observationnelles et théoriques sur le taux de perte de masse stellaire des étoiles de la séquence principale, ce qui permettrait de confirmer ou d'infirmer une telle option. De plus, la courbe en pointillés bleus sur la Figure 1.36 montre une régression linéaire faite sur l'ensemble des points, en considérant qu'un changement de régime n'aurait pas lieu. La qualité de la régression est alors moins bonne, mais le faible nombre de données ne permet pas d'exclure totalement cette hypothèse.

Dans ce contexte, il est crucial de comprendre ce qui se passe exactement pour la perte de moment angulaire des étoiles en rotation lente. De récentes observations suggèrent que les deux scénarios discutés pourraient agir conjointement afin d'affaiblir significativement le freinage (Metcalfe et al., 2022). Dans les Chapitres 5 et 6, nous tenterons d'apporter des contraintes théoriques supplémentaires à partir d'une étude numérique multi-dimensionnelle, afin d'explorer la possibilité d'une transition de tendance magnétique pour les rotateurs lents. En particulier, nous nous demandons si l'affaiblissement du freinage suggéré pourrait résulter d'une transition vers un profil de rotation différentielle particulier, apparaissant dans les modèles globaux de convection stellaire en rotation lente : le profil de rotation *anti-solaire*.
<figure><img src="image_46.png" /></figure>
# Deuxième partie

Dynamique des étoiles de type solaire :
Convection turbulente, écoulements à grandes échelles et implications
observationnelles
<figure><img src="image_47.png" /></figure>
# CHAPITRE 2

## ÉCOULEMENTS GRANDES ÉCHELLES DANS LES ÉTOILES DE TYPE G ET K

### Sommaire

<table><tr><td>2.1</td><td>Contraintes observationnelles et avancées des modèles</td><td>62</td></tr><tr><td>2.1.1</td><td>Contexte observationnel</td><td>62</td></tr><tr><td>2.1.2</td><td>L'apport des simulations numériques : la rotation différentielle antisolaire</td><td>64</td></tr><tr><td>2.2</td><td>Simulations magnétohydrodynamiques globales de dynamos stellaires : Description de l'outil et de l'expérience numérique</td><td>65</td></tr><tr><td>2.2.1</td><td>Description des équations résolues</td><td>66</td></tr><tr><td>2.2.2</td><td>Conditions initiales et conditions limites</td><td>68</td></tr><tr><td>2.2.3</td><td>Espace des paramètres exploré</td><td>69</td></tr><tr><td>2.2.4</td><td>Illustration d'une dynamo convective 3D</td><td>72</td></tr><tr><td>2.3</td><td>Caractérisation des écoulements à grande échelle</td><td>72</td></tr><tr><td>2.3.1</td><td>Analyse globale</td><td>74</td></tr><tr><td>2.3.2</td><td>Différents régimes et tendances de la rotation différentielle</td><td>75</td></tr><tr><td>2.3.3</td><td>Circulation Méridienne</td><td>78</td></tr><tr><td>2.4</td><td>Transferts de moment cinétique et établissement du profil de rotation différentielle</td><td>81</td></tr><tr><td>2.5</td><td>Bilan énergétique et intensité des écoulements moyens</td><td>85</td></tr><tr><td>2.6</td><td>Discussion et perspectives</td><td>89</td></tr><tr><td>2.7</td><td>Article publié : Powering Stellar Magnetism : Energy Transfers in Cyclic Dynamos of Sun-like Stars</td><td>89</td></tr></table>

Nous avons vu en Section 1.2.1 qu'il existe une grande diversité de structures stellaires. Dans le cas des étoiles de type solaire, la profondeur relative de l'enveloppe convective diminue pour une masse considérée plus grande. Nous avons ensuite vu en Section 1.3.2 que chacune de ces étoiles rencontre différentes phases évolutives et rotationnelles au cours sa vie. En particulier, il semblerait que le taux de rotation et l'amplitude du champ magnétiques diminuent durant la séquence principale. On peut alors se demander quels impacts auront ces changements sur la dynamique interne d'une étoile.

Pour y répondre, il faut étudier différentes étoiles à différents instants de leur vie. Cependant, la quasi-totalité des étoiles est bien trop éloignée pour que nous puissions observer directement leur dynamique. Grâce à sa proximité, le Soleil offre alors une formidable opportunité pour comprendre les mécanismes physiques à l'œuvre au sein de ces astres. Mais bien que l'héliosismologie nous ait permis de sonder et de confirmer sa structure, il reste certaines incertitudes, comme la forme de la circulation méridienne (voir Section 1.2.2.2). De plus, nous assistons à un instant de la vie du Soleil bien court en comparaison des 10 milliards d'années de sa séquence principale, tout en sachant qu'il est un représentant particulier pour une masse donnée. Ainsi, il est important de considérer le Soleil dans son contexte stellaire et de développer des outils pour entrevoir les mécanismes d'autres étoiles, afin de mieux comprendre ceux à l'œuvre au sein de notre étoile et vice-versa.

Dans ce chapitre, nous présenterons en Section 2.1 le contexte actuel pour l'étude des écoulements grandes échelles dans les étoiles de type solaire. Nous décrirons ensuite en Section 2.2 l'outil numérique nous permettant de modéliser la dynamique de 15 étoiles de type solaire, à différents taux de rotation et pour des
<figure><img src="image_48.png" /><figcaption>FIGURE 2.1 – Profils de la rotation différentielle de HD 173701 (panneaux du haut) et HD 187160 (panneaux du bas) obtenus par inversions sismiques. Des coupes méridiennes de cette rotation sont proposées à gauche avec l'emplacement le plus probable de la tachocline, et des profils latitudinaux de surface sont proposés en panneaux de droite. Sur ces derniers, la couleur correspond à la densité de probabilité, soulignée par l'interval à 1σ en trais pointillés. Crédits : Benomar et al. (2018).</figcaption></figure>

masses variant entre 0.5 et 1.1 $M_\odot$. La Section 2.3 sera dédiée à la caractérisation des écoulements grandes échelles (rotation différentielle, circulation méridienne) prenant place dans ces milieux convectifs, turbulents, magnétisés et en rotation. Enfin, pour comprendre la mise en place de ces profils, nous quantifierons en Section 2.4 les transferts de moment cinétique au sein des enveloppes convectives, avant de finalement analyser en Section 2.5 les transferts d'énergie soutenant ces écoulements axisymétriques.

## 2.1 Contraintes observationnelles et avancées des modèles

### 2.1.1 Contexte observationnel

Dans la Section 1.2.1, nous avons présenté la structure du Soleil ainsi que les différents écoulements à grandes échelles qui ont été détectés en son sein. Plus précisément, nous avons vu que son enveloppe convective extérieure est en rotation différentielle (DR, Figure 1.16) et est le siège d'une circulation méridienne (MC, Figure 1.17). Ces résultats ont en grande partie été trouvés par l'observation de sa surface, et en particulier en étudiant les vibrations de cette dernière, nous renseignant alors directement sur la propagation d'ondes au sein du Soleil (hélio-sismologie). Qu'en est-il des autres étoiles ?

Si l'observation directe de la surface des autres étoiles est quasi-impossible à cause de la distance nous en séparant, il est cependant possible de retrouver des traces de ces vibrations de surface dans leur courbe de lumière (variation temporelle de luminosité), et ainsi de sonder indirectement l'intérieur de ces étoiles (astero-sismologie : voir García & Ballot 2019 pour une revue récente). Cependant, la mise en place de telles techniques reste difficile et réservée à certaines étoiles observées de manière prolongée. Ainsi, il devient possible depuis maintenant 10 ans de sonder l'intérieur profond d'étoiles évoluées telles que des géantes rouges (Beck et al., 2012; Deheuvels et al., 2012), et plus récemment des sous-géantes (Deheuvels et al., 2020). Des profils radiaux de rotation différentielle ont pu aussi être reconstruits sur des étoiles de type solaire dont le signal sismique est assez fort (Benomar et al., 2018; Bazot et al., 2019).
<figure><img src="image_49.png" /></figure>

FIGURE 2.2 – Tendances observationnelles de la rotation différentielle, représentée ici comme fonction de $\Omega_*$ sur le panneau de gauche, et $T_{\text{eff}}$ sur le panneau de droite. Les données proviennent de différentes études H95, DSB96, RS03, B05 et CC07, explicitées dans le corps du texte. Les flèches orientées vers le haut signifient que le point concerné ne possède pas de limite haute (et vice-versa pour les limites basses). Crédits : Barnes et al. (2005) et Collier Cameron (2007).

Nous illustrons en Figure 2.1 les premières reconstructions de rotations différentielles dans d'autres étoiles de type solaire. Elle montre des résultats représentatifs de l'étude Benomar et al. (2018) pour les étoiles HD 173701 et HD 187160. On y voit une estimation du profil latitudinal de la DR (à gauche), ainsi que la position la plus probable de la tachocline. Le profil latitudinal en surface de l'étoile est détaillé avec une densité de probabilité. Dans cette étude, Benomar et al. (2018) ont analysé 40 étoiles de type solaire à partir des données du satellite Kepler. Il détecte clairement un profil de rotation différentielle de type solaire (équateur plus rapide que les pôles) pour 13 d'entre eux, cependant l'écart de rotation entre l'équateur et les pôles est plus fort, en moyenne de 60%, contre seulement 30% dans le Soleil. Bazot et al. (2019) détectera par la suite des DR plus proches du régime solaire au sein de deux analogues 16 Cyg A and B. Cette différence dans le contraste latitudinal $\Delta\Omega$ peut venir du fait que la majorité des étoiles de Benomar et al. (2018) tournent plus rapidement que le Soleil. Dans ce cadre, il est important de comprendre les différentes tendances de la DR en fonction de différents paramètres stellaires, tels que $M_*, \Omega_*, Z$ ou encore l'âge.

Premièrement, il a été observé que ce contraste rotationnel peut dépendre du taux de rotation absolu, sous la forme d'une loi de puissance $\Delta\Omega \propto \Omega_*^n$. Ainsi, différentes observations sont illustrées sur le panneau gauche de la Figure 2.2, tiré de Barnes et al. (2005) (B05). Lors de suivis photométriques sur le long terme menés par Henry et al. (1995) (H95), des variations temporelles de la courbe de lumière ont pu être reliées à la présence et rotation de starspots (taches stellaires, similaires aux taches solaires illustrées en Figure 1.21). Différentes périodes de rotation peuvent être observées au sein d'une même courbe de lumière, indiquant alors la présence potentielle de cisaillement rotationnel. En effet, il est possible d'attribuer les différentes périodes à différentes taches localisées à différentes latitudes, témoignant ainsi de la DR $\Delta\Omega$ de surface. Leur étude tend à montrer que $n = 0.24 \pm 0.06$, soit une dépendance faible de ce cisaillement de surface à la rotation de l'étoile. Par la suite, des observations similaires ont pu être faites dans les bandes H & K d'émission du Ca II par Donahue et al. (1996) (DSB96), mais pointant une tendance significativement plus élevée, telle que $\Delta\Omega \propto \Omega_*^{0.7\pm0.1}$. Cette valeur a par la suite été consolidée par Reiners & Schmitt (2003) (RS03), ayant trouvé $n = 0.66 \pm 0.26$ par l'intermédiaire d'une méthode de déconvolution spectrale<fnref n="1" />. Il est important de noter que ces mesures de cisaillement sont des limites basses. Les mesures ne donnent en effet que très peu d'indication sur la latitude de la tache à laquelle la modulation est rattachée. Les $\Delta\Omega$ ainsi déduits sont donc seulement rattachés à la gamme de latitudes explorées par les régions actives lors du cycle d'activité stellaire (voir Section 1.2.2 pour le cas du

<footnote n="1">1. Plus exactement, Least Square Deconvolution</footnote>
Soleil). Un fit global de l'ensemble des données citées dans ce paragraphe amène B05 à *n* = 0.15 ± 0.1, ce qui montre que cet indice est très sensible au degré de diversité spectrale de l'échantillon de cibles considéré, comme pressenti par RS03. En effet, Balona & Abedigamba (2016) montrent que cet indice pourrait dépendre du type spectral de l'étoile, comme suggéré par B05 qui relevait des ΔΩ plus forts pour des types spectraux plus chauds (ici G et F). Finalement, les études basées sur l'astéro-sismologie trouvent quant à elle une valeur pour *n* comprise entre 0.3 et 0.45, et soulignent l'importance de la dépendance de ΔΩ à la température effective *T*<sub>eff</sub> (voir (Reinhold & Gizon, 2015) et les références s'y trouvant).

La dépendance de la DR au type spectral a été d'abord caractérisée par Collier Cameron (2007) (CC07) sur des
rotateurs rapides M, K, G et F de la séquence principale. Leurs mesures sont aussi basées sur la méthode LSD
et sont présentées en panneau droit de la Figure 2.2. Elles montrent une forte dépendance à la température
effective, compatible avec $\Delta\Omega \propto T_{\text{eff}}^{8.6}$, aussi relevée par B05. De grandes barres d'erreurs sont indiquées pour les
types spectraux les plus chauds, cependant les modélisations de Kueker & Ruediger (2011) semblent confirmer
une forte dépendance, avec différents régimes en fonction de la $T_{\text{eff}}$ observée (voir (Reinhold & Gizon, 2015)
pour une mise en contexte détaillée).

Pour prendre en compte cet aspect, Saar (2010) proposa d'étudier cette dépendance en fonction du Rossby
$Ro_s = \tau_{\text{conv}}/\Omega$ (voir Section 1.2.2.4). En effet, il est possible d'évaluer le temps de retournement convectif
$\tau_{\text{conv}}$ en fonction de $T_{\text{eff}}$ (Cranmer & Saar, 2011), qui sera alors pris en compte dans ces dépendances. Il trouve
en particulier que $\Delta\Omega \propto Ro_s^{-1}$ pour les rotateurs non-saturés ($\Omega \le 12\Omega_\odot$, voir Section 1.3), suggérant alors
$\Delta\Omega \propto \Omega_*$ à $T_{\text{eff}}$ fixée pour ces derniers. Nous reviendrons sur ces tendances en Section 2.3.2.

Ainsi, on voit que la rotation différentielle est une caractéristique des écoulements en enveloppe convective stellaire, partagée par la majorité des étoiles de type solaire observées. Néanmoins, la quantification exacte du contraste rotationnel de surface reste difficile, car elle nécessite un haut degré de précision des instruments utilisés et une période d'acquisition assez longue pour chaque étoile. Les nouvelles observations qui seront apportées par la future mission PLATO devraient permettre des avancées majeures en ce sens. Dans le même temps, la modélisation théorique de ces écoulements s'avère capitale en tant que support aux observations outre la compréhension fondamentale du transport de moment cinétique dans les étoiles, afin de mieux comprendre celles déjà effectuées, et de mieux guider celles à venir.

2.1.2 L'apport des simulations numériques : la rotation différentielle antisolaire

En parallèle des observations stellaires, les simulations numériques ont permis d'apporter une meilleure compréhension
de la dynamique de la convection et des écoulements à grandes échelles dans les systèmes sphériques en rota-
tion. Cela a notamment commencé avec les travaux pionniers de Gilman & Glatzmaier, dont un historique non
exhaustif est énuméré dans Glatzmaier & Gilman (1982). Ils ont pu mettre en évidence que différents régimes
de rotation différentielle peuvent être atteints en fonction de la contrainte appliquée par la rotation. En par-
ticulier, un profil de rotation caractéristique a été relevé lorsque le transport advectif devient aussi important
que celui apporté par la force de Coriolis (Ro = 1) : le profil anti-solaire.

Nous avons vu en Figure 1.16 que le profil de rotation différentielle du Soleil est caractérisé par une rotation équatoriale *prograde*, c'est-à-dire tournant plus rapidement que la rotation globale, ainsi qu'une rotation polaire rétrograde, à l'inverse moins rapide que le taux de rotation moyen de l'étoile. Les profils anti-solaires sont caractérisés par une morphologie inverse à celle du Soleil. Il possède ainsi des pôles rapides et un équateur lent. Nous illustrons un profil anti-solaire idéalisé sur le panneau gauche de la Figure 2.3. Les détails de construction de ce profil sont donnés en Section 5.3.1.2.

On voit donc ici les écoulements polaires progrades en rouge, ainsi que l'écoulement équatorial rétrograde
en bleu. Le panneau de droite illustre $\alpha_e$ en fonction du Rossby convectif $Ro_c$, pour différentes simulations,
<figure><img src="image_50.png" /></figure>

FIGURE 2.3 – Gauche : Coupe méridienne d'un profil idéalisé de rotation différentielle anti-solaire. La barre de couleur représente le taux de rotation normalisé au taux de rotation de la zone radiative interne. Crédits : Noraz et al. (2022a). Droite : Graphique représentant le contraste de rotation $\alpha_e$ en fonction du Rossby convectif $Ro_c$, pour différentes simulations numériques. Le régime II où les valeurs d'$\alpha_e$ sont négatives représente des modèles avec une rotation différentielle anti-solaire. Crédits : Gastine et al. (2014).

calculé avec différents codes (Gastine et al., 2014). $\alpha_e = (\Omega_{Eq} - \Omega_0)/\Omega_0$ est la DR relative entre le taux de rotation équatoriale en surface $\Omega_{Eq}$ et le taux de rotation global de la simulation $\Omega_0$. Une valeur positive représente un équateur prograde, tandis qu'une valeur négative représente un équateur rétrograde. On observe une transition claire se faisant entre les rotations solaires (régime I) et les rotations anti-solaires (régime II) à une valeur autour $Ro_c = 1$, i.e. lorsque la flottaison, moteur de la convection, devient aussi importante que la force de Coriolis. De plus, la différence de rotation relative $|\partial\Omega/\Omega|$ est en moyenne plus forte pour les rotations anti-solaires ($\alpha_e \sim -0.4$ contre $\sim 0.1$ pour les rotations solaires). Finalement, on voit que la rotation différentielle décroît dans le régime II si on continue d'augmenter le Rossby, suggérant l'existence possible d'un autre régime où la puissance de la convection (non-axisymétrique) supprime la rotation différentielle (Gastine et al., 2013).

Afin d'étudier précisément la mise en place et le maintien de ces différents régimes dans un contexte stellaire, nous proposons dans ce chapitre une expérience numérique. Dans les prochaines sections, nous présenterons ainsi une étude paramétrique réalisée à l'aide de simulations MHD, globales sphériques et 3D ; nous permettant de modéliser différents régimes de DR, étudiés en fonction du taux de rotation et de la masse pour différentes étoiles.

## 2.2 Simulations magnétohydrodynamiques globales de dynamos stellaires : Description de l'outil et de l'expérience numérique

Dans cette section, nous allons utiliser le code ASH, validé par le benchmark pour les dynamos convectives globales de Jones et al. (2011), afin d'étudier quinze étoiles de types solaires de masses et rotations variées, afin de comprendre l'influence de ces paramètres sur les écoulements aux grandes échelles.
## 2.2.1 Description des équations résolues

Afin de résoudre les équations régissant la physique des étoiles, nous allons utiliser le code ASH (Anelastic Spherical Harmonic). Ce code permet une simulation à grande échelle, en coordonnées sphériques $(r, \theta, \phi)$ et globale, d'un fluide en rotation. Inspiré du travail de Glatzmaier (1984), le code était à l'origine développé pour la modélisation hydrodynamique des enveloppes convectives de type solaire (Clune et al. 1999, Miesch et al. 2000), puis a été étendu à la résolution complète des équations de magnétohydrodynamique (MHD) (Brun et al., 2004), avec possibilité d'une zone radiative stable. Ici l'approximation MHD est justifiée pour 3 raisons dans le cas d'un intérieur stellaire de type solaire : les mouvements résolus sont non relativistes ($v \ll c$)<fnref n="2" />, se font sur des échelles de distances plus grandes que la longueur de Debye $\lambda_{Debye}$<fnref n="3" />, et de temps plus longues que celle reliée à la fréquence plasma $(1/\omega_p)^4$.

Le code ASH supporte de hautes définitions spatiales et un fonctionnement via MPI sur des machines de calcul massivement parallèles. Le code ASH utilise une méthode de résolution pseudo-spectrale : les équations sont résolues par projection sur les harmoniques sphériques $Y_{l,m}(\theta, \phi)$ le long des directions $\theta$ et $\phi$; pour la résolution selon $r$, la possibilité est soit de projeter spectralement sur les polynômes de Chebyshev $T_n(r)$, soit d'utiliser une méthode de différences finies. La résolution temporelle est assurée par une méthode semi-implicite de Crank-Nicholson pour les termes linéaires, et par une méthode explicite d'Adams-Bashford d'ordre 2 pour les termes non-linéaires ainsi que le terme de Coriolis. Cette résolution se fait entièrement dans l'approximation anélastique (Gough, 1969).

Dans cette approximation, si les équations sur le champ de vitesse $v = (v_r, v_\theta, v_\phi)$ et le champ magnétique $B = (B_r, B_\theta, B_\phi)$ restent entièrement non linéaires, les variables thermodynamiques sont quant à elles linéarisées autour d'un état de référence à symétrie sphérique (état moyen), ayant une densité $\bar{\rho}(r, t)$, une pression $\bar{P}(r, t)$, une température $\bar{T}(r, t)$ et une entropie spécifique $\bar{s}(r, t)$. Les fluctuations autour de cet état de référence sont désignées par $\rho$, $P$, $T$ et $S$. L'approximation anélastique permet de lever l'hypothèse d'incompressibilité proposée par des méthodes comme Boussinesq, sans avoir à considérer les ondes sonores (ce qui implique $\partial\rho/\partial t = 0$). Cela permet alors de conserver les effets importants causés par la stratification de densité dans les intérieurs stellaires, tout en maintenant un critère CFL acceptable. En effet, le critère de stabilité de Courant est alors basé sur les mouvements convectifs plutôt que sur la vitesse du son, permettant des pas de temps beaucoup plus grands quand le nombre de Mach est subsonique, et donc un temps de calcul réduit. Dans le contexte de la MHD, l'approximation anélastique élimine les ondes magnéto-soniques rapides, tout en conservant les ondes d'Alfvén et les ondes magnéto-soniques lentes.

Les équations résolues dans cette approximation sont alors (Clune et al., 1999) :

$$ \nabla \cdot (\bar{\rho}v) = 0, \qquad (2.1) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (2.2) $$

$$ \bar{\rho} \left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} + 2\Omega_* \times \mathbf{v} \right) = -\nabla P + \rho \mathbf{g} + \frac{1}{4\pi} (\nabla \times \mathbf{B}) \times \mathbf{B} + \nabla \cdot \mathbf{D}, \qquad (2.3) $$

$$ \bar{\rho} \bar{T} \frac{\partial s}{\partial t} = - \bar{\rho} \bar{T} \mathbf{v} \cdot \nabla (\bar{s} + s) - \nabla \cdot \mathbf{q} + \Phi_d + \bar{\rho} \epsilon, \qquad (2.4) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times [\mathbf{v} \times \mathbf{B} - \eta \nabla \times \mathbf{B}]. \qquad (2.5) $$

<footnote n="2">2. $c$ est la vitesse de la lumière.</footnote>

<footnote n="3">3. $\lambda_{Debye} = \sqrt{\frac{\epsilon_0 k_B T}{m_e e^2}}$ avec $\epsilon_0$ la permittivité du vide, $k_B$ la constante de Boltzmann, $m_e$ la densité électronique dépendante du nombre de charge, et $e$ la charge élémentaire.</footnote>

<footnote n="4">4. $\omega_p = \sqrt{\frac{n_e e^2}{m_e \epsilon_0}}$ avec $m_e$ la masse de l'électron (voir note précédente)</footnote>
On retrouve premièrement l'équation de conservation de la masse (1.28), dans laquelle la dérivée locale $\partial_t \rho$ est nulle par définition de l'approximation anélastique. L'équation suivante est celle de Maxwell-flux, aussi appelée Maxwell-Thompson, nous indiquant que le flux du champ magnétique se conserve et qu'il n'existe pas de monopôle magnétique. Ensuite viennent respectivement l'équation de Navier-Stokes, l'équation de conservation de l'énergie interne (exprimée ici sur l'entropie), ainsi que l'équation d'induction. On y note $\Omega_* = \Omega_* \hat{e}_z$ la vitesse angulaire dans le repère tournant, $\hat{e}_z$ le vecteur unitaire orienté selon l'axe de rotation. $\bar{\rho}\epsilon$ est un terme de chauffage volumique, représentant l'apport d'énergie par les réactions nucléaires au cœur de l'étoile. Les réactions sont modélisées de façon simple en supposant l'expression $\epsilon = \epsilon_0 T^{n_c}$, où $\epsilon_0$ et $n_c$ sont déterminés en correspondance avec la luminosité de surface de l'étoile. $\mathcal{D}$ est le tenseur des contraintes visqueuses, $\Phi_d$ le terme de dissipation, somme de la contribution visqueuse et ohmique, et $\mathbf{q}$ un flux d'énergie, s'exprimant comme

$$D_{ij} = 2\bar{\rho}\nu \left[ e_{ij} - \frac{1}{3} \nabla \cdot \mathbf{v} \delta_{ij} \right], \qquad (2.6)$$

$$\Phi_d = 2\bar{\rho}\nu \left[ e_{ij}e_{ij} - \frac{1}{3}(\nabla \cdot \mathbf{v})^2 \right] + \frac{4\pi\eta}{c^2} \mathbf{J}^2, \qquad (2.7)$$

$$\mathbf{q} = \kappa_r \bar{\rho} c_p \nabla (\bar{T} + T) + \kappa \bar{\rho} \bar{T} \nabla s + \kappa_0 \bar{\rho} \bar{T} \frac{\partial \bar{s}}{\partial r} \hat{\mathbf{e}}_r. \qquad (2.8)$$

$e_{ij}$ le tenseur d'effort égal à $1/2(\partial_j v_i + \partial_i v_j)$, $\delta_{i,j}$ le symbole de Kronecker, et $\mathbf{J} = c/4\pi \nabla \times \mathbf{B}$ le vecteur densité de courant (exprimé ici par l'équation de Maxwell-Ampère). $c_P$ est la capacité calorifique massique à pression constante, représentant la quantité d'énergie qu'il faut transférer au système pour augmenter sa température. Le flux d'énergie $\mathbf{q}$ s'exprime, respectivement de gauche à droite, comme la somme de trois termes, qui sont le flux radiatif (transport par le rayonnement), le flux conductif (transport par les collisions microscopiques) et le flux irrésolu. Ce dernier flux est une paramétrisation du mouvement des plus petites échelles de convection, présentes à la surface des étoiles que nous modélisons, et qui sont spatialement plus petites que la résolution des simulations. Ainsi, nous pouvons maintenant présenter les coefficients diffusifs $\nu$, $\kappa$ et $\eta$, respectivement la viscosité cinématique, la diffusion thermique et la diffusivité magnétique. Ils représentent (dans ce même ordre) les transports diffusifs de quantité de mouvement, chaleur et champ magnétique, par les interactions microscopiques qui ne sont pas résolues par la simulation. On note aussi en ce sens $\kappa_r$, le coefficient de diffusion radiative, et $\kappa_0$ le coefficient de diffusion effective (ne s'appliquant que sur la composante à symétrie sphérique du gradient d'entropie, $l = 0$). Comme pour $\kappa$, $\kappa_0$ paramétre le transport d'énergie par les processus non résolus, et assure le transport de la luminosité à travers la cocquille sphérique jusqu'à la surface. Cela paramétrise notamment les cellules convectives plus petites que la maille de la grille modélisée, qui apparaissent en particulier près de la surface du modèle (Gilman & Glatzmaier, 1981; Clune et al., 1999; Brun & Toomre, 2002). On supposera ces grandeurs stationnaires et à symétrie sphérique. Toutefois, elles peuvent posséder une dépendance radiale, dont le profil est prescrit pour $\nu$ tel que

$$\nu(r) = \nu_{\text{bot}} + \nu_{\text{top}} f_{\text{step}}(r), \qquad (2.9)$$

$$f_{\text{step}}(r) = (\bar{\rho}/\bar{\rho}_{\text{top}})^{\alpha}[1 - \beta]f(r), \\ f(r) = 0.5(\tanh((r - r_t)/\sigma_t) + 1), \\ \beta = \nu_{\text{bot}}/\nu_{\text{top}} = 10^{-3},$$

et illustré en Figure 2.4. On prend $\alpha = -0.5$, et $\sigma_t$ et $r_t$ varient en fonction de la simulation (précisés par la Table 7 de Brun et al. 2017). Le profil de $\kappa$ peut ensuite être déduit par l'intermédiaire du nombre de Prandtl $Pr = \nu/\kappa$ que l'on garde constant, au sein d'une simulation et à travers l'ensemble proposé dans ce chapitre, à 0.25. $\eta$ peut être déduit de façon similaire, par l'intermédiaire du Prandtl magnétique $Pm = \nu/\eta$ que l'on ajuste soit à 1, soit à 2, selon la simulation considérée (voir Table 2.2). Une transition croissante est placée depuis la zone radiative stable, vers la zone convective turbulente, afin de prendre en compte le transport diffusif effectif supplémentaire, apporté par le brassage des mouvements convectifs. Le profil diffusif $\kappa_0$ est quant à lui
<figure><img src="image_51.png" /><figcaption>FIGURE 2.4 – Profil radial typique de la viscosité cinématique $\nu$ utilisé dans cette étude. Les profils de $\kappa$ et $\eta$ ont une forme identique, mais avec une amplitude dépendant respectivement des nombres de Prandtl $Pr$ et de Prandtl magnétique $Pm$ choisis (voir Table 2.2).</figcaption></figure>

localisé à la surface, et ajusté de manière à transporter le flux d'énergie hors du domaine. Il décroît de manière exponentielle lorsque l'on s'éloigne de la surface, à mesure que le transport enthalpique se fait par des cellules convectives à plus grande échelle.

Afin d'assurer le caractère sans divergence du flux de masse $\bar{\rho}v$ et du champ magnétique $B$ (Équations (2.1 et (2.2)), nous les exprimons selon des composantes poloïdales et toroïdales, tel que

$$ \bar{\rho}v = \nabla \times \nabla \times (W\hat{e}_r) + \nabla \times (Z\hat{e}_r), \qquad (2.10) $$

$$ B = \nabla \times \nabla \times (C\hat{e}_r) + \nabla \times (A\hat{e}_r). \qquad (2.11) $$

avec avec $W$, $Z$, $C$ et $A$ des champ scalaires. Pour fermer le système d'équations, nous considérons une équation d'état du type Gaz Parfait, i.e. $\bar{P} = R\bar{\rho}\bar{T}$ où $R = c_P(\gamma-1)/\gamma$ et $\gamma$ est l'exposant adiabatique. Ainsi, la linéarisation induite par l'approximation anélastique donne

$$ \frac{\rho}{\bar{\rho}} = \frac{P}{\bar{P}} - \frac{T}{\bar{T}} = \frac{P}{\gamma\bar{P}} - \frac{S}{c_P}. \qquad (2.12) $$

### 2.2.2 Conditions initiales et conditions limites

Afin de déterminer l'état moyen unidimensionnel en initialisation du modèle, nous construisons d'abord les profils de gradient d'entropie $d\bar{s}/dr$, de gravité $g$ et de diffusivité radiative $\kappa_{rad}$ à partir du code d'évolution stellaire CESAM (Morel, 1997; Brun et al., 2002). Après avoir obtenu le profil d'entropie spécifique $\bar{s}$, on construit le profil de densité $\bar{\rho}$ satisfaisant l'équilibre hydrostatique comme solution de l'équation suivante

$$ \frac{d\bar{\rho}}{dr} + \frac{g}{\gamma}\bar{\rho}^2 - \gamma e^{-\gamma \frac{\bar{s}}{c_P}} + \frac{\bar{\rho}}{c_P} \frac{d\bar{s}}{dr} = 0, \qquad (2.13) $$

grâce à une méthode de Newton-Raphson. On peut alors déduire le profil de pression $\bar{P}^5$, et ainsi celui de température $\bar{T}$ via l'équation d'état.

<footnote n="5">5. Sachant que pour un gaz parfait, $\bar{S} = C_V \ln \bar{P} - C_P \ln \bar{\rho}$, avec les capacités calorifiques $C_V$ et $C_P$, respectivement à volume et pression constante.</footnote>
Enfin, il nous faut les conditions aux limites afin de bien définir l'ensemble des équations MHD anélastiques présentées. Afin de correctement évaluer et conserver le moment cinétique, on définit ces conditions comme — impénétrables, au sommet (top) et à la base (bot) de la simulation :

$$v_r(r = r_{\text{top}}) = v_r(r = r_{\text{bot}}) = 0, \qquad (2.14)$$

— sans frottements/couple :

$$\frac{\partial}{\partial r} \left( \frac{v_{\theta}}{r} \right) = \frac{\partial}{\partial r} \left( \frac{v_{\phi}}{r} \right) = 0 \bigg|_{r=r_{\text{top}}, r_{\text{bot}}}, \qquad (2.15)$$

— conducteur parfait en limite interne :

$$B_r|_{r_{\text{bot}}} = \left. \frac{\partial}{\partial r} \left( \frac{B_{\theta}}{r} \right) \right|_{r_{\text{bot}}} = \left. \frac{\partial}{\partial r} \left( \frac{B_{\phi}}{r} \right) \right|_{r_{\text{bot}}} = 0, \qquad (2.16)$$

— raccordé à un champ magnétique potentiel en surface :

$$\mathbf{B}|_{r_{\text{top}}} = \nabla \Psi \quad \text{où} \quad \Delta \Psi = 0, \qquad (2.17)$$

— et on impose un flux d'énergie, via le gradient d'entropie :

$$\frac{\partial \bar{s}}{\partial r}(r = r_{\text{bot}}) = a, \quad \frac{\partial \bar{s}}{\partial r}(r = r_{\text{top}}) = b. \qquad (2.18)$$

Cette dernière condition implique de plus que le gradient d'entropie fluctuante $\partial s/\partial r$ est conservé nul à la base et en surface. Les quantités $r_{\text{bot}}$, $r_{\text{top}}$, $a$ et $b$ sont des valeurs précisées à l'initialisation et dépendant de chaque modèle<fnref n="6" />.

### 2.2.3 Espace des paramètres exploré

#### 2.2.3.1 Description des modèles de structure stellaire de référence

Dans ce chapitre, nous nous interessons plus précisemment aux étoiles G et K, ayant une $T_{\text{eff}}$ allant de 4000 à 6000K. Nous considérerons des étoiles ayant la même composition que le Soleil (i.e. même métallicité), et échantillonnons l'espace des paramètres avec les 4 masses 0.5, 0.7, 0.9 et 1.1 $M_{\odot}$. On détermine ensuite leurs différentes quantités physiques au cours de la séquence principale<fnref n="7" /> à l'aide du code d'évolution stellaire CESAM, et les récapitulons en Table 2.1. On y voit notamment qu'une masse $M_{\ast}$ plus grande correspond à une $T_{\text{eff}}$ plus chaude, témoignant alors d'une énergie transportée vers l'extérieur plus importante. En effet, on explore des luminosités allant de 4.6% à 179% de la luminosité solaire. On peut aussi noter la différence de stratification dans la CZ simulée, passant entre 0.5 et 1.1 $M_{\odot}$ d'un $\Delta_{cz}\bar{\rho} = \bar{\rho}_{cz}/\bar{\rho}_{\text{top}}$ de 42 à 81, et témoignant donc d'une stratification plus forte pour ce dernier.

Par la suite, nous avons vu à la Section 1.3.2 que la rotation des étoiles évolue durant leur séquence principale, ralentissant alors par perte de moment cinétique. Nous proposons donc ici de considérer différents taux de rotation pour chacun des modèles proposés, afin de représenter différents âges sur la séquence principale, depuis une rotation rapide caractérisant en majorité les jeunes étoiles, jusqu'à une rotation lente atteinte proche de la TAMS. Pour chaque étoile, nous modélisons des taux de rotation à 5, 3 et $1\Omega_{\odot}$. L'impact de la rotation nous amène à modifier les diffusivités pour chaque modèle. En effet, pour maintenir un haut degré de super-criticalité convective, nous diminuons les diffusivités à mesure que nous augmentons le taux de rotation en suivant

<footnote n="6">6. Voir les rayons en Table 2.2 et les valeurs a,b en Table 7 de Brun et al. (2017).</footnote>

<footnote n="7">7. En toute rigueur, il faut noter que ces valeurs peuvent changer significativement durant la séquence principale (Amard et al., 2019). Dans le cas du Soleil, ces variations peuvent être entre la TAMS et la ZAMS de 200K sur la $T_{\text{eff}}$, un facteur 2 sur $L_{\ast}$ ou encore un facteur 1.3 sur $R_{\ast}$. Nous nous intéresserons ici aux quantités obtenues sur ces modèles d'étoiles en milieu de MS.</footnote>
<table><caption>TABLE 2.1 – Paramètres globaux des 4 étoiles de types solaires en séquence principale, utilisées dans nos modèles dynamo. Les valeurs listées ont été calculées avec le code d'évolution stellaire CESAM (Morel, 1997). Depuis la gauche, nous listons la masse de l'étoile, son rayon, sa luminosité, la température effective, son type spectral, le rayon à la base de la zone convective (BCZ), la température à la BCZ, la densité à la BCZ, le ratio de densité sur la CZ puis sur tout le domaine (f), ainsi que le rayon inférieur (bot) et supérieur (top) du domaine simulé. Crédits : Brun et al. (2022).</caption><thead><tr><th>Mass<br>(M<sub>☉</sub>)</th><th>Rayon<br>(R<sub>☉</sub>)</th><th>L<sub>*</sub><br>(L<sub>☉</sub>)</th><th>T<sub>eff</sub><br>(K)</th><th>Type Sp.</th><th>r<sub>bcz</sub><br>(R<sub>☉</sub>, R<sub>*</sub>)</th><th>T<sub>bcz</sub><br>(K)</th><th>ρ&#773;<sub>bcz</sub><br>(g.cm<sup>-3</sup>)</th><th>Δ<sub>cz</sub>ρ&#773;<br>-</th><th>Δ<sub>f</sub>ρ&#773;<br>-</th><th>r<sub>bot</sub><br>(R<sub>*</sub>)</th><th>r<sub>top</sub><br>(R<sub>*</sub>)</th></tr></thead><tbody><tr><td>0.5</td><td>0.44</td><td>0.046</td><td>4030</td><td>K7</td><td>0.25,0.56</td><td>4.3 × 10<sup>6</sup></td><td>14.0</td><td>42</td><td>193</td><td>0.13</td><td>0.95</td></tr><tr><td>0.7</td><td>0.64</td><td>0.15</td><td>4500</td><td>K4/K5</td><td>0.42,0.66</td><td>3.0 × 10<sup>6</sup></td><td>2.1</td><td>50</td><td>605</td><td>0.32</td><td>0.97</td></tr><tr><td>0.9</td><td>0.85</td><td>0.55</td><td>5390</td><td>G8</td><td>0.59,0.69</td><td>2.6 × 10<sup>6</sup></td><td>0.51</td><td>67</td><td>1013</td><td>0.38</td><td>0.97</td></tr><tr><td>1.1</td><td>1.23</td><td>1.79</td><td>6030</td><td>G0</td><td>0.92,0.75</td><td>1.6 × 10<sup>6</sup></td><td>0.048</td><td>81</td><td>830</td><td>0.5</td><td>0.97</td></tr></tbody></table>

$\nu \propto 1/\Omega_*$, Il s'agit d'un compromis entre le fait de garder la diffusivité constante, rendant alors potentiellement le régime trop laminaire, et la proportionnalité exacte de $\nu \propto 1/\Omega_*^2$ (voir Section 2.2 de Brown et al. 2008), permettant de conserver le degré de super-criticalité mais impliquant alors des niveaux de turbulence trop ambitieux pour les ressources numériques à notre disposition. Ce niveau de super-criticalité sera vérifié a posteriori avec la quantification de $Ra/Ra_{crit}$ pour chaque modèle (voir Table 2.2 et Takehiro et al. 2020).

Á partir d'ici, nous nommerons les modèles comme M**R°m, afin de s'y référer selon sa masse '**', son taux de rotation '°' et la présence de champ magnétique notifié par l'indice 'm' (en comparaison aux progéniteurs hydrodynamiques présentés dans Brun et al. (2017)). Pour les modèles à 0.5, 0.7 et 0.9$M_☉$, nous considérons de plus un taux de rotation encore plus faible, respectivement 1/8, 0.3 et 0.5, que nous noterons M**Sm. Ces derniers taux sont ajoutés pour amener chacune de ces étoiles dans un régime rotationellement peu contraint, où l'émergence d'un profil de DR anti-solaire est attendu (ce qui est déjà le cas pour le modèle M11R1m à $1.1M_☉$).

### 2.2.3.2 Nombres de Rossby attendus

Afin de caractériser au mieux le régime que nous comptons atteindre au travers des différentes simulations que nous proposons dans ce chapitre, nous choisissons de considérer dans un premier temps le Rossby fluide $Ro_f$. Sa définition est rappelée en Section 1.2.2.4 et comparée au Rossby convectif $Ro_c$ utilisé par Gastine et al. (2014) en Section 2.1.2 (voir aussi Aurnou et al. 2007; Anders et al. 2019). Pour mieux comprendre comment celui-ci évolue dans un contexte stellaire, il est possible de le développer de la manière suivante.

Nous avons vu en Section 1.2.1 que les étoiles de type solaire se caractérisent par une structure constituée d'une zone radiative en profondeur et convective en enveloppe de surface. En ce sens, le flux d'enthalpie $F_{en}$ transporté par les cellules convectives est considéré comme l'une des sources principales de transport de luminosité stellaire à travers la zone convective de ces étoiles. Ce transport se fait donc par mouvements de matière, que nous pouvons donc relier à une vitesse convective, afin d'approximer ce transport comme un flux cinétique $F_{en} \sim \rho v^3$. En se plaçant à la base de la zone convective (BCZ), nous avons, par conséquent, $L_* \propto \rho_{BCZ} v_{conv}^3 R_*^2$, avec $\rho_{BCZ}$ la densité en cette position. À partir de sa définition (1.23), le nombre de Rossby fluide $Ro_f$ peut alors s'écrire comme

$$Ro_f \sim \frac{v_{conv}}{2\Omega_* R_*} \propto \left( \frac{L_*}{\rho_{BCZ} R_*^5} \right)^{1/3} \Omega^{-1}. \qquad (2.19)$$
<figure><img src="image_52.png" /><figcaption>FIGURE 2.5 – Position des modèles dans l'espace des paramètres ($M_*$ vs. $P_{\text{rot}}$). La barre de couleur quantifie $Ro_f/Ro_f^\odot$, déduit à partir de considérations théoriques et de modèles d'évolution stellaire (voir Équation (2.22). La place du Soleil est représentée par le signe $\odot$.</figcaption></figure>

Comme les étoiles que nous étudions sont considérées durant leur séquence principale, nous utilisons les relations d'homologie bien connues de la physique stellaire (Kippenhahn et al., 2012), telles que

$$L_* \propto M_*^m \quad \text{et} \quad R_* \propto M_*^n. \qquad (2.20)$$

En paramétrant $\rho_{BCZ}$ pendant la MS, comme $\rho_{BCZ} \sim M_*^p$, nous obtenons alors finalement

$$Ro_f \propto M_* \frac{\frac{m-5n-p}{3}}{\Omega^{-1} \propto M_* \frac{m-5n-p}{3} P_{\text{rot}}}, \qquad (2.21)$$

où $m, n > 0$ et $p < 0$ à priori. On obtient donc une expression du nombre de Rossby fluide en fonction de $M_*$ et de la période de rotation $P_{\text{rot}}$, dont la constante de proportionnalité restera à calibrer.

Dans un second temps, il est maintenant intéressant d'essayer de prévoir, avec cette formule, la valeur qui sera atteinte pour chacune des étoiles simulées entre 0.5 et 1.1 $M_\odot$. Pour cela, des régressions linéaires sur les valeurs de la Table 2.1 permettent de déduire les valeurs de $m = 4.6$, $n = 1.3$ et $p = -6.9$, donnant alors $Ro_f \propto M_*^{1.7} P_{\text{rot}}$ (Brun et al., 2017). On voit donc que le Rossby augmente avec la masse $M_*$ de l'étoile, ce qui implique des mouvements convectifs plus intenses en comparaison à la rotation. En effet, une masse plus grande nécessite une pression thermique plus élevée pour lutter contre l'effondrement gravitationnel du système. Une pression plus grande amène alors à une température centrale plus élevée, ce qui signifie plus d'énergie libérée par les réactions de fusion<fnref n="8" /> et donc une luminosité transportée jusqu'à la surface plus intense. Afin d'atténuer la dépendance aux différentes définitions du nombre de Rossby, nous le normaliserons à la valeur solaire, donnant

$$\frac{Ro_f}{Ro_f^\odot} = \left(\frac{M_*}{M_\odot}\right)^{1.7} \left(\frac{P_{\text{rot}}}{P_{\text{rot},\odot}}\right), \qquad (2.22)$$

qui est donc valable dans l'intervalle de masse utilisé pour faire les régressions linéaires ($0.5 < M_* < 1.1 M_\odot$) et à la métallicité solaire. Cela nous permet de tracer la Figure 2.5 et d'y placer chacun des modèles considérés.

<footnote n="8">8. Voir définition de $\epsilon$ Section 2.2.1</footnote>
On y représente en couleur le Rossby fluide $Ro_f$, normalisé à la valeur solaire et obtenu en Équation (2.22), en fonction de $M_*$ et $P_{\text{rot}}$. Nous indiquons par des cercles blancs la position des 15 modèles dans cet espace des paramètres, ainsi que la position du Soleil avec un cercle noir. Cela nous indique donc que 11 des 15 modèles sont prédits avec un Rossby inférieur à la valeur solaire, jusqu'à 10 fois inférieur, tandis que 4 modèles sont au-dessus, un dépassant même le double. Nous verrons par la suite comment de telles différences de régimes peuvent impacter la dynamique interne des modèles considérés dans ce chapitre, pour lesquels j'ai contribué aux calculs de l'état dynamo.

### 2.2.4 Illustration d'une dynamo convective 3D

Les 15 simulations de dynamo que nous présentons ici ont été développées à partir des progéniteurs hydro-dynamiques convectifs présentés et analysés par Brun et al. (2017). Ces états hydrodynamiques initiaux ont été calculés suffisamment longtemps pour atteindre un état statistiquement stationnaire, et avoir développé un profil de rotation bien établi. Ces profils possèdent une tachocline construite de manière auto-consistante, définie comme la transition entre la zone radiative en rotation solide, et l'enveloppe convective externe en rotation différentielle, formant alors une région à fort cisaillements (Spiegel & Zahn, 1992).

En repartant de ces états, nous avons ajouté une graine magnétique dipolaire au sein de la zone convective, que nous avons laissé se développer. Après convergence des différents modèles MHD, nous obtenons donc des simulations globales 3D de dynamo convectives pour 15 étoiles de type solaire. Nous présentons en Figure 2.6 l'une d'entre elle, le modèle M09R3m, possédant des caractéristiques proches de notre Soleil. En effet, on voit tout d'abord que la zone équatoriale de basse latitude possède un écoulement prograde, caractéristique du profil de DR solaire. Ces écoulements convectifs prennent la forme de "cellules bananes", dont les extrémités possèdent un écoulement rétrograde à moyenne latitude. Ces motifs sont caractéristiques d'une dynamique significativement impactée par la force de Coriolis ($Ro_f < 1$), et sont cohérents avec les corrélations statistiques $\langle v_r v_\phi \rangle$ produisant un transport efficace du moment cinétique via le tenseur de Reynolds, comme nous le verrons en Section 2.4. Nous voyons ensuite une présence de motifs convectifs plus petits dans la région polaire, signe que la dynamique radiale de ces hautes latitudes ressent un impact moins fort de la rotation (Anders et al., 2019).

On propose ensuite une extrapolation potentielle du champ magnétique pris à cette même surface de la simulation, pour une surface source située à $r_{ss} = 1.8 R_{\text{top}}$ (voir l'extrapolation potentielle en Annexe C). Les parties pourpres représentent une composante magnétique radiale pointant à l'extérieur de l'étoile, tandis que les parties bleues correspondent à une composante radiale pointant vers son intérieur. On observe que la majorité des lignes de champ magnétique ouvertes se situe dans les régions polaires, cependant le modèle montre une topologie générale plus complexe, en particulier dans la région équatoriale. Cela vient du fait que l'instant de la simulation illustré est un maximum d'activité d'un cycle magnétique de type solaire. Notons tout de même ici que la topologie exacte des boucles de surface est fortement dépendante du $r_{ss}$ choisi, et que cette visualisation n'est qu'une illustration qualitative d'une topologie possible. Nous réservons l'étude et la caractérisation de l'activité magnétique de la dynamo pour le Chapitre 6, et allons tout d'abord étudier les écoulements grandes échelles établis au sein des différents modèles dans le cas MHD.

## 2.3 Caractérisation des écoulements à grande échelle

Nous présentons dans cette section les différents écoulements à grande échelle obtenus dans l'ensemble des 15 simulations MHD, et étudions comment ces derniers sont impactés par la masse, le taux de rotation et la présence d'un champ magnétique. Nous commençons par présenter les quantités importantes pour l'analyse
<figure><img src="image_53.png" /><figcaption>FIGURE 2.6 – Illustration de la convection et du champ magnétique du modèle M09R3m. On représente la convection à la surface de la simulation par un nivelage orangé de celle-ci. Les zones plus claires dont l'altitude est remontée représentent des $v_r > 0$ pouvant atteindre 168 cm.s⁻¹ pour les zones blanches, tandis que les creux plus sombres représentent des $v_r < 0$ pouvant descendre à -210 cm.s⁻¹. Les lignes de champ représentent le champ magnétique extrapolé à partir de la surface du modèle, où les parties bleues représentent des zones de $B_r < 0$, tandis que les parties pourpres représentent $B_r > 0$.</figcaption></figure>
des différents modèles, puis nous caractériserons les différents régimes de rotation différentielle ainsi que de circulations méridiennes établies dans les modèles dynamos considérés.

### 2.3.1 Analyse globale

En repartant de ces états, nous avons ajouté une graine magnétique dipolaire au sein de la zone convective, que nous avons laissé se développer. Les 15 modèles développent une dynamo qui crée et maintient un champ magnétique contre sa dissipation ohmique. Les modèles sont tous calculés sur une dizaine d'années au minimum (pouvant monter à 30 comme les cas M07R1m ou M09R3m), à l'exception du cas M11R1m ayant atteint rapidement un état quasi-stationnaire, et dont le calcul est toujours en cours pour confirmation sur le long terme. Une fois les différentes solutions dynamos établies, nous avons pu les analyser et construire la Table 2.2, récapitulant les principales quantités physiques d'intérêt, évaluées en milieu de zone convective.

<table><caption>TABLE 2.2 – Quantités et nombres adimensionnés caractérisant les modèles. Les nombres adimensionnés et $\tilde{v}_r$ sont moyennés radialement sur 0.01$R_*$, au milieu de l'enveloppe convective. Sont listés depuis la gauche : le taux de rotation, la viscosité cinématique au sommet du domaine, la moyenne quadratique de la vitesse radiale, le nombre de Reynolds $Re = \tilde{v}_r \Delta_{conv}/\nu$, le Prandtl magnétique $P_m = \nu/\eta$, le degré de super-criticalité $Ra^*/Ra_c$, défini comme le ratio entre le Rayleigh modifié et Rayleigh critique, calculés par Takehiro et al. (2020), et où $\Delta_{conv} = r_{top} - r_{bot}$. Vient ensuite les différents nombres de Rossby : fluide $Ro_f = |\omega|/2\Omega_*$, convectif $Ro_c = \sqrt{Ra/TaPr}$ et stellaire $Ro_s = P_{rot}/\tau_c^{CS}$. Le Rossby stellaire $Ro_s$ a été calculé en considérant le temps de retour-nement convectif $\tau_c$ empirique, proposé par Cranmer & Saar (2011) ($\tau_c^{CS} = 314.24\exp[-\frac{T_{eff}}{1952.5 K} - (\frac{T_{eff}}{6250 K})^{18}]$ + 0.002 jours), afin de comparer directement les simulations aux observations. Précisons de plus que nous trouvons une linéarité telle que $Ro_s \approx 2.26Ro_f$. Finalement, nous présentons le contraste de rotation différentielle, pris entre l'équateur et une latitude de 60°, pour les modèles avec et sans champ magnétique (respectivement MHD et HD). Adapté de Brun et al. (2017) et Brun et al. (2022)</caption><thead><tr><th></th><th>Ω<sub>*</sub><br>(Ω<sub>☉</sub>)</th><th>ν<sub>top</sub><br>(cm<sup>2</sup>.s<sup>-1</sup>)</th><th>v&#773;<sub>r</sub><br>(m/s)</th><th>Ra<sup>*</sup><br>(Ra<sub>c</sub>)</th><th>Re</th><th>P<sub>m</sub></th><th>Ro<sub>f</sub></th><th>Ro<sub>c</sub></th><th>Ro<sub>s</sub></th><th>ΔΩ<br>(nHz)<br>MHD (HD)</th></tr></thead><tbody><tr><td>M05Sm</td><td>1/8</td><td>10.5 × 10<sup>11</sup></td><td>13.52</td><td>67.64</td><td>56.34</td><td>1</td><td>1.74</td><td>0.38</td><td>3.96</td><td>-23 (-24)</td></tr><tr><td>M05R1m</td><td>1</td><td>4.47 × 10<sup>11</sup></td><td>7.27</td><td>16.83</td><td>71.71</td><td>2</td><td>0.33</td><td>0.31</td><td>0.70</td><td>112 (129)</td></tr><tr><td>M05R3m</td><td>3</td><td>2.58 × 10<sup>11</sup></td><td>6.21</td><td>8.58</td><td>106.08</td><td>1</td><td>0.15</td><td>0.17</td><td>0.23</td><td>200 (85)</td></tr><tr><td>M05R5m</td><td>5</td><td>2.00 × 10<sup>11</sup></td><td>6.95</td><td>5.69</td><td>118.72</td><td>1</td><td>0.07</td><td>0.10</td><td>0.14</td><td>9 (146)</td></tr><tr><td>M07Sm</td><td>1/4</td><td>5.31 × 10<sup>12</sup></td><td>25.44</td><td>73.32</td><td>26.16</td><td>1</td><td>1.24</td><td>0.72</td><td>2.98</td><td>-53 (-32)</td></tr><tr><td>M07R1m</td><td>1</td><td>2.91 × 10<sup>12</sup></td><td>16.34</td><td>28.75</td><td>30.69</td><td>2</td><td>0.42</td><td>0.50</td><td>0.89</td><td>111 (120)</td></tr><tr><td>M07R3m</td><td>3</td><td>1.68 × 10<sup>12</sup></td><td>14.74</td><td>14.78</td><td>47.94</td><td>2</td><td>0.16</td><td>0.23</td><td>0.30</td><td>68 (187)</td></tr><tr><td>M07R5m</td><td>5</td><td>1.30 × 10<sup>12</sup></td><td>13.42</td><td>8.91</td><td>56.36</td><td>1</td><td>0.09</td><td>0.17</td><td>0.18</td><td>-2 (223)</td></tr><tr><td>M09Sm</td><td>1/2</td><td>1.43 × 10<sup>13</sup></td><td>53.51</td><td>54.09</td><td>27.69</td><td>1</td><td>1.28</td><td>0.74</td><td>3.02</td><td>-36 (-25)</td></tr><tr><td>M09R1m</td><td>1</td><td>1.01 × 10<sup>13</sup></td><td>38.74</td><td>24.71</td><td>28.48</td><td>2</td><td>0.68</td><td>0.66</td><td>1.51</td><td>102 (108)</td></tr><tr><td>M09R3m</td><td>3</td><td>5.81 × 10<sup>12</sup></td><td>30.61</td><td>11.02</td><td>38.90</td><td>1</td><td>0.27</td><td>0.35</td><td>0.50</td><td>265 (288)</td></tr><tr><td>M09R5m</td><td>5</td><td>4.50 × 10<sup>12</sup></td><td>27.94</td><td>7.00</td><td>38.21</td><td>1</td><td>0.10</td><td>0.20</td><td>0.30</td><td>76 (338)</td></tr><tr><td>M11R1m</td><td>1</td><td>3.80 × 10<sup>13</sup></td><td>130.77</td><td>47.17</td><td>32.06</td><td>2</td><td>1.38</td><td>0.78</td><td>3.30</td><td>-102 (-131)</td></tr><tr><td>M11R3m</td><td>3</td><td>2.20 × 10<sup>13</sup></td><td>90.23</td><td>17.16</td><td>38.31</td><td>1</td><td>0.54</td><td>0.50</td><td>1.10</td><td>278 (291)</td></tr><tr><td>M11R5m</td><td>5</td><td>1.70 × 10<sup>13</sup></td><td>88.50</td><td>6.87</td><td>38.81</td><td>1</td><td>0.27</td><td>0.34</td><td>0.66</td><td>109 (435)</td></tr></tbody></table>

Tout d'abord, on confirme qu'une masse plus importante implique de plus grandes vitesses de convection (voir la section précédente). On voit en effet que la moyenne quadratique de la vitesse radiale $\tilde{v}_r$, augmentant par exemple de 7.27 à 130.77 m.s⁻¹, respectivement pour M05R1m et M11R1m. Cependant, on voit qu'elle diminue légèrement avec un taux de rotation plus important, passant par exemple de 25.44 à 13.42 m.s⁻¹ pour les modèles à 0.7$M_☉$. Ceci témoigne de la baisse du niveau de super-criticalité lorsque l'impact de la rotation augmente (Chandrasekhar, 1961), cette dernière inhibant les mouvements convectifs (ces deux tendances sont aussi vérifiées pour la composante poloïdale $\tilde{v}_\theta$, listée en Table 2 de la section 2.7). Ce niveau de super-criticalité est quantifié par le nombre de Rayleigh $Ra$, ici normalisé à sa valeur critique pour chaque simulation. Cette
valeur critique est la valeur minimum pour déclencher une instabilité convective, et donc le transport d'énergie par mouvement de matière dans le système considéré. Les nombres de Reynolds nous montrent ensuite que ces simulations ont été construites de manière à évoluer dans un régime de turbulence similaire (entre $Re = 26$ et 106), autant d'un point de vue fluide que magnétique (le Reynolds magnétique $Rm = Re * Pm$).

Ensuite, nous calculons les différents Rossby stellaires $Ro_s$, convectifs $Ro_c$ et fluides $Ro_f$ définis en section 1.2.2.4 des différents modèles. Concernant $Ro_f$, cela nous permet de les comparer à ceux prédits analytiquement en Section 2.2.3.2 et donc de vérifier la solidité des arguments théoriques proposés. Pour cela, nous procédons à une régression linéaire de $Ro_f$ sur les paramètres $M_*$ et $P_{\text{rot}}$, ce qui nous donne

$$Ro_f = (0.87 \pm 0.1) \times \left(\frac{\Omega_*}{\Omega_{\odot}}\right)^{-0.9 \pm 0.06} \times \left(\frac{M_*}{M_{\odot}}\right)^{1.63 \pm 0.24} \quad (2.23)$$

On note que la relation trouvée par l'expérience numérique est donc relativement proche de ce qui était attendu en équation (2.22). On note de plus que la constante de proportionnalité nous donne ici directement une valeur de $Ro_f = 0.87$ pour le régime solaire. Ensuite, il est intéressant de noter qu'une autre régression linéaire nous permet de trouver une relation linéaire $Ro_s \approx 2.26Ro_f$ entre le Rossby fluide $Ro_f$ tiré des simulations, et le Rossby stellaire $Ro_s$ résultant d'une définition généralement utilisée dans un contexte observationnel (Cranmer & Saar, 2011). Ceci permet de confirmer l'utilité et la pertinence de cet ensemble de simulations numériques, en tant que support aux observations. C'est pourquoi nous allons maintenant pouvoir interpréter par la suite les résultats de la dernière colonne, dans un contexte observationnel. Il s'agit du contraste de la rotation différentielle de surface, pris entre une latitude de 60° et l'équateur. Une valeur positive correspond à une rotation différentielle de type solaire, et une valeur négative à une rotation différentielle anti-solaire. Nous complétons ces valeurs, avec celles des progéniteurs hydrodynamiques entre parenthèses dans la dernière colonne de la Table 2.2.

Nous allons maintenant présenter les différents régimes observés dans l'ensemble des 15 simulations, puis nous illustrerons les tendances aux paramètres stellaires globaux, dans le contexte observationnel présenté en début de chapitre.

## 2.3.2 Différents régimes et tendances de la rotation différentielle

Dans l'ensemble des 15 simulations calculées, nous relevons 3 régimes principaux de rotation différentielle. Nous les résumons sur la partie inférieure de la Figure 2.7 avec des profils représentatifs, et ajoutons sur la partie supérieure les profils des progéniteurs hydrodynamiques (HD) correspondant pour une comparaison visuelle directe. Les barres de couleur représentent les vitesses d'écoulements zonaux, i.e. le long de la direction longitudinale. Elles s'étendent de la couleur bleu (écoulements rétrogrades) à la couleur rouge (écoulements progrades) et sont centrées en jaune pâle sur le taux de rotation global, qui est le taux de la rotation solide de l'intérieur radiatif. Dans un but de concision, seulement 3 modèles représentatifs sont illustrés sur cette figure. La totalité des profils établis en sont disponibles en Figure 5 de Brun et al. (2022) (attaché en Section 2.7) pour le régime MHD, et en Figure 11 de Brun et al. (2017) pour les progéniteurs hydrodynamiques correspondants.

- Le premier est le régime *anti-solaire*. Visible sur le panneau de gauche, il est caractérisé par un équateur lent et des pôles rapides. Il s'agit des cas M05Sm, M07Sm, M09Sm et M11R1m, pour lesquels $Ro_f \ge 1.24$. Ce sont des cas où la dynamique est peu contrainte par la rotation, au sens de la mécanique des fluides. Ces profils sont donc attendus pour de vieilles étoiles de la séquence principale en rotation lente, ou potentiellement plus évoluées. On note que la présence d'un champ magnétique a peu d'impact sur ces profils de rotation, en dehors du cas M11R1m pour lequel le contraste a diminué.
<figure><img src="image_54.png" /></figure>

FIGURE 2.7 – Profils de rotation différentielle, émergeant pour différents taux de rotation (0.5, 3 et 5 Ω⊙ de gauche à droite) et moyennés temporellement sur 10 temps de retournement convectif. Le panneau supérieur montre les profils de DR obtenus dans les modèles hydrodynamiques de Brun et al. (2017) (M09S, M09R3 et M09R5), puis les cas magnétiques correspondants sur le panneau inférieur (M09Sm, M09R3m et M09R5m). La barre de couleur quantifie le taux de rotation Ω(r, θ). L'équateur et la tachocline sont représentés par les lignes pointillées. Les panneaux de gauche représente des cas *anti-solaires*, correspondants à des nombres de Rossby fluides élevés > 1. Ceux du milieu sont des cas de *type solaire*, correspondant à des Rossby fluides intermédiaires [0.15; 0.9]. Enfin, à droite sont représentés des profils caractérisés par des Rossby les plus faibles, dit *cylindriques/Joviens* dans le cas hydrodynamique, ou quasi-uniformes ("gelés") dans le cas magnétique.

— Le second régime est appelé *type solaire*, en référence au profil de rotation de notre étoile. Visible sur le panneau du milieu, il est caractérisé, à l'inverse, par un équateur rapide et des pôles lents. On le retrouve dans les cas M05R1m, M05R3m, M07R1m, M07R3m, M09R1m, M09R3m et M11R3m, pour des Rossby fluides dans l'intervalle $0.68 \ge R_{of} \ge 0.15$. Ces profils sont donc attendus pour des $R_{of}$ du même ordre de grandeur que la valeur solaire, et à priori inférieurs à 1. On note que cette fois-ci, le contraste de rotation différentielle a diminué avec la présence d'un champ magnétique, pour la majorité de ces cas. Cet effet est d'autant plus important que le taux de rotation est élevé pour une masse donnée. Il s'agit de la rétroaction du champ magnétique sur le profil de rotation qui l'a créé, appelée *Ω-quenching*.

— Finalement, un troisième régime *cylindrique* apparaît pour les rotateurs les plus rapides à $\Omega_* = 5\Omega_⊙$. Les profils de rotation sont alors fortement contraints par la rotation, et se cylindrisent. Il s'agit de l'expression de la contrainte de Taylor-Proundman $\partial\langle v_\phi \rangle / \partial z = 0$, apparaissant dans les enveloppes convectives en rotation rapide (Gilman, 1975; Elliott et al., 2000; Miesch et al., 2006). On les trouve $R_{of} \le 0.27$.
<figure><img src="image_55.png" /></figure>

FIGURE 2.8 – Tendances de la rotation différentielle en fonction de différents paramètres stellaire. Pour les 3 panneaux, on représente les différents modèles avec un code couleur et différentes formes. Les masses 0.5, 0.7, 0.9 et $1.1M_\odot$ sont respectivement représentées en rouge, jaune, vert et bleu. Les rotations à taux sous-solaire, 1, 3 et $5\Omega_\odot$ sont respectivement représentés par des triangles, losanges, carrés et ronds. Les modèles MDH possédant un champ magnétique sont représentés par les grandes formes pleines, tandis que les modèles HD sans champ magnétique Brun et al. (2017) le sont par les petites formes évidées. *Gauche* : Rotation différentielle de surface $|\Delta\Omega|$ en fonction de la période de rotation $P_{\text{rot}}$. Les petits points noirs représentent l'échantillon de la mission *Kepler* extraits de Reinhold & Gizon (2015), tandis que les lignes en pointillés noirs illustrent les seuils de détection de cette campagne observationnelle. *Milieu* : $|\Delta\Omega|$ en fonction du taux de rotation $\Omega_\star$. Nous y représentons des régressions linéaires afin de déduire la tendance dans le cas MHD en noir, et HD en gris. *Droite* : Contraste de rotation différentielle, normalisé au taux de rotation, en fonction de l'inverse du Rossby fluide $Ro_f^{-1}$. Crédits : Brun et al. (2022)

Lorsque le profil cylindrique présente de plus une structure en bandes avec une alternance de jets cylindriques, on dit de plus qu'il est de type *jovien*, en référence à Jupiter. Ce comportement émerge timidement sur le profil de $M11R5m$, mais est plus flagrant lorsque l'influence du champ magnétique est faible (voir le cas HD M09R5 sur le profil en haut à droite). En effet, la rétroaction du champ magnétique dans ces rotateurs rapides est amplifiée par un processus dynamo plus puissant, rapprochant le profil de DR d'un profil de rotation solide qu'on désignera comme *quasi-uniforme* (ou "gelé", voir profil en bas à droite). Pour les 2 modèles au Rossby les plus faibles (M05R5m et M07R5m), les contrastes de DR sont devenus si faibles qu'ils passent en-dessous du seuil de détection de la mission *Kepler*, comme illustré en Figure 2.8.

Le panneau de gauche de la Figure 2.8 représente la valeur absolue du contraste de rotation différentielle de surface $|\Delta\Omega|$ en fonction de la période de rotation de l'étoile $P_{\text{rot}}$, comparée aux détections de l'échantillon *Kepler*. On voit qu'en dehors des deux modèles en rotation solide, les étoiles que nous modélisons entrent dans les seuils de détection de la mission. Cependant, on voit que 3 des 4 modèles anti-solaires en rotation lente sont dans une zone vide d'observations pour cet espace des paramètres. Ce sera l'objet de perspectives observationnelles présentées au Chapitre 3.

On trace ensuite sur le panneau du milieu cette même quantité $|\Delta\Omega|$, mais cette fois-ci en fonction du taux de rotation $\Omega_\star = 2\pi/P_{\text{rot}}$. Brun et al. (2017) avaient montré que les progéniteurs hydrodynamiques suivaient une tendance en $\Delta\Omega \propto \Omega_*^n$ avec l'indice $n = 0.66$. Nous montrons ici que cet indice prend la valeur $n = 0.46$ pour les cas MHD. Cette réduction de la dépendance témoigne directement de la retro-action du champ magnétique sur les écoulements globaux, d'autant plus forte que la rotation est rapide, et crée un champ magnétique intense. Bien que les tendances observationnelles soient encore très incertaines à l'heure actuelle (voir Section 2.1.1), la régression trouvée pour les modèles HD ainsi que celle que nous trouvons en MHD correspondent aux tendances générales. Il faut garder en tête que les intervalles de Rossby et masses étudiés ici sont relativement petits, cependant nous pouvons noter que la régression MHD est plus significative que celle des modèles HD pour $\Omega_* \le 3\Omega_\odot$. Pour les rotateurs rapides, la rétroaction du champ magnétique crée une chute du contraste de DR et ajoute une forte dispersion des points.
Un moyen de prendre en compte l'effet du taux de rotation global $\Omega_*$ est de l'utiliser pour normaliser le contraste discuté. Nous le représentons ainsi sur le panneau de droite de la Figure 2.8 en fonction de l'inverse du nombre de Rossby (aussi appelé nombre de Coriolis $Co = Ro^{-1}$). On remarque que nos données numériques sont en accord avec les tendances observationnelles de Saar (2010). Premièrement, on voit que les modèles avec un $Ro_f^{-1}$ compris entre 1 et 5 montre une tendance relativement constante pour $|\Delta\Omega/\Omega_*|$. Il s'agit de modèles possédant un profil de rotation de type solaire, sauf pour le modèle M11R5m qui est un peu à l'écart. On observe ensuite une baisse de ce contraste normalisé pour les rotateurs rapides à faible Rossby ($Ro_f^{-1} > 5$). Le nombre de modèles MHD dans ce régime est clairement trop faible pour caractériser précisément cette pente de type $|\Delta\Omega/\Omega_*| \propto Ro_f^d$, néanmoins il est possible de la situer telle que $d \in [2, 6]$, ce qui reste cohérent avec la tendance observationnelle $d = 2$ de Saar (2010). Il reste intéressant de noter qu'une diminution plus faible est trouvée dans le cas des modèles HD, ce qui confirme le contrôle de cette tendance en $d$ par la rétroaction magnétique dont une analyse énergétique sera détaillé en Section 2.5. Finalement, nos modèles montrent une possible remontée de $|\Delta\Omega/\Omega_*|$ pour les rotateurs lents avec un profil anti-solaire ($Ro_f^{-1} < 1$). Comme illustré sur le panneau de gauche, c'est un régime pour lequel nous manquons de contraintes observationnelles afin de confirmer un tel comportement. Ainsi, nous proposerons en Chapitre 3 une méthodologie pour une recherche d'étoiles à fort Rossby, et donc potentiellement hôtes de profils anti-solaires.

Cette comparaison de modèles avec et sans champ magnétique, puis leur confrontation aux données observationnelles, nous montre donc que les observations actuelles semblent bien témoigner de la présence d'un champ magnétique au sein des étoiles, mais surtout d'une potentielle rétroaction significative de ce champ sur le profil de rotation différentielle. En effet, on a vu en Figure 2.7 que cet impact est visible sur les profils axisymétriques dans le plan méridien pour les $Ro_f$ les plus faibles. Nous verrons de plus au Chapitre 6 que cette rétroaction magnétique à grande échelle peut se faire sur des échelles de temps plus courtes ("oscillations torsionnelles"), potentiellement à l'origine des cycles décennaux de type solaire. Puisqu'il existe des régimes dans lesquels la rétroaction magnétique peut contrôler la dynamique ces écoulements longitudinaux à grande échelle, on peut se demander si cette rétroaction aura aussi un impact sur la circulation méridienne, et si oui, quelle en sera l'ampleur. Poursuivant ce questionnement, nous proposons en section suivante d'illustrer et d'étudier les profils des écoulements méridiens, résultant de la dynamique grande échelle que nous venons de voir.

### 2.3.3 Circulation Méridienne

Comme nous l'avons vu en Section 1.2.2.2, le Soleil possède une circulation méridienne, i.e., un écoulement à grande échelle axisymétrique dans le plan méridien de l'étoile. Cet écoulement $v_m$ apparaît en réaction à la conservation du moment cinétique. On peut obtenir son équation d'évolution en prenant la moyenne longitudinale et temporelle de la composante $\phi$ dans l'équation de Navier-Stokes (2.3), telle que

$$ \bar{\rho} \frac{\partial \mathcal{L}}{\partial t} + \langle \bar{\rho} v_m \rangle \cdot \nabla \mathcal{L} = \mathcal{T}, \qquad (2.24) $$

où $\mathcal{L} = r \sin\theta(v_\phi + r \sin\theta\Omega_*)$ est le moment cinétique spécifique et $\Omega_*$ le taux de rotation global. Le second terme représente le transport de moment cinétique par la circulation méridienne, tandis que le terme de droite $\mathcal{T}$ contient le reste des termes appliquant un couple : le transport turbulent, la viscosité, et le magnétisme. Dans un état statistiquement stationnaire ($\partial\mathcal{L}/\partial t$), on peut ainsi considérer la circulation méridienne comme une circulation émergeant en réaction à un déséquilibre de ces derniers termes ($\mathcal{T} \neq 0$) construisant la rotation différentielle. Dans l'exemple du Soleil où $\nabla \mathcal{L}$ est dirigé vers l'extérieur au niveau de l'équateur, un couple $\mathcal{T} > 0$ produira une circulation méridienne s'éloignant de l'axe de rotation, tandis qu'un couple $\mathcal{T} < 0$ produira un écoulement en direction de cet axe. Cet effet est appelé "pompage gyroscopique" (McIntyre, 1999; Miesch et al., 2008).
Afin de représenter cette circulation $\langle \bar{\rho}v_m \rangle = \langle \bar{\rho}v_r \rangle e_r + \langle \bar{\rho}v_\theta \rangle e_\theta$, nous calculons la fonction de courant $\Psi$ définie dans Miesch et al. (2000) comme

$$r \sin \theta \langle \bar{\rho}v_r \rangle = -\frac{1}{r} \frac{\partial \Psi}{\partial \theta} \quad \text{et} \quad r \sin \theta \langle \bar{\rho}v_\theta \rangle = \frac{\partial \Psi}{\partial r} \qquad (2.25)$$

Cette définition nous permet ensuite de représenter en Figure 2.9 la fonction de courant de chacun de nos modèles MHD. Pour comparaison directe avec les progéniteurs hydrodynamiques sans champ magnétique, une Figure similaire est disponible dans Brun et al. (2017). Nous retrouvons de manière similaire à la section précédente 3 grands régimes, correspondant aux différentes morphologies des fonctions de courant $\Psi$ :

- Pour les modèles à haut Rossby ($R_{0f} \ge 1.24$) où les profils de rotation anti-solaire sont observés, la circulation méridienne est peu impactée par la présence d'un champ magnétique. Chaque hémisphère possède une cellule de circulation bien définie, dont l'écoulement se fait vers les pôles en surface, puis revient en direction de l'équateur au fond de la zone convective. Dans ce cas, l'écart à l'équilibre géostrophique se fait par la baisse significative de l'influence de la force de Coriolis dans la balance globale, créant alors de fortes circulations.

- Pour les modèles de Rossby intermédiaire dans l'intervalle $0.68 \ge R_{0f} \ge 0.15$, les profils restent aussi assez similaires aux progéniteurs hydrodynamiques, et se complexifient par rapport aux modèles à rotation anti-solaire, tout en diminuant d'un ordre de grandeur en amplitude par rapport à ces derniers. Il s'agit des cas M05R3m, M07R1m, M09R1m, M09R3m, M11R3m et M11R5m, possédant des profils multi-cellulaires dans chaque hémisphère. Cette différence par rapport à au régime précédant est en accord avec les résultats trouvés par l'étude de Featherstone & Miesch (2015), ne considérant ni zone stable interne, ni champ magnétique. On retrouve une symétrie équatoriale des profils de circulation, avec des cellules allongées le long de l'axe de rotation lorsqu'elles sont en dehors du cylindre tangent 9. Dans ces cas-ci, l'écart à l'équilibre géostrophique est principalement dû au transport de moment cinétique par le tenseur de Reynolds, résultant des corrélations au sein de la turbulence convective. On peut cependant noter deux exceptions, M05R1m et M07R3m, pour lesquelles la symétrie équatoriale n'est pas respectée, et dont le profil méridien se rapproche du prochain régime.

- Finalement, un troisième régime apparaît pour les nombres de Rossby les plus faibles ($R_{0f} < 0.16$), ainsi que pour le modèle M05R1m. On observe alors les plus fortes différences morphologiques par rapport aux progéniteurs hydrodynamiques, dues à l'impact du champ magnétique sur l'établissement de la rotation différentielle. En particulier, des cellules de circulation trans-équatoriale apparaissent au niveau de l'équateur, et brisent la symétrie équatoriale de la circulation présente dans les autres modèles. Notons que l'amplitude des écoulements méridiens de ce régime reste toutefois du même ordre de grandeur que ceux du régime de Rossby intermédiaire.

On conclut ainsi que la présence d'un champ magnétique modifie la forme de la circulation méridienne des modèles simulés, en particulier pour les rotations les plus rapides. Cependant, une analyse des transferts énergétiques nous montrera en Section 2.5 que les contraintes magnétiques ne s'appliquent pas d'une manière directe sur cette circulation. En effet, ce sont les changements dans la rotation différentielle qui entrainent la modification des profils méridiens par conservation du moment cinétique. Pour mieux comprendre la mise en place de ces changements de DR, nous proposons dans la prochaine section d'étudier les transferts de moment cinétique des différents termes de la balance (2.24) établissant les écoulements moyens.

<footnote n="9">9. Cylindre dont le rayon cylindrique est défini par la base équatoriale de la zone convective.</footnote>
<figure><img src="image_56.png" /><figcaption>FIGURE 2.9 – Profils des circulations méridiennes résultant des cas MHD, moyennés temporellement et longitudinalement sur 10 temps de retournement convectif. Les contours sont normalisés et indiquent des circulations dans le sens horaire (rouge) et anti-horaire (bleu).</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>FIGURE 2.10 – *Haut* : Profils radiaux de la fréquence angulaire Ω/2π pour les modèles M05R1m (bleu), M07R1m (jaune), M09R1m (vert) et M11R1m (rouge), aux latitudes 0, 30, 45, et 60°. On ajoute de plus les données solaires obtenues par GONG en noir. Les courbes pointillées de chaque panneau représentent la vitesse angulaire d'une parcelle de fluide se déplaçant radialement au sein du modèle indiqué par la couleur, et pour laquelle son moment cinétique spécifique est conservé. La ligne en pointillés grise indique l'origine des ordonnées. (voir aussi DeRosa 2001) *Bas* : Profils méridiens de la rotation différentielle pour les 4 modèles considérés.</figcaption></figure>

## 2.4 Transferts de moment cinétique et établissement du profil de rotation différentielle

On commence par illustrer en Figure 2.10 la vitesse angulaire au sein de différents régimes. Les 4 modèles au taux de rotation solaire sont ainsi comparés aux observations du Soleil obtenues par le réseau GONG (Global Oscillations Network Group). Les courbes en pointillées correspondantes représentent la vitesse angulaire qu'aurait une parcelle de fluide isolée, conservant son moment cinétique spécifique *L*, et se déplaçant dans une direction purement radiale. On aurait alors $\Omega \propto \omega^{-2} \approx 10^2$ avec la constante de proportionnalité dépendant du modèle. On voit qu'en dehors de certaines régions spécifiques (la partie équatoriale profonde du modèle anti-solaire ou encore une profondeur moyenne de ses latitudes médianes), le fluide n'a pas tendance à conserver son moment cinétique spécifique localement. Il y a donc une redistribution du moment cinétique, dont nous proposons d'étudier les transferts dans cette partie. Cela nous permettra de comprendre le rôle de chaque processus physique dans la construction du profil de DR pour les différents régimes de rotation identifiés dans la section précédente.

<footnote n="10">10. $\omega = r \sin(\theta)$</footnote>
On considère de nouveau la conservation du moment cinétique en équation 2.24, tout en incluant cette fois-ci le transport de la circulation méridienne dans le membre de droite, tel que

$$ \bar{\rho} \frac{\partial \mathcal{L}}{\partial t} = \mathcal{T}. \qquad (2.26) $$

$\mathcal{T}$ est donc maintenant le couple global appliqué sur l'enveloppe convective. Les conditions aux limites (2.14), (2.15), (2.16), (2.17) et (2.18) que nous avons choisies ont la propriété de n'appliquer aucun couple extérieur sur les modèles. Si de plus on suppose un état statistiquement stationnaire, la conservation globale du moment cinétique est assurée et $\mathcal{T} = 0$. Ensuite, nous pouvons exprimer le couple global comme la divergence d'un flux tel que (Brun et al., 2004)

$$ \mathcal{T} = \nabla \cdot \mathcal{F} = \nabla \cdot (\mathcal{F}_r e_r + \mathcal{F}_\theta e_\theta) = \frac{1}{r^2} \frac{\partial(r^2 \mathcal{F}_r)}{\partial r} + \frac{1}{r \sin \theta} \frac{\partial(\sin \theta \mathcal{F}_\theta)}{\partial \theta} = 0. \qquad (2.27) $$

Les composantes radiale $\mathcal{F}_r$ et latitudinale $\mathcal{F}_\theta$ du flux de moment cinétique moyen se décomposent alors telles que

$$ \mathcal{F}_i = \mathcal{F}_i^{\text{VD}} + \mathcal{F}_i^{\text{RS}} + \mathcal{F}_i^{\text{MC}} + \mathcal{F}_i^{\text{MS}} + \mathcal{F}_i^{\text{MT}}, \qquad (2.28) $$

où les différents termes représentent des transports opérés par différents processus, qui sont respectivement de gauche à droite : les frottements par diffusion visqueuse (VD), le tenseur de Reynolds (RS), la circulation méridienne (MC), le tenseur de Maxwell (MS), ainsi que le couple magnétique à grande échelle (MT pour Maxwell Torques). Le tenseur de Reynolds quantifie les corrélations de vitesses à petite échelle, venant ici principalement de la turbulence convective. Le tenseur de Maxwell quantifie quant à lui, les corrélations magnétiques à petite échelle créées par la turbulence magnétisée, tandis que le couple magnétique MT représente la rétroaction appliquée par le champ magnétique grandes échelles (“magnetic web”) sur les écoulements. On note que le transport par la circulation méridienne peut se décomposer comme la somme de deux transports $\mathcal{F}_i^{\text{MC}} = \hat{v}_i(\hat{v}_\phi + \Omega r \sin \theta)$, qui sont opérés respectivement par le couple imposé par la rotation différentielle à grande échelle, et par la force de Coriolis. Nous référons le lecteur aux Équations 10 et 11 de la Section 2.7 pour le développement exact des autres termes de l'Équation 2.28. Enfin, nous intégrons ces composantes $\mathcal{F}_r$ et $\mathcal{F}_\theta$ de manière à obtenir

$$ \Phi_r(r) = \int_0^\pi \mathcal{F}_r(r, \theta) r^2 \sin \theta \, d\theta \quad \text{et} \quad \Phi_\theta(\theta) = \int_{r_{\text{bot}}}^{r_{\text{top}}} \mathcal{F}_\theta(r, \theta) r \sin \theta \, dr. \qquad (2.29) $$

$\Phi_r(r)$ représente le flux de moment cinétique total passant à travers la sphère de rayon $r$, tandis que $\Phi_\theta(\theta)$ représente le flux de moment cinétique total passant à travers le cône défini par la co-latitude $\theta$. Cela nous permet maintenant d'évaluer la direction et l'amplitude du transport par chacun des processus physiques composants l'Équation (2.28). Nous illustrons ces différentes composantes normalisées par $R_*^2$ en Figure 2.11, pour 3 modèles à $0.7 M_\odot$ représentatifs des 3 régimes rotationnels discutés en Section 2.3.

Nous remarquons tout d'abord que le transport total du moment cinétique, a atteint un bon équilibre (somme totale proche de zéro) dans la plupart des modèles, sur la balance radiale (à gauche) et latitudinale (à droite). Ces équilibres correspondent, de haut en bas, aux modèles M07Sm, M07R1m et M07R5m, possédant respectivement des profils de rotation anti-solaire, de type solaire et quasi-uniforme :

- Le panneau du haut est donc caractéristique des modèles à haut Rossby ($R_0f \ge 1.24$). On voit sur l'équilibre radial que la contribution du tenseur de Reynolds transporte le moment cinétique vers l'intérieur de l'étoile, ce qui maintient la circulation équatoriale rétrograde de ces modèles. Cette contribution est équilibrée par celles de la diffusion visqueuse (courbes oranges) et de la circulation méridienne (courbes bleus), tandis que les composantes magnétiques jouent un rôle négligeable dans le transfert radial. On les retrouve tous aussi négligeables sur l'équilibre latitudinal. Sur cet équilibre, on voit que le tenseur de Reynolds est la composante principale qui transporte le moment cinétique en direction de l'équateur
<figure><img src="image_58.png" /></figure>

FIGURE 2.11 – Équilibre des transports de moment cinétique, pour les cas et M07Sm (en haut), M07R1m (au milieu) et M07R5m (en bas). Le bilan radial (resp. latitudinal) est représenté sur la colonne de gauche (resp. droite). On retrouve dans le transport global (en trait plein noir) les différentes contributions : du tenseur de Reynolds en vert, de la MC en bleu, de la friction visqueuse en orange, du tenseur de Maxwell en violet, et du couple magnétique à grande échelle en pointillés noirs. Crédits : Brun et al. (2022)
(cette caractéristique sera la seule commune à tous les régimes rotationnels). On voit aussi que le transport visqueux contribue aussi à ce transfert vers l'équateur. Cela était attendu puisque la direction de ce transfert est toujours orienté à l'opposé du gradient de rotation $\nabla\Omega$. Ainsi, avec les pôles rapides et l'équateur lent de ces cas anti-solaires, les contraintes visqueuses ont tendance à ralentir les pôles et à accélérer l'équateur. Cependant, ces transferts vers l'équateur sont équilibrés par une forte composante de la circulation méridienne, reflétant le profil majoritairement unicellulaire de cette circulation dans chaque hémisphère (illustrée en Figure 2.9).

— Sur le panneau du milieu, nous représentons les transferts de moment cinétique menant aux profils de DR de type solaire de nos modèles. On note premièrement que si le tenseur de Reynolds est toujours une des contributions majoritaires du transport radial, celui-ci transporte maintenant le moment cinétique vers l'extérieur de l'étoile, et participe donc à la mise en place d'une rotation équatoriale prograde. Une nouvelle fois, les contraintes visqueuses équilibrent ce transport, mais ont cette fois-ci changé de sens par inversion radiale du profil de rotation <fnref n="11" />. De plus, on note que le transport par la circulation méridienne a chuté en amplitude et n'est plus de même signe sur l'ensemble de l'enveloppe. Ceci témoigne d'une meilleure balance entre les transport turbulent RS et visqueux VD, créant alors un pompage gyroscopique plus faible et menant à une profil de circulation méridienne multi-cellulaire, illustré en Figure 2.9. Concernant l'équilibre latitudinal, c'est maintenant le transport visqueux <fnref n="12" /> et une faible contribution du couple magnétique grandes échelles qui compensent le transport vers l'équateur imposé par le tenseur de Reynolds. Ces équilibres sont en assez bon accord avec le schéma global proposé par Brun (2004). Enfin, on ne note toujours pas de contribution magnétique significative sur l'équilibre radial des modèles de ce régime, à l'exception de M07R3m. La rétroaction magnétique au sein de ce modèle devient non négligeable (voir le profile de rotation en Figure 5 de la Section 2.7), ce qui en fait un modèle en transition vers le régime suivant.

— Nous arrivons finalement au régime atteint par les plus faibles $R_o$, dont des équilibres représentatifs sont illustrées sur le panneau du bas. Dans ce régime, la contrainte rotationnelle est forte et à l'origine d'un champ magnétique rétroagissant fortement sur la dynamique du système. On peut le voir sur l'équilibre des flux, où c'est maintenant le tenseur de Maxwell qui représente entièrement le transport en direction du centre de l'étoile. Il est principalement compensé par le tenseur de Reynolds, qui lui transporte encore ici le moment cinétique en direction de l'équateur. La part transportée par la circulation méridienne est quant à elle de nouveau significative, sans devenir une des composantes majoritaires de l'équilibre radial. Les profils de circulation méridienne sont redevenus en effet relativement moins complexe que pour le régime précédent, avec des cellules plus grandes et mieux définies (voir Figure 2.9). En ce qui concerne l'équilibre latitudinal, la tendance devient moins claire. En effet, le transport par le tenseur de Reynolds se fait toujours vers l'équateur, mais les autres contributions fluctuent significativement le long des latitudes. Cela peut être pré-senti en regardant le profil rotation différentielle quasi-entièrement lissé par la rétroaction du champ magnétique. Cet effet d'Ω-quenching augmente avec une diminution du Rossby et entraine donc une diminution du transport par diffusion visqueuse à mesure que le contraste de rotation $\nabla\Omega$ est amorti (voir par exemple le modèle intermédiaire M07R3m de Brun et al. 2022, disponible en Figure 9 de la Section 2.7).

En résumé, le transport du moment cinétique $\mathcal{L}$ par le tenseur de Reynolds est compensé par une combinaison de circulation méridienne, de contraintes visqueuses et de composantes magnétiques, dans la plupart des cas étudiés. Le magnétisme influence particulièrement les transferts pour les rotateurs les plus rapides, confirmant ainsi les résultats et interprétations des sections précédentes, sur la dynamique des écoulements moyens. Dans la section suivante, le but sera maintenant d'estimer la puissance soutenant ces différents écoulements. Pour

<footnote n="11">11. On a bien $\nabla \mathcal{F}_i^{VD} \propto -\nabla \Omega$, le transport visqueux en sens opposé du gradient de $\Omega$.</footnote>

<footnote n="12">12. Ayant aussi changé de sens par inversion du cisaillement latitudinal du passage de la DR anti-solaire à la DR solaire.</footnote>
cela, nous procédons à une analyse détaillée, afin de quantifier les différents échanges au sein des différents réservoirs énergétiques.

## 2.5 Bilan énergétique et intensité des écoulements moyens

Après avoir analysé et décrit les transferts de moment cinétique au sein des modèles dynamo d'étoiles G et K, nous procédons cette fois-ci à une analyse des échanges au sein des différents réservoirs d'énergies. Dans un premier temps, nous définissons ces derniers et quantifions l'énergie contenue dans chacun d'eux, en nous intéressant ici principalement aux composantes cinétiques. Par la suite, nous analyserons de la même manière les différents transferts alimentant ou diminuant l'énergie, en fonction du régime rotationnel caractérisé par le Rossby fluide $Ro_f$. Nous déléguons cette même analyse pour les composantes magnétiques au Chapitre 6, où nous analyserons spécifiquement le champ magnétique créé par ces dynamos.

Tout d'abord, nous illustrons sur la Figure 2.12 les différents réservoirs d'énergie ainsi que leurs échanges au sein d'une étoile. On y représente premièrement, le réservoir d'énergie interne, représentant l'agitation microscopique des particules du système. Il est initialement construit par contraction gravitationnelle à la naissance de l'étoile, puis est maintenu tout au long de sa vie par les réactions nucléaires en son centre. Par la suite, les différences de températures, si elles sont assez grandes, vont pouvoir mettre en mouvement la matière (voir Section ??). Il s'agit de la convection dans les enveloppes externes des étoiles de type solaire, qui alimente alors le réservoir d'énergie cinétique. Ce réservoir va contenir l'énergie cinétique des différents écoulements au sein de l'étoile, comprenant aussi les écoulements moyens (rotation différentielle et la circulation méridienne). Les mouvements ascendants de matière vont pouvoir transformer une partie de l'énergie cinétique en énergie potentielle, qui sera ensuite redonnée au réservoir cinétique par la descente des bulles convectives froides. Au cours de ces mouvements, des frottements visqueux vont avoir lieu, amenant lieu à un chauffage et donc un transfert d'énergie de l'énergie cinétique vers l'énergie interne. On note de plus que lors de son ascension, la transformation n'est pas parfaitement adiabatique et la bulle va aussi perdre de l'énergie par conduction thermique et rayonnement. Cette énergie reste alors au sein du réservoir d'énergie interne, sauf au niveau de la surface où elle crée la luminosité de l'étoile en s'échappant du système. Enfin, les mouvements du plasma vont donner naissance à une génération de champ magnétique par effet dynamo. De la même manière que la diffusion visqueuse ralentit les écoulements, la dissipation ohmique va diminuer ce champ magnétique, menant là aussi à un chauffage, et donc à un transfert de l'énergie magnétique vers l'énergie interne.

On définit la densité d'énergie cinétique, ainsi que la densité d'énergie magnétique, telles que

$$ KE = \frac{1}{2} \rho_{\text{tot}} (v_r^2 + v_{\theta}^2 + v_{\phi}^2), \qquad (2.30) $$

$$ ME = \frac{1}{8\pi} (B_r^2 + B_{\theta}^2 + B_{\phi}^2). \qquad (2.31) $$

On rappelle que l'on travaille ici avec les unités du code ASH qui sont définies dans le système CGS, avec un champ magnétique qui est donc en Gauss. Dans ce chapitre, nous nous intéressons en particulier au réservoir d'énergie cinétique, que nous décomposons de la manière suivante :

$$ DRKE = \frac{1}{2} \rho_{\text{tot}} \langle v_{\phi} \rangle_{\phi}^2, \qquad (2.32) $$

$$ MCKE = \frac{1}{2} \rho_{\text{tot}} \left( \langle v_r \rangle_{\phi}^2 + \langle v_{\theta} \rangle_{\phi}^2 \right), \qquad (2.33) $$

$$ FKE = \frac{1}{2} \rho_{\text{tot}} \left[ \left( v_r - \langle v_r \rangle_{\phi} \right)^2 + \left( v_{\theta} - \langle v_{\theta} \rangle_{\phi} \right)^2 + \left( v_{\phi} - \langle v_{\phi} \rangle_{\phi} \right)^2 \right], \qquad (2.34) $$
<figure><img src="image_59.png" /><figcaption>FIGURE 2.12 – Schéma simplifié des réservoirs d'énergie représentés en bleu, ainsi que leurs principaux échanges symbolisés par des flèches, au sein d'une étoile magnétisée en rotation. Nous indiquons à côté des flèches les forces dont le travail provoque un transfert d'énergie, dont les principales manifestations sont indiquées à côté des réservoirs cinétiques et magnétiques. Inspiré de Starr & Gilman (1966); Brandenburg et al. (1996); DeRosa et al. (2002); Rempel (2006)</figcaption></figure>

<table><caption>TABLE 2.3 – Récapitulatif des densités d'énergies cinétiques (KE) et magnétique (ME) moyennes pour chaque étoile simulée. On décompose la composante cinétique en ses différentes sous composantes : le mouvement toroïdal grande échelle dans le sens de la rotation (DRKE), le mouvement poloidal dans le plan méridien de l'étoile (MCKE) ainsi que les mouvements fluctuants non-axisymétriques (FKE). Adapté de Brun et al. (2022).</caption><thead><tr><th></th><th>KE<br>×10<sup>6</sup></th><th>DRKE (%KE)<br>×10<sup>6</sup></th><th>MCKE (%KE)<br>×10<sup>4</sup></th><th>FKE (%KE)<br>×10<sup>6</sup></th><th>ME (%KE)<br>×10<sup>6</sup></th></tr></thead><tbody><tr><td>M05Sm</td><td>9.3</td><td>2.1 (22.8%)</td><td>17.1 (1.8%)</td><td>7.0 (75.3%)</td><td>1.7 (17.7%)</td></tr><tr><td>M05R1m</td><td>20.3</td><td>17.0 (83.9%)</td><td>2.3 (0.1%)</td><td>3.3 (16.0%)</td><td>0.8 (3.9%)</td></tr><tr><td>M05R3m</td><td>74.3</td><td>71.4 (96.1%)</td><td>0.9 (0.01%)</td><td>2.9 (3.9%)</td><td>1.1 (1.5%)</td></tr><tr><td>M05R5m</td><td>2.7</td><td>1.2 (44.7%)</td><td>0.3 (0.1%)</td><td>1.5 (55.2%)</td><td>2.7 (98.2%)</td></tr><tr><td>M07Sm</td><td>4.0</td><td>1.2 (30.3%)</td><td>6.1 (1.5%)</td><td>2.7 (68.2%)</td><td>0.3 (7.6%)</td></tr><tr><td>M07R1m</td><td>4.8</td><td>3.3 (68.5%)</td><td>1.0 (0.2%)</td><td>1.5 (31.3%)</td><td>0.9 (18.5%)</td></tr><tr><td>M07R3m</td><td>3.5</td><td>2.3 (67.1%)</td><td>0.4 (0.1%)</td><td>1.1 (32.8%)</td><td>2.0 (57.8%)</td></tr><tr><td>M07R5m</td><td>1.1</td><td>0.3 (24.6%)</td><td>0.2 (0.1%)</td><td>0.8 (75.3%)</td><td>1.3 (116.2%)</td></tr><tr><td>M09Sm</td><td>3.2</td><td>0.4 (13.3%)</td><td>3.0 (0.9%)</td><td>2.7 (85.8%)</td><td>0.04 (1.4%)</td></tr><tr><td>M09R1m</td><td>4.0</td><td>2.1 (51.2%)</td><td>1.7 (0.4%)</td><td>1.9 (48.4%)</td><td>0.2 (5.7%)</td></tr><tr><td>M09R3m</td><td>11.9</td><td>10.2 (85.5%)</td><td>0.8 (0.1%)</td><td>1.7 (14.4%)</td><td>0.5 (4.4%)</td></tr><tr><td>M09R5m</td><td>2.6</td><td>1.6 (62.4%)</td><td>0.3 (0.1%)</td><td>1.0 (37.5%)</td><td>1.3 (49.3%)</td></tr><tr><td>M11R1m</td><td>2.0</td><td>0.4 (20.5%)</td><td>2.0 (1.0%)</td><td>1.5 (78.4%)</td><td>0.2 (12.1%)</td></tr><tr><td>M11R3m</td><td>3.9</td><td>2.9 (74.0%)</td><td>0.6 (0.2%)</td><td>1.0 (25.9%)</td><td>0.3 (8.8%)</td></tr><tr><td>M11R5m</td><td>2.3</td><td>1.0 (43.3%)</td><td>8.9 (3.9%)</td><td>1.2 (52.7%)</td><td>4.4 (192.3%)</td></tr></tbody></table>

où DRKE est l'énergie cinétique des écoulements moyens longitudinaux, associée à la rotation différentielle. MCKE est l'énergie cinétique des écoulements moyens méridiens, et donc représentatifs de la circulation méridienne. Enfin, on définit par complémentarité la troisième énergie cinétique, représentant les écoulements fluctuants non-axisymétriques, principalement dus à la convection. Après avoir moyenné sur le volume simulé les différentes composantes, nous listons leurs valeurs pour les 15 modèles en Table 2.3.

Premièrement, on voit que KE diminue en augmentant la masse du modèle pour un même taux de rotation. Cela est dû à la densité moyenne de la zone convective $\rho_{CZ}$ de plus en plus faible, à mesure que l'épaisseur de cette enveloppe diminue (voir Table 2.1). Cela est cependant en partie compensé par l'augmentation de la luminosité avec la masse, augmentant alors les vitesses convectives (voir Section 2.2.3.2), ce qui maintient l'amplitude de KE autour de $10^6 - 10^7$ erg/cm$^3$. Si l'on considère maintenant des masses équivalentes, on voit
que KE augmente avec le taux de rotation, avant de diminuer de nouveau pour les rotateurs les plus rapides. Cette évolution est corrélée avec celle de DRKE, qui augmente avec le contraste du profil de rotation $|\Delta\Omega|$ (voir Table 2.2). DRKE peut alors représenter jusqu'à ~ 96% du réservoir total KE. La diminution qui s'ensuit lorsque l'on considère les rotateurs les plus rapides vient donc de la rétroaction du champ magnétique créé par la dynamo de ces modèles. Cet effet appelé Ω-quenching, qui lisse le profil de rotation des modèles aux plus faibles $R_{of}$, diminuant le contraste $|\Delta\Omega|$. Concernant MCKE, on note que cette énergie représente une part très minoritaire de KE dans la majorité des modèles, à l'exception de M05Sm qui présente une circulation méridienne hémisphérique particulièrement forte (voir Figure 2.9). En parallèle, la partie non-axisymétrique FKE domine proportionnellement le réservoir KE dans les cas à rotation lente. Cela s'explique par une faible influence de la rotation sur la convection de ces modèles. Pour tous les autres modèles, on voit que l'amplitude de FKE reste majoritairement autour de $10^6$ erg/cm$^3$. Pour les Rossby intermédiaires, la proportion de FKE dans le réservoir total KE diminue alors par augmentation de DRKE pour les profils de rotation de type solaire. Cette proportion de FKE/KE remonte ensuite pour les modèles où la rotation est la plus rapide et où les $R_{of}$ sont les plus faibles, puisque DRKE diminue dans ces modèles par action de l'Ω-quenching. FKE conservant son ordre de grandeur, la diminution d'amplitude de KE vient de la baisse significative de DRKE. On note cependant que FKE diminue progressivement en augmentant la rotation, notamment à cause d'une rétroaction de la force de Lorentz, et potentiellement une diminution du niveau de super-criticalité des modèles à rotation rapide dans l'étude (Takehiro et al., 2020). Finalement, nous indiquons de la même manière l'énergie magnétique globale ME dans la dernière colonne du tableau. Sa comparaison avec KE nous permet de quantifier l'amplitude de l'instabilité dynamo transformant KE en ME, et ainsi de voir où la rétroaction du champ magnétique sur les écoulements sera d'autant plus ressentie par la dynamique. Ainsi, on voit que pour les cas à rotation lente, ME représente une part minoritaire de KE. Lorsque nous diminuons $R_{of}$, cette proportion augmente pour atteindre un équilibre entre le terme d'advection et la force de Lorentz au sein de l'Équation (2.3), appelé régime d'équipartition. Cela se fait par une combinaison d'échanges non-linéaires entre ME et KE, qui augmentent tous les deux. De plus, il est intéressant de noter que certains des cas les plus rapides montrent des ME plus grandes que l'énergie KE des écoulements à l'origine de ce magnétisme. Ces modèles sont alors dans un état dit de super-équipartition globale, et nous montrent que l'énergie cinétique disponible dans l'enveloppe convective n'est pas la valeur maximale que la ME peut atteindre. La force de Lorentz n'est alors plus en équilibre avec le terme de transport advectif, mais avec celui de la force de Coriolis qui augmente avec le taux de rotation : c'est ce qu'on appelle l'état magnétostrophique (voir aussi Nataf et al. 2008; Brun et al. 2015). Les différentes composantes de ME seront étudiées dans une analyse des mécanismes dynamos au Chapitre 6.

Après avoir vu l'évolution des réservoirs d'énergie liés aux écoulements moyens DRKE et MCKE, on s'intéresse aux processus les soutenant. En suivant les développements proposés par Starr & Gilman (1966), Brandenburg et al. (1996), DeRosa et al. (2002) et Rempel (2006), on peut complexifier la Figure 2.12. Ainsi, on obtient la Figure 23 de l'article en Section 2.7, où le réservoir d'énergie cinétique est décomposé en DRKE et MCKE (FKE étant alors compté avec l'énergie potentielle et l'énergie interne). Les développements analytiques de l'annexe A du même article (voir les équations A.8 à A.23) nous permettent de déduire les budgets d'énergie des écoulements moyens (i.e. axisymétriques).

$$ \partial_t \text{DRKE} = Q_{\text{RS}}^{\text{DR}} - Q_{\Omega} - Q_C - Q_{\text{MS}}^{\text{DR}} - Q_{\nu}^{\text{DR}} + C^{\text{DR}} - S^{\text{DR}}, \quad (2.35) $$

est le bilan énergétique de la rotation différentielle, où les termes représentent respectivement le tenseur de Reynolds (RS), l'effet-Ω, la force de Coriolis (C), le tenseur de Maxwell (MS), le chauffage visqueux, le terme de courbure, ainsi que les effets de bords. Le développement de chacun de ces termes est détaillé en Annexe A de la Section 2.7. On développe ensuite le bilan énergétique de la circulation méridienne, tel que

$$ \partial_t \text{MCKE} = Q_{\text{RS}}^{\text{MC}} + Q_{\text{MS}}^{\text{MC}} + Q_{\nabla P} + Q_C + Q_{\text{TM}}^{\text{MC}} - Q_{\text{PM}}^{\text{MC}} - Q_b - Q_{\nu}^{\text{MC}} + C^{\text{MC}} - S^{\text{MC}}, \quad (2.36) $$

où les termes représentent respectivement le tenseur de Reynolds (RS), le tenseur de Maxwell (MS), le travail
<figure><img src="image_60.png" /><figcaption>FIGURE 2.13 – Équilibres énergétiques pour la rotation différentielle (gauche) et la circulation méridienne (droite), sur chacun des 15 modèles dynamo, en fonction du nombre de Rossby fluide. On normalise les transferts pour une comparaison directe à la luminosité de l'étoile $L_*$. Crédits : Brun et al. (2022)</figcaption></figure>

par le gradient de pression ($\nabla P$), la force de Coriolis (C), le cisaillement par cette circulation du champ magnétique toroïdal (TM) en champ magnétique poloidal (PM), la flottaison, le chauffage visqueux, le terme de courbure, ainsi que les effets de bords. Nous avons calculé chacun de ces termes et les avons normalisés à la luminosité stellaire $L_*$ pour les différents modèles, afin d'illustrer les bilans énergétiques en Figure 2.13. Les équilibres 2.35 et 2.36 ont été moyennées temporellement sur une centaine de temps de retournement convectif afin d'obtenir un bon équilibre de chacune d'elles.

Tout d'abord, on remarque sur l'équilibre énergétique de DRKE (panneau de gauche) que la contribution des différents processus est cohérente avec ce qui est observé dans l'analyse des transferts de moment cinétique en section précédente. Le profil de rotation différentielle est majoritairement entretenu par la composante radiale du tenseur de Reynolds dans tous les modèles. Pour les profils anti-solaires à $Ro_f > 1$, ce maintien est équilibré par un transfert négatif de la composante latitudinale du même tenseur, ainsi que de la diffusion visqueuse. Cette diffusion devient alors progressivement le principal transfert au détriment de DRKE, pour les $Ro_f$ intermédiaires. C'est dans ce régime intermédiaire qu'on observe le maintien le plus fort de la rotation différentielle, adoptant alors un profil de type solaire. Les puissances impliquées peuvent atteindre jusqu'à ~30% de la luminosité $L_*$ pour les modèles autour de $log_{10}Ro_f \sim 0.2$, créant alors de forts contrastes de rotation $|\Delta\Omega|$. Pour des Rossby plus élevés, la dynamique n'est plus assez influencée par la rotation pour maintenir aussi fortement DRKE. Pour des Rossby plus faibles, la rotation est si forte que les composantes magnétiques ne sont plus négligeables dans l'équilibre énergétique (composantes bleues et rouges). La composante bleue représente le cisaillement des lignes de champ magnétique poloïdales par la rotation différentielle, appelé effet-$\Omega$ (voir Section 1.1.6), ce qui extrait une partie de la DRKE. La composante rouge représente la rétroaction du tenseur de Maxwell, venant alors diminuer le contraste de DR $|\Delta\Omega|$ à grande échelle, et donc les cisaillements (c'est-à-dire l'$\Omega$-quenching). Enfin, pour les nombres de Rossby les plus faibles, la rétroaction du magnétisme sur les écoulements est si forte qu'elle impacte les mouvements convectifs à plus petite échelle, eux-mêmes à l'origine du tenseur de Reynolds (contributions oranges) alimentant la DRKE. On a donc une diminution des puissances fournies dans le maintien de DRKE.

Ensuite, le panneau de droite montre que l'équilibre énergétique MCKE de la circulation méridienne se fait généralement par un équilibre entre les contributions du gradient de pression, de la flottaison et de la force de Coriolis. On note que cette dernière contribution est généralement négative, signifiant que la circulation méridienne est déviée par la force de Coriolis, contribuant alors à la rotation différentielle. De plus, bien que sa proportion dans l'équilibre totale évolue, on note que cette contribution garde le même ordre de grandeur pour la majorité des modèles (quelques % de $L_*$). La puissance totale impliquée pour la mise en place de la circulation méridienne augmente avec le nombre de Rossby, ce qui semble cohérent avec les profils en Figure 2.9. Finalement, on voit que les contributions magnétiques restent toujours minoritaires dans l'équilibre MCKE,
voir totalement négligeables dans certains modèles. En ce sens, le champ magnétique n'a pas d'impact direct sur la circulation méridienne, même pour les cas en rotation rapide. Cette circulation s'adapte morphologiquement au profil de rotation, qui lui est bien impacté par le magnétisme.

## 2.6 Discussion et perspectives

Au cours de cette étude, nous avons simulé 15 étoiles, dont les masses varient de 0.5 à 1.1$M_\odot$, et pour des taux de rotation débutant à des valeurs inférieures au taux solaire, et allant jusqu'à $5\Omega_\odot$. En continuité des études Brun et al. (2017) et Varela et al. (2016), on cherche à comprendre l'impact de la rotation et de la structure interne sur les écoulements moyens de ces solutions dynamos.

Premièrement, nous avons vu que la force de Lorentz peut avoir une rétroaction forte sur la rotation différentielle. Cela limite l'exposant de la loi de puissance à une relation du type $\Delta\Omega \propto \Omega_*^{0.46}$. Les différents profils de rotation se regroupent en 3 morphologies principales, correspondant à différents régimes rotationnels caractérisés par le nombre de Rossby. Ces différents régimes se retrouvent aussi dans les profils de circulation méridienne établis dans ces simulations. Nous proposons par la suite une caractérisation de ces régimes par une analyse des transferts de moment cinétique, puis nous réalisons une décomposition énergétique, afin d'en analyser les différents transferts au sein des réservoirs d'énergie cinétique grandes échelles.

En résumé, nous avons vu que la rétroaction du champ magnétique pouvait être très significative, allant même jusqu'à lisser le profil de rotation différentielle pour les valeurs de Rossby les plus faibles. Pour un régime de Rossby intermédiaire, où le Soleil est très probablement situé, la rotation différentielle peut être alimentée par une puissance allant jusqu'à ~30% de la luminosité de l'étoile. Enfin, des profils de rotation anti-solaires apparaissent dans un régime à hauts nombres de Rossby.

Par la suite, nous verrons aux Chapitres 5 et 6 que les modèles anti-solaires ont la particularité de perdre l'aspect cyclique de leur dynamo, les rendant alors d'autant plus important pour la compréhension de l'évolution des dynamos de type solaire. Nous étudierons dans ce contexte le potentiel caractère cyclique de ces solutions, et proposerons au Chapitre 6 une analyse détaillée du rôle des composantes grandes échelles de ME dans cette activité magnétique.

Bien que des changements de comportement semblent avoir lieu pour le régime anti-solaire (voir Figure 2.8), le nombre de modèles ou d'observations couvrant ce régime n'est pas suffisant pour confirmer statistiquement cette tendance. Notre objectif pour le prochain chapitre sera donc de proposer une liste de candidats de type solaire sur la séquence principale susceptibles de posséder un profil de rotation différentielle anti-solaire, en se basant sur nos attentes théoriques du nombre de Rossby, afin d'étudier et valider les changements de comportements à la transition $R_o \sim 1$.

## 2.7 Article publié : Powering Stellar Magnetism : Energy Transfers in Cyclic Dynamos of Sun-like Stars
# Powering Stellar Magnetism: Energy Transfers in Cyclic Dynamos of Sun-like Stars

Allan Sacha Brun<fnref n="1" />, Antoine Strugarek<fnref n="1" />, Quentin Noraz<fnref n="1" />, Barbara Perri<fnref n="1" />,<fnref n="2" />, Jacobo Varela<fnref n="1" />,<fnref n="3" />, Kyle Augustson<fnref n="1" />, Paul Charbonneau<fnref n="4" />, and Juri Toomre<fnref n="5" />

Received 2021 June 30; revised 2021 November 30; accepted 2021 December 23; published 2022 February 9

## Abstract

We use the anelastic spherical harmonic code to model the convective dynamo of solar-type stars. Based on a series of 15 3D MHD simulations spanning four 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 \propto (|\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$ (antisolar-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 also find efficient energy transfer between reservoirs, leading to the conversion of several percent of the star's luminosity into magnetic energy that could provide enough free energy to sustain intense eruptive behavior at the star's surface. 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}$).

*Unified Astronomy Thesaurus concepts: Solar dynamo (2001); Solar magnetic fields (1503); Stellar magnetic fields (1610); Stellar rotation (1629); Solar differential rotation (1996); Magnetohydrodynamics (1964); Stellar convection envelopes (299); Magnetohydrodynamical simulations (1966); Solar analogs (1941); K stars (878); G stars (558)*

## 1. Introduction

Sun-like stars go through various magnetic activity phases in their lives. From young very active T Tauri 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 the source of key temporal variabilities like the Schwabe 11 yr or Gleissberg 90 yr magnetic cycles in the Sun and at the heart of a complex feedback loop between stellar magnetism and rotation via wind braking and the loss of mass and angular momentum by the star (Matt et al. 2015; Brun & Browning 2017; Vidotto 2021). It is also key in providing the free energy reservoir needed to power eruptive events such as flares or coronal mass ejections (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 yr 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 & Aallot 2019) to connect the spectral class 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 $\langle B_V \rangle \propto Ro_s^{-1.38 \pm 0.14}$ 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. (2019b) have revisited this trend and found a similar result with $\langle B_V \rangle \propto Ro_s^{-1.40 \pm 0.10}$. It was also proposed in Petit et al. (2008) and later by See et al. (2015) that the toroidal magnetic field dominates over the 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 is 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 et al. 2007; Hall 2008) can be found to be either irregular with no obvious cyclic activity or to possess

<figure><img src="image_61.png" /></figure>

<footnote n="1">¹Département d'Astrophysique/AIM, CEA/IRFU, CNRS/INSU, Univ. Paris-Saclay & Univ. de Paris, F-91191 Gif-sur-Yvette, France; sacha.brun@cea.fr</footnote>

<footnote n="2">²Centre for Mathematical Plasma Astrophysics, KU Leuven, Celestijnenlaan 200b-box 2400, B-3001 Leuven, Belgium</footnote>

<footnote n="3">³Universidad Carlos III de Madrid, Leganés, E-28911, Spain</footnote>

<footnote n="4">⁴Département de physique, Université de Montréal, C.P. 6128 Succ. Centre-Ville, Montréal, QC H3C-3J7, Canada</footnote>

<footnote n="5">⁵JILA and Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309-0440, USA</footnote>
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 $\Omega(t) \propto t^{-0.5}$ and gyrochronology trend (Skumanich 1972; Barnes 2003, 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. 2018, 2019; do Nascimento et al. 2020) while others do not (Metcalfe et al. 2016; 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 versus long-term monitoring of individual stars) used, as each has rotation rate and age determination that sometimes differ significantly (Lorenzo-Oliveira et al. 2016; do Nascimento et al. 2020). Another alternative would be that stars temporarily stop spinning down before starting again (Curtis et al. 2020; Spada & Lanzafame 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 (DR) realized 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 $\Omega$ effect (e.g., stretching of the poloidal magnetic field lines by large-scale shear). Hence, the role of 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 DR. The combination of all these observations on stellar rotation and magnetism helps constrain the trends linking rotation with

stellar DR and magnetic activity. Various analyses of stellar DR revealed different dependencies between DR and star's rotation ($\Delta\Omega \propto \Omega''$), with $n$ varying between 0.15 and 0.7 (Barnes et al. 2005; Reiners 2006; Reinhold et al. 2013). 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 DR 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 ($\Delta\Omega \propto T_{\text{eff}}^{8.92}$, Barnes et al. 2005; Reinhold et al. 2013; and $\Delta\Omega \propto T_{\text{eff}}^{8.6}$, 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 antisolar DR 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 DR and stellar magnetism in the convection zone of solar-like stars (Glatzmaier & Gilman 1982; Miesch et al. 2000, 2006; Brun et al. 2004, 2011; Brown et al. 2008, 2010; Ghizaru et al. 2010; Käpylä et al. 2011, 2014; Gastine 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; Brown et al. 2011; Racine et al. 2011; Augustson et al. 2013, 2015; Käpylä et al. 2013; Nelson et al. 2013; Beaudoin et al. 2016; Guerrero et al. 2016, 2019; Strugarek et al. 2017, 2018; Viviani et al. 2018, 2019; Warnecke 2018; 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. 2016, 2019; Kä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 Brun et al. (2017) with the MHD anelastic spherical harmonic code (ASH) (Brun et al. 2004). 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 DR 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 give a quick overview of one of the dynamo solutions emphasizing the main properties of a cyclic solution. In Section 4 we discuss the 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 DR contrast $\Delta\Omega$. In Section 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 $\alpha$-$\Omega$ classification. We further expand our data set with the 17 simulations published previously Strugarek et al. (2017, 2018) with the Eulag-MHD code (Smolarkiewicz & Charbonneau 2013), in order to improve the statistics. In Section 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 Section 7 we reflect on our findings in an astrophysical context, comparing our results with recent observational results and then conclude.

## 2. 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.

### 2.1. 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 magnetoacoustic 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_{lm}(\theta, \phi)$ for their horizontal structure and in Chebyshev polynomials $T_n(r)$ for their radial structure (Brun et al. 2004). The density ($\rho$), entropy, pressure, and temperature ($T$) are linearized about the spherically symmetric background values, denoted by the symbol ($^\circ$). The equations solved by ASH are (Brun et al. 2004)

$$ \nabla \cdot \hat{\rho}v = 0, \quad \nabla \cdot \mathbf{B} = 0 \qquad (1) $$

$$ \begin{aligned} \hat{\rho} \left( \frac{\partial v}{\partial t} + (v \cdot \nabla) v + 2\Omega_{*} \times v \right) \\ = -\nabla P + \rho g + \frac{1}{4\pi} (\nabla \times B) \times B + \nabla \cdot \mathbf{D} \end{aligned} \qquad (2) $$

$$ \hat{\rho} \hat{T} \frac{\partial S}{\partial t} = -\hat{\rho} \hat{T} v \cdot \nabla (\hat{S} + S) - \nabla \cdot q + \Phi_d + \hat{\rho} \epsilon \qquad (3) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times [\mathbf{v} \times \mathbf{B} - \eta \nabla \times \mathbf{B}], \qquad (4) $$

with the velocity field $v = v_r \hat{e}_r + v_\theta \hat{e}_\theta + v_\varphi \hat{e}_\varphi$, the magnetic field $B = B_r \hat{e}_r + B_\theta \hat{e}_\theta + B_\varphi \hat{e}_\varphi$, the angular velocity in the rotation frame $\Omega_*$ = $\Omega_* \hat{e}_z$, $\hat{e}_z$ the unit vector along the rotation axis, and $g$ the magnitude of the gravitational acceleration. A volumetric heating term $\hat{\rho}\epsilon$ is also taken into account to approximate the generation of energy by nuclear reactions in the stellar core. The nuclear reactions are modeled very simply by assuming that $\epsilon = \epsilon_0 \hat{T}^{n_c}$. By enforcing that the integrated luminosity of the star matches its known surface value, we can determine $\epsilon_0$ and $n_c$ as listed in Table 7 of Brun et al. (2017). Note that only the low-mass star series of models (e.g., 0.5 and 0.7 $M_\odot$) require that heating source term, because their computational domain includes a portion of the nuclear energy generation core.

The diffusion tensor $D$ and the dissipative term $\Phi_d$ are defined as

$$ \begin{aligned} D_{ij} &= 2\hat{\rho}\nu \left[ e_{ij} - \frac{1}{3} \nabla \cdot \mathbf{v} \delta_{ij} \right], \\ \Phi_d &= 2\hat{\rho}\nu \left[ e_{ij} e_{ij} - \frac{1}{3} (\nabla \cdot \mathbf{v})^2 \right] + \frac{4\pi\eta}{c^2} \mathbf{J}^2, \end{aligned} $$

with $e_{ij}$ the stress tensor and $\mathbf{J} = c/4\pi \nabla \times \mathbf{B}$ the current density. The energy flux $\mathbf{q}$ is the sum of a radiation flux and of a turbulent entropy diffusion flux:

$$ \mathbf{q} = \kappa_r \hat{\rho} c_p \nabla (\hat{T} + T) + \kappa \hat{\rho} \hat{T} \nabla S + \kappa_0 \hat{\rho} \hat{T} \frac{\partial \hat{S}}{\partial r} \hat{e}_r, $$
<figure><img src="image_62.png" /><figcaption>**Figure 1.** Typical radial profile of kinematic viscosity $\nu$ used in this study, here for case M11R3m. Profiles of $\kappa$ and $\eta$ are the same, but their amplitude depends respectively for each case on the chosen Prandtl and magnetic Prandtl numbers (see Table 3).</figcaption></figure>

with $\nu$, $\kappa$, and $\eta$ the effective eddy diffusivities of the momentum, heat, and magnetic field transport, $\kappa_r$, the atomic radiation diffusion coefficient, $\kappa_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 fields still exist 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 $\nu$, $\kappa$, and $\eta$ represent momentum, heat, and magnetic field transport by motions that 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, $\nu$, $\kappa$, and $\eta$ have the following profile:

$$ \nu(r) = \nu_{\text{bot}} + \nu_{\text{top}} f_{\text{step}}(r), $$

where

$$ \begin{aligned} f_{\text{step}}(r) &= (\hat{\rho}/\hat{\rho}_{\text{top}})^{\alpha}[1 - \beta]f(r), \\ f(r) &= 0.5(\tanh((r - r_t)/\sigma_t) + 1), \\ \beta &= \nu_{\text{bot}}/\nu_{\text{top}} = 10^{-3}, \end{aligned} $$

and with $\nu$, $\nu_{\text{bot}}$ and $\nu_{\text{top}}$ in cm$^2$ s$^{-1}$ and $r_t$ and $\sigma_t$ in cm, as provided in Table 7 of Brun et al. (2017), and $\alpha$ is -0.5 for all cases. All models assumed a Prandtl number $Pr = \nu/\kappa$ of 0.25, so that $\kappa$ can be directly obtained from the amplitude and profile of $\nu$. The magnetic Prandtl number $P_m = \nu/\eta$ is equal to 1 or 2 depending on the case considered (see Table 3), so that $\eta$ can as well be deduced from $\nu$. These tapered profiles are chosen in order to take into account the much smaller subgrid-scale transport expected in the stably stratified radiative interior. A representative profile is shown in Figure 1. Their amplitudes are adapted for each rotation rate and stellar mass 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 $\kappa_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 stream function formalism (Brun et al. 2004):

$$ \hat{\rho}\nu = \nabla \times \nabla \times (W\hat{e}_r) + \nabla \times (Z\hat{e}_r), \quad (5) $$

$$ \mathbf{B} = \nabla \times \nabla \times (C\hat{e}_r) + \nabla \times (A\hat{e}_r). \quad (6) $$

A perfect ideal gas equation is used for the mean state and the fluctuations are linearized:

$$ \begin{aligned} \hat{P} &= (\gamma - 1)c_p \hat{\rho} \hat{T} / \gamma \\ \rho / \hat{\rho} &= P / \hat{P} - T / \hat{T} = P / \gamma \hat{P} - S / c_p \end{aligned} $$

with $\gamma = 5/3$ the adiabatic exponent.

The anelastic MHD system of equations requires 12 boundary conditions (BCs). We use an impenetrable and stress-free BCs at the top and bottom of the domain, i.e.,

$$ v_r = \frac{\partial}{\partial r} \left( \frac{v_{\theta}}{r} \right) = \frac{\partial}{\partial r} \left( \frac{v_{\phi}}{r} \right) = 0. $$

Magnetic BCs are perfectly conducting at the lower radial boundary and the magnetic field matches a potential field in the upper boundary: $B_r|_{r_{\text{bot}}} = \frac{\partial}{\partial r}\left(\frac{B_{\theta}}{r}\right)|_{r_{\text{bot}}} = \frac{\partial}{\partial r}\left(\frac{B_{\phi}}{r}\right)|_{r_{\text{bot}}} = 0$ and $B_r|_{r_{\text{top}}} = \nabla\Psi \Rightarrow \Delta\Psi = 0$, with $r_{\text{top}}$ and $r_{\text{bot}}$, respectively, the radius of the top and bottom of the numerical domain and $r_{\text{bcz}}$ that of the base of the convective layer (see Table 1).

Finally, we maintain the entropy flux at the top and bottom. Keeping the values of $d\hat{S}/dr|_{r_{\text{top}},r_{\text{bot}}}$ fixed at all times in the simulations further implies that the fluctuating $dS/dr$ is set to zero at both BCs.

## 2.2. 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 ($\hat{S}$), density ($\hat{\rho}$), temperature ($\hat{T}$) derived from a one-dimensional solar-structure model CESAM (Morel 1997; Brun et al. 2002). Our starting point is the G- and K-star rotating convective 3D models published in Brun et al. (2017) (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 four mass bins and four rotation rates. We
<figure><img src="image_63.png" /><figcaption>**Figure 2.** (a) Radial dependence of the mean entropy gradient for case M11R3m. The region of the tachocline is shown in the inset figure. (b) Time and horizontally averaged radial energy fluxes as luminosities (normalized to the star luminosity) for case M11R3m. The solid line is the total flux, the long-dashed line the radiative flux, the dashed-triple-dotted line the enthalpy flux, the dotted line the conductive entropy flux, the thick dashed-dotted line the kinetic energy, the dashed line the viscous diffusion flux, and the thick dashed line the Poynting flux.</figcaption></figure>

have models for stellar masses 0.5, 0.7, 0.9, and 1.1 $M_\odot$ and rotation rates ranging from 1/8 to 5 times the solar rotation rates. In keeping with the naming nomenclature of Brun et al. (2017), 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/antisolar ($x=s$) and prograde ($x=R$) DR models (except model M11R1m, which is also antisolar) 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_\theta, N_\phi$) 769 × 256 × 512 except for few cases in the M09m and M11m series rotating at $\Omega_* = 3$ or $5\Omega_\odot$, where $N_\theta$ is 512 and $N_\phi$ is 1024.

In Figure 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 flux balance 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. Because we are in this study mostly interested in the magnetic state of our simulation, we refer the reader to the following multidimensional studies of internal gravity wave generation in solar-like stars (Rogers & Glatzmaier 2006; Brun et al. 2011; Alvan et al. 2014, 2015). Turning to the radial flux balance, we note that the enthalpy flux (dashed-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 dashed) and unresolved near the top (dotted line)) carry the stellar luminosity at each end of the domain. We note an inward kinetic energy flux (dashed-dotted 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, which is compensated by a local increase of the radiative flux, such as to reach a satisfactory radial energy balance and thermal equilibrium.

## 2.3. G- and K-star 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 genuinely established tachocline, defined as the transition between DR 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 solar-like stars as reported in simulations performed by several authors (Glatzmaier 1985b; Browning et al. 2006; 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_{\hat{\rho}_{bcz}} = \ln(\hat{\rho}_{out}/\hat{\rho}_{bcz})$ and $N_{\hat{\rho}_{tot}} = \ln(\hat{\rho}_{out}/\hat{\rho}_{in})$. For the M05 model, $N_{\hat{\rho}_{bcz}} = 3.25$ and $N_{\hat{\rho}_{tot}} = 4.70$; M07 model, $N_{\hat{\rho}_{bcz}} = 3.48$ and $N_{\hat{\rho}_{tot}} = 5.78$; M09 model, $N_{\hat{\rho}_{bcz}} = 3.31$ and $N_{\hat{\rho}_{tot}} = 5.99$; and M11 model, $N_{\hat{\rho}_{bcz}} = 3.28$ and $N_{\hat{\rho}_{tot}} = 5.60$. The convective flows at the middle of the convective envelope vary from $5 \text{ m s}^{-1}$ up to about $300 \text{ m s}^{-1}$ in our sample, and the convective turnover time from 7 (case M11R1m) to 222 (case M05R3m) days. The surface DR between the equator and latitude $60^\circ$ varies from $-102$ to $+278$ nHz in our sample, and we will study its maintenance in detail in Section 4.1. In the middle of the convection zone, the rms magnetic field typically varies from 1 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
<table><caption>**Table 1**
Global Properties on the Main Sequence of the Four Stars Used in Our ASH Dynamo Models</caption><thead><tr><th>Mass<br>(M<sub>☉</sub>)</th><th>Radius<br>(R<sub>☉</sub>)</th><th>L<sub>*</sub><br>(L<sub>☉</sub>)</th><th>T<sub>eff</sub><br>(K)</th><th>Sp. T.</th><th>M<sub>bcz</sub><br>(M<sub>☉</sub>, M<sub>*</sub>)</th><th>r<sub>bcz</sub><br>(R<sub>☉</sub>, R<sub>*</sub>)</th><th>Î<sub>bcz</sub><br>(K)</th><th>ρ<sub>bcz</sub><br>(g cm<sup>-3</sup>)</th><th>Δ<sub>cz</sub>ρ<br>...</th><th>Δ<sub>f</sub>ρ<br>...</th><th>r<sub>bot</sub><br>(R<sub>*</sub>)</th><th>r<sub>top</sub><br>(R<sub>*</sub>)</th></tr></thead><tbody><tr><td>0.5</td><td>0.44</td><td>0.046</td><td>4030</td><td>K7</td><td>0.18, 0.36</td><td>0.25,0.56</td><td>4.3 × 10<sup>6</sup></td><td>14.0</td><td>42</td><td>193</td><td>0.13</td><td>0.95</td></tr><tr><td>0.7</td><td>0.64</td><td>0.15</td><td>4500</td><td>K4/K5</td><td>0.079, 0.11</td><td>0.42,0.66</td><td>3.0 × 10<sup>6</sup></td><td>2.1</td><td>50</td><td>605</td><td>0.32</td><td>0.97</td></tr><tr><td>0.9</td><td>0.85</td><td>0.55</td><td>5390</td><td>G8</td><td>0.042, 0.046</td><td>0.59,0.69</td><td>2.6 × 10<sup>6</sup></td><td>0.51</td><td>67</td><td>1013</td><td>0.38</td><td>0.97</td></tr><tr><td>1.1</td><td>1.23</td><td>1.79</td><td>6030</td><td>G0</td><td>0.011, 0.010</td><td>0.92,0.75</td><td>1.6 × 10<sup>6</sup></td><td>0.048</td><td>81</td><td>830</td><td>0.5</td><td>0.97</td></tr></tbody></table>

Note. All the listed values were computed with the CESAM stellar evolution code (Morel 1997). We adopt $M_{☉} = 1.989 \times 10^{33}$ g, $R_{☉} = 6.9599 \times 10^{10}$ cm, and $L_{☉} = 3.846 \times 10^{33}$ erg s$^{-1}$. The density ratios $\Delta_{cz}\hat{\rho}$ and $\Delta_f\hat{\rho}$ 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><caption>**Table 2**
Models Dimensional Characteristics</caption><thead><tr><th></th><th>Ω<sub>*</sub><br>(Ω<sub>☉</sub>)</th><th>Ṽ<sub>r</sub><br>(m s<sup>-1</sup>)</th><th>Ṽ<sub>θ</sub><br>(m s<sup>-1</sup>)</th><th>Ṽ<sub>φ</sub><br>(m s<sup>-1</sup>)</th><th>ΔΩ<br>(nHz)</th><th>B̃<sub>r</sub><br>(G)</th><th>B̃<sub>θ</sub><br>(G)</th><th>B̃<sub>φ</sub><br>(G)</th><th>τ<sub>c</sub><br>(days)</th><th>τ<sub>ν</sub><br>(yr)</th><th>τ<sub>κ</sub><br>(yr)</th><th>τ<sub>η</sub><br>(yr)</th></tr><tr><td></td><td></td><td></td><td></td><td></td><td colspan="6">MHD (HD)</td></tr></thead><tbody><tr><td>M05Sm</td><td>1/8</td><td>13.52</td><td>12.19</td><td>29.35</td><td>-23 (-24)</td><td>13.54</td><td>15.01</td><td>24.87</td><td>102.23</td><td>15.78</td><td>3.94</td><td>15.78</td></tr><tr><td>M05R1m</td><td>1</td><td>7.27</td><td>7.39</td><td>29.72</td><td>112 (129)</td><td>15.92</td><td>15.66</td><td>40.92</td><td>190.03</td><td>37.33</td><td>9.33</td><td>74.65</td></tr><tr><td>M05R3m</td><td>3</td><td>6.21</td><td>6.80</td><td>56.85</td><td>200 (85)</td><td>10.32</td><td>9.34</td><td>25.91</td><td>222.49</td><td>64.65</td><td>16.16</td><td>64.65</td></tr><tr><td>M05R5m</td><td>5</td><td>6.95</td><td>4.69</td><td>6.59</td><td>9 (146)</td><td>39.36</td><td>49.61</td><td>70.71</td><td>198.79</td><td>64.65</td><td>16.16</td><td>64.65</td></tr><tr><td>M07Sm</td><td>1/4</td><td>25.44</td><td>18.14</td><td>30.35</td><td>-53 (-32)</td><td>11.85</td><td>11.05</td><td>17.72</td><td>62.82</td><td>4.50</td><td>1.13</td><td>4.50</td></tr><tr><td>M07R1m</td><td>1</td><td>16.34</td><td>14.48</td><td>44.72</td><td>111 (120)</td><td>5.28</td><td>5.00</td><td>8.12</td><td>97.82</td><td>8.22</td><td>2.06</td><td>16.45</td></tr><tr><td>M07R3m</td><td>3</td><td>14.74</td><td>11.21</td><td>38.41</td><td>68 (187)</td><td>29.62</td><td>33.21</td><td>67.96</td><td>108.46</td><td>14.24</td><td>3.56</td><td>28.49</td></tr><tr><td>M07R5m</td><td>5</td><td>13.42</td><td>11.55</td><td>14.34</td><td>-2 (223)</td><td>35.85</td><td>42.33</td><td>54.71</td><td>119.11</td><td>18.39</td><td>4.60</td><td>18.39</td></tr><tr><td>M09Sm</td><td>1/2</td><td>53.51</td><td>36.80</td><td>48.98</td><td>-36 (-25)</td><td>1.70</td><td>1.66</td><td>1.68</td><td>35.83</td><td>2.72</td><td>0.68</td><td>2.72</td></tr><tr><td>M09R1m</td><td>1</td><td>38.74</td><td>35.32</td><td>68.55</td><td>102 (108)</td><td>2.32</td><td>2.44</td><td>3.17</td><td>49.50</td><td>3.86</td><td>0.97</td><td>7.72</td></tr><tr><td>M09R3m</td><td>3</td><td>30.61</td><td>32.42</td><td>148.70</td><td>265 (288)</td><td>1.07</td><td>1.07</td><td>1.93</td><td>62.64</td><td>6.67</td><td>1.67</td><td>6.67</td></tr><tr><td>M09R5m</td><td>5</td><td>27.94</td><td>19.74</td><td>56.43</td><td>76 (338)</td><td>20.33</td><td>19.82</td><td>47.02</td><td>68.62</td><td>7.18</td><td>1.80</td><td>7.18</td></tr><tr><td>M11R1m</td><td>1</td><td>130.77</td><td>93.56</td><td>140.61</td><td>-102 (-131)</td><td>12.69</td><td>11.81</td><td>13.06</td><td>16.67</td><td>1.46</td><td>0.37</td><td>2.93</td></tr><tr><td>M11R3m</td><td>3</td><td>90.23</td><td>81.57</td><td>272.53</td><td>278 (291)</td><td>4.49</td><td>4.66</td><td>6.83</td><td>24.17</td><td>2.54</td><td>0.63</td><td>2.54</td></tr><tr><td>M11R5m</td><td>5</td><td>88.50</td><td>61.74</td><td>88.73</td><td>109 (435)</td><td>18.63</td><td>18.48</td><td>32.93</td><td>24.63</td><td>2.62</td><td>0.65</td><td>3.27</td></tr></tbody></table>

Note. Characteristic velocities, differential rotation, magnetic fields, and timescales are listed, using averages over a small interval of 0.01R<sub>*</sub> at the middle of the convective envelopes (unless stated otherwise). The differential rotation is taken between latitude 60° and the equator at the surface of the models (see Section 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 Section 5.4). The rms velocity and magnetic field are $\bar{v} = (\bar{V}_r^2 + \bar{V}_\theta^2 + \bar{V}_\phi^2)$ and $\bar{B} = (\bar{B}_r^2 + \bar{B}_\theta^2 + \bar{B}_\phi^2)$. Here, $\tau_c = D/\bar{v}$ is the overturning convection time, and the dissipation timescales are defined as $\tau_x = D^2/x$ with $x \in [\nu, \kappa, \eta]$, where $D = r_{top} - r_{bcz}$ the thickness of the convective layer that differs for each mass bin.

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, compared 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 times 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 DR as well as 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 Brun et al. (2017) 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 antisolar DR regimes (the transition is nearly at $Ro_f \simeq 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 into a more detailed analysis of our dynamo simulations ensemble in Sections 4–6.

### 3. 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.
<table><caption>**Table 3**
Model Characteristics Nondimensional Numbers</caption><thead><tr><th></th><th>Re</th><th>R<sub>m</sub></th><th>Pe</th><th>P<sub>m</sub></th><th>Ra<br>[10<sup>6</sup>]</th><th>Ra*/Ra<sub>c</sub></th><th>Ta<br>[10<sup>6</sup>]</th><th>Ro<sub>f</sub></th><th>Ro<sub>c</sub></th><th>Ro<sub>s</sub></th><th>Ek<br>[10<sup>-3</sup>]</th><th>Λ<br>[10<sup>-3</sup>]</th></tr></thead><tbody><tr><td>M05Sm</td><td>56.34</td><td>56.34</td><td>14.08</td><td>1</td><td>0.01</td><td>67.64</td><td>0.21</td><td>1.74</td><td>0.38</td><td>3.96</td><td>4.37</td><td>4.57</td></tr><tr><td>M05R1m</td><td>71.71</td><td>143.42</td><td>17.93</td><td>2</td><td>0.89</td><td>16.83</td><td>37.49</td><td>0.33</td><td>0.31</td><td>0.70</td><td>0.33</td><td>1.88</td></tr><tr><td>M05R3m</td><td>106.08</td><td>106.08</td><td>26.52</td><td>1</td><td>7.15</td><td>8.58</td><td>1012.05</td><td>0.15</td><td>0.17</td><td>0.23</td><td>0.06</td><td>0.14</td></tr><tr><td>M05R5m</td><td>106.08</td><td>106.08</td><td>29.68</td><td>1</td><td>7.60</td><td>5.69</td><td>2811.24</td><td>0.07</td><td>0.10</td><td>0.14</td><td>0.04</td><td>4.59</td></tr><tr><td>M07Sm</td><td>26.16</td><td>26.16</td><td>6.54</td><td>1</td><td>0.01</td><td>73.32</td><td>0.05</td><td>1.24</td><td>0.72</td><td>2.98</td><td>9.03</td><td>7.18</td></tr><tr><td>M07R1m</td><td>30.69</td><td>61.38</td><td>7.67</td><td>2</td><td>0.11</td><td>28.75</td><td>1.82</td><td>0.42</td><td>0.50</td><td>0.89</td><td>1.48</td><td>0.39</td></tr><tr><td>M07R3m</td><td>47.94</td><td>95.88</td><td>11.98</td><td>2</td><td>0.64</td><td>14.78</td><td>49.12</td><td>0.16</td><td>0.23</td><td>0.30</td><td>0.29</td><td>8.40</td></tr><tr><td>M07R5m</td><td>56.36</td><td>56.36</td><td>14.09</td><td>1</td><td>1.59</td><td>8.91</td><td>227.44</td><td>0.09</td><td>0.17</td><td>0.18</td><td>0.13</td><td>8.68</td></tr><tr><td>M09Sm</td><td>27.69</td><td>27.69</td><td>6.92</td><td>1</td><td>0.01</td><td>54.09</td><td>0.05</td><td>1.28</td><td>0.74</td><td>3.02</td><td>8.98</td><td>0.12</td></tr><tr><td>M09R1m</td><td>28.48</td><td>56.96</td><td>7.12</td><td>2</td><td>0.04</td><td>24.71</td><td>0.40</td><td>0.68</td><td>0.66</td><td>1.51</td><td>3.16</td><td>0.14</td></tr><tr><td>M09R3m</td><td>38.90</td><td>38.90</td><td>9.73</td><td>1</td><td>0.33</td><td>11.02</td><td>10.79</td><td>0.27</td><td>0.35</td><td>0.50</td><td>0.61</td><td>0.01</td></tr><tr><td>M09R5m</td><td>38.21</td><td>38.21</td><td>9.55</td><td>1</td><td>0.36</td><td>7.00</td><td>34.70</td><td>0.10</td><td>0.20</td><td>0.30</td><td>0.34</td><td>5.33</td></tr><tr><td>M11R1m</td><td>32.06</td><td>64.12</td><td>8.02</td><td>2</td><td>0.01</td><td>47.17</td><td>0.06</td><td>1.38</td><td>0.78</td><td>3.30</td><td>8.33</td><td>11.36</td></tr><tr><td>M11R3m</td><td>38.31</td><td>38.31</td><td>9.58</td><td>1</td><td>0.10</td><td>17.16</td><td>1.56</td><td>0.54</td><td>0.50</td><td>1.10</td><td>1.60</td><td>0.51</td></tr><tr><td>M11R5m</td><td>38.81</td><td>38.81</td><td>9.70</td><td>1</td><td>0.13</td><td>6.87</td><td>4.62</td><td>0.27</td><td>0.34</td><td>0.66</td><td>0.93</td><td>13.06</td></tr></tbody></table>

Note. All quantities listed have been averaged over a small interval of 0.01R$_*$ at the middle of the convective envelopes. Re = $\bar{v}D/\nu$ is the Reynolds number. The Prandtl number Pr = $\nu/\kappa = 1/4$ in all cases. P$_m$ = $\nu/\eta$ is the magnetic Prandtl number, R$_m$ = ReP$_m$ is the magnetic Reynolds number, and Pe = RePr is the Péclet number. Ra = $(-\partial\hat{p}/\partial S)\Delta SgD^3/\hat{p}\nu\kappa$ is the Rayleigh number, and Ra*/Ra$_c$ is the modified Rayleigh number as computed by Takehiro et al. (2020). Ta = $4\Omega_c^2 D^4/\nu^2$ is the Taylor number. We also list three Rossby numbers: the fluid Rossby number Ro$_f = \bar{\omega}/2\Omega_*$, the convective Rossby number Ro$_c = \sqrt{Ra/TaPr}$, and the stellar Rossby number Ro$_s = P_{tot}/\tau_c^{CS}$. The latter is useful for comparison with observationally derived Rossby numbers. For Ro$_s$ we have therefore considered the empirical convective turnover time derived by Cranmer & Saar (2011), which is $\tau_c^{CS} = 314.24 \exp[-\frac{T_{eff}}{1952.5 \text{ K}} - (\frac{T_{eff}}{6250 \text{ K}})^{18}] + 0.002$ days. We note that it correlates well with our fluid Rossby number and find Ro$_s \approx 2.26Ro_f$. The Ekman number is defined as Ek = $\nu/(\Omega_c D^2)$ and the Elsässer number as $\Lambda = \bar{B}^2/8\pi\hat{p}D\bar{v}\Omega_*$.

<figure><img src="image_64.png" /><figcaption>**Figure 3.** Temporal evolution in case M09R3m of the 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 yr period. Case M09R3m is indeed one of our cyclic cases (see also Figure 4).</figcaption></figure>

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 ($Ro_f = 0.27$) and therefore lies in the middle of our sample of models. The temporal evolution of kinetic and magnetic energies is shown in Figure 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, and fluctuating) oscillate in phase in this model. The mean toroidal kinetic energy also presents oscillation of the same amplitude, albeit anticorrelated 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 Figure 3) and points toward a similar dynamo mechanism involving a strong feedback of the magnetic field on the DR within the convective envelope. We will perform a detailed analysis of this mechanism in Section 5. Here, we first illustrate the dynamics of the dynamo achieved in model M09R3m in Figure 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 the 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_\phi$ 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 entrophy $|\nabla \times \mathbf{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 of up to $1.5 \times 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 100 G and exhibits a complex topology, mixing dipolar and quadrupolar
<figure><img src="image_65.png" /><figcaption>**Figure 4.** Temporal evolution of a magnetic cycle for case M09R3m, taken at three different instances during one reversal. The upper row shows a 3D potential extrapolation of the modeled magnetic field, with blue lines denoting field lines oriented outward and black lines oriented inward. Behind a semitransparent representation, the radial velocity is 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 entrophy at the same depth.</figcaption></figure>
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) show little imprint on the convective flows themselves and mainly act on the mean flows (see Section 4.1). The specific entropy fluctuations have two distinctive features. First, 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 DR (Brun & Toomre 2002) and can be generally related to the pressure field required to drive the observed meridional flow. Second, patterns in $S'$ are imprinted by the nonaxisymmetric 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 (Section 4), magnetic properties (Section 5), and energetic balances (Section 6) achieved in the 15 models. The reader mostly interested in the astrophysical consequences of our study may consider going directly to Section 7 for a summary.

## 4. Large-scale Flows in the Models

In this section, we analyze the DR 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 DR profiles of the hydrodynamical and MHD models. We confirm our preliminary results (Varela et al. 2016) and those 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 DR 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.

### 4.1. Differential Rotation Profiles as a Function of Rossby Number

We analyze the DR 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 DR averaged over 10 overturning convective times, defined as $\tau_c = \int_{r_{bcz}}^{r_{top}} dr / \tilde{v}_r$ (see Table 2).

We observe that for the simulations M05Sm, M07Sm, M09Sm, and M11R1m, there is an antisolar DR, with the poles rotating faster than the equator, like their hydrodynamical counterparts (see Figure 6 and also Brun et al. 2017). The cases rotating at an intermediate rotation rate show a solar-like DR. Finally, the cases rotating the fastest (R5 series) show almost no DR (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 a magnetic field. This effect is observed to be stronger as the rotation rate increases. Conversely, the effect of the 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 DR. As one may expect, the radial gradient of the DR near the tachocline is generally weaker in all MHD simulations compared with hydrodynamic progenitors. This points to a magnetic feedback of the dynamo field on the DR itself, a feedback that is observed to strengthen as the rotation rate increases.

We have calculated the surface latitudinal DR $\Delta\Omega$ for each model, defined as the difference between the equator and 60° latitude. A positive value thus denotes a solar-like DR and a negative value an antisolar DR. We report these values for the magnetic cases as well as the hydrodynamic progenitors in Table 2 (fourth column).

The DR of our sample spans a range between -102 and +278 nHz, with some fast-rotating models presenting an extremely weak DR, like M05R5m with $\Delta\Omega = 9$ nHz. We find that the absolute DR generally weakens in MHD models compared to their hydrodynamic progenitors, as expected from the radial profiles shown in Figure 6. This is particularly striking for fast rotators such as M09R5m, which goes from 338 nHz in hydro to 76 nHz in MHD.

We investigate in Figure 7 the DR trends with respect to the rotation rate (left panel) and rotation period (right panel). The DR of the hydrodynamic progenitors and of the MHD cases are respectively shown in small semitransparent and large opaque symbols. The shape of the symbol labels the rotation of the model, and the color the mass of the modeled star, as indicated in the legend. In the right panel, we compare the model DR to the DR 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 DRs agrees well with the observed values. In addition, the DR range in our sample increases as the rotation period decreases, like what is observed in the Kepler satellite sample. Several of our models nevertheless lie outside the observed values: the three antisolar DRs on the right (triangles) and two of our fast-rotating models. Several reasons can explain this discrepancy. Slowly rotating stars could produce very few starspots or even no starspots at all (see, for instance, van Saders et al. 2019), making their DR 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 DRs due to magnetic feedback, which are outside the detection limits of Kepler (dotted black lines).

The left panel of Figure 7 shows the DR trend with the rotation rate. Using only the hydrodynamic progenitors, we
<figure><img src="image_66.png" /><figcaption>**Figure 5.** Temporal and longitudinal averages of the angular velocity profiles over 10 convective overturning times (10 $\tau_c$) in our suite of models. Prograde flows are in reddish tones and retrograde ones in bluish tones. In each panel, the dashed semicircle represents the base of the convective envelope and the dashed horizontal line the equator.</figcaption></figure>

previously showed that the DR scales as $\Omega^{0.66}$ (Brun et al. 2017). Blindly trying to fit such a power law to the MHD sample, we find that the exponent reduces to $\Omega^{0.46}$. This weaker dependency is expected due to the magnetic feedback on the DR 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 left panel of Figure 7, it clearly appears that a power-law fit is a poor representation of the DR in our sample. Rather, we see that $\Delta\Omega$ increases with $\Omega$ for slow rotators, while it dramatically drops for fast rotators due to the magnetic feedback. Following Saar (2011), we recast in Figure 8 the DR trend in terms of relative DR $\Delta\Omega/\Omega$ with respect to the fluid Rossby number $R_o$ (for the different definitions of Rossby number used in this work, see the caption of Table 3 or the Appendix A of Brun et al. 2017). We

find a trend that is very similar to the observational trend reported by Saar (2011) (shown by the dashed line in Figure 8): $\Delta\Omega/\Omega$ is roughly constant for inverse Rossby numbers lower than a certain threshold (here $R_o^{-1} \lesssim 5$), and it drops for fast rotators as $\Delta\Omega/\Omega \propto R_o^p$. Saar (2011) proposed that $p=2$, but here our sample agrees with a somewhat large range $p \in [2, 6]$. Additional models with even higher turbulence level are required to confirm the exact amplitude of the drop in DR contrast found in the fast-rotating cases. Finally, our sample also shows some hint of an increase of $\Delta\Omega/\Omega$ 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 antisolar rotation states.

In Brun et al. (2017), we have proposed that the DR could follow two power laws with respect to the Rossby number and the stellar mass. Here, we find that the DR is weakened at high Rossby number, and therefore we do not recover a simple
<figure><img src="image_67.png" /><figcaption>Figure 6. Radial cuts of the temporal- and longitudinal-averaged angular velocity from the equator to 75° latitude every 15° (black lines are for the hydrodynamical cases and the blue lines the MHD/dynamo cases). The solid vertical line in each panel shows the bottom of the convective layer.</figcaption></figure>

power-law trend, as we saw in Figure 7. We can nevertheless attempt to fit such a combined power law on a subsample of our models, excluding the fast-rotating case but retaining the slow rotators. We obtain in this way

$$ \Delta\Omega = 107Ro_{f}^{-0.73\pm0.13} \frac{M_{*}^{1.93\pm0.42}}{M_{\odot}} \text{nHz (HD)}, \quad (7) $$

$$ \Delta\Omega = 84Ro_{f}^{-0.40\pm0.20} \frac{M_{*}^{0.78\pm0.62}}{M_{\odot}} \text{nHz (MHD)}. \quad (8) $$

In the MHD case, we find again that the DR 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
<figure><img src="image_68.png" /><figcaption>**Figure 7.** Absolute value of the differential rotation between the equator and 60° latitude vs. rotation rate (left panel) and rotation period (right panel). The symbols denote the rotation rate of the model and the color the mass of the modeled star, as shown in the legend. MHD models are shown by the large plain symbols, and the hydrodynamic progenitors by the smaller open ones. On the right 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.</figcaption></figure>

nevertheless included the results of the fit here to compare with the purely hydrodynamic case (Brun et al. 2017). We can conclude here that the clear trend in stellar mass and effective temperature found in the hydrodynamic study (Brun et al. 2017) 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.

The MHD simulations therefore show that the magnetic field changes the angular momentum redistribution, especially for fast-rotating stars. In the next section, we perform a detailed analysis of this balance for four representative models.

## 4.2. Angular Momentum Transfer

We can better understand how the DR 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 BCs at the top and stress-free and perfect conductor BCs 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 ($\mathcal{F}_r$) and latitudinal ($\mathcal{F}_\theta$) angular momentum fluxes, applying the procedure used in Brun et al. (2004). Starting from the $\phi$ component of the momentum equation expressed in conservative form and averaged in time and longitude:

$$ \frac{1}{r^2} \frac{\partial(r^2 \mathcal{F}_r)}{\partial r} + \frac{1}{r \sin \theta} \frac{\partial(\sin \theta \mathcal{F}_\theta)}{\partial \theta} = 0, \quad (9) $$

involving the mean radial angular momentum flux

$$ \mathcal{F}_r = \hat{\rho} r \sin \theta \left[ -\nu r \frac{\partial}{\partial r} \left( \frac{\hat{v}}{r} \right) + \widehat{v'_r v'_\phi} + \hat{v}_r (\hat{v}_\phi + \Omega r \sin \theta) - \frac{1}{4\pi\hat{\rho}} \widehat{B'_r B'_\phi} - \frac{1}{4\pi\hat{\rho}} \hat{B}_r \hat{B}_\phi \right], \quad (10) $$

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /><figcaption>**Figure 8.** Relative differential rotation between the equator and 60° latitude as a function of the fluid Rossby number. The symbols shape and color are the same as in Figure 7. The trend found in the MHD sample is highlighted by the gray area, and the observational trend reported by Saar (2011) is shown by the dashed black line.</figcaption></figure>

and the mean latitudinal angular momentum flux

$$ \mathcal{F}_\theta = \hat{\rho} r \sin \theta \left[ -\nu \frac{\sin \theta}{r} \frac{\partial}{\partial \theta} \left( \frac{\hat{v}_\phi}{\sin \theta} \right) + \widehat{v'_\theta v'_\phi} + \hat{v}_\theta (\hat{v}_\phi + \Omega r \sin \theta) - \frac{1}{4\pi\hat{\rho}} \widehat{B'_\theta B'_\phi} - \frac{1}{4\pi\hat{\rho}} \hat{B}_\theta \hat{B}_\phi \right]. \quad (11) $$

In these equations, the terms on the right-hand side represent, for both fluxes, contributions respectively from viscous diffusion (which we denote as $\mathcal{F}_r^{\text{VD}}$ and $\mathcal{F}_\theta^{\text{VD}}$), Reynolds stresses ($\mathcal{F}_r^{\text{RS}}$ and $\mathcal{F}_\theta^{\text{RS}}$), meridional circulation ($\mathcal{F}_r^{\text{MC}}$ and $\mathcal{F}_\theta^{\text{MC}}$), Maxwell stresses ($\mathcal{F}_r^{\text{MS}}$ and $\mathcal{F}_\theta^{\text{MS}}$), and large-scale magnetic torques ($\mathcal{F}_r^{\text{MT}}$ and $\mathcal{F}_\theta^{\text{MT}}$). The Reynolds stresses are
<figure><img src="image_71.png" /><figcaption>Figure 9. Angular momentum transport in the 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, the radial (left) and latitudinal (right) angular momentum balance are shown, with Reynolds stress contribution shown as green dashed-dotted lines, viscous stresses as dashed-triple-dotted yellow lines, meridional circulation as dashed cyan lines, Maxwell stresses as long dashed magenta lines, the large-scale magnetic torques as a black dotted line and the sum of all contributions as a 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.</figcaption></figure>

linked to correlations of the fluctuating velocity components coming from organized tilts within the convective structures, especially in the downflow plumes. Likewise, the Maxwell stresses are associated with correlations of the fluctuating magnetic field components due to the twist and tilt of the dynamo-generated magnetic structures.

In Figure 9 we show the components of $\mathcal{F}_r$ and $\mathcal{F}_\theta$ for the M07 case series, having integrated over colatitude and radius as follows:

$$ 
\begin{aligned}
\Phi_r(r) &= \int_0^\pi \mathcal{F}_r(r, \theta) r^2 \sin \theta \, d\theta, \\
\Phi_\theta(\theta) &= \int_{r_{\text{bot}}}^{r_{\text{top}}} \mathcal{F}_\theta(r, \theta) r \sin \theta \, dr.
\end{aligned} 
\qquad (12) 
$$

Thus, $\Phi_r$ represents the net angular momentum flux through horizontal shells at different radii and $\Phi_\theta$ 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 $\mathcal{F}_r$ and $\mathcal{F}_\theta$.

For each of the four cases, we display $\Phi_r$ on the left panel and $\Phi_\theta$ 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 dashed-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 Section 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 dotted line) have little influence on the overall radial angular momentum balance.

Considering now $\Phi_\theta$, we can assess the balance of latitudinal angular momentum transport. We first notice that the Reynolds stresses (green curves) are systematically equatorward in both
<figure><img src="image_72.png" /><figcaption>**Figure 10.** Meridional flows in M07m case series. The contours are normalized and denote clockwise (red) and anticlockwise (blue) circulations.</figcaption></figure>

hemisphere (positive in the northern hemisphere and negative in the southern one). Because 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 the 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 DR 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.

### 4.3. 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 Section 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 DR (see, e.g., Passos & Charbonneau 2014).

We illustrate the meridional flow pattern achieved in the M07m set of simulations in Figure 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 counterrotating 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 DR is strongly quenched by magnetic feedback as seen in the previous section. In particular, the radial shear of DR vanishes at the equator as seen in Figure 6. As a result, gyroscoping pumping (Miesch et al. 2006; 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 Figure 10. Having presented the large-scale flows achieved in the simulations, we now turn to discussing their magnetic properties.

## 5. Magnetic Properties

In this section, we discuss in more detail various aspects of our dynamo simulations, such as their type, their temporal
<figure><img src="image_73.png" /><figcaption>**Figure 11.** Various dynamo states achieved in our sample by our models, as illustrated by M07R5m (top panels), M09R3m (middle panels), and M09R1m (lower panels). In the top panel, we show the time–latitude, time–radius, and instantaneous meridional plane of $B_r$ (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 the latitude–time (at the base of the convection zone) and radius–time (at midlatitude) 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.</figcaption></figure>

variability, the amount of magnetic flux they generate, and the distribution in space and size of their magnetic fields.

## 5.1. Properties of Dynamo Solutions: Long Cycles, Short Cycles, and Steady Dynamos

We find three dynamo states in our sample of 15 MHD models: long (decadal) magnetic cycles, short (yearly) magnetic cycles, and stable magnetic wreaths concentrated close to the bottom of the convection zones. These three states are illustrated in Figure 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 Figure 11). In this model,

we find that the global magnetic field of the star reverses with a period of 10 yr (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 polar and equatorial branches. The magnetic field is generally consistent with dipolar symmetry, with the azimuthal field of inverse polarities in each hemisphere. We also see some departures from hemispheric symmetry (for instance around $t=42$ yr). The azimuthal field is found to be concentrated at the base of the convective envelope and in the tachocline, where the radial shear of $\Omega$ 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