© The Authors 2024

1. Introduction

Since the pioneering work of Skumanich (1972), followed by Barnes (2003), there has seemed to be a link between the age of a solar-type star and its rotation rate. Older stars usually rotate slower than their younger equivalent; this is the well-known ‘gyrochronology’ proposed by Barnes (2003). This in turn seems to influence the degree of magnetic activity that a given solar-type star harbours. Young stars are usually much more magnetically active than older ones like the Sun. Vidotto et al. (2014) propose that this link between stellar age and magnetic activity be called ‘magnetochronology’ (see also Mathur et al. 2023).

Recently, some authors (van Saders et al. 2016; Hall et al. 2021; Metcalfe et al. 2022) have questioned this link between rotation, magnetic activity, and age, arguing that after a certain age (about the age of the Sun for solar twins, i.e 4.5 Gyr), such a link is broken, with stellar rotation being ‘stalled’. Others have found that it still holds (Lorenzo-Oliveira et al. 2016, 2019; do Nascimento et al. 2020), with possibly only a temporary pause (Curtis et al. 2020). The difficulty comes from the large uncertainty of the age determination and the observational method used. Supporters of the ‘stalling’ scenario are usually basing their analysis on asteroseismically determined ages and rotation period calibration using Kepler data. So it is important to also consider theoretical aspects when discussing the relation between age, rotation, and magnetic activity levels in the so-called ‘Sun in time’ approach (Ayres 1997; Guinan et al. 2002; Ribas et al. 2005; Güdel 2007; Ahuir et al. 2020; Lorenzo-Oliveira et al. 2020; Johnstone et al. 2021). For instance, one can turn to numerical simulations of solar-like star dynamos to assess trends between various stellar parameters characterising mass, age, rotation rates, or magnetic states. Since the work of Durney & Latour (1977), it has been quite clear that the dynamo number D, characterising the state of the dynamo, can be directly linked to the Rossby number, Ro, of the star such that D ∝ 1/Ro², in classical α − ω dynamos. The Rossby number is a key non-dimensional number that is widely used in the study of stellar evolution and activity and can be used to bridge observations and numerical simulations (Brun & Browning 2017; Käpylä et al. 2023). Simply put, it allows us to quantify if the turbulence and internal magnetohydrodynamics in rotating stars is strongly influenced (small Rossby numbers) or not (large Rossby numbers) by rotational effects. The most classical definition of the Rossby number is the so-called ‘stellar Rossby number’, and is defined as the ratio between a measure of the convective turnover time, τ_c, at a given depth of the convective envelope and the stellar rotation period, P_rot; that is, Ro_s = P_rot/τ_c (see Landin et al. 2010; Brun et al. 2017 for further discussions on the many definitions of this number). When Ro_s is small, rotation influences the convection dynamics, tending towards the so-called magnetostrophic state, when Lorentz and Coriolis forces (and horizontal pressure gradients) balance one another out, for very small Ros values (Davidson 2014; Augustson et al. 2019). Thanks to an extensive dynamo study published in Brun et al. (2017, 2022) with the ASH code along with a similar study with the Eulag-MHD code (Strugarek et al. 2017, 2018), we now have a database of more than 30 fully 3D magnetohydrodynamic (MHD) rotating convective dynamos of solar-type stars, spanning several mass and rotation bins; hence, Rossby numbers. In the present paper, we wish to study how the properties of the surface magnetism change as we vary the Rossby number. We are helped by an equivalent systematic observational study of solar-like star magnetism performed by See et al. (2015, 2019b,a; see also Reiners et al. 2022) in the context of the Bcool consortium (Marsden et al. 2014) using Zeeman-Doppler imaging (ZDI) techniques.

In Sect. 2, we briefly present the set of 15 ASH simulations published in Brun et al. (2022) that we will analyse further in Sect. 3. We conclude in Sect. 4 by proposing a plausible magneto-rotational scenario over secular timescales for solar-type stars.

2. Magnetic properties of magnetohydrodynamic models of stellar dynamos

2.1. Brief overview of stellar dynamo simulation database

In this paper, we focus on the ASH dynamo simulations published in Brun et al. (2017, 2022). The properties of the various large-scale flows and magnetic states, as well as their associated energy and non-linear angular momentum transports in purely HD and MHD conditions, have been characterised in detail in these two publications and will not be repeated here. Instead, we wish to focus the analysis in this paper on the global properties of their magnetic field with respect to various stellar parameters.

The simulations represent solar-type stars in a mass range of 0.5–1.1 M_⊙, with rotation rates of 1/8–5 times Ω_⊙ at solar metallicity. The effective temperature range considered for the simulations lies approximately between 4030 and 6030 K, covering mostly G and K-type stars on the main sequence.

In all the simulations, the rotating convective envelope of the stars is modelled from the base of their convection zone up to about 0.97 R_*. The range of values covered by the stellar radius is between about 0.44–1.23 R_⊙ and the luminosity ranges from about 0.04–1.8 times the solar luminosity, L_⊙. In the simulations computed with the ASH code, the models also include a stable radiative layer about the same thickness as the convective envelope above, and hence possess a tachocline at the base of their convective zone (in the middle of the computational domain, approximately). The diffusivity profiles (viscous, thermal, and magnetic) are adapted (tapered) such that they maintain an almost constant Reynolds number throughout the simulations. All these simulations develop a genuine multi-scale convective dynamo, and most have been numerically integrated for several decades of physical time over many years.

From a numerical point of view, the ASH code is a semi-implicit pseudo-spectral method of solving the anelastic MHD equations in a frame rotating at Ω_* (Clune et al. 1999; Brun et al. 2004). Each simulation has a significant stratification, the level of which depends on the spectral type considered. The radial density contrast varies from about 40–80 in the convective envelope and from about 200 to 1000 when including the stable layer. The numerical resolution is Nr = 769 in radius and N_θ = 512 or 1024 in latitude, with N_ϕ = 2N_θ (a higher horizontal resolution for the low Rossby cases that develop smaller scale dynamics; see Takehiro et al. 2020 for more details on the critical convection mode excitation). We now briefly summarise some of their key magnetohydrodynamical properties, which were first analysed in Brun et al. (2022).

2.2. Rotation profiles and their link to magnetic properties: The role of the Rossby number

In Fig. 1 we display the evolutionary tracks of four stars that have stellar masses ranging from 0.5 to 1.1 M_⊙, starting from the PMS all the way to the TAMS; in other words, similar to the mass range used in the 3D MHD dynamo solutions used in this study. To do so, we plotted their evolution in a normalised Rossby number versus age diagram, using 1D stellar structure and evolution models computed with the Starevol code (Amard & Matt 2020). We superimposed the 15 ASH models on the plot, to show how our parameter space study can cover several temporal phases of evolution (see also Emeriau-Viard & Brun 2017 for a detailed specific study of the PMS phase). We used the following definition of the Rossby number:

