© The Authors 2023

1. Introduction

Magnetic fields in stars have been studied both theoretically and through observations, particularly magnetic fields of solar-type main-sequence stars (e.g., Brun & Browning 2017, and references therein). M dwarfs are low-mass main-sequence stars whose structure undergoes a transition from fully convective for masses up to 0.35 M_⊙ to solar-like (radiative core and convective envelope) for higher masses (Chabrier & Baraffe 1997). These stars are found to be magnetically active, as shown by Saar & Linsky (1985), who confirmed surface magnetic activity for M dwarfs with infrared measurements. Today there is considerable observational evidence of magnetic activity in M dwarfs that shows magnetic field strengths reaching up to a few kG (see Kochukhov 2021, and references therein). Because of the lack of a tachocline, the shear layer between the radiative and convective zones, fully convective M dwarfs are quite interesting from the point of view of dynamo theory and can help us to understand whether a tachocline has a strong impact on the dynamo itself. In this context, Wright & Drake (2016) report that the X-ray emission of fully and partially convective stars follows a similar trend to the Rossby number Ro = P_rot/τ, which is the ratio of the rotation period to convective turnover time, and which measures the rotational influence on convective flows. It was found that the X-ray emission increases with decreasing Ro until Ro ≈ 0.1, and for smaller Ro the X-ray luminosity saturates. Newton et al. (2017) find a similar trend, a saturated relation between the chromospheric Hα emission and Ro for rapidly rotating M dwarfs and a power-law decay of the Hα emission with increasing Ro for slowly rotating stars. The transition occurs near Ro = 0.2. In addition, Doppler and Zeeman-Doppler inversions have revealed that fully convective M dwarfs often show large-scale magnetic fields, and that for rapid enough rotation both dipolar and multipolar fields are possible (e.g., Morin et al. 2010; Kochukhov 2021). Furthermore, Klein et al. (2021) found that the fully convective star Proxima Centauri has a seven-year activity cycle.

Numerical simulations of stars are performed to achieve a better understanding of their magnetic fields, dynamos, and convection as functions of stellar parameters and dimensionless quantities, such as the magnetic Prandtl number, which is an intrinsic property of the fluid defined by the ratio of kinematic viscosity ν to resistivity η of the plasma. Some authors have performed magnetohydrodynamic (MHD) simulations of fully convective M dwarfs, which are particularly interesting for comparison with solar dynamo models due to the lack of a tachocline. The first simulations of fully convective M dwarfs were presented by Dobler et al. (2006), who used a star-in-a-box model to study dynamos as a function of rotation. They found predominantly quasi-static large-scale magnetic fields and typically weak or anti-solar differential rotation with faster poles and a slower equator. These simulations had relatively modest fluid and magnetic Reynolds numbers as well as low density-stratification. Browning (2008) presented simulations of fully convective M dwarfs using the anelastic magnetohydrodynamic equations, considering a spherical domain extending from 0.08 to 0.96 stellar radius, finding magnetic fields with significant axisymmetric components. In simulations without magnetic fields, the differential rotation is strong and solar-like with a fast equator and slow poles; instead, in magnetic simulations it is reduced, and tends to a solid-body rotation in the most turbulent magnetohydrodynamical simulations. A similar numerical approach was taken in the studies of Yadav et al. (2015, 2016) who used strongly stratified anelastic simulations to study the coexistence of dipolar and multipolar dynamos and cyclic solutions at relatively slow rotation corresponding to a parameter regime similar to that of Proxima Centauri. More recently, Brown et al. (2020) performed simulations of fully convective M dwarfs in spherical coordinates, finding cyclic hemispheric dynamos in their models.

The rotation period of the star P_rot is a key factor that determines the nature of the dynamo. This is evidenced by observational studies of M dwarfs, which demonstrate that with decreasing P_rot the magnetic field strength increases (e.g., Wright et al. 2018; Reiners et al. 2022). This has also been shown numerically by Käpylä (2021), among others, who used a star-in-a-box model for fully convective stars and found increasing magnetic field strength with decreasing rotation period. Furthermore, different dynamo modes were found as a function of rotation in that work. Slowly rotating stars have mostly axisymmetric and quasi-steady large-scale magnetic fields; for intermediate rotation the large-scale field is mostly axisymmetric and cyclic, and for rapid rotation the large-scale magnetic fields are predominantly non-axisymmetric with a dominant m = 1 mode. As demonstrated by Käpylä (2021), the large-scale dynamo is sustained even in the absence of a tachocline. In this sense, the work by Bice & Toomre (2020) using simulations of early M dwarfs supports the hypothesis that the tachocline is not necessary for producing strong toroidal magnetic fields, although it may generate stronger fields in faster rotators.

In this paper we present three-dimensional MHD simulations of fully convective M dwarfs with the star-in-a-box set-up described in Käpylä (2021; see also Dobler et al. 2006). Our main goal is to explore the dependence on dimensionless parameters, in particular the magnetic Prandtl and Reynolds numbers Pr_M and Re_M, which are crucial ingredients for dynamos and plasmas in general. High and low values of Pr_M and Re_M lead to very different dynamo scenarios; at low Pr_M the magnetic energy is dissipated in the inertial range of the flow and small-scale dynamo action requires a much higher Re_M to be excited (e.g., Schekochihin et al. 2007; Käpylä et al. 2018). On the other hand, stars typically have Pr_M ≪ 1 and Re_M ≫ 1 (e.g., Augustson et al. 2019; Jermyn et al. 2022). Our simulations were performed for a set of rotation periods P_rot ranging from 43 to 90 days, the latter being the rotation period of Proxima Centauri, and for values of Pr_M and Re_M ranging from 0.1 to 10 and 21 to over 1400, respectively, which is the numerically feasible range for this type of simulation. The methods and model are described in Sect. 2, while the description and analysis of the results is provided in Sect. 3. We discuss the conclusions in Sect. 4.

