On Convective Turnover Times and Dynamos in Low-mass Stars

We compiled a catalog of measured ratios of X-ray to bolometric luminosity as well as rotation periods from the literature (see Table 1 for a summary). These data are available alongside our MESA models on Zenodo.¹¹ The majority of this sample comprises solar-type stars with masses larger than the fully convective boundary of ∼0.35 M_⊙ retrieved from the work of N. J. Wright et al. (2011). In that study, the X-ray luminosities were computed using the ROSAT bandpass ranging from 0.1 to 2.4 keV, which we also adopt as the reference X-ray luminosity bandpass for the study in hand. When necessary, we converted X-ray luminosities of data from other studies to the same energy band using the webPIMMS tool,¹² assuming a plasma/APEC model with K and N_H = 10²⁰ cm⁻² for the conversion. If unavailable, we assumed a 20% uncertainty for L_X/L_bol, consistent with the mean uncertainties reported in other studies (e.g., B. Stelzer et al. 2016; E. González-Álvarez et al. 2019; E. Magaudda et al. 2020, 2022).

This base sample was extended by new and archival data provided by E. Magaudda et al. (2020), who homogenized various different catalogs from the literature (B. Stelzer et al. 2016; N. J. Wright & J. J. Drake 2016; N. J. Wright et al. 2018; E. González-Álvarez et al. 2019) to provide a uniform sample of X-ray activity and rotation properties of M dwarfs. We additionally included data for field stars from Y. Shan et al. (2024, CARMENES/ROSAT), K. G. Stassun & M. Kounkel (2024, TESS/eROSITA), E. Magaudda et al. (2022, TESS/eROSITA), and D. Pizzocaro et al. (2019, Kepler/XMM-Newton), as well as data for the open clusters M34, M35, Praesepe, and Hyades from P. Gondoin (2012), P. Gondoin (2013), and A. Núñez et al. (2022, ROSAT/Chandra/Swift observatory/XMM-Newton/K2), respectively.

This resulted in a total of 9344 sources in the base catalog, out of which, only 1451 sources remained in our “clean catalog” after applying various quality cuts to these data. Below, we describe each of the quality cuts, improvements, and considerations we included for our sample.

Crossmatch with Gaia DR3. We performed a simple crossmatch using a 3^″ search radius in the Tool for Operations on Catalogues and Tables (Topcat; M. B. Taylor 2005). We found that 98% (i.e., all but 189 stars) of our sample had a match. This allowed us to update the distances in our sample using the most accurate and precise measurements available, and apply the improvements and quality cuts described below.

Remove duplicates. Using the respective Gaia ID, we identified in total 190 stars that had more than one X-ray measurement in our sample. From these, we only removed 10 duplicated sources that had identical measurements of X-ray luminosity, meaning that the same measurement was being included more than once. The duplicated stars with varying values of X-ray luminosity (due to either intrinsic variability of X-ray emission or uncertainty) were kept to improve the calibration of the scatter within the rotation–activity relation.

Remove possible binaries. We subsequently removed stars with renormalized unit weight error (RUWE) > 1.4, which is a parameter provided by the Gaia survey, and allows to filter out stars for which the single-star model does not provide a good fit to the astrometric observations.¹³ Although extremely useful, the RUWE parameter does not identify tight unresolved binaries. To improve this, we added the cut ipd_frac_multi_peak > 1, which removes stars that are visually resolved double stars. In addition, we used the flag rv_amplitude_robust, which indicates the difference between the largest and the smallest radial velocity measured by Gaia. This number depends on color and magnitude, becoming larger for fainter stars. Therefore, we classified stars as binaries when rv_amplitude_robust > 9 and G < 11. In total, our binary flag removed 2216 stars.