$$\begin{array}{c}\hfill {\mathrm{Ro}}_{\mathrm{f}}=|\omega |/2{\mathrm{\Omega }}_{\ast }\end{array}$$(1)

Fig. 1. Stellar evolutionary tracks in a Rossby vs. age diagram. Shown on the figure are, 1D evolutionary track from Amard & Matt (2020) as a colour shaded envelope taking into account the initial rotational spread. The particular case of the solar track is shown in cyan colour. Symbols represent the ASH simulations and are changed according to the internal differential rotation profile achieved in the convective envelope of the simulations, either banded/quenched, solar like (fast equator slow poles) and anti-solar (slow equator fast poles). The coloured vertical bars on the left represent the magnetic states either short, decadal and no cycles.

with |ω| being the mean vorticity of the convective flows, taken at the middle of the convection zone in the ASH simulation and stellar evolution tracks, D the thickness of the convective envelope, and Ω_* the model rotation rate. This definition corresponds to the ‘fluid’ Rossby number, a measure of the influence of the Coriolis force on the non-linear advection term in the Navier-Stokes equation. The fluid and more usually stellar Rossby numbers can be related to one another, as is shown in Appendix B (see also Brun et al. 2017 for an overview of the different definitions). In Fig. 1 we normalised it to the value of the Sun’s Ro_⊙, here chosen to be 0.9. Before getting into the details of the figure, we wish to quickly recall how the Rossby number characterises the dynamics. For low values of the Rossby number, the rotational effects are dominant and force the dynamics to be aligned along the rotation axis (the so-called Taylor-Proudman constraint; Pedlosky 1987; Brun & Toomre 2002; Miesch et al. 2006; Busse 1983; Busse & Simitev 2006). This usually results in an internal cylindrical DR profile. For intermediate values, thermal effects via baroclinic torques can bend the iso-contours of Ω to be more conical at mid-latitudes (Miesch et al. 2006), as in the Sun and its helioseismically inferred angular velocity with a fast equator and slow poles. For large values of the Rossby number, the rotational effects are weaker and the local angular momentum conservation can lead to anti-solar DR profiles, with a slow equator and fast poles (Gastine et al. 2014). These different states of internal angular velocity profiles are represented by the symbols (cross, square, and circle) in Fig. 1.

We can notice several key pieces of information about our set of 15 ASH simulations that are plotted as symbols of various shapes and colours. In this stellar evolutionary diagram, we first see a clear trend of stars evolving from the bottom left towards the upper right. Indeed, as stars age they tend to slow down (these 1D stellar evolutionary tracks do not consider the possible stalling of the spin-down advocated by some authors, as is discussed in the introduction and Sects. 2.2 and 3). We note that the large range of Rossby numbers of the 3D MHD dynamo study covers a significant part of the stellar evolutionary track of solar-type stars. Some stars are in the low Rossby number regime (square symbols), whereas others are in the slow rotation regime (cross symbols). We also note that there is a continuous change in the associated DR regimes, going from banded or quenched in the early stages towards becoming anti-solar for a long-enough secular evolution. As of today, it is difficult to say if solar-type stars will become anti-solar before turning into sub-giant or red giant stars. In Noraz et al. (2022a), we conducted a systematic study of the ‘Kepler’ sample published in Santos et al. (2021) and found 22 possible candidates that would be worth observing further in order to put more constraints on slowly rotating stars. Nevertheless, the continuous transition of states is interesting in and of itself and could sustain a Sun in time scenario that describes a magneto-rotational dynamical evolution of stars like the Sun over secular time frames. Indeed, we also added in Fig. 1, the magnetic dynamo states of the simulations that were found in Brun et al. (2022). For low Rossby numbers, most models harbour short cycle period dynamo actions. Periods of the order of half a year to two years are found in the models. For intermediate, more solar-like Rossby number values, corresponding to most of the main sequence of those stars, we find decadal cycle periods as in the Sun or 18Sco (Do Nascimento et al. 2023). Finally, evolved old stars may lose their cyclic magnetic behaviour and display instead statistically stationary magnetic states, with very stable polarity in each hemisphere over secular time frames before turning into more evolved stages out of the main sequence, for which our 3D MHD study is not designed.

2.3. Magnetic butterfly diagrams versus stellar dynamo types

In order to illustrate a bit more the dynamo states achieved in the solar-type stars modelled in the study of Brun et al. (2022), we represent in Fig. 2 three typical magnetic butterfly diagrams found in this 3D MHD parameter study. The butterfly diagrams were formed by azimuthally averaging the toroidal field of the simulation at any depth in the simulation (usually either in the tachocline at the base of the convection zone or near the surface) and by stacking these latitudinal bands in time to form time–latitude contour plots. Each of these diagrams covers several decades of evolution and is strikingly different at first sight.

Fig. 2. Typical magnetic butterfly diagrams (i.e., time–latitude diagram of the axisymmetric toroidal magnetic field) achieved in the simulations published in Brun et al. (2022). On the top panel, the butterfly diagram is shown as a colour contour plot in Gauss for the rapidly rotating cases (low Rossby numbers), in the middle panel for intermediate value (Sun-like) and in the bottom panel for statistically steady solutions for slowly rotating stars with large Rossby numbers.

In the top panel, we show a representative butterfly diagram near the surface for rapidly rotating stars, those with small Rossby numbers. We clearly see the small red-blue alternating colour bands at mid-latitudes, illustrating here local polarity inversions. The short cycle period is of the order of six months in the case illustrated. The bands propagate polewards, and the dynamo wave follows the Parker-Yoshimura rule for α − Ω dynamo types (as was demonstrated in Brun et al. 2022).

In the middle panel, we show the butterfly diagram near the base of the convection zone for the case with an intermediate Rossby number. We clearly see the long cyclic behaviour, with three consecutive cycles with a typical global polarity reversal, as is seen in the Sun; in other words, the magnetic features have reversed polarity in each hemisphere and the polarity swap signs from one cycle to the next, as Hale et al. (1919) reported from the Sun in his seminal paper. This is due to the dominance of the dipolar-antisymmetric dynamo mode, although some desynchronisation between the northern and southern hemispheres can be seen, which is due to a non-negligible quadrupolar-symmetric mode that leads to a more independent hemispherical magnetic response (Gallet & Pétrélis 2009; DeRosa et al. 2012). Another key feature of the middle panel butterfly diagram is the mid-latitude equatorward branch, and a high latitude polar branch. Unlike the fast rotating dynamo cases at a low Rossby number, these long period cycle dynamos do not follow the classical Parker-Yoshimura rule. Both α and Ω effects do play a role, but the dynamo loop leading to a cyclic behaviour involves the non-linear retroaction of the large-scale toroidal field on the large-scale shear (Strugarek et al. 2017). A new cycle starts with the reversed polarity, when the Maxwell stress modifies locally the sign of the gradient, ∂Ω/∂θ. The longer cycle period comes from the time it takes for the field to alter the angular velocity shear, as this only occurs above a certain field strength. Indeed, the Ω effect is a linear field stretching mechanism that takes, in the specific case illustrated in Fig. 2, about ten years to act. We also find that this dynamo operates much deeper, straddling the base of the convective envelope, where significant energy transfers (up to several % of the solar luminosity) allow global polarity reversals of the large-scale magnetic field in a prey-predator-type mechanism, generating torsional oscillations (see Brun et al. 2022 for their analysis).