2. Methods

We used the star-in-a-box model described in Käpylä (2021), which is based on the set-up of Dobler et al. (2006). The model allows dynamo simulations of entire stars. In the present scenario we used a sphere of radius R that is enclosed in a cube with side 2.2 R. We solved the induction, continuity, momentum, and energy conservation equations:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}}{\partial t}& =\mathit{u}\times \mathit{B}-\eta {\mu }_0\mathit{J},\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}ln\rho }{\mathrm{D}t}& =-\mathbf{\nabla }·\mathit{u},\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}& =-\mathbf{\nabla }\mathrm{\Phi }-\frac{1}{\rho }(\mathbf{\nabla }p-\mathbf{\nabla }·2\nu \rho \mathsf{S}+\mathit{J}\times \mathit{B})-2\mathbf{\Omega }\times \mathit{u}+{\mathit{f}}_{\mathrm{d}},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill T\frac{\mathrm{D}s}{\mathrm{D}t}& =\frac{1}{\rho }[\mathcal{H}-\mathcal{C}-\mathbf{\nabla }·({\mathit{F}}_{\mathrm{rad}}+{\mathit{F}}_{\mathrm{SGS}})]+2\nu {\mathsf{S}}^2+\frac{{\mu }_0\eta {\mathit{J}}^2}{\rho }.\hfill \end{array}$$(4)

Here A is the magnetic vector potential; u is the velocity field; B = ∇ × A is the magnetic field; μ₀ is the magnetic permeability of vacuum; η is the magnetic diffusivity; ρ is the density of the fluid; D/Dt = ∂/∂t + u·∇ is the advective derivative; T is the temperature; Φ is the gravitational potential; p is the pressure; ν is the kinematic viscosity; s is the specific entropy; J = ∇ × B/μ₀ is the current density; $\mathbf{\Omega }={\mathrm{\Omega }}_0\widehat{z}$ is the rotation vector (with Ω₀ being the mean angular velocity of the star and $\widehat{\mathit{z}}$ the unit vector along the rotation axis); 𝖲 is the traceless rate-of-strain tensor

$$\begin{array}{c}\hfill {\mathsf{S}}_{\mathit{ij}}=\frac{1}{2}(u_{i,j}+u_{j,i})-\frac{1}{3}{\delta }_{\mathit{ij}}\mathbf{\nabla }·\mathit{u},\end{array}$$(5)

where the commas denote differentiation; ℋ and 𝒞 describe heating and cooling; and f_d describes the damping of flows outside the star (see Käpylä 2021, for more details). The radiative flux is given by

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{rad}}=-K\mathbf{\nabla }T,\end{array}$$(6)

where K corresponds to the Kramers opacity law, where its power law exponents are the same as in Käpylä (2021). The subgrid-scale (SGS) entropy flux F_SGS damps fluctuations near the grid scale, but contributes only negligibly to the net energy transport. It is given by

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{SGS}}=-{\chi }_{\mathrm{SGS}}\rho \mathbf{\nabla }{s}^{\prime },\end{array}$$(7)

where χ_SGS is the SGS diffusion coefficient, $s^{\prime }=s-{\overline{s}}_{\mathrm{t}}$ is the entropy fluctuation, and ${\overline{s}}_{\mathrm{t}}$ is a running temporal mean of the entropy. We note that the SGS flux used here does not include the temperature T. This form of the SGS flux is appropriate if the entropy equation has been solved, whereas the T factor appears in the SGS term if the corresponding energy equation has been solved (Rogachevskii & Kleeorin 2015).

The simulations were run with the PENCIL CODE¹ (Brandenburg 2021), which is a high-order finite-difference code for solving partial differential equations with primary applications in compressible astrophysical magnetohydrodynamics.

2.1. Dimensionless parameters

Each simulation is characterised by various dimensionless numbers. These parameters are usually order of magnitude ratios of various terms in the MHD equations or of the corresponding timescales.

The effect of rotation relative to viscosity is measured by the Taylor number, given by

$$\begin{array}{c}\hfill \mathrm{Ta}=\frac{4{\mathrm{\Omega }}_0^2{R}^4}{{\nu }^2}.\end{array}$$(8)

The Coriolis number is a measure of the influence of rotation on the flow

$$\begin{array}{c}\hfill \mathrm{Co}=\frac{2{\mathrm{\Omega }}_0}{u_{\mathrm{rms}}{k}_R},\end{array}$$(9)

where u_rms is the volume-averaged root mean square velocity and k_R = 2π/R is the scale of the largest convective eddies. Another definition of the Coriolis number used in other studies (e.g., Brown et al. 2020; Käpylä 2021) is based on the vorticity and considers the local length scale. This is defined by

$$\begin{array}{c}\hfill {\mathrm{Co}}_{\omega }=\frac{2{\mathrm{\Omega }}_0}{{\omega }_{\mathrm{rms}}},\end{array}$$(10)

where ω_rms is the volume averaged rms vorticity, with ω = ∇ × u. The fluid and magnetic Reynolds numbers, SGS and magnetic Prandtl numbers, and SGS Péclet numbers are defined as

