Rossby and Magnetic Prandtl Number Scaling of Stellar Dynamos

Magnetic fields influence both the dynamics and evolution of stars. Detecting such magnetic fields can be difficult given that the bulk of their energy resides below the surface. Such magnetic fields either are generated through convective dynamos or are remnant fields left over from a star's formation situated in stable radiative zones. At some point in their evolution, nearly all stars have convective regions capable of generating magnetic fields. In particular, during the pre-main-sequence, most stars are fully convective and rapidly rotating, which likely leads to the development of magnetic fields that are then frozen into the forming radiative regions of these stars. This implies that there may be a strong connection between the dynamical magnetic fields in convection zones and the secularly evolving magnetic fields in radiative regions, unless significant magnetic boundary layers form to prevent an active codevelopment of the two regions (e.g., Gough & McIntyre 1998; Wood & Brummell 2018).

From the perspective of stellar evolution, the presence of a magnetic field in stable regions can lead to instabilities and angular momentum transport. Therefore, it is of keen interest to estimate the magnitude of magnetic fields generated in convective regions that may be radiatively stable in later stages of evolution (Mestel & Weiss 1987; Spruit 2002; Zahn et al. 2007; Mathis 2013). Doing so is quite challenging owing to the inherently nonlinear character of convective dynamos. However, such an undertaking has been attempted in Emeriau-Viard & Brun (2017). Alternatively, there exist several models of convection and its scaling with rotation rate and magnetic field, which can provide an estimate of the kinetic energy (KE) in convective regions. There are also several regimes for which scaling relationships for the ratio of magnetic energy (ME) versus KE may permit the estimation of the magnitude of the rms magnetic fields, where they are distinguished by the assumptions about which processes dominate the energetic and force balances.

The effects of astrophysical dynamos can be detected at the surface and in the environment of magnetically active objects such as stars (e.g., Christensen et al. 2009; Donati & Landstreet 2009; Donati 2011; Brun et al. 2015). One direct way to approximate the dynamics occurring within such an object is to conduct laboratory experiments with fluids that have some equivalent global properties, while observing their response to controllable parameters, such as the strength of a thermal forcing or a rotation rate. In those cases, all the observable scales of the system can be accounted for, from the global or driving scale to the dissipation scales. In practice, this has proven to be quite difficult when trying to mimic geophysical or astrophysical dynamos, but they are still fruitful endeavors (e.g., Gailitis et al. 2001; Laguerre et al. 2008; Spence et al. 2009). However, recent experiments with liquid gallium have shown that magnetostrophic states, where the Coriolis force balances the Lorentz force, seem to be optimal for heat transport (King & Aurnou 2015), which is interesting given the strong likelihood that many astrophysical dynamos are in such a state. Another method is to simulate a portion of those experiments, but these numerical simulations are limited in the scales they can capture: either an attempt is made to resolve a portion of the scales in the inertial range down to the physical dissipation scale (e.g., Mininni & Pouquet 2009; Mininni et al. 2009; Brandenburg 2014), or an attempt is made to approximate the equations of motion for the global scales while modeling the effects of the unresolved dynamical scales (e.g., Gilman 1983; Brun et al. 2004; Christensen & Aubert 2006; Strugarek et al. 2016; Varela et al. 2017).

These varying approaches to gathering data about the inner workings of convective dynamos provide a touchstone for thought experiments. Further, one can attempt to identify a few regimes for which some global-scale aspects of those dynamos might be estimated with only a knowledge of the basic parameters of the system (e.g., Christensen 2010; Jones 2011; Davidson 2013; Oruba & Dormy 2014) and in quasi-geostrophic systems as discussed in Calkins et al. (2016, 2017) and Calkins (2018). The following questions are examples of such parametric dependencies: how may the ME contained in the system change with a modified level of turbulence (or stronger driving)? How does the ratio of the dissipative length scales impact that energy, or how does rotation influence it? Establishing the global-parameter scalings of convective dynamos, particularly with stellar mass and rotation rate, is useful given that they provide an order-of-magnitude approximation of the magnetic field strengths generated within the convection zones of stars as they evolve from the pre-main-sequence to a terminal phase. This, in turn, permits the placement of better constraints on transport processes, such as those for chemical elements and angular momentum, most of which occur over evolutionary timescales.

In the hunt for a simple set of algebraic equations to describe the basic convective and dynamo processes at work in stellar convection zones, it is useful to consider the following set of MHD equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}=-{\boldsymbol{\nabla }}\cdot \left(\rho {\boldsymbol{v}}\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\rho \displaystyle \frac{\partial {\boldsymbol{v}}}{\partial t} & = & -\rho \left({\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}\right){\boldsymbol{v}}-2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{v}}-{\boldsymbol{\nabla }}P-\rho {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}\\ & & +\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,\sigma ,\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times \left[{\boldsymbol{v}}\times {\boldsymbol{B}}-\displaystyle \frac{4\pi \eta }{c}{\boldsymbol{J}}\right],\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{B}}=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}=-{\boldsymbol{\nabla }}\cdot \left[\left(E+P-\sigma \right){\boldsymbol{v}}+{\boldsymbol{q}}\right]+H,\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{v}}$ is the velocity, ${\boldsymbol{B}}$ is the magnetic field, ρ is the density, P is the pressure, and s is the entropy per unit mass. Moreover, the following variables are also defined: ${{\rm{\Phi }}}_{\mathrm{eff}}={\rm{\Phi }}+{\lambda }^{2}{{\rm{\Omega }}}^{2}/2$, where Φ is the gravitational potential, Ω is the rotation rate of the frame, and λ is the distance from the axis of rotation; the current is ${\boldsymbol{J}}=c{\boldsymbol{\nabla }}\times {\boldsymbol{B}}/4\pi$, with c being the speed of light; ${\sigma }_{{ij}}=2\rho \nu \left({{\boldsymbol{e}}}_{{ij}}-1/3{\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}{\delta }_{{ij}}\right)$ is the viscous stress tensor, where $e=1/2\left({\boldsymbol{\nabla }}{\bf{v}}+{\boldsymbol{\nabla }}{{\bf{v}}}^{T}\right)$ and ν is the kinematic viscosity; ${\boldsymbol{q}}\,=-\kappa {\rm{\nabla }}T$; κ is the thermal diffusivity; η is the magnetic diffusivity; and H is an internal heating rate per unit mass that is due to some prescribed exoergic process (e.g., chemical, nuclear, or otherwise). The total energy is $E=\rho {\text{}}{v}^{2}/2\,+{B}^{2}/8\pi +\rho {{\rm{\Phi }}}_{\mathrm{eff}}+\rho e$, where e is the internal energy per unit mass. Since only the most basic scaling behavior of the stellar system is sought, the following assumptions are made: the total energy of the system is conserved, and the system is in a nonlinearly saturated statistically stationary equilibrium.

In considering stellar dynamos, the choice of an appropriate magnetic Prandtl number arises: should one consider those defined by atomic diffusivities or those defined by turbulent diffusion mechanisms. If the former, then the typical Prandtl number can be quite small. If the latter, then the Prandtl number is typically near unity. Here it will be considered that the ordering of timescales is preserved even if the turbulent diffusivities are invoked; therefore, the thermal and magnetic Prandtl numbers will remain fixed at their atomic values, thereby preserving the hierarchy of length and timescales that they set, which is not without precedence (Schekochihin et al. 2004; Brandenburg & Subramanian 2005). On the other hand, with turbulent diffusivities, it would seem that all stellar dynamos would be in the moderate magnetic Prandtl number regime, where $\mathrm{Pm}\approx 1$.

If the timescale of the flows at a scale ℓ is given by ${\tau }_{{\ell }}={\ell }/{\text{}}{v}_{{\ell }}$, where ${\text{}}{v}_{{\ell }}$ is the characteristic velocity at a length scale of ℓ and the gradient operators for a magnetic and velocity field of scales ${{\ell }}_{B}$ and ℓ, respectively, yield ${\rm{\nabla }}{\boldsymbol{B}}\approx {B}_{{{\ell }}_{B}}/{{\ell }}_{B}$ and ${\rm{\nabla }}{\boldsymbol{v}}\approx {\text{}}{v}_{{\ell }}/{\ell }$, then the magnitude of each term in Equation (3) is considered as