In the bottom panel, we also display the butterfly diagram near the base of the convection zone for a typical large Rossby number case (Ro_f > 1). This type of dynamo possesses an anti-solar DR. As is shown in Noraz et al. (2022b; see also Karak et al. 2020), such reverse angular velocity profiles (with respect to the Sun) often yield stationary dynamos that do not show clear and systematic polarity reversals. Temporal variability is still present, with the large-scale toroidal magnetic wreaths exhibiting amplitude variations, sometimes not going all the way around the 360° longitudes of the star (see Nelson et al. 2013). Here, we do not find any sign of polarity reversal in the non-axisymmetric components of the field that were found by Viviani et al. (2018).

Our simulations, despite their sophistication, may not faithfully reproduce every nuance of real stellar surface dynamics. Meridional circulation and torsional oscillations amplitudes observed in the Sun, for instance, may not be precisely mirrored in our simulations (see Hotta et al. 2023 for a review). Nevertheless, we do not aim to reproduce the precise solar case or every detail of solar surface magnetism in this parametric study, but instead wish to unveil broader trends as a function of the Rossby number over evolutionary timescales. We believe that the trends that will be discussed in Sects. 3 and 4 are indicative of genuine energy and force balances occurring within solar-type star convective envelopes (as is discussed in Davidson 2013; Aubert et al. 2017; Augustson et al. 2019, and references therein).

In summary, one can thus imagine that, as a solar-type star ages, it will respectively go through these three magnetic and rotation states. In order to further verify if such a stellar magneto-rotational scenario is plausible, we wish to compare other magnetic proxies with recent observations of stellar magnetism. To this end, we now turn to search in our dynamo database for various trends with respect to some global stellar parameters by splitting the field into various components (toroidal, poloidal, axisymmetric, dipolar, multipolar, etc.). In the following section, two simulations out of the set of 15 simulations presented above (namely M07R3m and M11R5m, see Brun et al. 2022 for naming nomenclature) will not be considered because of a gap in the data needed at the time of the present study (spatial and temporal gaps, respectively). The top of the numerical domain is r_top = 0.95 R_* for M = 5 M_⊙ models and r_top = 0.97 R_* for all the others. Values referred to as near the surface in the rest of the paper are evaluated at r = 0.9997 r_top for M = 5 M_⊙ models, r = 0.9993 r_top for M11R3m, and r = 0.9998 r_top for all the others.

3. Magnetic dynamo trends with stellar parameters

Having recalled the main broad properties of the set of dynamo solutions considered in the present paper (see Brun et al. 2022 for more details), we now wish to look systematically at various trends regarding their magnetic properties with respect to key stellar parameters (such as the Rossby number, stellar mass, or field geometry). In doing so, we intend to assess how well the set of stellar dynamo simulations can further confirm our Sun in time scenario, by directly comparing our results to those published in the observational studies of See et al. (2015, 2019b,a). To that end, we use similar layouts for the figures to ease the direct comparison and plot observational scaling laws (fits) when available in the publications. We account for the difference in the Rossby number definition used in observational studies (stellar) and the present paper (fluid) when showing observational trends as a function of Ro_f ≃ Ro_s/2.26 (see Appendix B for more details).

3.1. Poloidal versus toroidal magnetic field properties

In Fig. 3 we display the magnetic ratio, defined here as the relation between the toroidal and the poloidal field amplitude squared. We show this relation near the top of the dynamo solution (left panel) and at the middle of the convective envelope (right panel) for each stellar spectral type. The definitions of the poloidal, B_pol, and toroidal magnetic field, B_tor, can be found in Appendix A. To ease direct comparisons to observations, we truncated the spherical harmonics decomposition, conserving only ℓ ≤ 5 in the computation of both quantities (see Appendix A, and also Vidotto et al. 2016; See et al. 2019b).

Fig. 3. Left: magnetic ratio between the squared toroidal, Btor, and poloidal, Bpol, components near the surface. Two observational fits (dashed black for M > 0.5 M⊙ and dotted lines for M < 0.5 M⊙) from See et al. (2015) are indicated, as well as the fit of the simulations (purple line) and their respective error bars. Right: same ratio within the convective envelope. The coloured symbols have the following meanings: red, yellow, green, and blue represent stars with 0.5, 0.7, 0.9, and 1.1 M⊙, respectively, and triangle, diamond, square, and filled circle shapes the rotation rates for, respectively, slow (the exact value depends on mass see Brun et al. 2022), one, three, and five times the solar rotation rate, Ω⊙. The error bars indicate the minimum and maximum values reached by the model during the integrated time (one or several cycle periods) used to compute the mean value indicated by the symbol.

Looking first at the result near the surface, we see that the amplitude of the toroidal field is smaller in the set of simulations. We note that the simulation trend, $B_{\text{tor}}^2\propto {(B_{\text{pol}}^2)}^n$ (purple fit), lies in between the two values proposed by See et al. (2015), and hence captures the mean trend found in the observational studies. The lower amplitude may be explained by the potential field surface magnetic boundary condition of the simulation, enforcing a zero longitudinally averaged magnetic toroidal field. Even though we probed that interdependence, a few mesh points below the top of the numerical domain (so not quite where the zero value is enforced numerically), this is likely to still have an influence on the ratio between the two field geometries. In order to quantify this effect, we now show in the right panel of Fig. 3 this relation in the middle of the convection zone, where this influence from the boundary condition does not hold any more. We indeed note that the simulations better match the data in terms of the amplitude reaching a value of $B_{\text{tor}}^2$ around 10⁴G² and above. The simulation trend matches the higher value of the exponent, n (i.e., 1.18 vs. 1.25); in other words, the simulations tend to qualitatively agree with the observations, which is pretty encouraging. If we further believe the scenario that starspots are created by the surface emergence of flux ropes coming from deeper in the convection zone (as a whole entity or subparts of it), then the observed relation may well represent the toroidal geometries of the deeper interior, as is seen here (see also Finley et al. 2024 for further investigations into the M11R3m model of this set). Since in the simulation we do not have yet the formation of large compact magnetic features due to the lack of local resolution or near-surface dynamics, and the choice of potential-field magnetic boundary conditions, the near-surface ratio in the left panel may be biased in the simulation to have much lower near-surface toroidal fields. Hence, assessing this ratio in the middle of the convection zone of the simulations is less affected by the potential field boundary conditions than at the surface, which is confirmed by comparing the two panels of Fig. 5. Work is in progress to add a realistic atmosphere on top of current global dynamos to have much improved surface magnetic field boundary conditions (Warnecke et al. 2016; Perri et al. 2021; Delorme et al. 2022; Kaneko et al. 2022).