Remove young stars. As we will discuss in more detail in Section 4.2, the convective turnover time depends among other properties on the stellar age. We therefore removed stars younger than 300 Myr for masses <0.6 M_⊙, younger than 200 Myr for masses in the range 0.6−0.8 M_⊙, and younger than 100 Myr for masses ≥0.9 M_⊙. This cut allows us to keep only stars for which the convective turnover time has stabilized, which does not always agree with the convergence into the main sequence. K. G. Stassun & M. Kounkel (2024) included the ages of the stars in their catalog, which we used to apply the cut described above, and for the rest of the sample, we used the Bayesian Analysis for Nearby Young AssociatioNs Σ tool (BANYAN Σ; J. Gagné & J. K. Faherty 2018) to estimate the probability that each star belongs to a known moving group from its position, proper motion, parallax, and radial velocity from Gaia DR3. We removed all the stars that had a probability larger than 0.9 of belonging to a group that would put them in the “young” category according to the cut described above. The cut in age removed 6739 stars, out of which 120 stars were removed using BANYAN Σ. We note that we removed most of the stars in the catalog with this age cut. In addition, we used the position in the color–magnitude diagram (CMD) to identify and remove 106 stars that belonged to the giant branch by visual inspection.

Quality cuts. We applied quality cuts to keep stars with a parallax over error > 20 (equivalent to 5% uncertainty) and also an uncertainty of smaller than 0.3 dex.

Metallicity. Although we do not have metallicity measurements for our sample, we can take into account several factors that justify the sample being mostly made of stars with solar-like metallicity. R. Kiman et al. (2019) showed that low-metallicity M dwarfs have a locus below the main sequence in the Gaia CMD. When compared to the CMD of our sample, none of the M dwarfs in our sample seem to be low metallicity. In addition, we crossmatched our clean sample with the APOGEE survey (S. R. Majewski et al. 2017; V. V. Smith et al. 2021) and found 245 stars in common (27%). We found that the values for [M/H] are close to solar (between −0.3 and 0.3 dex), which supports the conclusion we obtained from the CMD position of the stars.

Improve mass estimation. We found that the calculations of stellar masses were inconsistent among the different catalogs. In particular, there are known problems with the mass estimation from N. J. Wright et al. (2011; see W.-C. Jao et al. 2022). Therefore, we decided to reestimate the stellar masses for the complete sample by interpolating within the relations provided by M. J. Pecaut & E. E. Mamajek (2013)¹⁴ to estimate masses from the absolute G magnitude (M_G). This is a precise mass–luminosity relation that is valid for the widest mass range, in contrast to other calibrations that do not cover a wide range of masses (e.g., A. W. Mann et al. 2019; M. R. Giovinazzi & C. H. Blake 2022).

Our final clean sample comprises 1451 stars. We show a summary of the compiled sample in Table 1 and the X-ray luminosities as a function of the rotation period for the clean sample of stars color-coded by mass in the left panel of Figure 1.

Zoom In Zoom Out Reset image size

The Rossby number description of stellar activity requires the stellar rotation period in combination with the convective turnover time. The turnover time can be calculated empirically by finding the parameter that reduces the scatter in the activity–rotation relation (e.g., R. W. Noyes et al. 1984; K. Stepien 1994; N. J. Wright et al. 2011, 2018). In this section, we describe our calculation of the empirical calculation of the turnover time, τ_cE, following the methods outlined in N. J. Wright et al. (2011, 2018). This method allows us to determine τ_cE as a function of mass under the assumption that all stars follow the same two-step piecewise relationship of L_X/L_bol as a function of their Rossby number, Ro (defined as Ro = P_rot/τ_cE, where P_rot is the stellar rotation period),

where C is a normalization constant, and Ro_sat is the Rossby number at which the transition from saturated to unsaturated regimes occurs. The implication of Equation (1) is that all masses share the same saturated value of , and all masses become unsaturated at the same value of Ro_sat and, for higher values of Ro, follow the same power-law dependence of L_X/L_bol with Ro.