$$\begin{eqnarray}&&\displaystyle \frac{{\text{}}{v}_{{\ell }}{B}_{{{\ell }}_{B}}}{{\ell }}=2\displaystyle \frac{{\text{}}{v}_{{\ell }}{B}_{{{\ell }}_{B}}}{{\ell }}+\displaystyle \frac{{\text{}}{v}_{{\ell }}{B}_{{{\ell }}_{B}}}{{{\ell }}_{B}}+\displaystyle \frac{\eta {B}_{{{\ell }}_{B}}}{{{\ell }}_{B}^{2}}.\end{eqnarray} \tag{ 6 }$$

Therefore, assuming that the hierarchy of scales is preserved along the turbulent inertial range and defining the scale Reynolds number as ${\mathrm{Re}}_{{\ell }}\,=\,{\text{}}{v}_{{\ell }}{\ell }/\nu$, a reformulation of the above equation yields

$$\begin{eqnarray}&&\displaystyle \frac{{{\ell }}_{B}^{2}}{{{\ell }}^{2}}{\mathrm{PmRe}}_{{\ell }}\left(1+\displaystyle \frac{{\ell }}{{{\ell }}_{B}}\right)\approx 1,\end{eqnarray} \tag{ 7 }$$

where $\mathrm{Pm}=\nu /\eta$ is the magnetic Prandtl number, which implies that ${{\ell }}_{B}/{\ell }\approx {\mathrm{Pm}}^{-1/2}$ as seen in both the kinematic and nonlinear regimes of time-dependent ABC-flow dynamos (Arnold & Korkina 1983; Brummell et al. 2001). Therefore, at large magnetic Prandtl number, the length scales of the magnetic field are smaller than those of the velocity field. At low $\mathrm{Pm}$, the expectation should be the opposite, that the generated magnetic fields roughly tend to be of larger scale than the velocity field. This does not preclude smaller-scale intermittent fields for the low-$\mathrm{Pm}$ regime or vice versa; rather, it is more a statement about which scales of the magnetic field contain the most energy: the large scales for low-$\mathrm{Pm}$ dynamos and the small scales for high-$\mathrm{Pm}$ dynamos as is sketched in Figure 1 and seen in nonlinear 3D dynamo simulations (Brummell et al. 2001; Schekochihin et al. 2004; Brandenburg & Subramanian 2005).

Zoom In Zoom Out Reset image size

Using the Braginskii definitions of the viscous and magnetic diffusivities (see Appendix B) and turning to stellar models, one can determine the ranges of magnetic Prandtl numbers within a star. For dynamos in stars with convective envelopes, the low-$\mathrm{Pm}$ regime is of greatest interest. Using the Sun as an example, the atomic $\mathrm{Pm}$ ranges from about 10⁻¹ in the core to roughly 10⁻⁶ at the photosphere, decreasing exponentially toward it. Other low-mass stars are roughly similar, with a low $\mathrm{Pm}$ throughout the bulk of their interiors. High-mass stars, in contrast, have relatively large atomic values of $\mathrm{Pm}$ ranging from about 10 throughout the convective core and remaining nearly unity until the near-surface region is reached, at which point it drops to about 10⁻².

The prescription for the atomic values of the Prandtl numbers given in Appendix A has been applied to models of stars near the zero-age main sequence (ZAMS) with the solar metallicity in the mass range between 0.03 and 100 M_⊙ computed with MESA (Paxton et al. 2011, 2013, 2015, 2018). The resulting convective-zone-averaged Prandtl numbers are shown in Figure 2, where it is clear that all stars possess convective regions with a low thermal Prandtl number ($\Pr =\nu /\kappa$, as seen in Figure 2(a)), whereas one can find two regimes of magnetic Prandtl number (Figure 2(b)). The existence of these two regimes is directly related to where the convective region is located. For massive stars with convective cores, the temperature and density averaged over the convective volume are quite high when compared to lower-mass stars with a convective envelope. Such a high temperature leads to a large magnetic Prandtl number. The dichotomy in magnetic Prandtl number implies that there may be two fundamentally different kinds of convective dynamo action in low-mass versus high-mass stars.

Zoom In Zoom Out Reset image size

Of the convective dynamo topics, it is the geodynamo that has received the most focused theoretical attention (e.g., Buffett 2000, 2009; Roberts & Glatzmaier 2000), but it also has a unique niche of its own not shared by most stellar dynamos. The Rossby number of mantle flows is estimated to be approximately 10⁻⁶ (Finlay & Amit 2011; Aurnou et al. 2015), which implies that the convective overturning times are of the order of several thousand years since Earth rotates in 1 day. Therefore, the estimated convective velocities are of order $1\mathrm{km}\,\ {\mathrm{yr}}^{-1}$ in the mantle. Thus, Coriolis forces take center stage, and one may be justified in ignoring the nonlinear advection of vorticity though not of heat or composition. Primordial magnetism is ruled out by ohmic losses that could dissipate core magnetism on timescales of order ${10}^{4}\mathrm{yr}$ and magnetic reversals in the recent past of about ${10}^{5}\mathrm{yr}$.

The recent spherical shell simulations of the geodynamo have focused on trying to minimize the Ekman number $\mathrm{Ek}$ (the ratio of viscous to Coriolis forces) to as small values as computationally feasible and to examine the resulting force balances achieved and the scaling relationships that may result. These are discussed in overview for two of the groups of modeling efforts (e.g., Christensen & Aubert 2006; King & Buffett 2013; Oruba & Dormy 2014; Schrinner et al. 2014; Yadav et al. 2016). The guiding notion is that convective motions in the real geodynamo should be strongly influenced by the resulting magnetic fields, within a so-called Magneto-Archimedean-Coriolis (MAC) balance involving magnetic or Lorentz forces, buoyancy as the Archimedean force, and the rotational effects of the Coriolis force. However, the current hero simulations attaining $\mathrm{Ek}$ of order 10⁻⁷ still do not attain a clear MAC balance, as discussed by Hughes & Cattaneo (2016) and Cattaneo & Hughes (2017) in offering other limiting equations that may be used to explore weak and strong magnetic field states (e.g., Dormy 2016).