3.2. Trends with the Rossby number

We now turn to considering various trends of the magnetic field and its components with the Rossby number, Ro_f.

3.2.1. Poloidal and toroidal decomposition

In Fig. 4 we show how the near-surface poloidal magnetic field squared amplitude depends on the Rossby number. At low Rossby numbers (Ro_f < 1) where the observational data are concentrated, the agreement with the observations is quantitatively good regarding the tendency. Another interesting property can be seen in Fig. 4 for large Rossby number values. We see that the trend is opposite in sign, with $B_{\text{pol}}^2$ now increasing with Ro_f rather than decreasing as for the more rapidly rotating (low Ro_f) dynamo cases. This is due to a sharp transition in the DR in the model, going from solar to anti-solar dynamo (see Matt et al. 2011; Gastine et al. 2014; Brun & Browning 2017; Brun et al. 2022). This is a very interesting property that needs to be studied further, as this possible V-shape trend could explain a weaker temporary wind braking, due to a minimum in the poloidal field strength as a function of rotation (Rossby number). Brandenburg & Giampapa (2018) and Lehmann et al. (2023) seem also to find a reverse trend for the magnetic flux amplitude of slowly rotating stars in their observational study. Clearly, studying the high Rossby number states is becoming very timely, and observational investigations already started (Noraz et al. 2022a; Donati et al. 2023; Cristofari et al. 2023).

Fig. 4. Poloidal field component squared amplitude as a function of the Rossby number, Rof. At low Rof, the agreement between the simulations (purple line) and observation (dashed black line) is quantitatively good, with both showing a similar decreasing trend. At a higher Rossby number, Rof > 1 (yellow range), not covered by the observational database, an inverse trend is indicated by the simulations, and hence suggests a minimum in poloidal field strength near Rof ∼ 1, which could have interesting consequences for stellar spin-down via wind braking. An indicative trend proportional to ${\mathrm{Ro}}_{\mathrm{f}}^{10}$ is shown using a dotted grey line. We note the V shape that the two trends (purple and dotted lines) form, with the minimum being near Rof ∼ 1. The coloured symbols have the same meanings as in Fig. 3.

Turning now to the toroidal component, we show in Fig. 5 the dependency of $B_{\text{tor}}^2$ on Ro_f near the surface (left panel) and deeper inside the convection zone (right panel). Near the surface, as for Fig. 2 (left panel), the simulations only broadly match the observed properties. The field values are a bit too low and have only in the large the same declining trend for low Ro_f. The situation improves significantly when forming the same figure deeper down in the simulation and focusing on the low Ro_f part of the plot. We see that the simulation trend (purple fit) agrees better with the observations. The field amplitude is a bit too high, indicating that the toroidal field at the stellar surface is likely weaker than deep in the stellar dynamo, but not as weak as the potential field boundary condition imposes (see left panel). We also notice the inverse trend for high Ro_f, which we will need to look into in the near future (Noraz et al. 2022a). Very little observations are yet reporting toroidal field measurements for high Rossby number stars, except for a few M dwarfs in Donati et al. (2023) and Lehmann et al. (2023), not quite directly comparable with our study, which is more focused on G and K-type dwarfs.

Fig. 5. Left: toroidal field component squared amplitude as a function of the Rossby number, Rof, probed near the surface. At low Rof, the agreement between the simulations (purple line) and observational (dashed black line) fits is qualitatively good, with both simulations and observations showing a decreasing trend, although the field amplitude is too low near the surface (left panel). This is likely due to our choice of top magnetic boundary conditions. Right: same quantity probed in the middle of the convection zone, where the agreement is quantitatively good, with both fits being close in terms of the power law index. A possible inverse trend is indicated by the simulations at higher Rossby numbers, Rof > 1 (not covered by the observational database), and hence suggests a minimum in the toroidal field strength near Rof ∼ 1, which could have interesting consequences for stellar spin-down via wind braking.

Indeed, the errors from the simulations fits are not as small as in the observations due to the moderate number of dynamo models compared to the number of observed stars, but overall the fits agree. It should be noted that running these 15 3D MHD global convective dynamo simulations over many decades of physical time is already a challenge that took several years on massively parallel supercomputers, and that nothing in particular was done or tuned in the simulations to get this comprehensive match with the observations. This is reassuring and reinforces our confidence in the set of simulations published in Brun et al. (2022) and further discussed here.

3.2.2. Multipolar decomposition

It is interesting to further study the properties of the magnetic field of the dynamo simulations by considering the behaviour of single low ℓ spherical harmonics degree magnetic field components such as the dipole, quadrupole, and octupole (ℓ = 1, 2, 3), which can be observed in most ZDI studies of magnetic stars (Petit et al. 2008; Marsden et al. 2014; See et al. 2015, 2019b). We do so in Fig. 6, using the definition listed in Appendix A.

Fig. 6. Multipolar decomposition of the magnetic field near the surface (with m components summed quadratically), showing the dipole, ℓ = 1 (left), quadrupole, ℓ = 2 (middle), and octupole, ℓ = 3 (right), as a function of the Rossby number, Rof. At low Rof, the agreement between the simulations (purple line) and observational (dashed black line) fits is quantitatively good for all low-degree dynamo modes, with all modes showing a decreasing trend and power laws comparable to observations. The possible inverse trend for high Rossby numbers is also found in all three multipoles. Symbols have the same meanings as in Fig. 5.

The dipole is the dominant magnetic ingredient for efficient wind braking (then the quadrupole; Réville et al. 2015; Finley & Matt 2017), and thus being able to predict its amplitude is key when trying to understand the magneto-rotational evolution of solar-type stars. Looking at the leftmost panel of Fig. 6, we report a good quantitative agreement for the low Rossby number values in terms of trend and amplitude. Again, for large values we see that the dipolar field increases for slowly rotating stars. This confirms that such a V-shaped dip in the trend could play as a minimum in wind braking efficiency at this Rossby number transition, explaining a possible stalling or weakening of solar-type star spin-downs at an intermediate age (van Saders et al. 2016; Curtis et al. 2020).