$$\begin{array}{cc}\hfill \mathrm{Re}& =\frac{u_{\mathrm{rms}}}{\nu k_R},{\mathrm{Re}}_{\mathrm{M}}=\frac{u_{\mathrm{rms}}}{\eta k_R},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\mathrm{Pr}}_{\mathrm{SGS}}& =\frac{\nu }{{\chi }_{\mathrm{SGS}}},{\mathrm{Pr}}_{\mathrm{M}}=\frac{\nu }{\eta },\mathrm{Pe}=\frac{u_{\mathrm{rms}}}{{\chi }_{\mathrm{SGS}}{k}_R}.\hfill \end{array}$$(12)

2.2. Physical units and non-dimensional quantities

We modelled a main-sequence dwarf star (M5) with the same stellar parameters as in Dobler et al. (2006) and Käpylä (2021). The mass, radius, and luminosity of the star are M_⋆ = 0.21 M_⊙, R_⋆ = 0.27 R_⊙, and L_⋆ = 0.008 L_⊙, respectively. We used an enhanced luminosity approach (Käpylä et al. 2020) to reduce the gap between the thermal and dynamical timescales, such that fully compressible simulations are feasible. This implies that the results need to be scaled suitably for comparison with real stars. The conversion factor between the rotation rate, length, time, velocity, and magnetic fields in the simulation and in physical units are the same as those used by Käpylä (2021; see their Appendix A). Non-dimensional quantities were obtained by using the stellar radius as the unit of length [x]=R. Time is given in terms of the free-fall time ${\tau }_{\mathrm{ff}}=\sqrt{R^3/GM}$, the unit of velocity is [U]=R/τ_ff, and magnetic fields are given in terms of the equipartition field strength $B_{\mathrm{eq}}=⟨\sqrt{{\mu }_0\rho {\mathit{U}}^2}⟩$, where ⟨.⟩ stands for time and volume averaging.

3. Results

We present a set of 3D MHD simulations in the slow to intermediate rotation regime with global Coriolis number Co ranging between 3.1 and 12.9 (see Table 1). The rotation rates are $\stackrel{\sim }{\mathrm{\Omega }}=1.0$, 0.7, and 0.5 (which correspond to P_rot = 43, 61, and 90 days) in sets A, B, and C, respectively. The magnetic Prandtl number Pr_M varies between 0.1 and 10 in set A and between 0.5 and 10 in sets B and C.

3.1. Flow properties

3.1.1. Differential rotation and meridional circulation

The averaged rotation rate in cylindrical coordinates is given by

$$\begin{array}{c}\hfill \overline{\mathrm{\Omega }}(\varpi ,z)={\mathrm{\Omega }}_0+{\overline{U}}_{\varphi }(\varpi ,z)/\varpi ,\end{array}$$(13)

where ϖ = r sin θ is the cylindrical radius and where the overbar denotes azimuthal averaging. The averaged meridional flow is given by

$$\begin{array}{c}\hfill {\overline{U}}_{\mathrm{mer}}(\varpi ,z)=({\overline{U}}_{\varpi },0,{\overline{U}}_z).\end{array}$$(14)

The angular velocity varies with depth and also with latitude. A way to quantify this is by measuring the amplitude of the radial and latitudinal differential rotation with

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}=\frac{{\overline{\mathrm{\Omega }}}_{\mathrm{top},\mathrm{eq}}-{\overline{\mathrm{\Omega }}}_{\mathrm{bot},\mathrm{eq}}}{{\overline{\mathrm{\Omega }}}_{\mathrm{top},\mathrm{eq}}},\hfill \end{array}$$(15)

$$\begin{array}{cc}& {\mathrm{\Delta }}_{\mathrm{\Omega }}^{(\overline{\theta })}=\frac{{\overline{\mathrm{\Omega }}}_{\mathrm{top},\mathrm{eq}}-{\overline{\mathrm{\Omega }}}_{\mathrm{top},\overline{\theta }}}{{\overline{\mathrm{\Omega }}}_{\mathrm{top},\mathrm{eq}}},\hfill \end{array}$$(16)

where the subscripts top, bot, eq, and $\overline{\theta }$ respectively correspond to R = 0.9R, r = 0.1R, θ = 0°, and an average of $\overline{\mathrm{\Omega }}$ for latitudes +θ and −θ in spherical coordinates.

The values of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ and ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(\overline{\theta })}$ are listed in Table 2. Positive values of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ indicate solar-like differential rotation. We find that ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ has a tendency to decrease with increasing Re_M, which is equivalent to an increasing Pr_M (second column of Table 2). Figure 1 shows ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ as a function of Re_M for sets A (circles), B (squares), and C (triangles) confirming the decreasing trend as a function of Re_M. The red circles show simulations without dynamos (A8 and A12) and low magnetic Reynolds numbers, where ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ is higher. This reduction or quenching of the differential rotation by magnetic fields especially at high Re_M has been shown by various simulations, for example in Brun et al. (2004), Schrinner et al. (2012), and Käpylä et al. (2017). The differential rotation profiles for three representative simulations A1, A16, and A8 are shown in Fig. 2. The rotation profile in run A1 is solar-like; the profile is similar in the rest of the simulations with Pr_M ≤ 2. In run A16 with Pr_M = 10 the amplitude of the latitudinal differential rotation is positive, whereas the amplitude of the radial differential rotation is negative since the angular velocity does not change considerably with depth at the equator. The middle panel of Fig. 2 shows that the rotation rate at the equator is in fact higher than average almost everywhere and the negative value of ${\mathrm{\Delta }}_{\mathrm{\Omega }}^{(r)}$ is due to the slower-than-average rotation only very near the surface. Therefore, the differential rotation is solar-like. Another method for classifying the rotation profile (solar-like or anti-solar) is to use the mean rotation profile at the equator, which, as indicated in Käpylä (2023), can help prevent erroneous conclusions. Furthermore, the profiles in Fig. 2 are symmetric with respect to the equator, as in the other simulations performed in this work.