In the spirit of the MAC balance, one approach to building a scaling relationship for the ratio of ME to KE considers the balances established in generating entropy, KE, and ME, as well as the force balance. To begin, note that the evolution of the ME is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }\right) & = & \displaystyle \frac{{\boldsymbol{B}}}{4\pi }\,{\boldsymbol{\cdot }}\,{\boldsymbol{\nabla }}\times \left[{\boldsymbol{v}}\times {\boldsymbol{B}}-\displaystyle \frac{4\pi \eta }{c}{\boldsymbol{J}}\right],\\ & = & -\displaystyle \frac{1}{c}\left[{\boldsymbol{\nabla }}\cdot \left(\eta {\boldsymbol{J}}\times {\boldsymbol{B}}\right)+\displaystyle \frac{4\pi \eta }{c}{{\boldsymbol{J}}}^{2}+{\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,\left({\boldsymbol{J}}\times {\boldsymbol{B}}\right)\right].\end{array}\end{eqnarray} \tag{ 8 }$$

If this equation is averaged over many dynamical times τ, when it is in a quasi-steady state, and if it is integrated over the volume of the convective region, it yields

$$\begin{eqnarray}&&\int \displaystyle \frac{{dt}}{\tau }{dV}\left[\displaystyle \frac{4\pi \eta }{c}{{\boldsymbol{J}}}^{2}+{\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,\left({\boldsymbol{J}}\times {\boldsymbol{B}}\right)\right]=-\int \displaystyle \frac{{dt}}{\tau }d{\boldsymbol{S}}\,{\boldsymbol{\cdot }}\,\left(\eta {\boldsymbol{J}}\times {\boldsymbol{B}}\right).\end{eqnarray} \tag{ 9 }$$

The Poynting flux of the right-hand side of Equation (9) can vanish at the boundaries of the convective region for an appropriate choice of boundary conditions. As an example, if the magnetic field satisfies a potential field boundary condition, then it is zero. If the field is force-free (e.g., ${\boldsymbol{J}}\propto {\boldsymbol{B}}$), then it is also zero. Supposing that this is the case, then one has that the average Lorentz work ($\int {dt}/\tau {dV}{\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,\left({\boldsymbol{J}}\times {\boldsymbol{B}}\right)$) is equal to the average Joule heating ${H}_{\eta }=4\pi \int {dt}/\tau {dV}\eta {{\boldsymbol{J}}}^{2}/c$. This is an important point, for it shows that the nature of the convection and magnetic field structures are directly impacted by the form of the resistive dissipation. Hence, the use of numerical dissipation schemes could yield unexpected results, and the boundary conditions may leak ME through a Poynting flux unless they are chosen carefully.

In a fashion similar to that used for the ME evolution above, one can find that the KE evolves as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{2}\displaystyle \frac{\partial \rho {{\boldsymbol{v}}}^{2}}{\partial t} & = & -{\boldsymbol{\nabla }}\cdot \left[\left(\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+P-\sigma \right){\boldsymbol{v}}\right]\\ & & +P{\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}-\sigma :{\boldsymbol{\nabla }}{\boldsymbol{v}}-\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}+\displaystyle \frac{{\boldsymbol{v}}}{c}\,{\boldsymbol{\cdot }}\,\left({\boldsymbol{J}}\times {\boldsymbol{B}}\right).\end{array}\end{eqnarray} \tag{ 10 }$$

To eliminate the compressibility term, the total internal energy must also be added to the system as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial t}\left[\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\displaystyle \frac{{{\boldsymbol{B}}}^{2}}{8\pi }+\rho e\right]\,\\ \quad =\,-{\boldsymbol{\nabla }}\cdot \left[\left(\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\rho e+P-\sigma \right){\boldsymbol{v}}-\displaystyle \frac{\eta }{c}{\boldsymbol{J}}\times {\boldsymbol{B}}\right]\\ \,-\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}-\displaystyle \frac{4\pi \eta }{{c}^{2}}{{\boldsymbol{J}}}^{2}-\sigma :{\boldsymbol{\nabla }}{\boldsymbol{v}}+H-{\boldsymbol{\nabla }}\cdot {\boldsymbol{q}}.\end{array}\end{eqnarray} \tag{ 11 }$$

One can also consider the time-averaged and volume-integrated evolution equation for this energy equation, which yields

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int \displaystyle \frac{{dt}}{\tau }{dV}\left[H-{\boldsymbol{\nabla }}\cdot {\boldsymbol{q}}-\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}-\displaystyle \frac{4\pi \eta }{{c}^{2}}{{\boldsymbol{J}}}^{2}-\sigma :{\boldsymbol{\nabla }}{\boldsymbol{v}}\right]\\ \quad =\displaystyle \int \displaystyle \frac{{dt}}{\tau }d{\boldsymbol{S}}\,{\boldsymbol{\cdot }}\,\left[\left(\displaystyle \frac{1}{2}\rho {{\boldsymbol{v}}}^{2}+\rho e+P-\sigma \right){\boldsymbol{v}}-\displaystyle \frac{\eta }{c}{\boldsymbol{J}}\times {\boldsymbol{B}}\right].\end{array}\end{eqnarray} \tag{ 12 }$$

The surface integral is zero if there are no outflows or net torque from the Lorentz force at the domain boundaries, implying the following:

$$\begin{eqnarray}&&\begin{array}{l}\int \displaystyle \frac{{dt}}{\tau }{dV}\left[\Space{0ex}{2.9ex}{0ex}H-{\boldsymbol{\nabla }}\cdot {\boldsymbol{q}}-\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}\right.\\ \,\left.-\,\displaystyle \frac{4\pi \eta }{{c}^{2}}{{\boldsymbol{J}}}^{2}-\sigma :{\boldsymbol{\nabla }}{\boldsymbol{v}}\right]=0.\end{array}\end{eqnarray} \tag{ 13 }$$

Note that $\int {dVH}=L(r)$, where L is the total luminosity of the star at a given radius r for a spherically symmetric heating. Likewise, the radiative luminosity of the star is given by $\int {dV}{\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,{\boldsymbol{q}}={L}_{r}(r)=-4\pi {r}^{2}\kappa \partial T/\partial r$ for the spherically symmetric component of the temperature field, which should be dominant. For the case of stars that are on the main sequence, there are three configurations of their primary regions of convection: either a convective core for high-mass stars, a convective envelope for lower-mass stars, or both for F-type stars (radiative luminosity of the star is given; see Figure 2). In these cases, one can assume that the region of integration is over the entire convection zone, and so the volume-integrated luminosity will be the nuclear luminosity, or the current total luminosity ${L}_{* }=L({R}_{* })$. Furthermore, the radiative luminosity will be nearly, but not exactly, equal in magnitude to L_*. The reason for this is that the thermal evolution of the system is a largely passive response to the changes in the nuclear burning rates. Because the nuclear luminosity is slowly increasing along the main sequence, ${L}_{R}={L}_{r}({R}_{* })$ will always lag behind L_* owing to the time required for thermal diffusion to modify the thermal gradient. From Equation (13), one can see that

$$\begin{eqnarray}&&{L}_{* }-{L}_{R}+{H}_{\nu }+{H}_{\eta }=-\int \displaystyle \frac{{dt}}{\tau }{dV}\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}={W}_{B},\end{eqnarray} \tag{ 14 }$$

where H_ν and H_η are the positive-definite, time-averaged, volume-integrated dissipation rates due to viscosity and resistivity, respectively. Thus, the rate of buoyancy work W_B is directly proportional to the mismatch of the two luminosities and the rates of viscous and resistive dissipation, implying that the latter result from the former.

Returning to the time-averaged curl of the momentum equation, though neglecting the viscous and inertial terms under the assumption that the Rossby number is sufficiently small and that the viscous terms play no role, one can find

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times \left[2\rho {\boldsymbol{v}}\times {\boldsymbol{\Omega }}+\displaystyle \frac{1}{c}{\boldsymbol{J}}\times {\boldsymbol{B}}\right]+{\boldsymbol{\nabla }}\rho \times {{\boldsymbol{g}}}_{\mathrm{eff}}=0,\end{eqnarray} \tag{ 15 }$$

where ${{\boldsymbol{g}}}_{\mathrm{eff}}=-{\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}$. This provides the basis of finding the length scales of the flows and the magnetic fields. Neglecting the inertial terms does mean that the domain of applicable Rossby numbers is restricted to being below unity, whereas neglecting the viscous terms is consistent with the assumption of a low magnetic Prandtl number as discussed in Davidson (2013). Specifically, the Rossby number is assumed to be small enough that the flow becomes roughly columnar and moderately aligned with the rotation axis. In such a case, there are two integral length scales of the flow: one parallel to the rotation axis, ${{\ell }}_{\parallel }$, and another perpendicular to it, ${{\ell }}_{\perp }$, with ${{\ell }}_{\perp }\lt {{\ell }}_{\parallel }$. These can be defined as ${{\ell }}_{\parallel }={\left(\langle {\text{}}{v}^{2}\rangle /\langle {\omega }_{z}^{2}\rangle \right)}^{1/2}$ and ${{\ell }}_{\perp }\,={\left(\langle {\text{}}{v}^{2}\rangle /\langle {\omega }^{2}-{\omega }_{z}^{2}\rangle \right)}^{1/2}$, where ${\boldsymbol{\omega }}={\boldsymbol{\nabla }}\times {\boldsymbol{v}}$, ${\omega }_{z}=\hat{z}\cdot {\boldsymbol{\nabla }}\times {\boldsymbol{v}}$, and $\langle \rangle$ denotes a volume average. Here, unlike Davidson (2013), the density stratification is retained. Therefore, for comparison purposes, note that the integral scale of the flow is then ${{\ell }}_{0}\,={\left(\langle {\text{}}{v}^{2}\rangle /\langle {\omega }^{2}\rangle \right)}^{1/2}={\left[{{\ell }}_{\perp }^{2}{{\ell }}_{\parallel }^{2}/\left({{\ell }}_{\perp }^{2}+{{\ell }}_{\parallel }^{2}\right)\right]}^{1/2}$. Moreover, the density perturbations rather than temperature perturbations are contained in the buoyancy work integral and the force balance below, which permits the treatment of the baroclinic flows present in the convection zones of most stars. The full density can be retained in the integral and in the scaling given that the gradient of the mean density is parallel to ${{\boldsymbol{g}}}_{\mathrm{eff}}$. Hence, their cross product vanishes, leaving only the product of the velocity and gradients of the density perturbations. Assuming further that the ME density per unit mass scales only with ${{\ell }}_{\parallel }$, since ${{\ell }}_{\perp }\lt {{\ell }}_{\parallel }$, and w_B, unit consistency requires that ${B}^{2}/(4\pi \rho )\,\propto F\left({{\ell }}_{\parallel },{w}_{B}\right)\approx {{\ell }}_{\parallel }^{2/3}{w}_{B}^{2/3}$, where w_B is the rate of buoyancy work per unit mass ${W}_{B}/{M}_{\mathrm{CZ}}$ and ${M}_{\mathrm{CZ}}$ is the mass contained in the convection zone. Therefore, given Equation (15) and assuming that each term is of the same order of magnitude, the basic proportionality is