Turning to the middle and right panels of Fig. 6, we see that both the quadrupolar and octupolar magnetic field components are also in good quantitative agreement with the observational trends. Both purple fit indexes match the observational trends of See et al. (2019a) at low Rossby numbers. They show an inverse trend for large Ro_f values too, reinforcing the dipolar trend discussed above and explaining why it is also clearly seen in the poloidal component in Fig. 4.

3.2.3. Large-scale versus total magnetic field relationships

We now wish to consider how the total magnetic field, B_I, and the large-scale field, B_V, obtained by different techniques based on the Zeeman effect, behave with respect to one another. We already showed in Brun et al. (2022) that filtering the surface magnetic field of the simulations was necessary for a meaningful direct comparison with observations. Indeed, standard equipartition dynamo scaling (Davidson 2014; Augustson et al. 2019), considering the bulk magnetic field, differs from the one derived using the large-scale surface magnetic field. Comparing the observed B_V obtained with ZDI techniques with respect to B_I obtained by Zeeman spectroscopy can help us disentangle the contribution of the large-scale and smaller-scale magnetic fields (down to observational limit of distant stars) to the overall dynamo mechanism occurring inside solar-type stars. We defined our proxy Stokes V magnetic field, B_V, as the filtered low ℓ degree magnetic field (up to ℓ_max = 5, see Appendix A and See et al. 2019b). We defined our proxy Zeeman total magnetic field, B_I, as the normalised near-surface integrated magnetic field, keeping all degrees ℓ (see Appendix A). The distinction between these two magnetic field definitions helps to characterise the relative sensitivities of the large-scale field and the total field to stellar parameter changes. In Fig. 7 we present the dependency of B_V on Ro_f (left panel), B_V vs. B_I (middle panel), and B_V/B_I vs. Ro_f in the rightmost panel.

Fig. 7. Various trends for our theoretical proxies of both the Stokes V magnetic field measurement, BV, and the total field measurement via Zeeman spectroscopy, BI. Left: BV dependency on the Rossby number is well captured, with both fits (purple for the simulations and dashed black lines for the observations) agreeing quantitatively well in terms of the power law index. Middle: magnetic relation between BV and BI showing, as in Fig. 3, an expected positive trend. Right: BV/BI ratio as a function of Rof, with a fit (purple line) indicating a weak dependency. Symbols have the same meanings as in Fig. 8.

In the left panel of Fig. 7 we also report a good quantitative agreement between the simulation and observational fits. We find a decreasing trend for the low Rossby number region of the plot, which matches the observational power law index of See et al. (2019b) with Ro_f. We also note the reverse trend for high Ro_f, similar to the ones shown for low-degree components in Fig. 6, as could be anticipated, given their close relationship. In the middle panel, we plot the total field B_I as a function of B_V, in a similar way to what we did in Fig. 3 for the poloidal and toroidal components. The trend is clear; both fields are positively correlated and not quite linearly related. We also note that B_V is larger and B_I smaller than the typical amplitude expected from observations, both near the surface and deeper in the convection zone (the latter is not shown here). This is likely to be a limitation of the large eddy simulation (LES) approach adopted in these simulations, whereby the dynamics of the smallest scales is not resolved. Nevertheless, the trends we find are robust, implying that the magnetic field energy distribution (spectrum) between large- and small-scale fields is expected to evolve as the star ages. In that context, we finally show in the rightmost panel how the ratio B_V/B_I behaves with respect to the Rossby number. This is to verify if, as the star rotate slower and slower, there is a tendency to build less large-scale magnetic fields. There is a weak trend, stating that indeed the large fields tend to diminish slightly faster than the total field, implying a busier and smaller-scale magnetic field near stellar surfaces. Such a broad tendency can be seen in the observational data of See et al. (2019b, see their Fig. 2) for M > 0.5 M_⊙; however, the number of stars considered is also small, and any stronger conclusion will need further investigations on both theoretical and observational sides. Again, when Ro_f goes over 1 and the DR of the star flips direction, now harbouring a slow equator and fast poles, the situation reverses and stronger large-scale magnetic fields are generated by the dynamo. Hence, there is still a peculiar region both in the field geometry and the amplitude near Ro_f ∼ 1, which we will investigate in the near future both observationally and theoretically.

In Fig. 8, we display how the large-scale magnetic field, B_V, scales as a function of the stellar mass. In the study, we cover four mass bins that can help us search for possible trends. We see that indeed there is a decreasing field amplitude with stellar mass, with the fit indicating a high negative exponent, −2, but with a relatively large error bar due to some rotational spreading. This result could seem counterintuitive, as more massive stars are more luminous (recall that for solar-type stars $L_{*}\propto M_{*}^4$; Hansen et al. 1995). However, this effect is compensated for by the fast narrowing of the convective envelope with stellar mass, resulting in a much lower averaged density in the convective envelope, and hence a lower kinetic energy reservoir (Brun et al. 2017, 2022). Observations by Johns-Krull & Valenti (2000) also found that the equilibrium field, B_eq, in M dwarfs possesses a higher amplitude when compared to that of G or F stars. They used the definition of the equilibrium field as satisfying $P_{\mathrm{mag}}\sim B_{\mathrm{eq}}^2=P_{\mathrm{gaz}}\sim \overline{\rho }{T}_{\mathrm{eff}}$ at the surface (Brun et al. 2015) and they found a relatively good agreement with their measurements for the same reason – that T_eff varies less than the mean density, $\overline{\rho }$, when going from stellar spectral M to F.

Fig. 8. Theoretical proxy for Stokes V magnetic field, BV, vs. stellar mass in the simulation’s dataset. Symbols have the same meanings as in Fig. 3. We note a rather strong decreasing trend with stellar mass (purple fit exponent around −2), the rotation rate leading to some spread, and the relatively large error bar on the fit.

3.2.4. Axi- versus non-axisymmetric magnetic field trends

In Fig. 9 we show how the axisymmetric (m = 0; left panel) and non-axisymmetric (m ≠ = 0; right panel) magnetic field components behave with respect to the Rossby number. We wish to see if the field tends to be less regular and axisymmetric under some conditions.

Fig. 9. Axisymmetric (m = 0, left) and non-axisymmetric (m ≠ 0, right) magnetic field decomposition variations with the Rossby number. No clear differences are seen between the two field components; the purple fit exponents in the low Rossby number regions being close and negative, with at best a tendency for the field to be more axisymmetric for larger Rossby number values. An inverse trend is also seen for large Rossby numbers, as was expected from the previous analysis. Symbols have the same meanings as in Fig. 8.

