Nature of Intense Magnetism and Differential Rotation in Convective Dynamos of M-dwarf Stars with Tachoclines

Many of the M-dwarf stars, though they are tiny and dim, are observed to possess strong surface magnetic fields and exhibit remarkably intense flaring. Such magnetism may severely impact habitability on the exoplanets now discovered nearby. The origin of the magnetism must rest with dynamo action achieved by turbulent convection coupled to rotation within the M-dwarfs. To further explore the nature and diversity of the magnetism that can result, we turn here to an extensive set of 45 global MHD simulations to explore dynamos operating within deep convective envelopes of rapidly rotating M2 (0.4 M ⊙) stars. We observe a wide range of cycle periods present in the convection zones, whose durations we find to scale with the Rossby number as Ro−1.66±0.07 in concurrence with scalings identified in simulations of more massive stars. We find a unifying relationship between the ratio of magnetic to convective kinetic energy (ME/CKE) and the degree to which the differential rotation is quenched by magnetic fields. We show that the presence of a tachocline in these model stars enhances their axisymmetric magnetic field components, leading to a surface dipole fraction on average 78% greater than an equivalent star with only a CZ, potentially shedding light on the nature of the tachocline divide through resultant effects on the spin-down rate.


Introduction
M-dwarfs, the smallest, coolest, and dimmest stars on the main sequence, have come to be known for their ubiquitous and striking magnetic activity. While less than 20% of more massive stars demonstrate chromospheric markers for magnetic activity, the active fraction of fully convective (FC) late-type M-dwarfs appears to be on the order of 90% (e.g., West et al. 2008West et al. , 2015. This behavioral transition occurs over a narrow range of stellar masses centered on 0.35 M e . Below this mass, main-sequence stars lose their radiative zones (RZs) and become FC; thus, the transition has come to be called the tachocline divide. The relatively laminar shearing flows of the solar tachocline are often presented as crucial to the dynamo action of Sun-like stars (e.g., Gilman & Fox 1997;Tobias 2005). Therefore the preclusion of a tachocline in FC stars is speculated to demand that their distinctly vigorous magnetism is a consequence of fundamentally different dynamo processes. However, this line of thought is in tension with recent measurements of activity on slowly rotating FC stars (Wright et al. 2018), which have found no significant deviations in their rotation-activity relation from that of more massive M-stars. Furthermore, recent work by Reiners et al. (2022) recovered this continuity across the tachocline divide through direct measurements of surface magnetic field strengths. These findings would suggest that the underlying dynamo mechanism may in fact be shared, and thus that the influence of the tachocline on surface magnetism may be marginal.

Solar and Stellar Differential Rotation
The advent of helioseismology allowed for the first maps of the differential rotation occurring in the interior of the Sun (e.g., Dziembowski et al. 1989;Schou et al. 1998). The discovery that the differential rotation established in the solar convection zone (CZ) did not imprint into the deep interior, but instead shifted to near solid-body rotation over a thin tachocline came as a real surprise. It remains one of the great open questions in solar physics how the shear of the tachocline is maintained, given that diffusive effects should have imprinted it deeply into the RZ over over the lifetime of the Sun. One prominent scenario by which the tachocline may be maintained hydrodynamically was proposed by Spiegel & Zahn (1992), and invokes rapid horizontal redistribution of angular momentum through anisotropic turbulence. Gough & McIntyre (1998) claimed instead that the tachocline is prevented from spreading by a primordial magnetic field trapped in the RZ. This theory has been tested in multiple computational experiments (e.g., Strugarek et al. 2011), though reconnections with the CZ magnetic fields inevitably lead angular momentum to be pumped into the RZ in accordance with Ferraro's law of isorotation. It has also been proposed that confinement of the tachocline is achieved via downward imprinting of magnetic fields from a cycling CZ dynamo, which is restricted to a plasma skin depth (Forgács-Dajka & Petrovay 2001). Recent simulations by Matilsky et al. (2022) suggest that the thinness of the tachocline may be a consequence of nonaxisymmetric dynamo action in the supposedly quiescent RZ by Rossby waves and shear instabilities.
While helioseismology has allowed us to become intimately familiar with the differential rotation of the Sun, asteroseismic inversions must be undertaken with far less granularity in their data (Gizon & Solanki 2004). It has been shown to perform well for red giant stars (Chaplin & Miglio 2013), and can be managed for Sun-like stars under suitable circumstances (e.g., Nielsen et al. 2015;Benomar et al. 2018). So far, however, asteroseismology has not been able to provide well-constrained estimates of the differential rotation in low-mass M-dwarfs. Instead, photometric estimates for the latitudinal differential rotation on a small number of M-dwarfs using inversions based on the relative shear between two starspots (e.g., Vida et al. 2014;Davenport et al. 2015;Zaleski et al. 2020) have repeatedly shown that rapidly rotating M-dwarfs have very little differential rotation. The requirement of multiple persistent starspots at well-constrained latitudes renders this methodology more opportunistic in nature, and it cannot be employed in all circumstances, particularly for slow rotators that may demonstrate less activity. There remain serious gaps and uncertainties in our ability to measure or predict the differential rotation of these stars. As we have limited capacity to observe the flows and magnetic fields beneath stellar photospheres, the study of their differential rotation and the associated dynamo action is perhaps best approached through theory and computation.

Convective Dynamo Simulations
Solar and stellar convection and dynamo theory has seen substantial development through the use of 3D global fluid simulations of the stellar interior carried out on massively parallel machines (e.g., Brun & Browning 2017;Charbonneau 2020). Early work on the solar dynamo problem found that rotationally influenced turbulent convection could simultaneously establish a solar-like differential rotation with fast equators and slow poles, and build strong magnetic fields within the CZ (e.g., Brun et al. 2004). With the discovery through helioseismology of a tachocline at the base of the solar CZ, suggestions arose that an interface dynamo operating off the shear there may be better able to explain the emergence of sunspots and the cyclic reversals of their polarities. Subsequent solar-like dynamo models that explicitly included the tachocline (e.g., Browning et al. 2006;Ghizaru et al. 2010;Augustson et al. 2013Augustson et al. , 2015Passos & Charbonneau 2014;Guerrero et al. 2016;Strugarek et al. 2018;Matilsky & Toomre 2019) found support for the idea, building strong, wellorganized fields that tended toward longer cycle periods. Separately, it was discovered that with increased turbulence and rotational constraint, the CZ itself could build and sustain strong toroidal wreaths of magnetic field with orderly cycles (e.g., Glatzmaier 1985;Brown et al. 2010Brown et al. , 2011Matilsky & Toomre 2020).
Relative to the solar context, few global convective dynamo studies have been conducted for M-dwarf stars. Browning (2008) considered lower-mass FC M-dwarfs, finding that very strong nonaxisymmetric fields could be built in the deep CZ, which suppressed the star's differential rotation. Subsequent simulations with great levels of turbulence yielded a variety of interesting results.  found strong, axisymmetric fields that did not cycle but reproduced many characteristics of observed M-dwarf surface fields. With a slower rotating model , it was found that the concentration of magnetic flux by merging downflow lanes could cause the formation of large, persistent starspots at high latitudes in these stars. A still slower rotator (Yadav et al. 2016) built weaker but largely axisymmetric magnetic fields, which cycled regularly and did not quench the star's differential rotation. Brown et al. (2020) explored the first simulation of a stratified, rotating FC star whose computational domain reached r = 0, finding an abundance of single-hemisphere dynamo configurations. A set of similar models spanning a range of rotation rates were subsequently presented in Käpylä (2021). They found that at their slowest rotation rates, differential rotation was antisolar with fast poles and a slow equator, and produced dipole-dominant dynamo action. Models with intermediate rotation rates achieved solar-like differential rotation, while their fastest rotators produced strong nonaxisymmetric magnetic fields, which quenched the differential rotation.
In Bice & Toomre (2020) (BT20), we explored the influence exerted by a tachocline in shell-convecting M2 (0.4 M e ) stars as a potential contributing factor to the differences in observed activity across the tachocline divide. We found a variety of possible field configurations, identifying a preference for singlehemisphere states. Nearly all models presented in BT20 had fields that strongly quenched their differential rotation, but surprisingly still exhibited substantial induction through the mean field Ω-effect. We found that the inclusion of a tachocline of shear in our models tended to increase their cycle periods, and induced changes in their surface fields, without significantly altering the overall character of the dynamo operating in the CZ. In BT20, models with tachoclines produced surface fields that were both stronger, and organized onto larger spatial scales, conditions which favor more rapid spin-down via the stars' magnetized winds (e.g., Parker 1965;Matt et al. 2012). Three of the models that appear in our current work were previously reported on in Bice & Toomre (2022, hereafter BT22). They exhibited a previously unexplored dynamo mechanism associated with a persistent traveling nest of enhanced convection, which resulted in locally accelerated flux emergence and the formation of surface structures akin to active magnetic longitudes.
Our work here is motivated by two main lines of inquiry. How does the presence of a tachocline in these stars impact the overall amplitudes of magnetic fields at both large and small scales? How are the characteristics and durations of magnetic cycles affected? Does the presence of a tachocline meaningfully alter emergent poloidal fields in a way that could contribute to the tachocline divide through enhanced spindown? We have previously approached these questions in BT20 with a modest number of simulations. We revisit them here to extend those analyses to more challenging and nonlinear regions of parameter space, particularly with respect to the Rayleigh and magnetic Prandtl numbers, also enabling more robust statistical analysis. Second, we seek to understand the balance between magnetism and differential rotation in these stars. What determines the extent to which magnetic fields eliminate the differential rotation? By what mechanism is their dynamo action maintained in the near-total absence of large-scale shear? By focusing on the magnetic Prandtl number, we complement existing work on the subject that predominantly considers only the influence of rotation rate, as encoded by the Rossby number.

Formulating the Problem
To investigate the dynamo action and consequences thereof in shell-convecting M-dwarfs, we consider 3D global MHD simulations of the interiors of M2 (0.4 M e ) stars run using the open-source Rayleigh code (Featherstone & Hindman 2016). Rayleigh evolves the anelastic compressible MHD equations in rotating spherical shells using a pseudospectral numerical scheme. The anelastic equations are fully nonlinear while filtering out compressive sound waves, the very short dynamical timescales of which can dramatically inflate computational expense. This provides an appropriate framework for the exploration of subsonic convection in the interior of the star, but requires that we truncate our computational domain below the photosphere, where flow speeds may become supersonic, and nonlocal radiative transfer becomes an important effect. In Rayleigh, the thermodynamic variables are linearized against a 1D, time-steady background state defined by the density r , pressure P, temperature T , and entropy S¯. Fully dimensional deviations from this background are written without overbars. As with all large eddy simulations, the viscosity ν, conductivity κ, and resistivity η are inflated by many orders of magnitude as a parameterization of mixing by subgrid scale turbulence. With the velocity vector as v and the magnetic field vector B, the form of the anelastic equations solved by Rayleigh is as follows: .
Here, Q is the volumetric heating function,  is the viscous stress tensor, and Φ represents the viscous heating, which are defined as with e ij as the strain-rate tensor. Additionally, both the mass flux and magnetic field are divergenceless: v B 0. r  = = · (¯) · Closure is achieved with a linearized equation of state: The thermodynamic background states against which our simulations were calculated were derived using the stellar evolution code MESA (Paxton et al. 2010(Paxton et al. , 2013(Paxton et al. , 2015. In particular, we consider ZAMS 0.4 M e stars with solar metallicity, a luminosity of 9.478 × 10 31 erg s −1 (0.025L e ), and a range of rotation rates. The masses of these stars place them just above the transition to full-sphere convection; thus, they contain deep RZs and potentially tachoclines. In order to compare the effects induced by the inclusion of a tachocline, we consider both models with lower boundaries at the base of the CZ with R i = 0.42R * = 1.087 × 10 10 cm, and those that extend deep into the RZ with R i = 0.1R * = 2.588 × 10 9 cm. In both cases, the outer boundary is located at R o = 0.92R * = 2.382 × 10 10 cm, giving a density stratification in the CZ of roughly N ρ = 3.5 scale heights.
The parameters of the 45 simulations that make up the data set considered in this work are summarized in Table 1. The naming convention we employ to distinguish the models hinges on four dimensions that vary across the data set. First, the effective rotation rate enters as fast (F) or slow (S). Next, the three explored levels of convective supercriticality enter as relatively low (a), medium (b), or high (c). Third, we give the magnetic Prandtl number Pm = ν/η or an "H" in the case of hydrodynamical models. An optional final character of "t" or "c" distinguishes models with tachoclines and CZ-only models with "perfect conductor" lower magnetic boundary conditions, described later in this section. For example, the fast-rotating, low turbulence, Pm = 2 model with a tachocline, which we have previously examined in BT22, becomes case "Fa2t" in this data set.
We employ a diffusive profile in the CZ that has seen frequent use by similar studies (e.g., Brun et al. 2004;Brown et al. 2011;Augustson et al. 2013;Fan & Fang 2014), proportional to 1 2 r -, with the top values for each diffusivity varying across simulations. In models containing the RZ, an additional radius-dependent factor dramatically reduces the diffusivities in the radiative interior in accordance with our expectation of limited turbulent mixing there. Additionally, the mean (l = 0) entropy field sees a separate conductivity that drops to a floor value a short distance below the outer boundary. This serves to discourage thermal conduction as a means of energy transport in the bulk of the CZ, and consequently forces the convective motions to carry nearly the full luminosity of the star. The exact functional forms of the diffusion profiles employed here are presented in BT20. The thermal conductivity and electrical resistivity follow the same profiles as the viscosity, yielding fixed values of the Prandtl number Pr = ν/κ = 1/4 and the magnetic Prandtl number Pm, which varies between simulations.
The upper and lower boundary conditions are impenetrable and stress free, The lower boundary is thermally insulating, and the top boundary extracts the star's luminosity through a fixed conductive gradient, with dS dr With no conductive input, energy balance is instead maintained through the volumetric heating function Q, which is adapted from the radiative flux reported by MESA. In most of our models, the magnetic field matches onto an external potential field at both boundaries, as CZ-only models of Sun-like stars that employ a "perfect conductor" lower magnetic boundary (e.g., Brown et al. 2011) tend to produce stronger, more axisymmetric toroidal magnetic fields than those with potential field lower boundaries. It has been suggested that such a boundary condition may be a more appropriate analog for the tachocline lying just outside the computational domain. To test this claim, we computed a small number of supplemental CZ-only models with the perfect conductor lower magnetic boundary condition: Hydrodynamical models were evolved first in the absence of magnetic fields. They were allowed to reach a state that is in equilibrium, except for the tachocline, which erodes viscously on a timescale of dozens to hundreds of years, depending on the diffusivity of the model. Once the described state was achieved, the solutions were forked into a collection of dynamo models distinguished by their value of Pm. Magnetic fields were introduced and allowed to self-consistently reshape the flows as they grew to their mature amplitudes. In general, the saturated magnetic models were allowed to run for several magnetic diffusion times to ensure that the captured dynamics are robust. However, very short near-surface Alfvén timescales  ¢. The Rossby number Ro = ω rms /2Ω * is calculated in terms of the vorticity ω = ∇ × v and is also averaged over the full CZ. The column LMBC denotes the lower magnetic boundary condition for each model, with RZ, PF, and PC denoting a model containing a radiative zone, a CZ-only model with a potential field lower boundary, and a CZ-only model with a perfect conductor lower boundary, respectively.
in some of the particularly turbulent fast rotators made this standard numerically impractical to achieve in all cases.

Diversity of Dynamo Action
To characterize the variety of magnetic responses realized in our survey of dynamo action in these stars, we will first examine in detail four models that together present a cross section of the behaviors found. The first model we selected is Sa2t, which rotates slowly, has relatively low turbulence, Pm of 2, and contains a tachocline. The next model we consider is Sc05, which again rotates slowly, but is highly turbulent, and its Pm of 0.5 gives it a very similar magnetic Reynolds number Rm = v rms L/η to Sa2t. We also examine Fb025, which rotates rapidly, has an intermediate level of turbulence, and a Pm of 0.25. The final model we present here is Fa2t, which appeared also as case C in BT22, and shares nearly all parameters with Sa2t, but rotates four times as rapidly, and too, has a tachocline.

Magnetic Structures
In Figure 1, we give instantaneous snapshots of the magnetic structures realized in each of our selected models. In model Sa2t, there is a clear signature of strong, axisymmetric dynamo action in the tachocline, which reaches mean toroidal field amplitudes 〈B f 〉 in excess of 12 kG. At the same time, the CZ fields strongly favor the southern hemisphere, where they form only a single choppy toroidal wreath of slightly lower amplitude. Away from this wreath, toroidal fields in the CZ of Sa2t are diffuse and incoherent. The poloidal fields of model Sa2t are dominantly octupolar at large scale, with two bands in each hemisphere imprinting outward from the tachocline. They are disrupted in the vicinity of the CZ wreath, but remain evident all the way to the surface. Near the surface, the radial magnetic fields B r become concentrated within downflow lanes by converging flows, allowing them to reach local peaks on the order of 1.5 kG.
Relative to model Sa2t, Sc05 has a nearly identical value of Rm due to a combination of greater fluid turbulence and lower Pm, and has a slightly higher Rossby number Ro. In Sc05, we observe toroidal fields that are prominent in both hemispheres, exhibiting clearly antisymmetric parity across the equator at mid-depths. Near the bottom boundary, one 〈B f 〉 structure of opposite polarity to the CZ is evident in the southern hemisphere, while two structures are present in the north at the time shown. A broad cap of radial field covers much of the surface in the southern hemisphere, whereas deeper in the CZ, poloidal fields are largely antisymmetric. Overall, the coherent and yet finely structured fields of model Sc05 appear to be weaker on average than those of Sa2t. While the two models share the same value of Rm, there are a number of other factors that could be contributing to this difference. The Ro is marginally higher in model Sc05, which is likely to lead to diminished field amplitudes, but not at the level of a factor of 2 or more. There may be additional scalings with Ra or Pm beyond their contribution to Rm. Additionally, the discrepancy could be explained by the presence of a tachocline in Sa2t, or by its hemispheric configuration. We will revisit this question of field scaling in Section 5.1.
Model Fb025 rotates more rapidly than Sa2t or Sc05, with a mid-CZ fluid Rossby number Ro = 0.163, which would put it within or near the saturation regime of the rotation-activity relation. As a result, it shows stronger magnetic fields than were realized in the CZs of either of the previous two, despite its Rm being lower by a factor of nearly 6. In model Fb025, we found highly axisymmetric wreaths of toroidal fields with amplitudes on the order of 15 kG. These wreaths appear isolated to the southern hemisphere, with almost no significant toroidal fields existing outside of them. At the center of panels (i) and (k) of Figure 1, a disruption of the toroidal fields and enhancement of the emergent poloidal fields is apparent. This is a consequence of accelerated turbulent induction within a traveling nest of magnified convection, as we explored in detail in BT22. This phenomenon occurs to some extent in all of the fast-rotating models presented in this work, with models containing stronger nonaxisymmetric fields tending to flatten the nest and reduce the size of the effect.
Relative to model Fb025, Fa2t has a very similar Rossby number, but a value of Rm that is approximately four times greater. The fields of Fa2t are highly nonaxisymmetric, with mid-CZ peak B f amplitudes exceeding 20 kG, and surface B r in places greater than 5 kG. The dynamo action in the CZ of model Fa2t occurs primarily in the southern hemisphere, which at the time shown is dominated by a single wreath of toroidal field. Unlike in Sa2t, the CZ dynamo produces fields of comparable amplitude to those generated in the tachocline. In the southern hemisphere, the fields in the tachocline of model Fa2t are overridden and diminished by the CZ fields. In the northern hemisphere, however, the fields generated in the tachocline are able to spread upward uninhibited, resulting in a relatively weak but almost completely axisymmetric cap of field.
The magnetic fields as they have been presented here are only momentary glimpses at the richly time-dependent behavior realized in each of these four models, as well as in the data set as a whole. We next explore how the mean magnetic fields of the selected models evolve in time.

Evolving Magnetic Fields
In Figure 2, we present butterfly-like time-latitude diagrams of 〈B f 〉 for a portion of the evolutionary history of each of our four example models. Figure 2(a) shows the mid-CZ evolution of model Sa2t, which cycles erratically with an average period computed via Lomb-Scargle periodogram to be 3.02 yr. Over the course of a cycle, magnetic fields emerge at the equator before migrating poleward, breaking up beyond a boundary latitude of ±50°. The hemispheric mode of Figure 1 is a recurring feature, which appears to randomly flip between hemispheres between cycles, on average persisting in one hemisphere for 2.5 cycles before flipping. In the relatively inactive hemisphere, cycles are not well defined; however, magnetic fields may bleed across the equator and loosely mirror the cycle on the other side. There are brief intervals during which the mean fields of Sa2t are not hemispheric, for instance from 164-168 yr in the panel shown. During these intervals, the mean fields exist in roughly equal measure in each hemisphere, though there is no sign of a coherent cycle during those times. Figure 2(b) shows 〈B f 〉 in the tachocline of Sa2t during the same period of its history. There, the fields maintain their antisymmetry through multiple cycles, which take 34.2 yr on average to play out. Rather than poleward migration, the fields in the tachocline of Sa2t cycle in place, fully dying out before being replaced with new fields of the opposite sense.
Figures 2(c) and (d) give the evolution of 〈B f 〉 near mid-CZ and at the bottom boundary of model Sc05, respectively. Three distinct cycles with different morphologies and timings are evident. At low latitudes near mid-CZ, the fields undergo a slow cycle with an average period of 12.5 yr, and are antisymmetric across the equator. Reversals are initiated by new fields forming at the equator and pushing the existing mean fields poleward; however, most of the duration of each cycle is characterized by no field propagation, and the fields break up at latitudes beyond ±40°. Near the bottom boundary, a separate, faster cycle exists at mid-latitudes. There, mean fields form at latitudes of ±50°and migrate equatorward with a mean cycle period of 2.73 yr. Although this cycle occurs in both hemispheres, their reversals are not synchronized, allowing the parity to drift between symmetric and antisymmetric states. The signature of this cycle imprints upward and is visible at mid-latitudes near mid-CZ, where it appears weakened and reversed in sign relative to the deeper fields. Finally, a still faster cycle is present in the polar regions. There, fields cycle in place with an average period of 1.61 yr. While this cycle is regular in the deep parts of the CZ, and in the northern hemisphere at mid-CZ, the southern hemisphere at that depth shows a much longer reversal time. Although the amplitude of the fields there fluctuates in time with the deeper reversals, they tend to maintain their sense for two or more cycles before flipping.
In Figures 2(e) and (f), we present the evolution of 〈B f 〉 in model Fb025 at mid-depth and near the bottom boundary, respectively. There, the highly axisymmetric fields cycle extremely regularly with a period of 5.08 yr. Like in the CZ of Sa2t, the cycle is poleward and hemispheric, traveling from the equator to −50°latitude. Unlike the irregular cycle of Sa2t, however, there is no intermittency in the cycle of model Fb025. Each cycle achieves very similar field amplitudes, and there is no evidence of the active hemisphere flipping across the equator.
Finally, in Figures 2(g) and (h) we show the evolutionary history of model Fa2t at mid-CZ and in the tachocline. Like in Sa2t and Fb025, the magnetic cycle in the CZ of Fa2t is largely restricted to the southern hemisphere, where it propagates poleward with an average period of 11.1 yr. Unlike the other models, the cycle of Fa2t is biased, with negative 〈B f 〉 wreaths tending to cling to the equator for longer before breaking poleward, and thus occupying a longer fraction of the cycle than positive wreaths. This is likely due to interactions with the fields in the tachocline, where a single wreath of negative 〈B f 〉 persists throughout the full duration of the simulation. This tachocline wreath fluctuates in amplitude in response to the reversals occurring in the CZ, but never reverses its sign.

Hemispheric Dynamo Action
The prevalence of dynamo states that are strongly localized within one hemisphere is a striking feature of our computational survey. Indeed, it appears to be a widespread phenomenon in convective dynamo models of deep-CZ and FC stars (e.g., BT20, Brown et al. 2020). In BT20, we speculated that such a field configuration could be achieved as a result of quasi-linear superposition of dynamo modes with even and odd parity with respect to the equator. This assessment aligns with work in the geodynamo literature, which finds that when the dipole and quadrupole are the two most critical dynamo modes and have comparable growth rates, weakly nonlinear interactions between them have the potential to lead to hemispheric dynamo states (Gallet & Pétrélis 2009). In particular, hemispheric states can become dominant when the symmetry of flow fields across the equator is broken. The degree of symmetry breaking required to instigate a hemispheric dynamo was found to depend on the difference in growth rates between the dipolar and quadrupolar dynamo modes. For the planetary geometries considered in Gallet & Pétrélis (2009), flow symmetry deviations of only a few percent, which occurred frequently via convective fluctuations, were sufficient to induce hemispheric dynamo states.
Although the deep CZs of our simulated M-dwarfs are reminiscent of the typical aspect ratios in geodynamo models, their thermodynamic and electrical properties exist in disparate regions of parameter space. Without detailed calculations of the mean field eigenmodes supported in our computational domain, it is difficult to conclude whether the models we present here satisfy the prescribed growth rate conditions. It is abundantly clear, however, that in the regime of strong magnetism, feedback from an asymmetric magnetic field can exaggerate and perpetuate hemispheric asymmetries in the flow field, providing a form of hysteresis. This feedback may explain why models Fa2t and Fb025 remain locked into a single hemisphere, whereas Sa2t chaotically flips between hemispheres, and the dynamo of Sc05 does not exhibit hemispheric preference at all.

Quenched Differential Rotation
Next, we will consider how the diverse magnetic configurations and amplitudes achieved across our data set impact the large-scale differential rotation (DR) in the simulated stars, again turning to our four selected models. for our four ambassador models Sa2t, Sc05, Fb025, and Fa2t, as well as for their hydrodynamical (HD) progenitors SaHt, ScH, FbH, and FaHt. From Figure 3(a), it is apparent that the solar-like DR with fast equators, slow poles, and radially inclined contours established in SaHt is largely intact in Sa2t, albeit at a diminished amplitude. In particular, the contrast at the surface from the equator to ±60°latitude

Characteristics of Magnetized DR
drops from 0.426R * to 0.378R * , a reduction of 11.2%. Additionally, the unbalanced diffusive spread into the RZ in the HD model has been sharply reduced and lifted toward R bcz in Sa2t. As shown in the prior Figure 1(d), the magnetic fields of Sa2t are strongest at the interface between the CZ and the RZ, where they are built up by the radial shear of the DR. These fields are responsible for driving the RZ toward solidbody rotation and enhancing the tachocline in Sa2t; however, the deep interior is not in equilibrium. Due to the simulation design with long diffusion times in the RZ, there remains a minute but unbalanced diffusive flux toward the inner boundary of the computational domain, which did not reach saturation within the allowed simulation time. Thus, although the sharpening of the tachocline by dynamo action is reminiscent of the magnetic confinement scenario observed in solar-like simulations by Matilsky et al. (2022), it is unclear whether it would persist after a full equilibration of the diffusive profile. In Figure 3(b), we can see that the rotational profile of ScH and Sc05 appears solar-like at low-and midlatitudes, but features a prograde vortex at each pole. In simulations of more massive stars, values of the Rossby number Ro = ω rms /2Ω * with ω the vorticity ω = ∇ × v approaching unity often signal a shift to antisolar rotation, with fast poles and slow equators. Thus, it is tempting to interpret the polar vortices seen in ScH and Sc05 as a sign that these models are occupying an intermediate space between solar and antisolar rotation. However, this rotational structure is also preferred in the fast-rotating sister model FcH, suggesting that it may be a consequence of the level of turbulence in the "c" models more than their rotational constraint. Further study is necessary to determine what implications this may have for the rotation of real stars. The relatively weak magnetic fields of Sc05 have very little effect on the rotational profile, and ΔΩ 60 drops by just 3.1% from 0.422Ω * to 0.409Ω * . In the more quickly rotating models Fb025 and Fa2t, Figures 3(c) and (d), respectively, the stronger magnetic fields lead to dramatic reductions of the overall DR. From FbH to Fb025, ΔΩ 60 drops from 0.127Ω * to 0.050R * , a decrease of 61%. Additionally, the concentration of magnetic fields in the southern hemisphere leads to a markedly asymmetric DR in Fb025. From FaHt to Fa2t, the DR is nearly entirely quenched. ΔΩ 60 drops by 74% from 0.133Ω * to 0.035, with what DR survives being almost exclusively in the northern hemisphere. As in Sa2t, the strong fields generated in the tachocline have offset the preexisting diffusive spread of DR into the RZ, but do not necessarily represent an equilibrium state for the shear there.

Angular Momentum Balances
In rotating spherical systems such as these, the primary mechanism driving the establishment of DR is the Reynolds stress (e.g., Tuominen & Rüdiger 1989;Rüdiger et al. 1998). Rotationally induced correlations between velocity components lead to equatorward transport of angular momentum in fastrotating systems, and poleward transport in particularly slow ones. In HD models, this transport is balanced partly by turbulent diffusion, and partly through a thermal wind meridional circulation that appears as a consequence of the temperature gradients associated with DR in a process known as gyroscopic pumping (e.g., Garaud & Bodenheimer 2010;Miesch & Hindman 2011). Our magnetic context allows for three additional primary mechanisms by which the DR can achieve balance. In what is known in the literature as the Malkus-Proctor effect (Malkus & Proctor 1975), the Lorentz force associated with the mean magnetic field may act directly upon the axisymmetric component of v f , which is precisely the DR. Similarly to the Reynolds stress, correlations in the turbulent magnetic field may also transport angular momentum through the Maxwell stress. Finally, magnetic feedback on the velocity fields at small scales can lead to a reduction of the Reynolds stress in a process known as Λ-quenching (Kitchatinov et al. 1994). The resulting balance can be understood in terms of five fluxes, representing angular momentum transport by Reynolds stresses ( RS  ), meridional circulations ( MC  ), viscous diffusion ( n  ), mean magnetic torques ( MT  ), and Maxwell stresses ( MS  ), the precise definitions of which can be found in BT20.
The latitudinal components of these fluxes are presented in the left column of Figure 4 for each of models Sa2t, Sc05, Fb025, and Fa2t. In Sa2t, which showed a modest reduction to its HD differential rotation, the equatorward RS  is balanced in roughly equal parts by n  and MS  . In Sc05, the magnetic torques are small, with the highly structured RS  balanced by the combination of MC  and n  . In Fb025, n  and MS  are comparable in amplitude; however, both are surpassed by MC  and MT  . In particular, while, Figure 4(e) would suggest that all three major magnetic feedback are taking place to inhibit the DR in the southern hemisphere of Fb025, the nearly symmetric quenching of the DR in the northern hemisphere where the magnetic fields are weak is entirely due to an asymmetric meridional circulation. Finally, in model Fa2t, the balance against RS  is achieved similarly to that in Fb025, but Maxwell stresses MS  have supplanted the mean magnetic torque in the active hemisphere.
In the right column of Figure 4, we present for these four models the evolving deviations of angular velocity from its mean state , 1 3 overlaid with contours of 〈B f 〉 f to indicate the phases and locations of cycling magnetic fields. In models Sa2t, Sc05, and Fa2t, there is no sign that the DR is responding to the cycling fields. This is an expected result, as in each case, the angular momentum transport by the mean magnetic field MT á ñ f  was not a significant factor in the overall balance. Figure 4(f) shows the fluctuations of angular velocity for model Fb025 with its markedly hemispheric dynamo action, and a regular, symmetric torsional oscillation is present. The formation of each subsequent wreath of magnetic field at the equator induces a local slowdown of the rotation. As the wreath propagates toward the south pole, flux transport by the meridional circulation mirrors its effects in the northern hemisphere. The mechanism by which Fb025 achieves near symmetry in its DR despite a strong influence from highly asymmetric magnetic fields is quite intriguing. From the flux balance, it is clear that the symmetrization is achieved via the thermal wind, and so we present in Figure 5 the aspherical temperature perturbations and corresponding meridional circulations in model Fb025 and its hydrodynamic progenitor FbH. In particular, we decompose the temperature in terms of a pressure term and an entropy term The pressure component of temperature T P is largely the same in both models due to their similar convective structures, featuring above-average temperatures at the poles and in a shallow layer at the equator, and cooler temperatures everywhere else. In FbH, the contrast of T S follows largely the same distribution as T P . The resulting meridional circulation features two strong cells in each hemisphere, which are symmetric across the equator. At the surface, a low-latitude cell drives a diverging flow away from the equator, and at the tangent cylinder, another strong cell reflects flows that diverge from the equator at the bottom boundary. Both of these circulations represent flows between relatively warm regions (shallow equator and poles) and cool ones (mid-latitudes and deep equator).
Unlike in FbH, the entropy component of the temperature T S in Fb025 is not symmetric across the equator. Instead, the southern hemisphere is substantially hotter than the north as a consequence of ohmic dissipation of the magnetic fields, which are strong there. As a result, the total temperature perturbation in Fb025 shows a southern hemisphere that is, in its entirety, warmer than the northern mid-latitudes. The meridional circulations of Fb025 accordingly connect the warm southern hemisphere to the cool northern hemisphere, and are not symmetric across the equator. This trans-equatorial circulation then serves to homogenize angular momentum between the two hemispheres, leading to the nearly symmetric DR we observe. In principle, this symmetrization mechanism for the DR occurs to some degree in all of our systems featuring significant magnetic asymmetry, and acts quickly enough that its signature is visible on the timescale of a single magnetic cycle.

Scalings of Magnetized DR
The same force balances explained in detail for the four models Sa2t, Sc05, Fb025, and Fa2t can be extended to describe every model in the data set. By considering the data set as a whole, we may examine the systematic variation of the DR and its sensitivities. In Figure 6(a), we present the ratio of three rms axial torques, averaged over time and the full CZ, and plotted against the ratio of magnetic energy density ME = B 2 /8π to the kinetic energy density of convection ) . The three torques correspond to the Reynolds stress RS , and the total Maxwell stress MS The ratio presented, τ MS /(τ RS + τ MC ), indicates the relative importance of magnetic and hydrodynamical forces in the maintenance of the DR. Two regimes are apparent: below approximately ME/ CKE = 0.7, the torque ratio climbs steadily with respect to the energy ratio, to which we fit a power-law relationship of ME CKE . Models Sa2t and Sc05 both belong to this regime, where the primary torque balance resembles that of the hydrodynamical models, despite the presence of magnetic fields. As the magnetic fields grow stronger, the ratio of the torques saturates to a value of about 0.5. Models Fb025 and Fa2t belong to this saturated regime, where the torques from the magnetic fields are equally important to the hydrodynamical torques. In the literature, this state is often referred to as magnetostrophic balance (e.g., Curtis & Ness 1986;King & Aurnou 2015;Augustson et al. 2019).
We next consider the differential rotation established in each of our models with respect to the same energy ratio ME/CKE. Measured rotational contrasts and energy densities for each model are reported in Table 2. Figure 6(b) shows the surface latitudinal contrasts ΔΩ 60 = Ω| eq − Ω| 60°o f each model, normalized against the contrast achieved in its hydrodynamical progenitor. Doing so, we find that all of our simulated M2 stars, regardless of Ra, Ro, Pm, or boundary conditions, obey the same relationship. Employing maximum likelihood estimation to the functional form with cross-validation to estimate the uncertainty, we find bestfit values of a = 1.079 ± 0.103 and b = 1.347 ± 0.133. Notably, this relationship does not describes the DR itself, but rather the degree to which the DR is quenched by magnetic activity in a star. As with the torque balance, we note two distinct regimes of behavior for the differential rotation. At relatively low magnetic amplitudes, with ME/CKE < 0.1, the quenching ratio converges to unity. This indicates that the magnetic fields in these models exerted virtually no influence over the large-scale DR, and that the balance there is effectively hydrodynamical. We observe that magnetic quenching of the differential rotation begins to set in for stars whose dynamo efficiency exceeds ME/CKE = 0.1, with a factor of 2 reduction achieved by ME/ CKE = 0.95, at the cusp of magnetostrophic balance.
In observational surveys of stellar DR (e.g., Barnes et al. 2005;Reinhold et al. 2013), rotational contrasts are typically found to be primarily sensitive to effective temperature, with only modest rotational dependence. Barnes et al. (2005) reported that the absolute surface shear scales as P rot 0.15 0.10 DW~ for stars in their sample. Reinhold et al. (2013) largely recovered this weak scaling in the massive Kepler catalog, but did not detect it among particularly fast rotators with P rot < 2 days. Individual observations on stars in that regime have repeatedly shown them to have very little differential rotation (e.g., Savanov 2012; Vida et al. 2014;Davenport et al. 2015). In the context of our findings, we interpret the bulk trend to describe the effectively hydrodynamical variation of ΔΩ, and the breakdown at P rot < 2 to reflect the possible onset of magnetostrophic balance. Alternatively, the apparently flat rotational dependence could come about from fitting across rising and falling regimes of ΔΩ centered on more moderate rotation rates.

Scalings of Stellar Variations
From what we have presented so far, it is clear that an abundance of different magnetic behaviors and dynamo configurations are present in our data set. Next, we seek to explore systematic trends in these variations with respect to our simulation parameters. To quantify these trends, we primarily utilize ANCOVA, a statistical model that generalizes linear regression to test the significance of categorical variables such as the treatment of the lower boundary of the CZ in our simulations, while controlling for continuous variation in response to other input parameters. Due to the long time series Note. The energy densities, defined in Section 5.1, are all given in units of 10 4 erg cm −3 and are averaged over the CZ for the mature lifetime of each simulation. The differential rotation contrast ΔΩ 60 /Ω * is measured at r = R o between the equator and ±60°latitude, and averaged over longitude and time.
The convection zone magnetic cycle periods P CZ are computed using Lomb-Scargle periodograms of 〈B r 〉 f and 〈B f 〉 f , and given in units of years. Figure 6. Scalings of axial torques and differential rotation across the data set. The colors red, green, and blue denote models of type "a," "b," and "c," respectively. Circular, square, and cross markers indicate the "PF," "PC," and "RZ" variations of "F" models, while triangles, diamonds, and stars denote the same for the "S" models. (a) The ratio of rms axial torques MS/(RS + MC), averaged over time and the full CZ and plotted against ME/CKE. Here, MS is the torque associated with the Maxwell stress, RS corresponds to the Reynolds stress, and MC is torque due to the meridional circulation. (b) The differential rotation contrast ΔΩ 60 /Ω * normalized against that of the hydrodynamical progenitor, and plotted against ME/CKE. that compose our averaged diagnostics, measurement errors are typically no more than 1%-2% for energy densities and magnetic field properties, and are less than 1% for all nondimensional numbers. Because these uncertainties are far exceeded by the unexplained variation in all of the fits we compute, we choose to employ a simple form of ANCOVA that neglects measurement error, calculated using the R statistical software (v4.1.2 R Core Team 2021).

Magnetic Field Amplitudes
We first consider the amplitude and configuration of magnetic fields across the simulated stars. In Table 2, we present averages over the CZ for five energy densities of relevance to the dynamo action taking place in each model. These are the convective kinetic energy (CKE), differential rotation kinetic energy (DRKE), poloidal magnetic energy (PME), toroidal magnetic energy (TME), and fluctuating magnetic energy (FME), precise definitions of which can be found in BT20.
Whereas the individual energy densities all show significant stratification with respect to the Rayleigh number of a simulation, we find that the energy ratio ME/CKE, where ME = PME + TME + FME, appears to have little dependence on Ra. Figure 7(a) shows the energy ratio plotted against the combination of predictors Rm/Ro 2 . The majority of simulations fall upon a unified line in log-log space, indicating a power-law relationship between the two quantities; however, there are eight notable exceptions. Models Sa05t, Sb05t, Sa1t, Sb1t, Fa05t, and Fa1t all yield smaller values of the CZ energy ratio than those models on the curve. This is unsurprising, given that the CZ-only pairs to these simulations (except Fa1) did not yield dynamos, and that the vast majority of their dynamo action occurs within the tachocline.
Two more exceptions exist in models Sc05 and Sc05c, which are not as easily explained. Given their slow rotation and small value of Pm, it is possible that the dynamos in these models are only marginally unstable, lending them similarity to the other tachocline driven dynamos. Alternatively, they possess the smallest values of Rm/Ro 2 of any CZ-only models in the data set, which may be sufficient to place them in an otherwise unsampled regime with different scalings. ANCOVA regressions were performed on a sample restricted to simulations with ME/CKE > 0.2, and with log-transformed continuous variables to fit power-law behaviors. Discriminating based on the Akaike information criterion (AIC; Akaike 1974), which establishes a logarithmic scale for the relative likelihood of a statistical model representing the true trend while penalizing overfitting, we identified the most likely model to be ) , which reduces to ME/KE ∼ PmRo −1 far from the onset of the dynamo instability. This Ro −1 scaling appears to hold true not just in theory, but in observations of real M-dwarfs as well (Reiners et al. 2022). By fitting coefficients for each term in the sum, we find that this trend is also upheld reasonably well in our simulations. Interestingly, we observe that although ME/ CKE shares the Ro −1 scaling with ME/KE, it differs in the amplitude of its Pm scaling. The most important difference between the two ratios is the inclusion of DRKE in KE, which, in our models, is less than CKE in only the most rotationally quenched CZs, and typically exceeds it by a factor of 3-10, leading it to dominate the KE in most situations. We posit that the reduced Pm dependence in ME/CKE then results from the greater ability for the magnetic fields to feedback upon the DR than upon the convection.
It is worth noting that despite the crashing differential rotation in the magnetostrophic regime, our detailed explorations in BT20 and BT22 show that the mean shear Ω-effect does not disappear as a significant contributor in the overall balance of magnetic energy generation there. To use the language of mean field dynamo theory, our less energetic dynamo models (ME/CKE 0.2) bear all of the hallmarks of classical αΩ dynamos. Their magnetic field generation is dominated by the rotational shear, and PME is many times smaller than TME. By contrast, our models approaching and exceeding equipartition in ME/CKE appear to have transitioned to α 2 Ω dynamos. The PME is a much more significant component of the global energy balance, ranging from 10%-50% of the TME, which reflects the shrinking role of differential rotation in the inductive balance. Nevertheless, that shear remains an important source term for all of the dynamos in our sample, indicating that they are not purely turbulent α 2 type dynamos.
In Figure 7(b) we present the fraction of the total magnetic energy contained within the axisymmetric components of the field (PME + TME)/ME. Unlike with the previous measures, the eight models whose fields are dominated by activity near BCZ appear to follow the global trend with their axisymmetric fraction, though the dispersion of that trend is greater. Fitting ANCOVA regression models once more, we find that the axisymmetric fraction in our data set is best described by PME TME ME It is important to note that there is no theoretical basis behind the functional form of this fit, and the resulting power laws should not be used to extrapolate to other parameter spaces. After controlling for this parametric variation, we find that the presence of an RZ and tachocline increases the axisymmetric fraction by 35% ± 13% on average relative to CZ-only models with the inert "potential field" lower magnetic boundary condition. We find this effect to be highly significant, with a statistical p-value of p = 0.0096. This effect is reproduced in our CZ-only models with "perfect conductor" lower magnetic boundaries, which showed an average increase of 35% ± 16% (p = 0.034).
In summary, we find that the energy ratio ME/CKE is very well predicted by a power law of the form Pm 0.5 Ro −1 . Furthermore, although the inclusion of a tachocline seems to have little influence over the amplitude of magnetic fields present in the CZ once sufficiently far from dynamo onset, we find that it does cause the distribution of magnetic energy to be significantly more axisymmetric than it might otherwise be. Additionally, the broad effect of a perfect conductor lower magnetic boundary on CZ field amplitudes and axisymmetry in our sample is very similar to that of a true RZ and tachocline, while being far less computationally expensive to model.

Field Evolution and Cycling
Next, we examine the periods of magnetic cycles for largescale fields in the CZs of our simulated stars. As with the four ambassador models, cycle times are computed via Lomb-Scargle periodograms of axisymmetric field components at fixed radius, and averaged over the band of latitudes in which they are prominent. The CZ cycle periods measured across the data set are shown in Figure 8 in units of the stellar rotation period. A wide range of scalings for the cycle period have been noted across different computational studies of dynamo action. Warnecke (2018) found in spherical wedge simulations of Sunlike stars that the cycle periods scaled as P/P rot ∼ Ro −0.98±0.04 . Star-in-a-box models by Käpylä (2022) found almost no sensitivity to the Rossby number P/P rot ∼ Ro −0.14±0.16 . In their work concerning spherical 3D models of solar-like stars with tachoclines, Guerrero et al. (2019) loosely constrained the Rossby scalings as a negative power for rapid rotators and positive for slow rotators. Past work by Strugarek et al. (2018) again considering Sun-like stars found that in regimes of moderate field amplitudes that do not significantly quench the DR, the cycle period scaled as P/P rot ∼ Ro −1.6±0.14 . For very strong fields, however, the dynamos were relocated to a much shallower depth in the stellar interior, and though they exhibited a similar scaling, their cycle periods were dramatically reduced in what they called the "short cycles" regime.
The models in our data set that achieve superequipartition field strengths do not appear to operate their dynamos in fundamentally different locations from those of the slower rotators. However, even our slower rotating models have values of Ro at the top of the domain that are just under 1. This would suggest that all of our models fall within their "short cycles" regime, thus precluding a break. Fitting for the variation of cycle times in our models, we find with additional minor dependencies on Ra, in distinct agreement with the scaling identified by Strugarek et al. (2018). Controlling for these variations, we find no significant effects on the CZ cycle period from the choice of magnetic boundary condition or the presence of a tachocline. The scalings that our M-dwarf models share with simulations of more massive stars in both their magnetic field saturation and their cycle periods lend further support to the emerging observational evidence that the dynamos operated within these stars are not fundamentally different.

Tachocline-dependent Spin-down
In BT20, we argued that the origin of the tachocline divide could lie in differences in the surface configurations and amplitudes of magnetic fields in stars containing tachoclines. We suggested that the abundance of highly active FC stars could result if the process of angular momentum loss through the magnetized stellar wind was less efficient in these stars, thus preserving their fast rotation and strong activity. The rate of angular momentum loss in a magnetized wind depends both on the magnitude of the fields and on the distance they extend from the surface of the star (Parker 1965), with field components of the lowest multipole orders, particularly the dipole, dominating the process. In BT20, we considered three pairs of simulated 0.4M e models with and without tachoclines, finding that in each case, the tachocline models yielded stronger average surface fields organized on larger scales, from which they argued that tachocline induced differences in spindown rates could qualitatively explain a tachocline divide. Here, we extend that analysis to a much larger and more diverse set of M-dwarf simulations.
In Figure 9(a), we present the time-and spherically averaged rms surface fields B surf for each simulation. As with the energy ratio, a power-law behavior is evident for most models, while the eight least-energetic dynamos fall far below the other amplitudes. Excluding these models and applying ANCOVA to B surf with AIC discrimination, we find the most likely scaling for this relationship to be  observationally in Reiners et al. (2022) for surface magnetic field amplitudes in the saturated and unsaturated regimes of the rotationactivity relation. Such a result is unsurprising, as the two target values for Ro in our survey were selected with the intent of straddling this transition. It is highly likely that sampling more rotation rates would reveal a more complicated structure in Rossby space than our fit here describes. Although the surface field strengths appear to be independent of the presence of a tachocline in our data set, the spin-down rate is in some ways even more dependent upon the configuration of these fields than their amplitudes. In Figure 9(b), we present the average dipole fraction of the surface fields f A A lm l m dip 1,0 2 , 2 = å where A l,m are the weights of the spherical harmonic decomposition of B r at the outer boundary of the simulation. Performing ANCOVA, we find that the continuous variation of f dip in our data set is best described by

( )
We stress again that a power-law form for the variation in this case is not derived from theory, but is instead just a useful descriptive framework. In addition to increasing the overall strength of fields, faster rotation further promotes greater organization of those fields. Under this model for the variation, we find that the presence of an RZ and tachocline in simulations increased their dipole fraction at the surface by an average of 78% (±20%, p = 0.0075). Despite finding that the choice of a perfect conductor lower magnetic boundary yielded similar results to a tachocline in the axisymmetry of fields throughout the CZ, such models show a much smaller effect on f dip that does not rise to the level of statistical significance. This likely results from the perfect conductor boundary promoting axisymmetry in toroidal fields more than poloidal fields, thus leading to little impact at the stellar surface.
To estimate the effect that this variation in f dip may have on an M-dwarf's spin-down, we turn to calculations by Matt et al. (2012) for the torque associated with a purely dipolar emergent field. They found that the wind torque scales as B w 0.87 * t~, among other stellar and fitting parameters, where B * is the strength of the dipole at the stellar surface. The effect we identified for the tachocline corresponds to a 33% increase in B * , and thus an average increase of 29% ± 6% to τ w under the assumption of dipole dominated spin-down.

Conclusions
We have presented the results of 45 global simulations of rotating M2-like stars under a variety of parameters and showcasing a rich diversity of behaviors. We have observed an abundance of asymmetric and single-hemisphere dynamos, some meandering between the northern and southern hemispheres while others remained fixed in place, some highly intermittent while others were quite regular. We suggest that these hemispheric dynamo states are a product of weakly nonlinear interactions between a dipolar and quadrupolar mode, which have nearly equal growth rates. Such interactions can lead to hemispheric dynamo states in the presence of equatorial symmetry breaking in the flow field, which we find can be reinforced by magnetic feedback when the magnetic fields are strong. Precise calculations of the mean field eigenmodes in our computational domain are still necessary to determine whether the proposed model for the hemispheric asymmetry is accurate.
We have observed models with independently time-steady or cycling fields in both the tachocline and the CZ, and in some cases multiple independent cycles in the CZ alone. Examining the durations of the dominant CZ cycles, we find that they were well predicted by the Rossby number, scaling as P CZ /P rot ∼ Ro −1.66±0.07 in concurrence with past trends identified in simulations of more massive stars.
We have also explored as an ensemble the variation of dynamo action through measures including the ratio of magnetic to convective kinetic energy, magnetic field axisymmetry, and the amplitudes and dipole fractions of surface magnetic fields. We find that the energy ratio ME/CKE in our models scales as Pm 0.52±0.06 Ro −1.02±0.06 , bearing similarities to a scaling for ME/KE proposed in theoretical work by Augustson et al. (2019), but exhibits weaker dependence on the magnetic Prandlt number than they find when including the energy of mean flows.
We find that when controlling for all other parameters, the presence of a tachocline confers no significant effect to the rms amplitudes of magnetic fields in the bulk of the CZ or at the stellar surface, nor to the durations of magnetic cycles in the CZ. However, we also find that including a tachocline causes the axisymmetric fraction of ME in a stellar CZ to increase by an average of 35% ± 13%, and the dipole fraction at the surface increases by 78% ± 20%. As the rms amplitudes of surface fields showed no dependence on the presence of a tachocline, this increase in f dip corresponds to 29% ± 6% greater stellar wind torques in M-dwarfs that contain tachoclines, assuming that spin-down is dominated by the dipole. This corroborates our earlier results in BT20, wherein we proposed that the absence of a tachocline in FC M-dwarfs could be extending their magnetic lifetimes by slowing their wind driven spindown. In concurrence with recent observations, this effect is achieved without requiring that M-dwarfs lacking tachoclines must have fundamentally different dynamo processes.
Magnetic quenching of the differential rotation is a prominent feature of our simulated M-dwarfs. We find that below the threshold value of about ME/CKE ∼ 0.3, the DR in the CZ is largely insensitive to magnetic activity. Above that threshold, Maxwell stresses begin to cause a significant reduction of the DR. We find that the balance between HD and magnetic torques was very well predicted by ME/CKE, as was the degree to which those torques quenched the DR, which indicated a clear distinction between weakly magnetized and magnetostrophic regimes. We observed this break to also coincide with a transition from αΩ to α 2 Ω type dynamos in our simulated stars both in this work and in BT20. Considering the variation of observed rotational contrasts on real stars, we find that the trends are qualitatively consistent with the quenching regimes we describe.
Although we have observed a number of tantalizing trends with this extensive data set, there are still important questions to which we continue to seek answers. Most pressingly: the comparisons we have made of spin-down rates are limited to M2 stars with and without tachoclines. Are the deep CZs of an M2 star an effective analog for FC stars, or does the changing shell geometry impose even greater effects than we have been able to observe in simulations of more massive stars? Explaining the origins of the tachocline divide with any degree of certainty will necessarily require exploring the internal dynamos of FC stars. We look forward to reporting on ongoing work regarding these small but fierce stars.
We thank Loren Matilsky, Ben Brown, and Sacha Brun for helpful discussions in the development of this work. We thank the reviewer for providing thorough comments, which improved this work. Computational resources for this project were provided by the NASA High End Computing (HEC) program through the Pleiades supercomputer at NASA Advanced Supercomputing (NAS) in the Ames Research Center, as well as local computational infrastructure. We thank Nick Featherstone for authoring the Rayleigh code, as well as the Computational Infrastructure for Geodynamics, which is funded by the National Science Foundation. This work was supported by NASA Astrophysical Theory Program grant No. NNX17AG22G and the FINESST grant No. 80NSSC20K1543. Additionally, the Heliophysics grant No. 80NSSC18K1127 provided partial shared support of the computational infrastructure needed to analyze these major dynamo simulations.