In order to determine the values of , Ro_sat, and β for our stellar sample, we fit Equation (1) to the observations, where . To compute Ro, we first used the convective turnover times provided by N. J. Wright et al. (2018). As we discuss in Section 3.2, the functional form of the empirical calibration from N. J. Wright et al. (2018) differs significantly from the models for high (>1 M_⊙) and low (<0.4 M_⊙) masses. Therefore, to avoid biasing our fit, we perform this first fit only to masses in the range 0.4−0.9 M_⊙. The rest of the analysis will be done using the complete sample. The ratios of L_X/L_bol as a function of Ro for the selected mass range of the sample are shown in the middle panel of Figure 1. As we cannot see evidence of supersaturation in our sample, we did not include this phenomenon in the fit; however, we cannot discard that there might be a significant effect for the fastest rotators (P_rot < 1 day, e.g., S. Randich 1998; D. J. James et al. 2000; C. Argiroffi et al. 2016).

We used the Markov Chain Monte Carlo sampler implemented in emcee (D. Foreman-Mackey et al. 2013) to determine β, , and Ro_sat. In addition, we included a parameter σ to fit the scatter of the relation and obtain a more precise estimation of uncertainties. We defined the likelihood of the model as

where is calculated using Equation (1), and is the sum of the uncertainty of the X-ray fractional luminosity squared and the characterization of the scatter of the relation squared (). To obtain the posterior, we combined the likelihood described above with flat priors for each of the parameters such that 0.01 < Ro_sat < 1, −5 < β < −1, , and 0 < σ < 5.

We report the median of the posterior as the parameter estimate, and the symmetric interval surrounding the median that contains 68.3% of the posterior distribution as the uncertainty. We find β = −1.96 ± 0.08, Ro_sat = 0.11 ± 0.01, and , and σ = 0.40 ± 0.01. The parametric fit to the data is also shown in the middle panel of Figure 1, where we show the median together with 2000 random samples from the posteriors of each parameter and the fit to the scatter of the relation. We note that our fit agrees within uncertainties with the results from N. J. Wright et al. (2018) shown in an orange dashed line in Figure 1 (, , and ).

With β and determined, we now fix the parameters β = −1.96 and . In addition, we fixed σ = 0.40. These assumptions allow us to determine the rotation period at which saturation occurs for groups of stars of approximately equal mass, independently of the convective turnover time. We divide our sample into 14 mass bins ranging from 0.1 to 1.2 M_⊙, with approximately equal numbers of stars in each bin, and fit the equation

using the grid sampling method to obtain the probability distribution of P_sat per mass bin. emcee is significantly slower when we only fit for one parameter, so we decided not to use it in this case. The mass-dependent constant can also be defined as using Equation (1), which now allows for the determination of τ_cE as a function of stellar mass. We show the individual fits to each mass bin in Figure 2. By visually inspecting each of the fits, we find that most of them agree with the data. There is a small difference between fit and data for the two highest-mass bins: 1.01 M_⊙ and 1.09 M_⊙. In these two cases, the slope seems to be less steep in the data. It is possible that the relationship between activity and rotation rate would change toward higher masses, as more massive stars have vanishingly thin CEs (disappearing perhaps after 1.3 M_⊙). So, it might be expected that the relation should change toward higher masses at some point, as the magnetic dynamo correspondingly becomes ineffective. However, we cannot confirm this difference because we do not have enough saturated stars to fully characterize the relation. We note that the difference between data and fit is not large enough to affect the discussions in the rest of the paper.

Zoom In Zoom Out Reset image size