In the low Rossby number range, we note that the field geometry tends to keep its axi- versus non-axisymmetric nature. There is a slightly larger exponent of the simulation fit for the non-axisymmetric magnetic field component, although much attention has to be paid to the relevance of such a power law fit in the case of the axisymmetric field near the surface for low Ro_f. This tends to indicate that the rotation near Ro_f ∼ 1 favours relatively more regular axisymmetric magnetic fields in the simulations. This trend is more pronounced for larger values of Ro_f. We clearly see that axisymmetric field is particularly enhanced for the faster rotating cases (5Ω_⊙, circle shapes). We also note that the triangle symbols, which again present the reverse trends, are higher in amplitude in the axisymmetric part, indicating a more symmetric magnetic field. This can easily be understood by returning to the bottom panel of Fig. 2, where we display a butterfly diagram for a typical slowly rotating case. We see that these types of models, possessing an anti-solar DR, usually develop large-scale, statistically steady magnetic wreaths in both hemispheres. Such wreaths are intricate and intertwined large-scale magnetic field ribbons that can be stable over relatively long periods of time (Brown et al. 2010). They are dominantly axisymmetric, even though they are known to have their own complex dynamics as the degree of turbulence of the simulation is increased (Brown et al. 2008; Nelson et al. 2013).

4. Discussion and conclusion

We have discussed how well a set of 15 3D MHD global dynamo solutions compare to stellar magnetism observations of solar-type stars. Overall, we find a very good quantitative agreement between our set of 15 simulations and the observational studies from the Bcool consortium (Marsden et al. 2014) published in See et al. (2015, 2019b,a). In this ensemble study, several interesting features have been found. For instance, we find that for low Rossby number values, Ro_f < 1, the various trends in the range Ro_f = [0.05, 1] imply a decreasing magnetic field amplitude with a relatively steep slope. The poloidal, toroidal, or multipolar decompositions all follow relatively clear decreasing trends with Ro_f, confirming that young, fast-rotating stars tend to have a larger field amplitude than their older counterparts. We also note that the near-surface toroidal magnetic field amplitude is too low in the simulations. We believe this is due to our choice of potential field magnetic boundary condition, which sets the toroidal field to zero at the top of the numerical domain. When comparing the toroidal magnetic field strength deeper down in the simulations, we recover a better agreement with the observations. This somewhat confirms that there is a clear link between the dynamo mechanism operating deep in the convective envelope and the subsequent emergence of magnetic fields on the surfaces of the observed solar-type stars.

An interesting result to keep in mind is a possible reversal of the magnetic trend when reaching higher Rossby numbers, Ro_f > 1, pointing towards an enhancement of the stellar magnetism as the rotational influence decreases. This change of trend suggests a minimum of the large-scale field amplitude around Ro_f ∼ 1, forming a V shape, which could have interesting consequences for stellar spin-down via wind braking. If such a behaviour is indeed occurring in stellar dynamos, then stars reaching this rotational state could potentially be in a minimum state of wind braking (due to the local weakening of the large-scale magnetic field (dipole) in the parameter space). This in turn could explain some observational claims advocating for a stalling scenario of stellar spin-down (van Saders et al. 2016; Curtis et al. 2020; Hall et al. 2021), deviating from classical gyrochronology rotational laws (Skumanich 1972; Barnes 2003). Some hints that a reverse trend of the magnetic flux for slow rotators may exist have been put forward in Brandenburg & Giampapa (2018). While the existence of a minimum in the field amplitude near Ro_f ∼ 1 is clear in our study, this does not lead to a full stop of the stellar spin-down process, but rather to a significant slowing down of the spin-down process that could then be revived once the star has slowly but surely crossed the Ro_f > 1 transition. Such a deceleration of the rotational evolution could be further amplified if a decrease in the mass-loss rate is also happening (Metcalfe et al. 2022). However, observational constraints on the high Rossby range, Ro_f > 1, are still limited as the sensitivity of current observational techniques decreases drastically for slow rotations (Donati et al. 2006; Benomar et al. 2018). The possible disappearance of starspots due to a change in the dynamo nature would further make observational characterisations of high Rossby targets difficult with photometric techniques. In that context, more study of stellar magnetism in stars near the Ro_f = 1 transition must be undertaken to confirm that a transition exists. In Noraz et al. (2022a), we propose a list of slowly rotating (probably anti-solar) stellar candidates, and we hope to be able to study their magnetic properties in the near future and confront it with our stellar magneto-rotational scenario. In that respect, the first large-scale magnetic field quantifications for high Rossby M-dwarfs were recently reported (Donati et al. 2023; Lehmann et al. 2023), and we expect to see them soon for the solar-like G-K stars modelled in the present paper.

While very encouraging, this study could be improved in several ways. For instance, the Reynolds number of the dynamo simulations presented in this study have low to intermediate values (see Brun et al. 2022, Table 2), and as such are only numerical experiments to try to understand the complex and highly non-linear nature of stellar dynamos. Each of the individual simulations could possibly be an approximate realisation of the real single star it is supposed to represent. Indeed, global simulations of rotating convection still struggle to perfectly reproduce the solar case. In particular, there is currently a mismatch between global convection simulations and helioseismic observations regarding the power contained in giant convection cells, known as the ‘convective conundrum’ (O’Mara et al. 2016; Hotta et al. 2023). This results in a slightly too large effective Rossby number in global convection simulations of the solar rotation rate. In this context, relative comparisons between models with different Rossby numbers can be done, but the absolute positioning of a given solar-type star should be considered with care. However, the comprehensive agreement found between Brun et al. (2022) and Strugarek et al. (2017), using intrinsically different numerical methods, makes us confident about the relevance of these dynamo solutions in exploring and discussing the physical nature of solar-like cyclic activity, and the robustness of the overall trends found in our study. It is now important to reiterate that our primary focus is identifying and elucidating overarching trends. We intentionally avoided fine-tuning our simulations so as to make them agree with every observational detail. We thus find it particularly encouraging that such overall trend agreements arise without ad hoc adjustments. The existence in the simulation data of a magnetic upsurge in the high Rossby number regime then becomes a promising avenue for future research. For instance, the possibility of a significant change in the dynamo process in this high Rossby number regime will need further investigation (Noraz et al. 2022a; Donati et al. 2023; Cristofari et al. 2023; Lehmann et al. 2023).

This relatively successful ab initio approach is encouraging for future studies that will introduce higher degrees of turbulence (Reynolds numbers) and more realistic surface boundary conditions (Perri et al. 2021; Delorme et al. 2022; Kaneko et al. 2022). To summarise, this study suggests that a coherent magneto-rotational scenario for solar-type stars over secular evolution time, as is summarised in Fig. 1 and in Sects. 2.2 and 3, is plausible. That is, the young stars start rotate relatively fast and are very active with short period magnetic cycles. Then, as they age and reach intermediate Rossby number values, their cycle period increases, reaching decade-long time spans due to a change in the dynamo operation. Near the Ro_f ∼ 1 transition, the stars may undergo a weakening of their wind braking, with what could appear to be a stall. As the stars age and continue to spin down more slowly, they may reach high Rossby number values and reverse their angular velocity profile. This new state of internal rotation modifies the stellar dynamo once more, leading possibly to a loss of its cyclic behaviour and the building of stronger large-scale magnetic fields, resulting in a revival of stellar spin-down.