$$\begin{eqnarray}&&\begin{array}{l}\rho {\boldsymbol{\Omega }}\cdot {\boldsymbol{\nabla }}{\boldsymbol{v}}\approx {\boldsymbol{\nabla }}\rho \times {{\boldsymbol{g}}}_{\mathrm{eff}}\approx \displaystyle \frac{1}{c}{\boldsymbol{\nabla }}\times \left({\boldsymbol{J}}\times {\boldsymbol{B}}\right)\\ \quad \Longrightarrow \ \displaystyle \frac{{\rm{\Omega }}{\text{}}{v}_{\mathrm{rms}}}{{{\ell }}_{\parallel }}\approx \displaystyle \frac{g}{{{\ell }}_{\perp }}\approx \displaystyle \frac{{B}^{2}}{4\pi \rho {{\ell }}_{\perp }^{2}},\end{array}\end{eqnarray} \tag{ 16 }$$

where ${\boldsymbol{\Omega }}\cdot {\boldsymbol{\nabla }}$ selects the parallel length scale ${{\ell }}_{\parallel }$, whereas ${\rm{\nabla }}\propto {{\ell }}_{\perp }^{-1}$ since ${{\ell }}_{\perp }\lt {{\ell }}_{\parallel }$. Thus, comparing the curl of the Lorentz force to the curl of the Coriolis force, one has

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Omega }}{\text{}}{v}_{\mathrm{rms}}}{{{\ell }}_{\parallel }}\approx \displaystyle \frac{{B}^{2}}{4\pi \rho {{\ell }}_{\perp }^{2}}\approx \displaystyle \frac{{{\ell }}_{\parallel }^{\tfrac{2}{3}}{w}_{B}^{\tfrac{2}{3}}}{{{\ell }}_{\perp }^{2}},\end{eqnarray} \tag{ 17 }$$

and moreover it can be shown that an estimate of the rms velocity is

$$\begin{eqnarray}&&{\text{}}{v}_{\mathrm{rms}}\approx \displaystyle \frac{{{\ell }}_{\parallel }^{\tfrac{5}{3}}{w}_{B}^{\tfrac{2}{3}}}{{{\ell }}_{\perp }^{2}{\rm{\Omega }}}.\end{eqnarray} \tag{ 18 }$$

Within the context of such estimates, the buoyancy term can be written as

$$\begin{eqnarray}&&{w}_{B}=\displaystyle \frac{\int \tfrac{{dt}}{\tau }{dV}\rho {\boldsymbol{v}}\,{\boldsymbol{\cdot }}\,{{\boldsymbol{g}}}_{\mathrm{eff}}}{\int \tfrac{{dt}}{\tau }{dV}\rho }\approx g{\text{}}{v}_{\mathrm{rms}},\end{eqnarray} \tag{ 19 }$$

which implies that the estimated magnitude for the curl of the buoyancy force in Equation (16) is $g/{{\ell }}_{\perp }\approx {w}_{B}/({\text{}}{v}_{\mathrm{rms}}{{\ell }}_{\perp })$. Then, it is easily seen that

$$\begin{eqnarray}&&\displaystyle \frac{{{\ell }}_{\perp }}{{{\ell }}_{\parallel }}\approx \displaystyle \frac{{w}_{B}}{{\rm{\Omega }}{\text{}}{v}_{\mathrm{rms}}^{2}}.\end{eqnarray} \tag{ 20 }$$

Hence, the ratio of integral length scales should vary as

$$\begin{eqnarray}&&\displaystyle \frac{{{\ell }}_{\perp }}{{{\ell }}_{\parallel }}\approx \displaystyle \frac{{w}_{B}^{\tfrac{1}{9}}}{{{\rm{\Omega }}}^{\tfrac{1}{3}}{{\ell }}_{\parallel }^{\tfrac{2}{9}}}={\left(\displaystyle \frac{{w}_{B}}{{{\rm{\Omega }}}^{3}{{\ell }}_{\parallel }^{2}}\right)}^{\tfrac{1}{9}}.\end{eqnarray} \tag{ 21 }$$

Likewise, the ratio of ME to KE then scales as

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\approx \displaystyle \frac{{B}^{2}}{8\pi \left(\tfrac{1}{2}\rho {v}_{\mathrm{rms}}^{2}\right)}\approx {\left(\displaystyle \frac{{w}_{B}}{{{\rm{\Omega }}}^{3}{{\ell }}_{\parallel }^{2}}\right)}^{-\tfrac{2}{9}}.\end{eqnarray} \tag{ 22 }$$

The Rossby number may then be defined as

$$\begin{eqnarray}&&\mathrm{Ro}=\displaystyle \frac{{\text{}}{v}_{\mathrm{rms}}}{{\rm{\Omega }}{{\ell }}_{\parallel }}\approx {\left(\displaystyle \frac{{w}_{B}}{{{\rm{\Omega }}}^{3}{{\ell }}_{\parallel }^{2}}\right)}^{\tfrac{4}{9}},\end{eqnarray} \tag{ 23 }$$

which implies that

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\propto {\mathrm{Ro}}^{-\tfrac{1}{2}}.\end{eqnarray} \tag{ 24 }$$

To make contact with stellar parameters, the buoyancy work may be estimated noting that ${W}_{B}\approx {{GM}}_{* }{M}_{\mathrm{CZ}}{\text{}}{v}_{\mathrm{rms}}/{R}_{\mathrm{CZ}}^{2}$, where ${R}_{\mathrm{CZ}}$ is the radius at which $\rho ({R}_{\mathrm{CZ}})={\rho }_{\mathrm{CZ}}$ and ${\rho }_{\mathrm{CZ}}$ is the average density of the convection zone. However, if one takes into account the scaling of the buoyancy work with Rossby number and magnetic Prandtl number, Equation (23) becomes an implicit relationship for the Rossby number that may be indeterminate for large magnetic Prandtl number (see Section 6 and Davidson 2013). This can be partially seen through Equation (18), which may be recast to solve for w_B as

$$\begin{eqnarray}&&{w}_{B}\approx {\left(\displaystyle \frac{{{\ell }}_{\perp }^{2}{\rm{\Omega }}{\text{}}{v}_{\mathrm{rms}}}{{{\ell }}_{\parallel }^{\tfrac{5}{3}}}\right)}^{\tfrac{3}{2}}={\text{}}{v}_{\mathrm{rms}}^{2}{{\rm{\Omega }}}^{\tfrac{3}{2}}\displaystyle \frac{\langle {\omega }_{z}^{2}{\rangle }^{\tfrac{5}{4}}}{\langle {\omega }^{2}-{\omega }_{z}^{2}{\rangle }^{\tfrac{3}{2}}}.\end{eqnarray} \tag{ 25 }$$

One simple method to estimate these quantities is to suppose that the energy-containing flows have roughly the same length scale as the depth of the convection zone and that the speed of the flows is directly related to the rate of energy injection (given here by the stellar luminosity) and inversely proportional to the density of the medium into which that energy is being injected (Augustson et al. 2012; Brun et al. 2017). This is encapsulated as ${\text{}}{v}_{\mathrm{rms}}\propto {\left[{L}_{* }/\left(2\pi {\rho }_{\mathrm{CZ}}{R}_{\mathrm{CZ}}^{2}\right)\right]}^{1/3}$, where ${\rho }_{\mathrm{CZ}}$ is the average density in the convection zone and ${R}_{\mathrm{CZ}}$ is the radius corresponding to this average density. Therefore, the buoyancy work scales as