We confirmed previous results that the rotation period at which saturation occurs clearly increases with decreasing stellar mass (e.g., N. Pizzolato et al. 2003; E. R. Newton et al. 2016; N. J. Wright et al. 2018; and references therein). As mentioned above, , which means that we calculate the convective turnover times using . This leaves a constant C to be determined, which will affect the normalization of our calculations of τ_cE, but not the relative values, as noted also in previous studies (N. Pizzolato et al. 2003; N. J. Wright et al. 2011). N. J. Wright et al. (2011, 2018) decided to choose this normalization, so that their solar value for τ_c agrees with the value in R. W. Noyes (1984), which was theoretically derived. Unfortunately, there is no precise independent estimation of the solar convective turnover time that could be used to normalize our values (W.-C. Jao et al. 2022); therefore, we need to choose between other calibrations. As we are interested in comparing our empirical values with our model values of convective turnover times (τ_cM, described in Section 2.2), and with the empirical values obtained by N. J. Wright et al. (2018), we defined C such that τ_cE agrees with our theoretical calculations for the mass bin corresponding to 0.96 M_⊙. This particular mass bin is the closest to the solar-mass bin that has the highest number of stars, which makes the fit to the X-ray versus rotation period more precise (see Figure 2). Furthermore, the values for τ_c from N. J. Wright et al. (2018) and our theoretical calculations (14.8 and 14.6 days, respectively) agree in this bin, which allows the comparison of our empirical τ_cE with both calibrations. To calculate the uncertainties of τ_cE, we did a Monte Carlo propagation of uncertainties of all parameters involved in the calculation (P_sat, β, and ). We note that, although β and were fixed to calculate P_sat, the uncertainty of the first fit was included in the calculation of P_sat given that we included the parameter σ to characterize the scatter of the relation. The masses and uncertainties for each bin were calculated as the median mass of the bin and the standard deviation, respectively. The resulting convective turnover times calculated with the empirical method are shown in Figure 3.

Zoom In Zoom Out Reset image size

In order to examine stellar structure model calculations of the convective turnover time, τ_cM, we used MESA, simulating stars in the mass range 0.1–1.3 M_⊙ with solar metallicity and a ratio of mixing length to pressure scale height α_MLT = 1.82. Our model profiles as well as full history output at 1, 5, and 14 Gyr, plus inlist and source files, are available on Zenodo.¹⁵ Many of our underlying physical assumptions are derived from the MISTv1.2 models (J. Choi et al. 2016), based on calibrations described in that text. The value of α_MLT = 1.82 can affect the convection zone size of stellar models, with larger values leading to deeper convection zones. While our choice of α_MLT follows from a calibration to solar helioseismic data (J. Choi et al. 2016), it is similar to the value of 2 chosen by R. W. Noyes (1984). There, it was found to minimize the scatter of the activity–rotation relationship and was also noted to reproduce solar values.

We summarize a number of the most salient parameter adoptions made for our models in Table 2 and refer to it for values of the parameters described below. As stated above, most of these values are unchanged and adopted from MIST (J. Choi et al. 2016), to which we refer the reader for further details. In our models, as we are using MESA version r11701; the standard amalgamation of equation of state tables is now augmented by the PTEH tables covering especially low density material (see B. Paxton et al. 2019). Similarly, we use a standard amalgamation of opacity tables.

Table 2. A Selection of Adopted Parameters in Our MESA Models

Ingredient	Adopted Parameters	References
Solar abundance scale	X⊙ = 0.7154Y⊙ = 0.2703, Z⊙ = 0.0142	(a)
Equation of state	OPAL+SCVH+MacDonald+HELM+PC+PTEH	(b, c, d, e, f, g)
Boundary conditions	M ≤ 0.5 M⊙: tau_100_tables, off table: gray and_kap	(h)
	M > 0.5 M⊙: simple_photosphere
Opacities	OPAL Type I for T ≳ 104 K; Ferguson for T ≲ 104 K	(i, j)
	OPAL Type I → Type II post-TAMS
Convection	αMLT = 1.82 (solar calibrated), ν = 1/3, y = 8	(k)
Overshoot	time-dependent, diffusive,	(l)
Semiconvective mixing	αsc = 0.1	(m)
Thermohaline mixing	αth = 666	(n, o)
Winds (mass loss)	Reimers (ηR = 0.1) for the main sequence and red giant branch	(q)
	Blöcker (ηB = 0.2) for the asymptotic giant branch	(r)

Note. Unless otherwise noted in the text, values are adopted from MIST (J. Choi et al. 2016), which may be consulted for further details.