Acknowledgments

All authors are thankful to L. Amard for providing the 1D stellar evolution models used to make Fig. 1, P. Petit for useful discussions regarding the Bcool consortium observations of stellar magnetism, V. See, T.S. Metcalfe, and the anonymous referee for his/her review. We acknowledge GENCI project 1623 for access to supercomputing infrastructures where most of these 15 simulations were computed. We acknowledge funding support by ERC Whole Sun #810218, INSU/PNST grant and Solar Orbiter and Plato CNES financial supports. This work benefited from discussions within international teams “The Solar and Stellar Wind Connection: Heating processes and angular momentum loss” and “Solar and Stellar Dynamos: a New Era” supported by the International Space Science Institute (ISSI). A.S.B and A.S. are thankful to S.H. Saar for fruitful discussions.

References

Appendix A: Magnetic field decomposition definitions

In order to compare our simulations with the published results and figures of See et al. (2015, 2019b,a), we used the formalism of (Donati et al. 2006; Vidotto et al. 2016; Folsom et al. 2018) (see also Rieutord 1987). In this formalism, we made use of the vectorial spherical harmonics basis defined as

$$\begin{array}{c}\hfill \{\begin{array}{cc}{\mathbf{R}}_l^m\hfill & =Y_l^m{\mathbf{e}}_r\\ {\mathbf{S}}_l^m\hfill & ={\mathbf{\nabla }}_{\perp }{Y}_l^m={\partial }_{\theta }{Y}_l^m{\mathbf{e}}_{\theta }+\frac{1}{sin\theta }{\partial }_{\phi }{Y}_l^m{\mathbf{e}}_{\phi }\\ {\mathbf{T}}_l^m\hfill & ={\mathbf{\nabla }}_{\perp }\times {\mathbf{R}}_l^m=\frac{1}{sin\theta }{\partial }_{\phi }{Y}_l^m{\mathbf{e}}_{\theta }-{\partial }_{\theta }{Y}_l^m{\mathbf{e}}_{\phi }\end{array},\end{array}$$(A.1)

where $Y_l^m$ are the orthonormalised classical spherical harmonics (defined as $\int Y_{l_1}^{m_1}{\left(Y_{l_2}^{m_2}\right)}^{*}\text{d}\mathrm{\Omega }={\delta }_{l,l_1}{\delta }_{m_1,m_2}$). The vector spherical harmonics basis A.1 can be used to decompose any vector, B, such that

$$\begin{array}{c}\hfill \mathbf{B}={\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell \le {\ell }_{\mathrm{cut}}\\ -\ell \le m\le \ell \end{array}}}{\mathcal{A}}_{l,m}{\mathbf{R}}_l^m+{\mathcal{B}}_{l,m}{\mathbf{S}}_l^m+{\mathcal{C}}_{l,m}{\mathbf{T}}_l^m.\end{array}$$(A.2)

Based on this decomposition and for a direct comparison with the observational results, we used the following definitions of the various magnetic field components in the present paper,

$$\begin{array}{c}\hfill B_I^2=B_{\mathrm{rms}}^2=\frac{1}{4\pi }{\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell \le {\ell }_{\max }\\ -\ell \le m\le \ell \end{array}}}|{\mathcal{A}}_{\ell ,m}{|}^2+\ell (\ell +1)(|{\mathcal{B}}_{\ell ,m}{|}^2+|{\mathcal{C}}_{\ell ,m}{|}^2),\end{array}$$(A.3)

$$\begin{array}{c}\hfill B_V^2=\frac{1}{4\pi }{\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell \le {\ell }_{\mathrm{cut}}\\ -\ell \le m\le \ell \end{array}}}|{\mathcal{A}}_{\ell ,m}{|}^2+\ell (\ell +1)(|{\mathcal{B}}_{\ell ,m}{|}^2+|{\mathcal{C}}_{\ell ,m}{|}^2),\end{array}$$(A.4)

$$\begin{array}{c}\hfill B_{\mathrm{pol}}^2=\frac{1}{4\pi }{\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell \le {\ell }_{\mathrm{cut}}\\ -\ell \le m\le \ell \end{array}}}|{\mathcal{A}}_{\ell ,m}{|}^2+\ell (\ell +1){|{\mathcal{B}}_{\ell ,m}|}^2,\end{array}$$(A.5)

$$\begin{array}{c}\hfill B_{\{\mathrm{dip};\mathrm{quad};\mathrm{oct}\}}^2=\frac{1}{4\pi }{\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell =\{1;2;3\}\\ -\ell \le m\le \ell \end{array}}}|{\mathcal{A}}_{\ell ,m}{|}^2+\ell (\ell +1){|{\mathcal{B}}_{\ell ,m}|}^2,\end{array}$$(A.6)

$$\begin{array}{c}\hfill B_{\mathrm{tor}}^2=\frac{1}{4\pi }{\displaystyle \sum _{\begin{array}{c}\ell ,m\\ \ell \le {\ell }_{\mathrm{cut}}\\ -\ell \le m\le \ell \end{array}}}\ell (\ell +1){|{\mathcal{C}}_{\ell ,m}|}^2,\end{array}$$(A.7)

where ℓ_cut = 5 and ℓ_max = 2N_θ/3. The numerical module, SHTns (Schaeffer 2013), was used here to compute the different decompositions.

Appendix B: Relation between the stellar and the fluid Rossby number

In the present paper, we decided to adopt the ‘fluid’ definition of the Rossby number, Ro_f, as it is the direct quantification of the advection term over the Coriolis term in the momentum equation of fluid dynamics, which we used in a previous study Brun et al. (2022) to characterise the different rotational and magnetic states of the set of simulations presented in Sect. 2. However, this number is not directly accessible from observations (see Noraz et al. 2022a), and the observational trends we refer to in the present paper were computed with the ‘stellar’ definition of the Rossby number, Ro_s (See et al. 2015, 2019b,a). In order to compare our simulations fairly to the observations, we illustrate in Fig. B.1 the relationship between both definitions in our set of models.

Fig. B.1. Comparison of the stellar Rossby number, Ros, as a function of the fluid one, Rof, in the different models (orange dots) from Brun et al. (2022). A linear regression is proposed with the solid orange line (Ro${}_{\mathrm{s}}=(2.26\pm 0.05)R{o}_{\mathrm{f}}^{1.05\pm 0.04}$), along with the direct proportionality, Ros = Rof, illustrated with the dashed black line.