Fig. 1. Amplitude of radial differential rotation as a function of magnetic Reynolds number for simulations of sets A (circles), B (squares), and C (triangles). Cyan and red is for simulations with dynamo and without dynamo, respectively.

The global Coriolis number (see Eq. (9)) in the current simulations ranges from 3.1 to 13. All of our runs show solar-like differential rotation, which is consistent with Käpylä (2021), where the shift from anti-solar to solar-like differential rotation occurs for Coriolis number between 0.7 and 2. This is also consistent with the simulations of spherical shell convection by Viviani et al. (2018), which show that the transition occurs at around Co = 3 (see Table 5 in their work). More recently, Käpylä (2023) found that the transition from anti-solar to solar-like differential rotation depends on the subgrid-scale Prandtl number (Pr_SGS), such that solar-like differential rotation is more difficult to obtain at high Pr_SGS than at Pr_SGS ≤ 1. In this work all the simulations have Pr_SGS ≤ 1.

Simulations A8, A12, and C2 do not have dynamos, and they are considered to be kinematic cases. The right panel of Fig. 2 displays the rotation profile for simulation A8. It demonstrates that a faster-than-average angular velocity spans a broader latitudinal range and a narrower radial range when compared to dynamo simulations. This depicts the influence of a magnetic field on differential rotation. In the regime Pr_M < 2 of our simulations, the meridional flow is composed of multiple small cells, while in the regime Pr_M ≥ 2 the pattern is composed of two to three large cells which are symmetric with respect to the equator. The maximum values of the normalised meridional flow amplitude, ${\stackrel{\sim }{\overline{U}}}_{\mathrm{mer}}^{\mathrm{rms}}={(GM/R)}^{-1/2}{({\overline{U}}_{\varpi }^2+{\overline{U}}_z^2)}^{1/2}$, in the cases shown in Fig. 2 correspond to ${\stackrel{\sim }{\overline{U}}}_{\mathrm{mer}}^{\max }=0.009$, 0.005, and 0.023 for simulations A1, A16, and A8, respectively. The rms value of the meridional velocity is given in the fifth column of Table 2. In the simulations with no dynamo, A8, A12, and C2 (A8 in right panel in Fig. 2), the meridional circulation also exhibits similar multiple patterns, which are also symmetric with respect to the equator.

Fig. 2. Normalised time-averaged mean angular velocity $\overline{\mathrm{\Omega }}/{\mathrm{\Omega }}_0$ for simulations A1 (left), A16 (middle), and A8 (right). The coloured contours indicate $\overline{\mathrm{\Omega }}(\varpi ,z)$. The streamlines indicate the mass flux due to meridional circulation. The amplitude of the meridional flow speed (${\stackrel{\sim }{\overline{U}}}_{\mathrm{mer}}^{\max }$) is indicated in the lower right corner of each panel. The surface is indicated by the grey line, while the equator is indicated by the grey dotted line.

3.1.2. Power spectra and kinetic helicity

To characterise the convective flows we calculated the normalised kinetic energy power spectra (e.g., Viviani et al. 2018; Navarrete et al. 2022) from

$$\begin{array}{c}\hfill P_{\mathrm{kin}}=\frac{E_{\mathrm{kin},\ell }}{{\sum }_{\ell }{E}_{\mathrm{kin},\ell }},\end{array}$$(17)

where E_kin, ℓ is the kinetic energy of the spherical harmonic degree ℓ, which is calculated from the decomposition of the radial velocity field at the surface into spherical harmonics. Figure 3 shows P_kin as a function of ℓ for selected simulations. For the simulations with lower rotation rates and large-scale dynamos, the convective power is slightly shifted towards lower ℓ, with peaks between 16 and 20 for set A, 12 and 15 for set B, and 7 and 15 for set C. In simulations with no dynamo, the peak is at considerably larger scales at ℓ = 4. This demonstrates the suppression of large-scale convective flows by magnetic fields. This is similar to results from recent solar-like simulations that suggest that suppression of large-scale convection may be important to maintain a solar-like rotation profile in the Sun (e.g., Hotta et al. 2022; Käpylä 2023). Furthermore, the large-scale convective amplitudes are in general higher in cases with slower rotation, in agreement with linear theory (Chandrasekhar 1961) and various earlier simulations (Featherstone & Hindman 2016; Viviani et al. 2018; Navarrete et al. 2022).

Fig. 3. Normalised convective power as a function of ℓ for simulations A10, A13, A16, B2, B3, B5, C2, C3, C4, and C6.

Kinetic helicity, defined as $\mathcal{H}=\overline{\mathit{\omega }=·\mathit{u}}$ is an important component in the operation of the dynamo. It is a proxy of the α-effect, which is responsible for producing poloidal fields from toroidal fields (and vice versa) by rising or descending and twisting convective eddies (Parker 1955; Steenbeck et al. 1966). In all of our simulations the kinetic helicity is negative (positive) in the northern (southern) hemisphere, as is shown in Fig. 4 for run A1. This, combined with a solar-like differential rotation, suggests that an αΩ dynamo is operating, in which case the direction of propagation of the dynamo waves is polewards (Parker 1955; Yoshimura 1975). This is consistent with our findings, which are discussed in more detail in Sect. 3.2.

Fig. 4. Azimuthally averaged normalised kinetic helicity, $\stackrel{\sim }{\mathcal{H}}=\mathcal{H}(\varpi ,z)$/urmsωrms, for simulation A1 with Co = 9.2 and PrM = 0.2.

3.1.3. Convective energy transport