$$\begin{eqnarray}&&{w}_{B}=f{\left(\frac{{L}_{\ast }{{\rm{\Omega }}}^{\frac{27}{10}}}{2\pi {\rho }_{{\rm{CZ}}}{R}_{{\rm{CZ}}}^{2} \langle {\omega }_{z}^{2}{ \rangle }^{\frac{1}{4}}}\right)}^{\frac{2}{3}},\end{eqnarray} \tag{ 26 }$$

where $f= \langle {\omega }_{z}^{2}{ \rangle }^{\frac{3}{2}}/ \langle {\omega }^{2}-{\omega }_{z}^{2}{ \rangle }^{\frac{3}{2}}$ is a flow asymmetry parameter. Thus, additional estimations of f, ω, and ω_z need to be made in order to fully approximate the buoyancy work in terms of basic stellar parameters. However, such a mixing-length velocity prescription only provides an order-of-magnitude estimate, as it too will depend on the level of turbulence and more fundamentally on the rotation rate (e.g., Featherstone & Hindman 2016; Brun et al. 2017; Augustson & Mathis 2018, 2019).

In addition to assessing the energetic balances of a dynamo, another way to assess the scaling behavior of the ME and KE in a dynamo is to find the balance of forces acting in the system when in a statistically steady but nonlinear regime. Before proceeding into the more delicate multiterm balances, and as a means of estimating the relative order of magnitude of each term that may contribute to the amplitude of the magnetic field, consider first the four pairwise Lorentz force balances in the nonlinearly saturated statistically stationary momentum equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}\approx \rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{\boldsymbol{v}}\Longrightarrow \displaystyle \frac{{B}_{{{\ell }}_{B}}^{2}}{4\pi {{\ell }}_{B}}\approx \displaystyle \frac{\rho {v}_{{\ell }}^{2}}{{\ell }},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}\approx 2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{v}}\Longrightarrow \displaystyle \frac{{B}_{{{\ell }}_{B}}^{2}}{4\pi {{\ell }}_{B}}\approx \rho {\rm{\Omega }}{v}_{{\ell }},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}\approx {\boldsymbol{\nabla }}P^{\prime} +\rho ^{\prime} {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}\Longrightarrow \displaystyle \frac{{B}_{{{\ell }}_{B}}^{2}}{4\pi {{\ell }}_{B}}\approx \epsilon \left(\displaystyle \frac{\overline{P}}{{\ell }}+\overline{\rho }{g}_{\mathrm{eff}}\right),\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}\approx {\boldsymbol{\nabla }}\,{\boldsymbol{\cdot }}\,\sigma \Longrightarrow \displaystyle \frac{{B}_{{{\ell }}_{B}}^{2}}{4\pi {{\ell }}_{B}}\approx \displaystyle \frac{\nu {\text{}}{v}_{{\ell }}}{{{\ell }}^{2}},\end{eqnarray} \tag{ 30 }$$

where the hydrostatic balance has been removed, and for a flow at the scale ℓ, for which there is an equivalent magnetic scale ${{\ell }}_{B}$ as discussed in Section 3.

Convective flows often possess distributions of length scales and speeds that are peaked near a single characteristic value. Thus, assuming a single characteristic scale for the magnetic and velocity fields, where ${{\ell }}_{B}={\ell }={{\ell }}_{0}$ for which ${\text{}}{v}_{{\ell }}={\text{}}{v}_{\mathrm{rms}}$ as in the previous section, for Equations (27)–30 under the anelastic approximation this respectively leads to

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\right|}_{I,\mathrm{rms}}\propto 1,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\right|}_{C,\mathrm{rms}}\propto {\mathrm{Ro}}_{\mathrm{rms}}^{-1},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\right|}_{P,\mathrm{rms}}\propto \displaystyle \frac{\epsilon }{{\mathrm{KE}}_{\mathrm{rms}}}\left(\overline{P}+\overline{\rho }{g}_{\mathrm{eff}}{\ell }\right),\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}}\right|}_{V,\mathrm{rms}}\propto {\mathrm{Re}}_{\mathrm{rms}}^{-1},\end{eqnarray} \tag{ 34 }$$

where the subscript I represents the inertial balance, C the magnetostrophic balance, P the buoyancy and pressure work balance, and V the viscous balance. In anelastic systems $\epsilon \approx {\mathrm{Ma}}^{2}$, where $\mathrm{Ma}$ is the Mach number, implying that the potential energy represented by $\overline{P}+\overline{\rho }{g}_{\mathrm{eff}}{\ell }$ is much smaller than the other possible balances for the bulk of most stellar convection zones since $\mathrm{Ma}\ll 1$ for most stars except near their photosphere. Since stars are often rotating fairly rapidly, taking, for instance, young low-mass stars and most intermediate- and high-mass stars, their dynamos may reach a quasi-magnetostrophic state wherein the Coriolis acceleration also plays a significant part in balancing the Lorentz force. Such a balance has been addressed and discussed at length in Christensen (2010) and Brun et al. (2015), for example.

To better characterize the force balance, consider Equation (2) while taking its curl, wherein one can see that

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{\omega }}}{\partial t}={\boldsymbol{\nabla }}\times \left[{\boldsymbol{v}}\times {{\boldsymbol{\omega }}}_{P}+\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma -{\boldsymbol{\nabla }}P\right)\right],\end{eqnarray} \tag{ 35 }$$

where ${\boldsymbol{\omega }}={\boldsymbol{\nabla }}\times {\boldsymbol{v}}$ and ${{\boldsymbol{\omega }}}_{P}=2{\boldsymbol{\Omega }}+{\boldsymbol{\omega }}$.

Taking the dot product of this equation with ${\boldsymbol{\omega }}$ gives rise to the equation for the evolution of the enstrophy. Integrating that equation over the volume of the convective domain and over a reasonable number of dynamical times such that the system is statistically steady yields

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int d{\boldsymbol{S}}\,{\boldsymbol{\cdot }}\,\left[{\boldsymbol{v}}\times {{\boldsymbol{\omega }}}_{P}+\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma -{\boldsymbol{\nabla }}P\right)\right]\times {\boldsymbol{\omega }}\\ \,+\displaystyle \int {dV}\left({\boldsymbol{\nabla }}\times {\boldsymbol{\omega }}\right)\,{\boldsymbol{\cdot }}\,\left[\Space{0ex}{3.2ex}{0ex}{\boldsymbol{v}}\times {{\boldsymbol{\omega }}}_{P}\right.\\ \,+\,\left.\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma -{\boldsymbol{\nabla }}P\right)\right]=0.\end{array}\end{eqnarray} \tag{ 36 }$$

If no enstrophy is lost through the boundaries of the convective domains, then the surface integral vanishes, leaving

$$\begin{eqnarray}&&\begin{array}{l}\int {dV}\left({\boldsymbol{\nabla }}\times {\boldsymbol{\omega }}\right)\,{\boldsymbol{\cdot }}\,\left[{\boldsymbol{v}}\times {{\boldsymbol{\omega }}}_{P}+\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma -{\boldsymbol{\nabla }}P\right)\right]\\ \quad =\ 0.\end{array}\end{eqnarray} \tag{ 37 }$$

This assumption effectively means that magnetic stellar winds will not be part of this scaling analysis. For timescales consistent with the dynamical timescales of the dynamos considered here, this is a reasonably valid assumption. Then, since ${\boldsymbol{\nabla }}\times \omega$ is not everywhere zero, the terms in square brackets must be zero, which implies that

$$\begin{eqnarray}&&\rho {\boldsymbol{v}}\times {{\boldsymbol{\omega }}}_{P}+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma -{\boldsymbol{\nabla }}P=0.\end{eqnarray} \tag{ 38 }$$