The fluid Rossby number, Ro_f, was computed following Eq. 1 and taken in the middle of the convection zone. The stellar Rossby number, ${\text{Ro}}_{\text{s}}=P_{\text{rot}}/{\tau }_c^{\text{CS}}$, was computed similarly to See et al. (2015, 2019b), considering the empirical expression of the convective turnover time derived by Cranmer & Saar (2011), which is ${\tau }_c^{\mathrm{CS}}=314.24exp[-\frac{T_{\mathrm{eff}}}{1952.5\mathrm{K}}-{(\frac{T_{\mathrm{eff}}}{6250\mathrm{K}})}^{18}]+0.002$ days. We see that the relationship between both definitions is close to being linear, with a proportionality factor, Ro_s ≃ 2.26Ro_f. An obvious second factor comes from their respective definitions, making the fluid Rossby number smaller, the remainder coming from small numerical differences. We use this calibration factor when plotting the observational trends as a function of the fluid Rossby number, Ro_f, in Figs. 4, 5, 6 and 7.

All Figures

Fig. 1. Stellar evolutionary tracks in a Rossby vs. age diagram. Shown on the figure are, 1D evolutionary track from Amard & Matt (2020) as a colour shaded envelope taking into account the initial rotational spread. The particular case of the solar track is shown in cyan colour. Symbols represent the ASH simulations and are changed according to the internal differential rotation profile achieved in the convective envelope of the simulations, either banded/quenched, solar like (fast equator slow poles) and anti-solar (slow equator fast poles). The coloured vertical bars on the left represent the magnetic states either short, decadal and no cycles.

In the text

Fig. 2. Typical magnetic butterfly diagrams (i.e., time–latitude diagram of the axisymmetric toroidal magnetic field) achieved in the simulations published in Brun et al. (2022). On the top panel, the butterfly diagram is shown as a colour contour plot in Gauss for the rapidly rotating cases (low Rossby numbers), in the middle panel for intermediate value (Sun-like) and in the bottom panel for statistically steady solutions for slowly rotating stars with large Rossby numbers.

In the text

Fig. 3. Left: magnetic ratio between the squared toroidal, Btor, and poloidal, Bpol, components near the surface. Two observational fits (dashed black for M > 0.5 M⊙ and dotted lines for M < 0.5 M⊙) from See et al. (2015) are indicated, as well as the fit of the simulations (purple line) and their respective error bars. Right: same ratio within the convective envelope. The coloured symbols have the following meanings: red, yellow, green, and blue represent stars with 0.5, 0.7, 0.9, and 1.1 M⊙, respectively, and triangle, diamond, square, and filled circle shapes the rotation rates for, respectively, slow (the exact value depends on mass see Brun et al. 2022), one, three, and five times the solar rotation rate, Ω⊙. The error bars indicate the minimum and maximum values reached by the model during the integrated time (one or several cycle periods) used to compute the mean value indicated by the symbol.

In the text

Fig. 4. Poloidal field component squared amplitude as a function of the Rossby number, Rof. At low Rof, the agreement between the simulations (purple line) and observation (dashed black line) is quantitatively good, with both showing a similar decreasing trend. At a higher Rossby number, Rof > 1 (yellow range), not covered by the observational database, an inverse trend is indicated by the simulations, and hence suggests a minimum in poloidal field strength near Rof ∼ 1, which could have interesting consequences for stellar spin-down via wind braking. An indicative trend proportional to ${\mathrm{Ro}}_{\mathrm{f}}^{10}$ is shown using a dotted grey line. We note the V shape that the two trends (purple and dotted lines) form, with the minimum being near Rof ∼ 1. The coloured symbols have the same meanings as in Fig. 3.

In the text

Fig. 5. Left: toroidal field component squared amplitude as a function of the Rossby number, Rof, probed near the surface. At low Rof, the agreement between the simulations (purple line) and observational (dashed black line) fits is qualitatively good, with both simulations and observations showing a decreasing trend, although the field amplitude is too low near the surface (left panel). This is likely due to our choice of top magnetic boundary conditions. Right: same quantity probed in the middle of the convection zone, where the agreement is quantitatively good, with both fits being close in terms of the power law index. A possible inverse trend is indicated by the simulations at higher Rossby numbers, Rof > 1 (not covered by the observational database), and hence suggests a minimum in the toroidal field strength near Rof ∼ 1, which could have interesting consequences for stellar spin-down via wind braking.

In the text

Fig. 6. Multipolar decomposition of the magnetic field near the surface (with m components summed quadratically), showing the dipole, ℓ = 1 (left), quadrupole, ℓ = 2 (middle), and octupole, ℓ = 3 (right), as a function of the Rossby number, Rof. At low Rof, the agreement between the simulations (purple line) and observational (dashed black line) fits is quantitatively good for all low-degree dynamo modes, with all modes showing a decreasing trend and power laws comparable to observations. The possible inverse trend for high Rossby numbers is also found in all three multipoles. Symbols have the same meanings as in Fig. 5.

In the text

Fig. 7. Various trends for our theoretical proxies of both the Stokes V magnetic field measurement, BV, and the total field measurement via Zeeman spectroscopy, BI. Left: BV dependency on the Rossby number is well captured, with both fits (purple for the simulations and dashed black lines for the observations) agreeing quantitatively well in terms of the power law index. Middle: magnetic relation between BV and BI showing, as in Fig. 3, an expected positive trend. Right: BV/BI ratio as a function of Rof, with a fit (purple line) indicating a weak dependency. Symbols have the same meanings as in Fig. 8.

In the text

	Fig. 8. Theoretical proxy for Stokes V magnetic field, BV, vs. stellar mass in the simulation’s dataset. Symbols have the same meanings as in Fig. 3. We note a rather strong decreasing trend with stellar mass (purple fit exponent around −2), the rotation rate leading to some spread, and the relatively large error bar on the fit.
In the text

Fig. 9. Axisymmetric (m = 0, left) and non-axisymmetric (m ≠ 0, right) magnetic field decomposition variations with the Rossby number. No clear differences are seen between the two field components; the purple fit exponents in the low Rossby number regions being close and negative, with at best a tendency for the field to be more axisymmetric for larger Rossby number values. An inverse trend is also seen for large Rossby numbers, as was expected from the previous analysis. Symbols have the same meanings as in Fig. 8.

In the text

	Fig. B.1. Comparison of the stellar Rossby number, Ros, as a function of the fluid one, Rof, in the different models (orange dots) from Brun et al. (2022). A linear regression is proposed with the solid orange line (Ro${}_{\mathrm{s}}=(2.26\pm 0.05)R{o}_{\mathrm{f}}^{1.05\pm 0.04}$), along with the direct proportionality, Ros = Rof, illustrated with the dashed black line.
In the text