The luminosities corresponding to radiative, enthalpy, kinetic energy, cooling, and heating fluxes according to Eqs. (31)–(36) of Käpylä (2021) are shown in Fig. 5 for run C4. The enthalpy and kinetic energy fluxes dominate almost everywhere, except near the surface where the cooling becomes important. This is similar to the results of Brown et al. (2020) and to the rotating runs of Käpylä (2021). The total flux reaches somewhat less than 90% of the luminosity from the heating near the surface. A possible reason for this discrepancy is a non-negligible contribution from the SGS flux.

Fig. 5. Luminosity profiles of kinetic energy (purple), cooling (blue), heating (orange), and enthalpy (red) fluxes of run C4.

3.2. Dynamo variation

As shown in Table 1, the main differences between the simulations are the input parameters Pr_M and the rotation rate. In this subsection we present the effects of varying these parameters on the large-scale magnetic field.

3.2.1. Dependence on rotation

We explored simulations with fixed Pr_M and varying P_rot with values between 43 and 90 days. These values were determined using the conversion method outlined in Appendix A of Käpylä (2021). In order to compare the large-scale magnetic field at the three rotation rates used here, we chose runs with comparable magnetic Reynolds numbers and different rotation rates.

Three representative runs with Re_M < 100 from each set are A6 with Re_M = 75, B1 with Re_M = 84, and C3 with Re_M = 67. Run A6 shows cycles in the azimuthally averaged toroidal magnetic field, ${\overline{B}}_{\varphi }(R,\theta ,t)$, as shown in the middle top panel of Fig. 6. The cycles were computed using the empirical mode decomposition with the libeemd library (Luukko et al. 2016), as in Käpylä (2022). To determine the periods we used ${\overline{B}}_{\varphi }(R,\theta ,t)$, from the range −60° < θ < 60°. The cycle was determined by taking the mode with the largest energy, and counting the period from the zero crossings of that mode.

Fig. 6. Azimuthally averaged toroidal magnetic field near the surface of the star as a function of time. The name of the simulation is indicated at the bottom right of each panel.

The left middle panel of Fig. 6 shows ${\overline{B}}_{\varphi }(R,\theta ,t)$ for run B1, which also exhibits cycles. The reversals are periodic for most of the run, and this run also shows longer term modulation in the northern hemisphere towards the end of the run. The left bottom panel of Fig. 6 is for run C3. Unlike the runs just mentioned, C3 does not exhibit cyclic reversals. However, it does reveal the presence of a dipolar field, with a positive (negative) polarity in the northern (southern) hemisphere. At similar values of magnetic Reynolds number, the third column of Table 1 indicates a slight reduction in B_rms at lower rotation rates.

Three representative runs with higher magnetic Reynolds number (Re_M ≈ 300) and different rotation rates are A15 with Re_M = 300, B3 with Re_M = 315, and C4 with Re_M = 337. The right top panel of Fig. 6 shows ${\overline{B}}_{\varphi }(R,\theta ,t)$ of run A15. This run has irregular reversals with the field mainly distributed from mid-latitudes (±45°) to the equator. Near the poles the field is quasi-stationary. The middle centre panel of Fig. 6 is for B3, where a dipole with a few random reversals is visible with a predominantly negative (positive) polarity at the northern (southern) hemisphere. The azimuthally averaged toroidal magnetic field of C4 is shown in the centre bottom panel of Fig. 6, where a predominantly positive (negative) polarity. Mollweide projections of the radial magnetic field at the surface of runs A15, B3, and C4 are shown in Fig. 7, where the field is less intense for the runs with lower rotation. In this sense, the B_rms decreases with decreasing rotation rate from A15 to B3, while B3 and C4 have similar values. We find that in general the saturation level of the magnetic field increases with Re_M. This behaviour is likely related to the presence of a small-scale dynamo that produces magnetic fields at spatial scales that are of the same order of magnitude as that of the turbulence. While this was not the focus of our current study, it remains an important area for future research.

Fig. 7. Mollweide projection of the radial magnetic field Br for runs A15, B3 and C4 with ReM ≈ 300 and different rotation rates.

In Fig. 8 we show the ratio of the rotation period to cycle period as a function of the global Coriolis number. We find that P_rot/P_cyc ∝ Co^β with β = −1.30 ± 0.26. When considering the data points on the right of the figure, we find that β = −2.98 ± 1.13. The uncertainty in the slope indicates that we need to take these results with caution. Nevertheless, earlier studies have also found β < 0 (e.g., Strugarek et al. 2017, 2018; Warnecke 2018; Viviani et al. 2018) with global simulations of solar-like stars. Even when the domain of these simulations differs from the one presented here, the similarity in the relationship between the cycle period and the Coriolis number implies a likeness in the dynamo processes of solar-like and fully convective stars. Nevertheless, the negative slope found here differs from the positive slopes for the inactive and active branches from observations (Brandenburg et al. 1998, 2017). However, some simulations show β ≳ 0 (Guerrero et al. 2019; Käpylä 2022), but the cause of this behaviour is currently unclear.

Fig. 8. Rotation period normalised to the cycle period as a function of the global Coriolis number. The green circles are for ReM ∼ 100, blue for 70 < ReM < 85, and yellow for ReM ≤ 55.

3.2.2. Dependence on magnetic Reynolds and Prandtl numbers