References. (a) M. Asplund et al. (2009), (b) F. J. Rogers & A. Nayfonov (2002), (c) D. Saumon et al. (1995), (d) J. MacDonald & D. J. Mullan (2012), (e) F. X. Timmes & F. D. Swesty (2000), (f) A. Y. Potekhin & G. Chabrier (2010), (g) O. R. Pols et al. (1995), (h) B. Paxton et al. (2011), (i) C. A. Iglesias & F. J. Rogers (1993, 1996), (j) J. W. Ferguson et al. (2005), (k) L. Henyey et al. (1965), (l) F. Herwig (2000), (m) N. Langer et al. (1983), (n) R. K. Ulrich (1972), (o) R. Kippenhahn et al. (1980), (p) A. Heger et al. (2000), (q) D. Reimers (1975), and (r) T. Blöecker (1995).

Download table as:  ASCIITypeset image

Our choices for the treatment of convection derive from the calibrated parameterizations utilized in J. Choi et al. (2016). Namely, our convective overshooting follows the time-dependent, diffusive, F. Herwig (2000) exponential overshooting scheme. In this scheme, overshoot mixing is parameterized through a diffusion coefficient scaled by a free parameter f_ov. This free parameter can take on different values in the core (), shell (f_ov,sh), or envelope (f_ov,env) of the model.

As we utilize the Ledoux criterion, regions that are convectively unstable according to the Schwarzschild criterion may be stabilized by a positive composition gradient, leading to semiconvective mixing. In the opposite case of a thermally stable region destabilized by a negative composition gradient, thermohaline mixing can arise. Both thermohaline and semiconvective mixing enter as time-dependent diffusive processes, parameterized by diffusion coefficients (as in R. K. Ulrich 1972; N. Langer et al. 1983), scaled by free parameters α_th and α_sc, respectively.

Mass loss for these models follows the same scheme as used in the MIST models of J. Choi et al. (2016). Most relevant for the mass range under study are the wind schemes used to describe mass loss on the main sequence and red giant branch (RGB; D. Reimers 1975) and on the asymptotic giant branch (AGB; T. Blöecker 1995). The mass-loss rates described by D. Reimers (1975) and T. Blöecker (1995) are scaled by free parameters ranging from 0 to 1 denoted as η_R and η_B in Table 2. The adopted values were tuned by J. Choi et al. (2016; see references therein) to match observational constraints that include the initial–final mass relation, the AGB luminosity function, and asteroseismic constraints from Kepler.

Lastly, compared to J. Choi et al. (2016), we adopt a simplified set of boundary conditions. We utilize the atmosphere tables referred to as simple_photosphere for masses above 0.5 M_⊙. This option calculates the temperature using the Eddington T(τ) relation, assuming τ = 2/3, and calculates the pressure similarly utilizing an expression for hydrostatic equilibrium. As noted in J. Choi et al. (2016), the simple_photosphere tables can be a poor choice in some cases.

We calculated models that use MESA’s photosphere_tables, which are fully modeled atmosphere tables as described in B. Paxton et al. (2011) as well in this mass range; ultimately, these alternate models made no difference to our results. For lower masses, we adopt a combination of MESA’s tau_100_tables and grey_and_kap tables for density and temperature values not covered by the former (see B. Paxton et al. 2011 for more on these tables). This choice is motivated similarly to the choice made in J. Choi et al. (2016), where anchoring the boundary condition to a point deeper in the atmosphere (at τ = 100 in this case) can provide a more realistic model for low-mass stars whose outer envelopes are heavily influenced by molecules not accounted for by MESA’s standard atmosphere tables.

These are 1D stellar models evolved to an age of 1 Gyr, by which time stars in this mass range have settled onto the main sequence. On the main sequence, the convective turnover time has mostly stabilized to a single value and only changes slightly with time. Our models are nonrotating, and while rotation can affect the stellar structure, and thus convective boundaries, we find that such effects are relatively small. Our models are of the solar metallicity (Z_⊙ = 0.0142 from M. Asplund et al. 2009), and while metallicity variations can substantially alter convection zone size, we find the majority of our observed stars are near solar metallicity (see Section 2.1.1); although, we aim to explore the effect of metallicity variations in greater detail in future work.