Taking the curl of this equation eliminates the pressure contribution and gives

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times \left[\rho {\boldsymbol{v}}\times {\boldsymbol{\omega }}+2\rho {\boldsymbol{v}}\times {\boldsymbol{\Omega }}+\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{c}+{\boldsymbol{\nabla }}\cdot \sigma \right]=0.\end{eqnarray} \tag{ 39 }$$

This is the primary force balance, being between inertial, Coriolis, Lorentz, and viscous forces. Here the buoyancy force vanishes through the first curl and the pressure force through the second one. However, another choice could be made where the density dependence is retained, leaving a baroclinic term that does not appear in earlier dynamo considerations, where the flows are typically barotropic. Taking fiducial values for the parameters and scaling the derivatives as the inverse of a characteristic velocity length scale ℓ, the scaling relationship for the above equation yields

$$\begin{eqnarray}&&\displaystyle \frac{\rho {{\rm{v}}}_{\mathrm{rms}}^{2}}{{{\ell }}^{2}}+\displaystyle \frac{2\rho {{\rm{v}}}_{{rms}}{\rm{\Omega }}}{{\ell }}+\displaystyle \frac{{B}^{2}}{4\pi {{\ell }}^{2}}+\displaystyle \frac{\rho \nu {{\rm{v}}}_{{rms}}}{{{\ell }}^{3}}\approx 0.\end{eqnarray} \tag{ 40 }$$

If it is assumed that the flows are self-similar at each scale, as in the Kolmogorov turbulence model (Kolmogorov 1941), then the ratio of the velocity length scale to the magnetic length scale follows from the argument in Section 3 as ${\ell }/{{\ell }}_{B}={\mathrm{Pm}}^{1/2}$. In this case, when divided through by $\rho {v}_{\mathrm{rms}}^{2}/{{\ell }}^{2}$, the previous equation simplifies to

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}\,\mathrm{Pm}}\propto 1+{\mathrm{Re}}^{-1}+{\mathrm{Ro}}^{-1}.\end{eqnarray} \tag{ 41 }$$

The Reynolds number used in the above equation is taken to be $\mathrm{Re}\,=\,{v}_{\mathrm{rms}}{\ell }/\nu$. Since the curl is taken, the approximation for the pressure gradient and buoyancy terms employed in Augustson et al. (2016) is eliminated from the force balance. However, the leading term of this scaling relationship is found to be less than unity, at least when assessed through simulations. Hence, it should be replaced with a parameter to account for dynamos that are subequipartition, leaving the following:

$$\begin{eqnarray}&&\displaystyle \frac{\mathrm{ME}}{\mathrm{KE}\,\mathrm{Pm}}\propto a+b{\mathrm{Ro}}^{-1},\end{eqnarray} \tag{ 42 }$$

where a and b are unknown a priori, as they depend on the intrinsic ability of the rotating system to generate magnetic fields, which in turn depends on the specific details of the system such as the boundary conditions and geometry of the convection zone. Instead, they must be fit for currently as is discussed below. As demonstrated in Augustson et al. (2016), Equation (42) may hold for a subset of convective dynamos, wherein the ratio of the total ME to the KE depends on the inverse Rossby number and a constant offset. The constant is sensitive to details of the dynamics, and in some circumstances it may also be influenced by the Rossby number. In any case, convective dynamos are sensitive to the degree of the rotational constraint on the convection, as it has a direct impact on the intrinsic ability of the convection to generate a sustained dynamo. Yet, even in the absence of rotation, there appears to be dynamo action that gives rise to a minimum ME state in the case of sufficient levels of turbulence. Hence, there is a bridge between two dynamo regimes: the equipartition slowly rotating dynamos and the rapidly rotating magnetostrophic regime, where $\mathrm{ME}/\mathrm{KE}\propto {\mathrm{Ro}}^{-1}$. For low Rossby numbers, or large rotation rates, it is even possible that the dynamo can reach superequipartition states where $\mathrm{ME}/\mathrm{KE}\gt 1$. Indeed, it may be much greater than unity, as is expected for Earth's dynamo (see Figure 6 of Roberts & King (2013)).

To assess which models provide an accurate description of the physics underlying convective dynamo action, first consider four increasingly general fits to the data with

$$\begin{eqnarray}&&{f}_{0}={a}_{0}+{b}_{0}{\mathrm{Ro}}^{-1},\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{f}_{1}={a}_{1}{\mathrm{Ro}}^{{b}_{1}}{\mathrm{Pm}}^{{c}_{1}},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{f}_{2}={a}_{2}+{b}_{2}{\mathrm{Ro}}^{{c}_{2}},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{f}_{3}={a}_{3}+{b}_{3}{\mathrm{Pm}}^{{c}_{3}}+{d}_{3}{\mathrm{Ro}}^{{e}_{3}}+{f}_{3}{\mathrm{Ro}}^{{g}_{3}}{\mathrm{Pm}}^{{h}_{3}},\end{eqnarray} \tag{ 46 }$$

where the latter function is effectively the first-order multivariate expansion of $\mathrm{ME}/\left(\mathrm{KE}\ \mathrm{Pm}\right)$.

Following Stelzer & Jackson (2013), the fitting procedure to determine the coefficients and their variance is the bounded nonlinear least-squares technique in log space, which provides one measure of the goodness of an individual fit as

$$\begin{eqnarray}&&{\chi }_{j}^{2}=\displaystyle \frac{\mathrm{ln}10}{n}\displaystyle \sum _{i=1}^{n}{\left(\displaystyle \frac{{\mathrm{log}}_{10}\left({y}_{i}\right)-{\mathrm{log}}_{10}\left({f}_{j}\left({\mathrm{Ro}}_{i},{\mathrm{Pm}}_{i}\right)\right)}{{\sigma }_{i}/{y}_{i}}\right)}^{2},\end{eqnarray} \tag{ 47 }$$

where y_i are the data $\mathrm{ME}/\left(\mathrm{KEPm}\right)$ for a given ${\mathrm{Ro}}_{i}$ and ${\mathrm{Pm}}_{i}$, and where ${\sigma }_{i}$ are the standard deviations of the temporal variance of the data, which can be determined for four sets of data from Augustson et al. (2016), Varela et al. (2016), Käpylä et al. (2017), and Strugarek et al. (2018). In order to ascertain the goodness of fit globally, the fitting procedure is applied to n subsets of the data, where each subset is constructed by removing the ith data point, following the cross-validation method. The cross-validation procedure yields a measure of the goodness of fit as

$$\begin{eqnarray}&&{\chi }_{\mathrm{CV},j}^{2}=\displaystyle \frac{\mathrm{ln}10}{n}\displaystyle \sum _{i=1}^{n}{\left(\displaystyle \frac{{\mathrm{log}}_{10}\left({y}_{i}\right)-{\mathrm{log}}_{10}\left({f}_{j}^{* }\left({\mathrm{Ro}}_{i},{\mathrm{Pm}}_{i}\right)\right)}{{\sigma }_{i}/{y}_{i}}\right)}^{2},\end{eqnarray} \tag{ 48 }$$

where the model function ${f}_{j}^{* }$ is constructed using the subset of data excluding the ith point. Finally, an additional measure of goodness of fit is the relative misfit defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{rel},j}=\sqrt{\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{\left(\displaystyle \frac{{y}_{i}-{f}_{j}\left({\mathrm{Ro}}_{i},{\mathrm{Pm}}_{i}\right)}{{y}_{i}}\right)}^{2}},\end{eqnarray} \tag{ 49 }$$

which again may be applied to each entire data set and within the cross-validation scheme as

$$\begin{eqnarray}&&{\chi }_{\mathrm{rel},{CV},j}=\sqrt{\displaystyle \frac{1}{n}\displaystyle \sum _{i=1}^{n}{\left(\displaystyle \frac{{y}_{i}-{f}_{j}^{* }\left({\mathrm{Ro}}_{i},{\mathrm{Pm}}_{i}\right)}{{y}_{i}}\right)}^{2}}.\end{eqnarray} \tag{ 50 }$$

Note further that, since each data subset for the cross-validation yields a different set of fitting parameters, one may construct a probability distribution for each of those parameters.