Magnetic Prandtl numbers from 0.1 to 10 were used in the simulations. For all the current runs, the magnetic field is predominantly axisymmetric. When converted to physical units, the azimuthally averaged toroidal magnetic field reaches strengths ranging from 10 to 16 kG in our models. These values are higher than those of the reported observations which are up to a few kG (e.g., Kochukhov 2021). Set A has cycles for Pr_M ≤ 2 with periods ranging from 309 to 471 free-fall times, which correspond to 6.3–9.6 yr, when considering the same time conversion factor used by Käpylä (2021). Run B1 also shows cycles with a period of 274 free-fall times. Table 3 lists the values of ν, η, the cycle periods (if applicable) together with the corresponding standard deviation for all the simulations presented here. We found that the calculated length of the cycle periods of the runs of set A has a very slight increase when increasing the magnetic Reynolds number as $P_{\mathrm{cyc}}\propto {\mathrm{Re}}_{\mathrm{M}}^{\alpha }$ with α = 0.25 ± 0.14. Additionally, when considering the runs with similar Re_M and different Pr_M, we found that the cycle period is virtually independent of Pr_M in the parameter regime explored here.

The azimuthally averaged toroidal magnetic field ${\overline{B}}_{\varphi }(R,\theta ,t)$ is shown in Fig. 6 for a set of representative runs. The top panels are for three runs of set A, which have the same rotation period and increasing Re_M from left to right. The top left panel is for run A1 with Re_M = 55, with P_cyc = 320 ± 10 free-fall times. The top centre panel is for run A6 with Re_M = 75 and P_cyc = 326 ± 11 free-fall times. In these cases the field is distributed in latitudes |θ|≲80°. Simulations with higher values of Pr_M and Re_M, such as run A15 with Re_M = 300, result in the loss of cycles and the emergence of irregular solutions. Similar irregularities of dynamo solutions have previously been observed in simulations with high Re_M (e.g., Käpylä et al. 2017), but the exact mechanism is still unknown. In this case the field is distributed at latitudes |θ|≲50 and also exhibits quasi-stationary solutions near the poles.

The polarity of the field changes from the surface to r = 0.5 R. Figure 9 shows ${\overline{B}}_{\varphi }$ at r = 0.5 R for runs A1 and A15. In simulations with cycles, such as A1, the cycles are visible throughout the convection zone. However, for runs with higher magnetic Prandtl number (Pr_M > 2), such as A15, the azimuthally averaged toroidal magnetic field changes with depth and shows less clear magnetic structures in the deeper layers.

Fig. 9. Azimutally averaged toroidal magnetic field at r = 0.5 for simulations A1 (top) and A15 (bottom) with ReM = 55 and ReM = 300, respectively.

The second row of Fig. 6 shows three runs from set B that have the same rotation rate and increasing Re_M. The left panel is for run B1 with Re_M = 84, which exhibits a cycle with P_cyc = 274 ± 38 free-fall times, as well as longer reversals or disappearing cycles towards the end of the run. In this run the field is distributed at latitudes |θ|≲80°. The centre panel is for run B3 with Re_M = 315, which exhibits an irregular solution with few polarity reversals and predominantly quasi-static fields. In this case the field spans slightly less latitudinally, distributed at latitudes |θ|≲75°. The right centre panel of Fig. 6 shows run B5, which has the highest Pr_M and Re_M in this set with Pr_M = 10 and Re_M = 1360. The field is more concentrated towards the equator (at latitudes |θ|≲50°) with seemingly irregular reversals. Similar to run A15, B5 also has a quasi-stationary solution near the poles. The third row of Fig. 6 are for runs of set C, which have the slowest rotation in the present work, with increasing Re_M from left to right. Runs C3 with Re_M = 67 and C4 with Re_M = 337 show a predominantly quasi-static dipolar field, which spans latitudes |θ|≲75°. A similar dipolar field was reported by Moutou et al. (2017) for the fully convective and slow-rotating M dwarf GJ 1289. However, in our models the toroidal magnetic energy is dominant (see Table 3), whereas the large-scale magnetic field of GJ 1289 is purely poloidal. Run C6 with Re_M = 1419 has a field concentrated near the equator and the large-scale structures are less clear than in C3 and C4.

The radial magnetic field, B_r, also varies as a function of Re_M. Mollweide projections of the radial magnetic field for runs B1, B2, and B5, with increasing Re_M from left to right are presented in Fig. 10. The main differences here are the structure and maximum values of the magnetic field strength. The size of the structures in runs B1 and B2 are similar, but the strength of the field is slightly higher in B2. Run B5 has smaller field structures than in the previous cases, and the magnetic field strength is higher.

Fig. 10. Mollweide projection of the radial magnetic field Br for runs with the same rotation rate and different ReM, B1 with ReM = 84, B2 with ReM = 168, and B5 with ReM = 1360.

Table 3 lists the energy densities of the simulations. The total magnetic energy E_mag is a significant fraction of the kinetic energy density in all of the runs with dynamos, sometimes also exceeding it. One may expect that E_mag grows with increasing Re_M (Pr_M), as found in other works (e.g., Käpylä et al. 2017). In this regard, there is no discernible trend in the simulations shown here in terms of the variation of E_mag. Since the kinetic energy density, E_kin, decreases with increasing Pr_M the ratio E_mag/E_kin grows. The decrease in the kinetic energy can be explained because at large Pr_M it is converted into magnetic energy more efficiently.

At low Pr_M the kinematic exponentially growing regime lasts longer than in the simulations with high Pr_M. Figure 11 shows the evolution of E_mag and E_kin in the kinematic and saturated regimes for runs A4 and A16. The kinematic regime of simulation A4 lasted about ten times longer than the kinematic regime of run A16. It can be seen in simulation A4 that E_mag is amplified by six orders of magnitude. In the saturated regime both energies are comparable, such that E_kin is about 1.5 times E_mag. However, in simulation A16 the kinetic energy density is slightly reduced, while the magnetic energy density is increased by a factor of roughly 1.8. The Pr_M at which E_mag overcomes E_kin occurs at Pr_M > 1 for sets A and B, and at Pr_M > 2 for set C. A similar behaviour of the kinetic and magnetic energy densities was reported by Browning (2008) for simulations of fully convective stars. In run Cm2 of that work with Pr_M = 5, E_mag/E_kin ≤ 1, while Cm with Pr_M = 8 has E_mag/E_kin = 1.2.