Stars less than approximately 1.3–1.4 M_⊙ possess CEs that start off extremely thin at the higher-mass end (A. C. Beyer & R. J. White 2024), and eventually extend all the way to the core at around M ≲ 0.35 M_⊙, as determined via the “Ledoux criterion” (P. Ledoux 1947) implemented in the MESA models. At 1 Gyr, stars in this mass range are also primarily on the main sequence, when interior conditions and parameters are mostly changing only very slowly (see, e.g., Y.-C. Kim & P. Demarque 1996; N. R. Landin et al. 2010), and are representative of all but the very lowest-mass stars in young open clusters as well as the Galactic disk.

The convective turnover time is a quantity typically understood through MLT (via L. Henyey et al. 1965; see M. Joyce & J. Tayar 2023 for a review) that quantifies the timescale of convective motion. It is a local quantity, calculated at a position r within a stellar model as

where v_c(r) is the local convective velocity, and H_P(r) is the local pressure scale height. Thus, τ_c(r) is the local ratio of the pressure scale height to convective velocity at some location r in the CE, taking on well-defined values in convection zones. It is common to scale H_P by some factor, α_MLT, the mixing length parameter, which typically ranges in value from 1 to 2, depending on the model. We neglect the scaling factor of α_MLT, and calculate τ_c as in Equation (4). To calculate the Rossby number (Ro = P_rot/τ_c), one must then decide on the appropriate position r to use when calculating τ_c(r).¹⁶

For example, as noted in Section 1, R. W. Noyes (1984) used the convective turnover time calculations of P. A. Gilman (1980). The models of P. A. Gilman (1980; similar to calculations made by B. R. Durney & J. Latour 1978) calculated τ_c with H_P(r) evaluated at the bottom of the CE (r = r_BCE), and v_c(r) at one pressure scale height above r_BCE (at r = r_BCE + H_P(r_BCE); see also R. L. Gilliland 1985). In the context of stellar dynamos, these choices were made under the consideration that solar and stellar (αΩ) dynamos are thought to reside near the tachocline (near r = r_BCE).

We take Equation (4) as the definition of τ_c and test the assumption that the dynamo may lie near r_BCE in our analysis. In calculating τ_c, there is some subtlety in calculating H_P, motivating a slightly different calculation method from P. A. Gilman (1980) to calculate τ_c in our work that we describe below. Throughout the text, we use τ_cM and τ_cE to refer to our model and empirical τ_c, respectively.

Here, we outline complications of calculating H_P in fully convective models and our handling of them. Classically, H_P is calculated as H_P(r) = P(r)/g(r)ρ(r); with P(r) being the local pressure, g(r) being the local gravitational acceleration, and ρ(r) being the local density, all functions of the radius r. In fully convective stars, r_BCE is simply the center of the star. Here, g(0) → 0, causing H_P(0) to diverge under the definition above for fully convective stars. This is problematic for the classical definition of the convective turnover time that evaluates H_P(r) at r = r_BCE (P. A. Gilman 1980; R. L. Gilliland 1985).

In MESA, the divergence of H_P is handled by switching to an alternate definition near the center of the star (P. P. Eggleton 1971), when (as described in B. Paxton et al. 2011, Section 5.1). While this avoids the numerical trap of infinite scale height at the stellar center, it still represents a discontinuity in the treatment of H_P(r) below the fully convective limit. In the fully convective regime, this would then be using the P. P. Eggleton (1971) formula for and the classical expression for H_P(r) in the partially convective regime. Additionally, as v_c(r) → 0 at r = r_BCE, P. A. Gilman (1980) calculated v_c(r) at a different position from H_P(r), i.e., at r = r_BCE + H_P(r_BCE), one pressure scale height from the bottom of the CE (R. L. Gilliland 1985).

We choose to evaluate both H_P(r) and v_c(r)—and thus τ_c(r)—at the same position, hereafter called . Conceptually, is similar to where R. L. Gilliland (1985; and P. A. Gilman 1980) calculated v_c(r) in their models. We utilize