The fits resulting from the cross-validation method applied to the data set with known variance are shown alongside the data and their variances in Figure 3, with the caveat that the Rossby number variance for the data is not given in Käpylä et al. (2017). For the full data set, which includes data from Schrinner et al. (2014), Augustson et al. (2016), Käpylä et al. (2017), Varela et al. (2017), Strugarek et al. (2018), and Viviani et al. (2018), corresponds to the data shown in Figure 4. The fit parameters and their uncertainties for both data sets are detailed in Table 1. Since the variance of the Rossby number and ME-to-KE ratio was not available for the full data set, it is set to be a constant equal to the average variance of the cases that do possess this measurement. Only the statistically significant fits are shown in Figures 3 and 4, with the determination of significance made from the goodness-of-fit parameters given in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To assess the applicability of the force-balance and energy-based scaling regimes discussed in the previous sections, the magnetic Prandtl number and Rossby number dependence of the ratio of the ME and KE has been assessed in seven sets of simulation data. Five sets of this data are from 3D spherical MHD stellar simulations given in Tables 1 and 2: in Augustson et al. (2016) for core dynamo simulations of massive stars, Varela et al. (2016) for low-mass stars with a stable region, Strugarek et al. (2018) for Sun-like stars with an outer convective layer, and Käpylä et al. (2017) and Viviani et al. (2018) for Sun-like stars. Although one may also qualitatively consider the simulations presented in Mabuchi et al. (2015), Guerrero et al. (2016), and Warnecke (2018) for differential rotation and dynamo action in Sun-like stars that cover a range of rotation rates and magnetic Prandtl numbers, the data provided in those papers are not sufficient to compare them quantitatively to the other simulations. Specifically, the trends of the ME-to-KE ratio seen in those papers are consistent with those determined here. Unpublished data that were used to construct the data shown in those papers were reanalyzed to determine the temporal variability of the data from Augustson et al. (2016) and Varela et al. (2016). To extend this analysis into the domain of geodynamo simulations, data from a suite of 3D spherical MHD geodynamo simulations are also included from Table C.1 of Schrinner et al. (2014) and the online data provided in Yadav et al. (2016). These simulations range from being Boussinesq to highly stratified with both convectively stable and unstable regions. The bulk of these simulations have magnetic Prandtl numbers that are close to unity, with $\mathrm{Pm}$ being typically between 0.25 and 10.

With these data in hand, one may compute the best fit to the data for each of the fitting functions described above, with the statistically significant fits for Equations (43) and (45) shown in Figures 3 and 4. Histograms of the distributions of these parameters are shown for the statistically relevant fits in Figures 5 and 6. In Figure 5 most parameters have monomodal distributions. However, it would seem that there are two populations for the fitting parameter c₂, yet this is most likely an artifact given that the distribution for c₂ in Figure 6 is no longer strongly bimodal. The fitting functions that include a magnetic Prandtl number dependence, with the coefficients given in Tables 1 and 2, are not well constrained given the scatter in the data with respect to both magnetic Reynolds number and magnetic Prandtl number. Hence, the best that can be said about the scaling of the ratio of the ME to KE is with respect to the general trend of the data. Yet there may still be statistically significant departures in terms of magnetic Prandtl number, but a greater range of $\mathrm{Pm}$ must be probed in global-scale dynamo simulations to assess this, especially given the sparsity of the data at high Rossby number. Nevertheless, none of the fits are consistent with the energy-based scaling arguments, with the exception of a subpopulation in Figures 3 and 5. Hence, when considering the available simulation data, the buoyancy work Rossby number scaling of $\mathrm{ME}/\mathrm{KE}\propto {\mathrm{Ro}}^{-1/2}$ derived in Section 4 and shown in Equation (24) may be safely ruled out as a description of the ME scaling.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Within the context of the simulations exhibited in the full data set or Figure 4, the rotation rates employed lead to nearly three decades of coverage in Rossby number. In that figure, the force-based scaling derived in Section 5 and given by Equation (42) is depicted by the blue dashed curve, which does a reasonable job of describing the nature of the superequipartition state for a given Rossby number. In particular, the force balance studied here predicts a high Rossby number asymptotic regime. This is captured in the data of Augustson et al. (2016), in which there is a plateau of nearly equipartition states at high Rossby number, which points to a balance between inertial and Lorentz forces in those simulations. This behavior is of note and should provide a motivation for further investigating higher Rossby number dynamos. Surprisingly, the scaling law derived in Section 4 does not capture the behavior of this set of dynamos very well, in contrast to the many dynamo simulations and data shown in Christensen et al. (2009) and Christensen (2010), for which it performs well. This could be related to the fact that for large $\mathrm{Pm}$ the buoyancy work potentially has an additional $\mathrm{Ro}$ dependence (Davidson 2013). In particular, following Brandenburg (2014), the ratio of the dissipation rates can be described as

$$\begin{eqnarray}&&\displaystyle \frac{{H}_{\nu }}{{H}_{\eta }}=k{\mathrm{Pm}}^{n},\end{eqnarray} \tag{ 51 }$$

where k and n could be determined from a suite of direct numerical simulations. In particular, when kinetic helicity is injected at the driving scale, Brandenburg (2014) found that $k=7/10$ and $n=2/3$. Note that their results have also considered rotating driven turbulence and found that this scaling is effectively independent of the rotation rate and thus the Rossby number. However, other studies indicate that there may be a stronger rotational influence (Plunian & Stepanov 2010). Since there is ambiguity in that scaling, let $k(\mathrm{Ro})$ be an unknown function of the Rossby number. Subsequently, assuming that ${L}_{* }\approx {L}_{r}$ in Equation (14), the buoyancy work w_B per unit mass can be described as

$$\begin{eqnarray}&&{w}_{B}=-\displaystyle \frac{\int \tfrac{{dt}}{\tau }{dV}\rho {\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{eff}}}{\int \tfrac{{dt}}{\tau }{dV}\rho }=\left(1+k(\mathrm{Ro}){\mathrm{Pm}}^{n}\right)\displaystyle \frac{{H}_{\eta }}{{M}_{\mathrm{CZ}}}.\end{eqnarray} \tag{ 52 }$$

Such a scaling may provide a link between the moderate to high magnetic Prandtl number regime examined here and the low magnetic Prandtl number regime assessed in Davidson (2013), which will be the focus of later work. It may also possibly extend the analysis of Davidson (2013) to baroclinic flows often found in stars rather than purely barotropic flows that are more applicable to the geodynamo.

The characteristics of two descriptions of stellar and planetary dynamos have been discussed, where those dynamos are distinguished by their magnetic Prandtl number and by the Rossby number of the convection that drives them. Scaling relationships for the level of the partitioning of ME and KE have been derived from two different principles. Particularly, there may be a shift in the kind of dynamo action taking place within stars that possess a convective core and those that possess an exterior convective envelope if the atomic values of the diffusivities provide the basis of the effective diffusivities, as demonstrated in Figure 2. The scaling relationship between the ME and KE of such convective dynamos, in turn, provides an estimate of the rms magnetic field strength in terms of the local rms velocity and density at a particular depth in a convective zone, which will be of particular use to the stellar evolution community dealing with magnetic-related instabilities and induced diffusivities.

Given the atomic values of the diffusivities relevant to stellar interiors shown in Figure 2, the atomic value of the magnetic Prandtl number shifts from a very low regime in low-mass stars to a moderate $\mathrm{Pm}\gt 1$ regime in higher-mass stars. Thus, two scaling relationships for the ratio of the ME to KE appear to be the most applicable for stellar and planetary dynamos: one for the high or moderate magnetic Prandtl number regime and another for the low magnetic Prandtl number regime. Within the context of the moderate to large magnetic Prandtl number simulations in the current data set, the ME of the system scales as the KE multiplied by an expression that depends on the inverse Rossby number plus an offset, which in turn depends on the details of the nonrotating system (e.g., on Reynolds, Rayleigh, and Prandtl numbers). On the other hand, for low magnetic Prandtl number and fairly rapidly rotating systems, such as the geodynamo and rapidly rotating stars, another scaling relationship that relies on an energetic balance of buoyancy work and magnetic dissipation (as well as a force balance between the buoyancy, Coriolis, and Lorentz forces) may be more applicable. When focused on in detail, this yields an ME that scales as the KE multiplied by the inverse square root of the convective Rossby number. Such a scaling relationship has been shown to be fairly robust (Davidson 2013) for some dynamo simulations. With currently available simulation data, there is no clear way to assess whether there is a shift between these dynamo scaling regimes.