Fig. 11. Time evolution of the kinetic (black) and magnetic (blue) energy densities in the kinematic (left) and saturated (right) regimes for simulations A16 with PrM = 10 (top) and A4 PrM = 0.5 (bottom). The energies are normalised by Ekin.

Table 3 also includes the energy densities of mean toroidal ($E_{\mathrm{mag}}^{\mathrm{tor}}$) and poloidal ($E_{\mathrm{mag}}^{\mathrm{pol}}$) magnetic fields (see Cols. 8 and 9). The mean toroidal $E_{\mathrm{mag}}^{\mathrm{tor}}$ accounts for up to 30% of total magnetic energy density and, in general, diminishes as Re_M increases. $E_{\mathrm{mag}}^{\mathrm{pol}}$ is less than 10% of E_mag for almost all simulations. In general, the ratio of the energy of the mean field to total energy decreases for high magnetic Reynolds numbers. Figure 12 shows the saturation level of the mean field as a function of Re_M for subsets of simulations from sets A and B. We do not find a clear trend in the saturation level of the mean energy as a function of the magnetic Reynolds number. A decrease in the mean energy with the inverse magnetic Reynolds number is usually associated with catastrophic quenching (e.g., Cattaneo & Vainshtein 1991; Brandenburg 2001). It can be interpreted as an outcome of magnetic helicity conservation, which becomes important as Re_M grows (e.g., Brandenburg & Subramanian 2005). Nevertheless, the boundary conditions in our simulations do allow magnetic helicity fluxes.

Fig. 12. Mean magnetic energy normalised by kinetic energy as a function of the magnetic Reynolds number. The blue markers are for the last four simulations of set A; the red markers are for the simulations of set B. Circles (triangles) are for runs with 2003 (5763) of resolution. The dotted line corresponds to a power law that shows what a decrease in mean energy with the inverse magnetic Reynolds number might look like.

Furthermore, the kinetic energy density of the differential rotation, $E_{\mathrm{kin}}^{\mathrm{DR}}$, and meridional circulation, $E_{\mathrm{kin}}^{\mathrm{MC}}$, are given in Table 3. For simulations with a dynamo $E_{\mathrm{kin}}^{\mathrm{DR}}$ decreases at higher Re_M, while for simulations A8 and A12 with no dynamo $E_{\mathrm{kin}}^{\mathrm{DR}}$ is significantly higher. More specifically, the runs without a dynamo in set A exhibit roughly five times higher $E_{\mathrm{kin}}^{\mathrm{DR}}$ than runs with a dynamo in the same set. This indicates magnetic quenching of differential rotation. In all of the simulations discussed here, $E_{\mathrm{kin}}^{\mathrm{MC}}$ is around 1–3% of E_kin, with the exception of C2, where $E_{\mathrm{kin}}^{\mathrm{MC}}\approx 0.12E_{\mathrm{kin}}$.

4. Summary and conclusions

We have performed a large number of simulations of fully convective M dwarfs using the star-in-a-box set-up presented in Käpylä (2021). We used the stellar parameters for an M5 dwarf with 0.21 M_⊙ at three rotation rates corresponding to P_rot = 43, 61, and 90 days, and varied the magnetic Prandtl number from 0.1 to 10. Our simulations explore the intermediate to slowly rotating regime. Consistent with previous work by Käpylä (2021), we find solar-like differential rotation in the simulations presented here.

We found different solutions for the large-scale magnetic field depending on the rotation period and the magnetic Prandtl number, which in our models fixes the magnetic Reynolds number. For the simulations with $\stackrel{\sim }{\mathrm{\Omega }}=1.0$ (set A) and Re_M ≤ 105, the large-scale magnetic field is cyclic, with P_cyc ranging from 309 to 471 free-fall times. In this set we found a slight increase in the length of the cycle period with increasing Re_M. For larger Re_M no clear cycles are found, and instead the behaviour of the magnetic field and its reversals become irregular. For the simulations with $\stackrel{\sim }{\mathrm{\Omega }}=0.7$ (set B) we found cycles for run B1 with Re_M = 84, while for higher values of Re_M the reversals are less regular, and instead a quasi-static configuration is found. For the highest Re_M the solutions become irregular. For the case with the lowest rotation rate ($\stackrel{\sim }{\mathrm{\Omega }}=0.5$) the field is mainly dipolar for Re_M ≤ 750. At higher magnetic Reynolds numbers the magnetic field is predominantly irregular and concentrated at mid-latitudes, with quasi-stationary fields near the poles. We note that in the three sets the large-scale field is irregular and concentrated near the equator for the highest Pr_M (Re_M). Additionally, the rms velocity increases for decreasing rotation for comparable Re_M. We also note that for a few of the simulations, particularly A8, A12, and C2, no dynamo was found because Re_M was below the critical value to drive a large-scale dynamo.

