Single-hemisphere dynamos in M-dwarf stars

M-dwarf stars below a certain mass are convective from their cores to their photospheres. These fully convective objects are extremely numerous, very magnetically active, and the likely hosts of many exoplanets. Here we study, for the first time, dynamo action in simulations of stratified, rotating fully convective M-dwarf stars. Importantly, we use new techniques to capture the correct full ball geometry down to the center of the star. We find surprising dynamo states in these systems, with the global-scale mean fields confined strongly to a single hemisphere, in contrast to prior stellar dynamo solutions. These hemispheric-dynamo stars are likely to have profoundly different interactions with their surroundings, with important implications for exoplanet habitability and stellar spindown.


INTRODUCTION
Tiny M-dwarfs are the most numerous stars in the universe and one of the most promising locations for finding habitable planets. In particular, they have many Earth-sized habitable planets, and the nearest habitable exoplanet, when found, is likely to orbit an M-dwarf (e.g., Dressing & Charbonneau 2015). M-dwarfs are intrinsically dim, and observations of them are hard, but our data set is growing rapidly and now includes longduration campaigns of slowly-rotating systems (Newton et al. 2018). The habitability of potential exo-Earths around other stars is at the mercy of the magnetic environment. M-dwarfs, in particular, display multikilogauss surface fields. Photometric monitoring infers incredibly large starspots (e.g., Newton et al. 2016) or many, many small spots (e.g., Jackson & Jeffries 2012), and the combination of strong fields and complex topologies leads exceptionally high levels of magnetic activity Corresponding author: Benjamin Brown bpbrown@colorado.edu (e.g., West et al. 2015). Mega-flares larger than any ever seen on the Sun (e.g., Schmidt et al. 2019) happen often enough to threaten the ability of planets to support life, or possibly even to retain atmospheres (e.g., Shkolnik & Barman 2014).
The interiors of type-M4 stars or later (masses of ∼ 0.35M or less) are convectively unstable from their nuclear burning cores to the photosphere. Despite the relative simplicity of their internal structure, we know very little about the internal processes that drive their magnetic dynamos. Only a limited number of simulations have studied global-scale dynamo processes in M-dwarfs (e.g., Browning 2008;Yadav et al. 2015aYadav et al. ,b, 2016. Worse, all of these have suffered a fundamental limitation: heretofore it has only been possible to simulate spherical shells, with the coordinate singularity at r = 0 avoided by making a "cut-out" in the middle of the simulation. The impact of this cut-out on the global solution has been unknown. Using new techniques in the open-source Dedalus pseudospectral framework, we compute the first spherical solutions that include r = 0 without any cut-out. These are the first dynamo models of truly, fully convective M-dwarf stars.
As in solar-type stars, our simulations of M-dwarfs have wreath-building convection zone dynamos that create global-scale fields. However, our M-dwarf simulations find one significant difference: We find an abundance of dynamo solutions localized to a single hemisphere. Single-hemisphere dynamo states have profound implications for stellar spin-down, exoplanet habitability, and the interpretation of observations of strong fields on stellar surfaces.

MODEL SYSTEM
The interiors of M-dwarfs are regions of low-Mach number convection with significant density stratification. To study convection and dynamo action in these stars, we solve the anelastic equations using the opensource, spectrally-accurate Dedalus framework (Burns et al. 2020). This work is enabled by new approaches to global-spherical domains including r = 0 Lecoanet et al. 2019).
The non-dimensional energy-conserving anelastic equations with MHD and momentum and thermal diffusion are (Brown et al. 2012;Vasil et al. 2013;Lecoanet et al. 2014): together with the Coulomb gauge ∇ · A = 0 and the anelastic ∇·u+u·∇ ln ρ = 0 constraints. S is a volume heating term and rotation is represented as Ω = Ωẑ. The magnetic field B = ∇×A, and the electrostatic potential φ and reduced pressure enforce the constraints. We use the ideal gas equation of state, and diffusive entropy heat flux Q = −κT ∇s (Lecoanet et al. 2014). The boundaries are constant entropy, impenetrable and stress free for velocity, and potential for magnetic field. The quantities ∇ ln ρ, ∇T , 1/ρ, ∇ ln T and S represent the internal M-dwarf structure. We non-dimensionalize this set of equations on a characteristic lengthscale L = R the radius of the domain (0.85R ), rotation timescale τ = 1/(2Ω), density ρ c , magnetic field B 2 = 16πρ c Ω 2 R 2 , and entropy S c σ 2 R 2 /(κ/ρ c ). The dynamic momentum, µ, and thermal diffusivities, κ, are constant. The magnetic diffusivity η ∝ ρ −1 . The global non-dimensional parameters are the Ekman number Ek, Rayleigh number Ra, Prandtl number Pr and magnetic Prandtl number Pm, where σ represents the radial extent of the nuclearburning region. The Rayleigh number has this dependance on the heating because mean entropy contrast across the domain scales with the internal heating to within an O(1) constant as ∆s ∼ ρ S c σ 2 /κ. The values are specified at the center r = 0 (denoted with a subscript c), which is the extremal value for Ek and Ra. In our simulations, the local Ekman Ek and Rayleigh Ra number depend on ρ(r)/ρ c . The Prandtl numbers Pr and Pm are constant.
We use the convective Rossby number as an input proxy for the degree of rotational constraint, This combination has a long history of usefulness for boundary-heated incompressible convection (Julien et al. 1996). In that context, the combination is independent of microphysical diffusivities and depends on the domain size. The situation here is less simple. The internal heating structure, S(r), and radiative diffusion κ(r), determine the superadiabatic entropy gradient. Even thought our definition produces a more complex parameter dependance, we find good scaling correlation between Ro C and the output measured local Rossby number; see Anders et al. (2019) for more details. Fully-convective M-dwarf stars are nearly adiabatic, and are well represented by polytropic solutions to the Lane-Emden equations. We simulate the inner 3 density scale heights (n ρ = 3) of an adiabatically stratified (m = 1.5) Lane-Emden polytrope, corresponding to the inner 85% of the star. Our background state is in excellent agreement with a stellar structure model computed with MESA (Paxton et al. 2011) for a 0.3M star of age 5 Gyrs, at a solar metallicity. 1 For S, we extract the volumetric heating and radiative cooling from the Top is the toroidal mangetic field B φ at mid-convection zone (r = 0.5). Bottom is the vector potential of the poloidal magnetic field A φ at the top of the domain (r = 1). Color gives the sense of polarity for both quantities. The system goes through many dynamo reversals of the global-scale field during the period simulated. The mean fields are predominantly in a single hemisphere (here the Northern). At all times the global poloidal fields have a dipolar component, with strong contributions from higher-order modes.
MESA model, and fit the radial dependence with (7) A nonlinear least squares fit gives σ ≈ 0.12 and Q 0 /Q 1 ≈ 11, with an L 2 error of about 3%. The simulations are normalized such that S(r = 0)σ 2 = 1.
In this letter we report on a solution with Ro C = 0.41, Ek = 2.5 × 10 −5 and Pr = Pm = 1. The system starts in an adiabatic state with small entropy noise perturbations. We initialize a weak poloidal magnetic field with spherical harmonic mode = 1, 2, m = 1. The simulation has L = 127 and N = 127, where L is the maximal degree of the spin-weighted spherical harmonics used in the spectral representation and N is the maximal degree of the Jacobi polynomials used in the radial expansion.
We computed many such dynamo simulations, but focus our attention on a single representative case.
The convection is a successful dynamo, building both small-and global-scale magnetic fields which change in strength and orientation over time. This simulation ran for t = 6250 Ω −1 , or about 0.3 diffusion times at the core r = 0 and about 6.3 diffusion times at r = 1. The internal heating from S drives a convectively unstable state. By roughly t ≈ 2000 Ω −1 , the potential, kinetic and magnetic energies all equilibrate; fluctuating around well-defined averages.
After time-averaging over the second half of the simulation, the RMS Reynolds number ≈ 240; the fluctuating Reynolds number ≈ 120; the RMS Rossby number Ro = |∇ × u|/2Ω ≈ 0.3. In this letter, mean quantities are the m = 0 components (denoted with angle brackets). Fluctuating quantities are those with the mean removed. The time-volume-averaged total magnetic energy is ≈ 25% or the total kinetic energy. Approximately 25% of the magnetic energy exists in mean toroidal fields, ≈ 68% exists in fluctuating fields, and the balance is in mean poloidal fields. The differential rotation contains ≈ 45% of the kinetic energy, with ≈ 55% in the fluctuating convection.
3. HEMISPHERIC DYNAMOS Figure 1 shows the surprising global-scale mean fields. The fields are organized in space and in time, and undergo regular periodic reversals of polarity in both the mean toroidal B φ and the mean poloidal field A φ . The fields share similarities with the "wreathy dynamo" states found in dynamo simulations of solar-type stars (e.g., Brown et al. 2010Brown et al. , 2011Nelson et al. 2013). In the fully convective dynamo simulations, the mean fields represent a significant fraction of the total magnetic field. We measure the global magnetic reversals using Lomb-Scargle periodograms of the mean magnetic field (e.g., Townsend 2010; VanderPlas 2018). The fields change with a period of P ≈ 190 Ω −1 , which corresponds to a polarity reversal roughly every 15 rotation periods. The temporal consistency of these cycles is high. The peak of the periodogram is clearly present in both B φ (normalized peak power ≈ 0.3, P ≈ 190.1 Ω −1 ) and A φ (peak power ≈ 0.4, P ≈ 191.8 Ω −1 ), while the aliasing Figure 2. Magnetic field topologies through a magnetic reversal. The left side shown the mean toroidal field B φ and poloidal vector potential A φ for the full domain, temporally averaged over one rotation period. The right side shows the radial magnetic field Br at the upper boundary r = 1. These are taken from periods of peak magnetic field strength during two successive dynamo cycles; the top row from t = 5950 Ω −1 , bottom row from t = 6025 Ω −1 .
peaks are at much lower amplitude ( 0.05 and 0.1 respectively). The surface-poloidal and depth-toroidal fields change with nearly the same phase.
It is surprising is that the mean fields of these fullyconvective dynamos generally occupy one hemisphere or the other, but not both. In rapidly rotating solar-type stars with convection-zone shells fields typically symmetrically fill both hemispheres (e.g., Brown et al. 2010Brown et al. , 2011Nelson et al. 2013;Augustson et al. 2013Augustson et al. , 2015. Importantly: all previous M-dwarf dynamo simulations in deep shells also find two-hemisphere solutions. Those simulations have either found mean fields symmetrically in both hemispheres (e.g., Yadav et al. 2015aYadav et al. , 2016 or have built mean fields that don't undergo very many reversals (Browning 2008).
The tendency towards mean fields in a single hemisphere is strong: the initial dynamo action (t 125 Ω −1 ) is symmetric in the two hemispheres, but strong southern-hemisphere fields rapidly replace the initial transients. At about t ≈ 1500 Ω −1 those strong fields shift to the northern hemisphere, where they remain for the period we have simulated. The other hemisphere is not devoid of field, but the mean fields there are much weaker than in the dominant hemisphere. We specu-late that the reversals would continue on some longer timescale if allowed. Figure 2 shows the detailed structure of the global magnetic structure during two phases of reversed polarity late in the simulation. The first happens around t ≈ 5950 Ω −1 averaged over one rotation period (∆t ≈ 2π Ω −1 ) when positive (red) polarity dominates the northern hemisphere. We see that a positive B φ wreath fills the star, reaching from the core to the surface. North of this structure and closer to the surface lie the weak remnant wreaths of the previous cycle. A distinctly weaker and non-dipolar poloidal field A φ accompanies the dominant toroidal structure. There is some global-dipole component, but it is much weaker than the higher-order modes. The near-surface cuts also show the poloidal-field structure; indicating a substantial quadrupolar component. As the wreath migrates up along the rotation axis, a new opposite-polarity wreath (blue) replaces it. This reversal is visible in the surface radial fields.
The convective flows fill the interior, as visible in Fig Figure 2; the differential rotation is time-averaged over about 10 rotation periods, centered on this time.
high-entropy plumes (red). The center at r = 0 is a region of vigorous dynamics, and plumes frequently cross the origin. Convective enthalpy transport (L h ) dominates the radial luminosity balance, except in the upper boundary layer where thermal diffusion takes over (L Q ). Though these simulations are stratified, the kinetic energy flux (L KE ) is small, similar to rotationallyconstrained simulations of solar-type stars (e.g., Brown et al. 2008). The sum of the luminosities closely matches the nuclear source derived from MESA (L S ) indicating that these models are in proper thermal equilibrium. This convection drives a robust differential rotation, shown by the mean angular velocity, Ω , with a fast equator and slow polar and core regions. The gradients of angular velocity are strong in both radius and latitude. These gradients likely build the mean toroidal field B φ via shear induction (a.k.a the "Ω-effect"). Correlations in the fluctuating flows and fields maintain the polar field. The latitudinal contrast in Ω is about 10%, and the mean magnetic fields coexist with this shearing flow. This is similar to results from rapidly rotating simulations of solar-type stars (Brown et al. 2011), but markedly different from the results of Browning (2008) or Yadav et al. (2015a) where the fields strongly quench the differential rotation. Whether differential rotation and fields can coexist likely depends on the Rossby number regimes of the dynamos, as discussed in Yadav et al. (2016).

IMPLICATIONS OF HEMISPHERIC DYNAMOS
This letter has focused on a single dynamo solution that shows striking hemispheric asymmetries. We find the same preference for single-hemisphere dynamos in every fully convective M-dwarf case computed to present; including a sweep of cases at Ek = 2.5 × 10 −5 and Ro C = 0.38-0.57, corresponding to RMS Reynolds numbers of 210-350 and Rossby numbers of 0.27-0.42. To check if single-hemisphere preference happens spuriously from a numerical source, we also rotated around thex rather than theẑ axis (for Ro C = 0.41, Ek = 2.5 × 10 −5 ). In this strange solution, cycling singlehemispheric dynamo action happened around the eastwest pole, relative to the numerical grid.
So far, we cannot determine why fully convective simulations so strongly prefer single-hemisphere dynamos; nor why these have not been found previously. Possibly, hemispheric dynamos are local to this region of parameter space. It is also possible that small core cutouts, as used in shell simulations, significantly change the allowed topologies of the global-scale fields even when those cutouts are small ( 10%R). Small solid objects are a relevant source of vorticity injection. In any case, we find that hemispheric dynamo states are robust and possibly ubiquitous in fully convective simulations, and this prompts us to consider their impacts if they exist in real stars.
Strong global-scale mean fields confined to a single hemisphere will couple dramatically differently to the outflowing stellar wind than a traditional dipole. The higher-order moments of field fall of much faster with distance and this likely produces weaker angular momentum transport and less efficient spindown in hemispheric-dynamo stars. The presence of meanfield in a single hemisphere could dramatically affect attempts to observe this field and lead to significant viewing angle effects; this could explain substantial differences in observed topologies for otherwise similar stars (e.g., Morin et al. 2010).
Lastly, we expect that the habitability around a hemispheric-dynamo star might significantly differ from one with a dipolar field. Hemispheric-dynamo stars would be less likely to turn nearby exoplanets into lava worlds via inductive heating (e.g., Kislyakova et al. 2017Kislyakova et al. , 2018, as the stellar field at the location of the planetary orbit will be much weaker. But the bad news is there might be much higher levels of magnetic activity, flares and coronal mass ejections; the mean-field topologies are more complicated than simple dipoles. Alternatively, though, space weather might show a strong directional preference. A planet with a low-inclination orbit might escape the worst. But even then we anticipate the preferred hemisphere to switch often enough to cause trouble for everyone. Altogether, fully convective stars have a rich landscape of dynamo states, and exploring this landscape shines new light on our understanding of the dynamos operating in solar-type stars. It seems many questions remain wide open. This work was supported by NASA LWS grant NNX16AC92G, NASA SSW grant 80NSSC19K0026, and NASA HTMS grant 80NSSC20K1280. Computations were conducted with support by the NASA High End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades with allocation GIDs s1647 and s2114.