However, there are some hints from dynamo simulations that there may be an influence of the magnetic Prandtl number. Specifically, Brun et al. (2005), Featherstone et al. (2009), and Augustson et al. (2016) find that the dynamos of intermediate- and high-mass stars tend to have weak mean-field components and comparatively large nonaxisymmetric ME, whereas the bulk of lower-mass star simulations appear to have stronger mean fields (e.g., Augustson et al. 2013, 2015; Brun et al. 2017), as well as many other simulations, such as those of Mabuchi et al. (2015), Guerrero et al. (2016), Varela et al. (2016), Strugarek et al. (2018), Viviani et al. (2018), and Warnecke (2018). Such behaviors comport with modern dynamo theory (e.g., Brandenburg & Subramanian 2005), where if the atomic values of the diffusivities as shown in Figure 2 influence the effective Pm, then it would be expected that low-mass stars should exhibit stronger mean magnetic fields, which appear to be characteristic of low-Pm dynamos and core convective dynamos of massive stars that exhibit small-scale, high-Pm dynamo characteristics.

Indeed, additional work is needed to establish more robust scaling relationships that cover a greater range in both magnetic Prandtl number and Rossby number, among other parameters. The current suite of simulations does not exhibit a statistically significant Pm dependence in the data, which is likely due to the spread in the supercriticality of these simulations and the restricted range of magnetic Prandtl numbers considered. Thus, numerical experiments should also explore a larger range of Reynolds number and levels of supercriticality. Indeed, as in Yadav et al. (2016), some authors have already attempted to examine such an increased range of parameters for the geodynamo. Nevertheless, to be more broadly applicable in stellar physics, there is a need to find scaling relationships that can bridge both the low and high magnetic Prandtl number regimes that are shown to exist within main-sequence stars.

The authors thank the anonymous referee for their helpful comments that improved the data analysis and strengthened the conclusions. The authors further thank A. Strugarek and J. Varela for the use of their data, as well as S. Mathis for extensive discussions and helpful advice. The authors also thank S. Tobias for his helpful comments. K.C.A. acknowledges support from the ERC SPIRE 647383 grant. A.S.B. acknowledges funding by ERC WHOLESUN 810218 grant, INSU/PNST, and CNES Solar Orbiter. J.T. was partly supported by NASA through grants NNX14AB56G and NNX17AG22G.

The mean value and standard deviation of the fitting coefficients, which are determined through the cross-validation analysis, are given below in Table 1 for the fits to the data shown in Figures 3 and 4. Figure 5 shows the histograms of the fitting coefficients for the functions f₀ and f₂ in Figure 3, which are the only functions of statistical relevance studied here. Likewise, Figure 6 shows the same for Figure 4.

The goodness-of-fit results for the various fits and measures as discussed above are provided in Table 2.

Table 2.  Goodness-of-fit Parameters for the Data of Figures 3 and 4

fj	${\chi }_{\mathrm{CV},j}^{2}$	${\chi }_{\mathrm{rel},\mathrm{CV},j}$	${\chi }_{j}^{2}$	${\chi }_{\mathrm{rel},j}$
Figure 3
f0	1.06	9.03	1.13	9.20
f1	46.9	11.5	2.07	13.8
f2	14.5	6.48	1.28	9.67
f3	157	12000	53.7	269000
Figure 4
f0	1.06	5.99	1.07	5.72
f1	6.81	23.23	4.38	12.1
f2	1.04	6.13	1.08	6.57
f3	4.57	60.1	5.40	33.5

Note. These goodness-of-fit parameters indicate that only the model functions f₀ and f₂ are statistically significant fits to the data shown in Figures 3 and 4.

Download table as:  ASCIITypeset image

It is useful to quantify the specific regimes in which different models of convective dynamos are most applicable. To do so, consider a fully resolved convective dynamo and its associated dynamics, wherein the dissipation of the energy injected into the system is governed by the atomic values of the diffusivities. Using Braginskii plasma diffusivities (Braginskii 1965) and a set of MESA stellar models with a solar-like metallicity to obtain the density and temperature profiles (Paxton et al. 2011), one can define the average expected atomic magnetic Prandtl number in either the convective core of a massive star or the convective envelope of a lower-mass star. In what follows, ν is the kinematic viscosity, η is the magnetic diffusivity, ${\kappa }_{\mathrm{cond}}$ is the electron thermal diffusivity, and ${\kappa }_{\mathrm{rad}}$ is the radiative thermal diffusivity. These quantities are computed under the assumptions that charge neutrality holds and that the temperatures are not excessively high, so that the Coulomb logarithm is well defined. In the charge-neutral regime, the magnetic diffusivity happens to be density independent because the electron collision time scales as the inverse power of ion density and the conductivity is proportional to the electron density times the electron collision time:

$$\begin{eqnarray}&&\nu =\left[2.21\times {10}^{-15}{Z}^{-3}{\mu }^{1/2}+2.77\times {10}^{-17}\right]{\lambda }^{-1}{\rho }^{-1}{T}^{5/2},\end{eqnarray} \tag{ 53 }$$

which has units of ${\mathrm{cm}}^{2}\,\ {{\rm{s}}}^{-1}$, and where λ is the Coulomb logarithm (defined below), μ is the mean atomic mass (and it is unitless since μ is the ratio of the ion mass to the proton mass), ρ is the density in ${\rm{g}}\,\ {\mathrm{cm}}^{-3}$, T is the temperature in K, and Z is the mean atomic charge (unitless), which is defined as $Z={\sum }_{s=1}^{{N}_{S}}{Z}_{s}(T,P){n}_{s}/{n}_{i}$ and where ${n}_{i}={\sum }_{s=1}^{{N}_{S}}{n}_{s}$.

$$\begin{eqnarray}&&\eta =1.02\times {10}^{12}(2.86Z-0.90)\lambda {T}^{-3/2},\end{eqnarray} \tag{ 54 }$$

which also has units of ${\mathrm{cm}}^{2}\,\ {{\rm{s}}}^{-1}$. The thermal diffusivity arising from thermal conductivity is

$$\begin{eqnarray}&&\begin{array}{l}{\kappa }_{\mathrm{cond}}=\left\{7.41\times {10}^{-7}{Z}^{-3}{\mu }^{-1/2}+7.20\right.\\ \quad \times \ \left.{10}^{-6}\tanh \left[0.13\left(1+Z\right)\right]\right\}{\lambda }^{-1}{T}^{5/2}.\end{array}\end{eqnarray} \tag{ 55 }$$

However, it is necessary to include the radiative contribution to the thermal conductivity:

$$\begin{eqnarray}&&\kappa ={\kappa }_{\mathrm{cond}}+{\kappa }_{\mathrm{rad}},\end{eqnarray} \tag{ 56 }$$

with

$$\begin{eqnarray}&&{\kappa }_{\mathrm{rad}}=\displaystyle \frac{4{{acT}}^{3}}{3{c}_{P}{\rho }^{2}\sigma },\end{eqnarray} \tag{ 57 }$$

where a is the radiation constant, c is the speed of light, c_P is the pressure specific heat, and σ is the opacity. With those in hand, one can then compute the thermal and magnetic Prandtl numbers and the magnetic Reynolds number as necessary. The Coulomb logarithm is quite complex, as it depends on the number of included species and on the details of their Landau scattering integrals. However, it only varies between roughly 1 and 30 within the context of stars. Hence, for now, set $\lambda =15$ (which is an average value for densities between 10⁻¹⁰ and 10¹⁰ ${\rm{g}}\ {\mathrm{cm}}^{-3}$ and temperatures between 10³ and 10¹⁰ ${\rm{K}}$). Therefore, the thermal Prandtl number $\Pr$ and magnetic Prandtl number $\mathrm{Pm}$ scale as

$$\begin{eqnarray}&&\Pr =\displaystyle \frac{\nu }{({\kappa }_{\mathrm{cond}}+{\kappa }_{\mathrm{rad}})},\qquad \mathrm{Pm}=\displaystyle \frac{\nu }{\eta }.\end{eqnarray} \tag{ 58 }$$