Furthermore, the ratio P_rot/P_cyc decreases with the Coriolis number, similar to the simulations of solar-like stars by Strugarek et al. (2017, 2018), Warnecke (2018), and Viviani et al. (2018). Our results confirm the important role of rotation and dimensionless parameters such as Re_M and Pr_M in determining the properties of fully convective dynamos. Depending on the parameters, the magnetic field can show a clear cyclic behaviour with the cycle period influenced by the rotation rate and dimensionless parameters such as Re_M and Pr_M. The large-scale magnetic field shows cycles for low and modest values of Re_M, but the cycles are lost for highest magnetic Reynolds numbers where irregular or quasi-static fields dominate. A similar loss of cyclic solutions was reported by Käpylä et al. (2017) who also increased Pr_M to increase Re_M. Whether the behaviour of the dynamo changes if Re_M is fixed and Pr_M is lowered is yet an open question. This is also closer to the parameter regime of late-type stars where Pr_M ≪ 1 and Re_M ≫ 1, but this parameter regime is extremely challenging numerically.

A very tentative comparison can be pursued with the Proxima Centauri system, where Klein et al. (2021) inferred a seven-year activity cycle. In principle, the activity cycle inferred in our simulations is in the range from five to nine years, and is thus consistent with the observed data. We note that this comparison is preliminary; even though the rotation rate we adopt here is similar, the magnetic Prandtl number is likely to be different, and even larger differences are expected when comparing the magnetic Reynolds number of the star with that of our simulations. It is well known that the solutions for the magnetic field depend on these parameters, and thus leads to uncertainty in the possible interpretation. It is nevertheless encouraging that the behaviour found in the simulations is quite similar.

Overall, the study presented here, to our knowledge, consists of the largest exploration of the parameter space for dynamo models of fully convective M dwarfs. Uncertainties remain, for instance regarding the role of the magnetic Reynolds number, which will still be much larger in realistic systems. While a clear signature of a small-scale dynamo is not found in our simulations, the expectation is that small-scale dynamos are present at larger magnetic Reynolds numbers, and interact with the large-scale dynamo, thereby changing the solution.

Acknowledgments

C.A.O.R., D.R.G.S. and J.P.H. thank for funding via Fondecyt Regular (project code 1201280). C.A.O.R., D.R.G.S. and R.E.M. gratefully acknowledge support by the ANID BASAL projects ACE210002 and FB210003. D.R.G.S. and R.E.M. gratefully acknowledge support by the FONDECYT Regular 1190621. D.R.G.S. thanks for funding via the Alexander von Humboldt – Foundation, Bonn, Germany. P.J.K. acknowledges finantial support from DFG Heisenberg programme grant No. KA 4825/4-1. F.H.N. acknowledges financial support by the DAAD (Deutscher Akademischer Austauschdienst; code 91723643). The simulations were made using the Kultrun cluster hosted at the Departamento de Astronomía, Universidad de Concepción, and on HLRN-IV under project grant hhp00052.

References

All Tables

All Figures

	Fig. 1. Amplitude of radial differential rotation as a function of magnetic Reynolds number for simulations of sets A (circles), B (squares), and C (triangles). Cyan and red is for simulations with dynamo and without dynamo, respectively.
In the text

Fig. 2. Normalised time-averaged mean angular velocity $\overline{\mathrm{\Omega }}/{\mathrm{\Omega }}_0$ for simulations A1 (left), A16 (middle), and A8 (right). The coloured contours indicate $\overline{\mathrm{\Omega }}(\varpi ,z)$. The streamlines indicate the mass flux due to meridional circulation. The amplitude of the meridional flow speed (${\stackrel{\sim }{\overline{U}}}_{\mathrm{mer}}^{\max }$) is indicated in the lower right corner of each panel. The surface is indicated by the grey line, while the equator is indicated by the grey dotted line.

In the text

	Fig. 3. Normalised convective power as a function of ℓ for simulations A10, A13, A16, B2, B3, B5, C2, C3, C4, and C6.
In the text

	Fig. 4. Azimuthally averaged normalised kinetic helicity, $\stackrel{\sim }{\mathcal{H}}=\mathcal{H}(\varpi ,z)$/urmsωrms, for simulation A1 with Co = 9.2 and PrM = 0.2.
In the text

	Fig. 5. Luminosity profiles of kinetic energy (purple), cooling (blue), heating (orange), and enthalpy (red) fluxes of run C4.
In the text

	Fig. 6. Azimuthally averaged toroidal magnetic field near the surface of the star as a function of time. The name of the simulation is indicated at the bottom right of each panel.
In the text

	Fig. 7. Mollweide projection of the radial magnetic field Br for runs A15, B3 and C4 with ReM ≈ 300 and different rotation rates.
In the text

	Fig. 8. Rotation period normalised to the cycle period as a function of the global Coriolis number. The green circles are for ReM ∼ 100, blue for 70 < ReM < 85, and yellow for ReM ≤ 55.
In the text

	Fig. 9. Azimutally averaged toroidal magnetic field at r = 0.5 for simulations A1 (top) and A15 (bottom) with ReM = 55 and ReM = 300, respectively.
In the text

	Fig. 10. Mollweide projection of the radial magnetic field Br for runs with the same rotation rate and different ReM, B1 with ReM = 84, B2 with ReM = 168, and B5 with ReM = 1360.
In the text

	Fig. 11. Time evolution of the kinetic (black) and magnetic (blue) energy densities in the kinematic (left) and saturated (right) regimes for simulations A16 with PrM = 10 (top) and A4 PrM = 0.5 (bottom). The energies are normalised by Ekin.
In the text

Fig. 12. Mean magnetic energy normalised by kinetic energy as a function of the magnetic Reynolds number. The blue markers are for the last four simulations of set A; the red markers are for the simulations of set B. Circles (triangles) are for runs with 2003 (5763) of resolution. The dotted line corresponds to a power law that shows what a decrease in mean energy with the inverse magnetic Reynolds number might look like.

In the text