which is the bottom of the CE plus half of a (local) pressure scale height. However, rather than setting r = r_BCE in Equation (5) (as in classical works), we solve for the position r to evaluate τ_c(r) as follows.

Starting from the surface of the 1D stellar model, we advance inward cell by cell. At each ith model cell, we evaluate

and stop iterating when Equation (6) is satisfied, taking r = r_i. If we rearrange Equation (6), one may see that we are solving

where Δr_i ≡ r_i − r_BCE is an extended distance above the bottom of the convection zone, i.e., the distance from the bottom of the convection zone to the ith cell above it. Equation (7) is satisfied when this distance becomes comparable to half the local pressure scale height in the CE. In our MESA models, this locale lies near the bottom of the CE, but always above where and therefore presents a consistent treatment of H_P, regardless of convection zone depth. We then use r = r_i with Equation (4) to estimate our model convective turnover times, τ_cM.

Going forward, we operate under the framework that the empirical convective turnover time (τ_cE, described in Section 2.1.2) is comparable to the theoretical turnover time. Accordingly, we consider that the empirical value correlates to τ_c(r) evaluated at a particular position (r) within the stellar model. The implication is that this position corresponds to the mean location of the magnetic dynamo, as revealed through the activity–Rossby relationship.

The Brunt–Väisälä or buoyancy frequency is a quantity describing the frequency of oscillation for a particle displaced in a stable medium. It also serves as an indicator of convective stability in stellar atmospheres. As defined in MESA (B. Paxton et al. 2013), it may be calculated as

where B is a quantity accounting for composition gradients, as described in Section 3.3 of that text. The terms χ_ρ and χ_T represent the partial derivatives and , respectively. Hence, the Brunt–Väisälä frequency is directly related to the Ledoux criterion. When N² takes on positive values, the frequency N is real and leads to oscillatory solutions of motion (stability). However, when N² is negative, N is imaginary and leads to exponentially growing solutions of motion, i.e., instability. In a stellar model, N² is positive in radiative zones and negative in convective zones.

We plot ∣N²∣ in the convective zones of our models in Figure 5 as the dark orange dashed lines. We find that a position in the stellar model () where ∣N²∣ > 10⁻¹⁴ Hz⁻² (dark orange points/vertical lines) correlates fairly well with the position inferred by τ_cE (purple points/vertical lines). In our higher-mass models, ∣N²∣ never gets this low, and this metric is simply placed at the minimum value in the convection zone, near 10⁻¹³ Hz⁻² for early M dwarfs and 10⁻¹² Hz⁻² for models representing G-type stars. In Figure 7, this location (dark orange squares and line) tracks well with r_cE.

This may offer some intuition as to what physical processes are at play in the dynamo at r_cE. This condition suggests that the Brunt–Väisälä growth rate for exponential motion (via unstable buoyant forces in the convection zone) should rise above a threshold before magnetic flux is efficiently transported by convection. Below this threshold, the buoyant force would grow at a rate such that it is overwhelmed by the Coriolis force, and magnetic flux rises with a trajectory emerging at relatively high latitudes. The values cited above for ∣N²∣ that are found to correlate with r_cE are close to those calculated by S. D’Silva (1995). For a solar-type star, S. D’Silva (1995) analytically calculated that ∣N²∣ should be greater than 10⁻¹² Hz⁻² for magnetic flux tubes to rise exponentially in an adiabatic medium and arise at surface latitudes in agreement with sunspots. Those authors derive a condition for threshold values of ∣N²∣ that would scale with mass and could potentially provide a viable alternative approximation for defining where τ_c(r) should be calculated, but is beyond the scope of this work.

The concept of buoyantly rising flux tubes has undergone significant evolution over the last several decades, as reviewed by Y. Fan (2021). Initially, these models tended to consider magnetic flux being stored and rising from a stable overshoot layer beneath the convection zone and rising to become amplified by rotational shear at the tachocline. As also reviewed by P. Charbonneau (2020) and R. H. Cameron & M. Schüssler (2023), such interface dynamo models face a number of challenges, and currently produce results in contention with solar observations. Our results, like those of N. J. Wright et al. (2018), suggest that the stellar dynamo may be similar in partially and fully convective stars, the latter of which the tachocline is absent in. If the dynamo is situated deep in the CE, as suggested by our findings, then our results lend some evidence to the possibility of a global dynamo operating in the convection zone itself. Our results would also roughly align with findings from C. P. Bice & J. Toomre (2023) that suggest that a global dynamo may form a reservoir of magnetic flux tubes deep in the convection zone of both partially and fully convective stars, from which these tubes may rise to form active regions on the stellar surface.

Such models are in the vein of those described by, e.g., N. J. Nelson et al. (2011, 2013, 2014), L. Jouve et al. (2013), and Y. Fan & F. Fang (2014) where shear via turbulent convection may produce magnetic flux tubes, without the need of shear at the tachocline. These simulations tend to find magnetic flux concentrated near the bottom of the CE, which rises to form active regions at the stellar surface. Such simulations are capable of reproducing many observed properties of the solar magnetic field and bypass the difficulties encountered by interface dynamo models. They are however computationally challenging, requiring high spatial resolution to investigate further, as achieved by, e.g., H. Hotta et al. (2016) and H. Hotta & H. Iijima (2020). Further advancements will be necessary to confirm whether these models can accurately reproduce the solar cycle, and better understand details such as where fields are generated.

Both r_cM and begin to settle at about 20% of the stellar radius in fully convective stars, suggesting the dynamo in these stars sits at a fixed fraction of the total stellar radius. Simulations by M. K. Browning (2008) suggest that convection operates relatively weakly in the centers of fully convective stars, perhaps consequently leading to relatively weak dynamo action there as well. Simulations by C. P. Bice & J. Toomre (2023) also show that the energy available to convection drops below about 0.4 R_⋆ due to the rise in heating associated with nuclear burning toward the core. Thus, downward advected magnetic fields tend to collect near where this drop off occurs. Motivated by this, we examined a third metric based on the nuclear energy generation rate, in this case due to the proton–proton (p–p) chain reaction (the dominant nuclear reaction in low-mass main-sequence stars).

The pink shaded region in Figure 5 shows where the specific nuclear energy generation rate (_pp) falls below 50% of its maximum value. In the outer convection zone of stars with M_⋆ > 0.3 M_⊙, _pp is negligible and far from this threshold; however, in fully convective stars with M_⋆≤ 0.3 M_⊙, one may see this boundary more clearly as the core comes into view, with _pp peaking at the center of the stellar model.

In Figure 7, we plot the metric , which is the point at which the distribution rises to 50% of its peak value, as the dotted pink line. The shaded region in this case represents the region that spans where the distribution falls from 70% (closer to the center, ) to 50% (closer to the surface, ) of its peak value. We find that the metric correlates fairly well with the location of the empirical turnover times in Figure 7, roughly agreeing within errors.

We note that the nuclear energy generation metric is flat in the radial coordinate for M < 0.3 M_⊙, similar to the pressure scale height metric. Conceivably, a metric based on the bottom of the convection zone and some number of pressure scale heights ≳1, as commonly used, could provide similar results for M < 0.3 M_⊙, in agreement with empirical values.

Since the α-effect in an elementary αΩ dynamo depends on convection, it might also be expected that dynamo action is not effective where convection is weak. We examined as a function of depth the ratio of convective to total luminosity, L_c/L, for the range of masses in our study. For masses M_⋆ > 0.3 M_⊙, the empirical turnover times consistently occur where the ratio of convective to total power is L_c/L ∼ 0.35, confirming that dynamo action does not appear to be associated with zones of weak convection. However, at the lowest masses, L_c/L does not reach as low as 0.35 throughout the envelope, rendering any metric for dynamo action based on L_c/L problematic to apply.