Rotational evolution of young-to-old stars with data-driven three-dimensional wind models

Solar-type stars form with a wide range of rotation rates. A wide range persists until a stellar age of 0.6 Gyr, after which solar-type stars exhibit Skumanich spin-down. Rotational evolution models incorporating polytropic stellar winds struggle to simultaneously reproduce these two regimes, namely the initially wide range and the Skumanich spin-down without imposing an a-priori cap on the wind mass-loss rate. We show that a three-dimensional wind model driven by Alfv\'en waves and observational data yields wind torques that agree with the observed age distribution of rotation rates. In our models of the Sun and twenty-seven open cluster stars aged from 0.04 to 0.6 Gyr that have observationally derived surface magnetic maps and rotation rates, we find evidence of exponential spin-down in young stars that are rapid rotators and Skumanich spin-down for slow rotators. The two spin-down regimes emerge naturally from our data-driven models. Our modelling suggests that the observed age distribution of stellar rotation rates arises as a consequence of magnetic field strength saturation in rapid rotators.


INTRODUCTION
Solar-type (FGK) stars enter the stellar spin-down era with a wide range of rotation rates Ω (Stauffer et al. 1985;Bouvier 1991;Edwards et al. 1993).The fastest rotators spin near the stellar break-up velocity and the slowest rotators are up to two orders of magnitude slower (Gallet & Bouvier 2013).Observations of stellar clusters show that this wide range of periods is eventually lost; some mechanism harmonises their rotation rates by the onset of the spin-down era at an age  s " 0.6 Gyr (Skumanich 1972), after which Ω 9  ´1{2 .To reproduce the observed Ω distribution, spin-down models have incorporated a 'knee' in the 9 ΩpΩq relation (e.g.Stauffer & Hartmann 1987;Chaboyer et al. 1995;Krishnamurthi et al. 1997), so that 9 Ω 9 Ω for rapid rotators and 9 Ω 9 Ω 3 for slow rotators as required to attain Skumanich spin-down (e.g.Durney & Latour 1978).
Solar-type stars spin-down by means of magnetised stellar winds (e.g.Parker 1958;Schatzman 1962) both before and during the Skumanich spin-down era.From one-dimensional models, we can write the wind torque as 9  9 9 Ω  a where 9  is the wind mass-loss rate and  a is the Alfvén radius (Weber & Davis 1967;Mestel 1968Mestel , 1984;;Kawaler 1988).Mass-loss rate saturation in rapid rotators (e.g.Cranmer & Saar 2011) can be imposed to produce the required knee in the 9 ΩpΩq relation.In polytropic magnetohydrodynamic wind models (Vidotto et al. 2009;Matt et al. 2012;Réville et al. 2015;Finley & Matt 2017;Finley et al. 2018), the wind mass loss rate is governed largely by the coronal temperature and wind density.By varying the wind density, and thus 9 , models are able to match spin-down 'gyrotracks' Ωpq found by integrating 9 Ω to observational data (Gallet & Bouvier 2013;Johnstone et al. 2015; Sadeghi ‹ E-mail: evensberget@strw.leidenuniv.nlArdestani et al. 2017).This mass loss rate saturation is, however, externally imposed on the wind model.In contrast to polytropic wind models, an Alfvén wave based wind model (Sokolov et al. 2013;van der Holst et al. 2014) can shift the wind model boundary downwards into the much colder and denser stellar chromosphere, and compute the coronal wind density and 9 consistently with 9  and  a (e.g.Garraffo et al. 2015;Alvarado-Gómez et al. 2016;Cohen 2017).
In this letter, we show that the resulting 9 Ω of the Sun and a set of twenty-seven young, solar-type cluster stars aged from "40 Myr to "600 Myr (Evensberget et al. 2021(Evensberget et al. , 2022(Evensberget et al. , 2023) ) appears compatible with the observed age distribution of stellar rotation rates Ω (Gallet & Bouvier 2013) for a solar mass star without imposing any a priori assumptions about mass loss saturation.In Section 2, we show this by comparing instantaneous wind torque 9  w of our models to a simple observation-based gyrotrack model based on the work of Gallet & Bouvier (2013) containing hundreds of rotation rates.The resulting wind torques do however need to be scaled by a factor of " 7 in order to yield the expected Skumanich spin-down rate at solar maximum conditions, as presented in Section 3.This scaling factor is the only free parameter in our calculated wind torques.In Section 4 we present our conclusions.

WIND TORQUE MODEL
In our recent work we have modelled the winds of thirty young, solartype stars (Evensberget et al. 2021(Evensberget et al. , 2022(Evensberget et al. , 2023) ) with well characterised ages, and whose surface magnetic fields have previously been mapped with Zeeman-Doppler imaging (ZDI) by Folsom et al. (2016Folsom et al. ( , 2018)).The modelled stars are part of the Hyades, Coma Berenices, Hercules-Lyra, Pleiades, AB Doradus, and Columba-Horologium-Tucana clusters1 .As cluster members, these stars' ages are known to within 5 % to 20 %.For reference the stars' ages, rotation rates, masses and radii are given in Table 1.This is the largest sample of stars with well characterised ages and ZDI surface magnetic fields derived with a consistent methodology.
For each star a wind model was obtained by letting the surface magnetic field drive a three-dimensional numerical wind model (Powell et al. 1999;Tóth et al. 2012) incorporating Alfvén wave coronal heating (Sokolov et al. 2013;van der Holst et al. 2014).The resulting wind solutions permit us to simultaneously compute consistent key quantities such as wind mass loss rate 9  w and wind angular momentum loss rate 9  w as described in Appendix A. Identical solar parameter values are used in all the wind models; the only varying parameters are the surface magnetic map, stellar rotation rate, mass, and radius as indicated in Table 2.

Rotational braking from wind torques
The stellar spin-down rates 9 Ω w are computed from (scaled) wind torques 9 w that we previously published (Evensberget et al. 2023).In this paper we apply a scaling factor of 7 to our stellar and solar 9  w values to match the expected Skumanich spin-down rate at solar maximum conditions.
To compute spin-down rates from wind torques, we use the stellar rotational inertia values  from Baraffe et al. (2015).As  " Ω we get (e.g.Johnstone et al. 2015;Matt et al. 2015) where 9  w is the wind torque and Ω 9  is the contraction torque (Gallet & Bouvier 2015).In the Baraffe et al. models  is essentially steady by  " 100 Myr so that the wind torque term dominates 9 Ω w .For convenience, the  and values used are given in Appendix B. In Fig. 1 the stellar rate of rotation is plotted against the stellar age (coloured circles).The stellar spin-down rates 9 Ω w derived from our wind models are shown as coloured arrow attached to the symbol of each star, and pointing in the direction of 9 Ω w .To guide the eye, the shaded region of Fig. 1 shows the range between the 25 th and 90 th percentile of observed Ω values (curves from Gallet & Bouvier 2013).The pre-Skumanich spin-down era may be identified as the time where there is significant difference between the two percentiles (up to  " 600 Myr).This era is followed by the Skumanich spin-down era where Ω 9  ´1{2 .From Fig. 1 we have good coverage of the Gallet & Bouvier (2013) sample around 0.1 Gyr and 0.6 Gyr, however our ZDI dataset has no stars between onset of the Skumanich relation and the solar models.Gallet & Bouvier (2013) fitted Ωpq curves, so-called 'gyrotracks', to the observed rotation rates of young and very young ( ă 30 Myr) clusters, incorporating mass loss saturation, core-envelope decoupling, and other effects.recently GAIA data has been used to infer the age distribution of rotation rates for a much larger sample of stars.The essential features of the revised gyrotracks are however retained (e.g.Godoy-Rivera et al. 2021, Fig. 13) in the age and mass range where we have wind model data.Here, we fit a very simple model on the form

Rotational braking from cluster spin rates
to the observation-based Gallet & Bouvier (2013) 25 th and 90 th percentile Ωpq gyrotrack curves, focusing on the spin-down era 30 Myr ă  ă 4.6 Gyr for a one solar mass star2 .We note that much of the signal for solar analogues is in stars less than 100 Myr old (Denissenkov et al. 2010;Gallet & Bouvier 2013, 2015).
Our simplified model (2) permits gyrotracks Ω obs pq to decay exponentially when Ω ą Ω sat , while adhering to the Skumanich (1972) relation Ωpq 9  ´1{2 for Ω ă Ω sat .We have parametrised the pre-Skumanich spin-down era in terms of an exponential decay constant  1 so that the rotational spin-down timescale in this era Ω{ 9 Ω obs "  ´1 1 " 120 Myr.The value of  2 is then found by matching the solar age and angular velocity:

From continuity of 9
Ω obs the saturation threshold is Ω sat " a  1 { 2 .This gives a reasonable Ω sat " 8.9 Ω d , given that this parameter can vary between 3 Ω d and 15 Ω d depending on model choices (Matt et al. 2015;Amard et al. 2016).In Fig. 1 the resulting 9 Ω obs values are shown as a vector field (little grey arrows).By design we have excellent agreement between the 25 th and 90 th percentile Ωpq gyrotrack curves (edges of the grey region) and the 9 Ω obs vector field (grey arrows), indicating that eq. ( 2) captures the dominant features of the Gallet & Bouvier (2013) spin-down era gyrotracks for a Solar-mass star.In addition to the vector field, we have also added small trailing gyrotrack segments Ω obs pq to each star.By comparing the slope of each stars' trailing trajectory segment Ω obs pq to its emerging arrow 9 Ω w it can be seen that there is good agreement between the wind torque based spin-down rates and the gyrotrack segments.
As mentioned in Section 2.1, we apply a scaling factor of 7 to the wind torques 9  w of Evensberget et al. (2023).In our dataset this corresponds to matching the wind torque 9  w at solar maximum to the Sun's gyrotrack predicted torque 9  obs .The scaling effect can be seen in Fig. 1 where the most downwards-pointing of the wind torque arrows extending from the Sun symbols is in good agreement with the gyrotrack trail.The wind torque scaling is discussed further in Section 3.

Saturation in rapid rotators
The good agreement between the wind torque based spin-down rates 9 Ω w (coloured arrows) and the gyrotracks 9 Ω obs pq (coloured trailing curve segments) in Fig. 1 suggests the presence of a saturation effect in the 9 Ω w values around Ω sat .This is confirmed in Fig. 2 where the top panel shows the calculated spin-down timescale Ω{ 9 Ω w plotted against the rotation period as coloured symbols, and Ω{ 9 Ω obs (see equation 2) as a grey curve consisting of two line segments.Vertical distance between the star symbols and the gray curve correspond to trail-arrow mismatch in Fig. 1.
The dashed black curve in each panel shows a fitted broken power law { sat " pΩ{Ω sat q (3) Ω w plotted against stellar age for the stellar wind models in Evensberget et al. (2021Evensberget et al. ( , 2022Evensberget et al. ( , 2023)).The angular velocity of each star is plotted against the star's age with a coloured symbol.The spin-down rate 9 Ω w is indicated with an arrow extending from the symbol.Larger values of 9  w and thus of 9 Ω w produce more downwards-pointing arrows in the plot.In order for the symbols not to overlap each other we have added a small pseudorandom shift to the stellar ages.As the Hyades and Coma Berenices stars are close together in both Ω and age, we have added an inset in the top right corner which gives a zoomed in view of this region.The Sun symbol and the two black arrows extending from it represent the solar minimum and solar maximum models from Evensberget et al. (2023).In the plot background we have added a shaded region corresponding to the region between the 25 th and 90 th percentile of stellar Ω obs pq gyrotracks as a function of age for a solar mass star from Gallet & Bouvier (2013).The plot background also contains a vector field (grey arrows) of 9 Ω obs values from the fitted model described by equation ( 2).The vertical dashed line denotes the onset of the Skumanich spin-down era, while the horizontal dash-dotted line separates the saturated and un-saturated spin-down regimes as noted in the top and right margins. in which Ω sat " 8.9 Ω d is fixed, while  sat and  are found by simultaneous fitting.The exponent  is permitted to take on a different value on each side of Ω sat .The  and  sat values are given in Table 3, and the  values are annotated in Fig. 2. From its top panel it can be seen that the fitted curve closely matches the observed Ω{ 9 Ω obs relation of equation ( 2) shown and annotated in grey.
The middle panel of Fig. 2 show the wind mass loss rate 9  w plotted against the rotation period.The wind mass loss rate exhibits a similar pattern of saturation (notwithstanding a slightly different scatter) around Ω sat as that of the spin-down timescale.The mass loss rate 9  appears to saturate for Ω ą Ω sat in a similar manner as Ω{ 9 Ω w .
The bottom panel of Fig. 2 shows the surface unsigned magnetic flux Φ " 4 2   plotted against the rotation period3 .Once again the fitted dashed curve and the scatter have a similar shape those of the preceding panels.Hence both the 9  w and Φ values exhibit similar trends as the top panel Ω{ 9 Ω w values, albeit with different slopes and intersections with the range of solar values.
Table 1.Relevant model parameters and results from the wind models (Evensberget et al. 2023) used in this work and stellar rotational inertia values.The 'case' column denotes the model case number from Evensberget et al. (2023).The 'type' column gives the star's spectral type.When the spectral type is calculated from  eff following Pecaut & Mamajek (2013) this is denoted by a dagger symbol (:).The 'age' column gives the stellar age, and the following columns give the stellar mass , radius , and rotational period  rot .The   value is the average unsigned radial magnetic field over the stellar surface from the Zeeman-Doppler maps of Folsom et al. (2016Folsom et al. ( , 2018)), from which the the previous columns values were also taken.The wind mass loss rate 9  w and wind torque 9  w are computed from the wind models.A scaling factor of 7 has been applied to the 9  w values as described in Section 3. In the main text this factor is incorporated into the definition of 9  w .The rotational inertia  values are taken from the models of Baraffe et al. (2015).The stellar age column can also be used to identify the parent cluster as follows: Columba-Horologium-Tucana: 42 ˘6 Myr, AB Doradus: 120 ˘10 Myr, Pleiades: 125 ˘8 Myr, Hercules-Lyra: 257 ˘46 Myr, Coma Berenices: 584 ˘10 Myr, Hyades: 625 ˘50 Myr.Poynting flux-to-field ratio 1.1 ˆ10 6 W m ´2 T

DISCUSSION
Multiple, possibly interacting physical phenomena have been proposed as causes for the torque saturation in rapid rotators.These  Ω obs (grey curve) from equation (2).The two Sun symbols connected by a black bar represents the solar maximum and minimum.The second panel shows the wind torque 9  w (1 erg " 10 ´7 N m) and the third panel shows the wind mass loss rate 9  w values (the commonly used solar mass loss rate value is 2 ˆ10 ´14  d yr ´1 " 1 ˆ10 9 kg s ´1).The bottom panels show the surface unsigned magnetic flux Φ which is determined directly from the magnetic maps driving the numerical models.The dashed black lines show the broken power law fit of equation ( 3).The vertical dash-dotted line shows Ω sat " 8.9 Ω d which separates the saturated (Ω ą Ω sat ) and un-saturated (Ω ă Ω sat ) spin-down regimes.

2015)
. Another cause could be a change in the nature of emerging photospheric magnetic field and/or the stellar corona (Vilhu 1984;Pizzolato et al. 2003;Wright et al. 2011;Reiners & Mohanty 2012;Marsden et al. 2014) around Ω sat .MacGregor & Brenner (1991), for example, suggested that the surface magnetic field strength  may saturate in rapid rotators.Alternatively, the knee could be a consequence of a change in the magnetic field complexity (Garraffo et al. 2015(Garraffo et al. , 2018;;Finley & Matt 2018).See et al. ( 2019), however, found that that field complexity plays a modest role for reasonable mass-loss rates.
In Fig. 2 we see that the plotted quantities Ω{ 9 Ω w , 9  w , 9  w and Φ exhibit a similar structure, in which the quantity saturates for Ω Á Ω sat .This suggests that the saturation of the wind torque and wind mass loss rate is a consequence of the saturation of the surface unsigned magnetic flux Φ in our models.The result is consistent with the findings of Evensberget et al. (2023) where no break was found when considering 9  w and 9  as a function of surface unsigned magnetic flux or field complexity, and the residual spread was attributed to stellar magnetic cycles.We thus attribute the saturation in Fig. 2 to a saturation in the surface unsigned magnetic flux (and surface magnetic field strength since Φ " 4 2   ) for Ω Á Ω sat .The mass loss rate 9 saturation arises naturally in our models, while mass loss rate saturation at high rates of rotation is externally imposed in many models of spin-down (Gallet & Bouvier 2013, 2015;Johnstone et al. 2015;Sadeghi Ardestani et al. 2017;Ahuir et al. 2020).Such mass loss rate saturation may (or may not) have observational support (Vidotto 2021;Wood et al. 2021).In the Alfvén wave based model used in this work there is an approximate power-law relationship between 9  and Φ (Evensberget et al. 2023) which means that the mass loss rate saturation is a consequence of magnetic saturation.Similar mass loss rate saturation may be seen also in the one-dimensional Alfvén wave based model of Cranmer & Saar (2011), and in the models of Suzuki (2011); Suzuki et al. (2013); Shoda et al. (2020).Stellar magnetic activity is known to saturate with increasing rotation (Noyes et al. 1984;Pizzolato et al. 2003) and there are signs of saturation of the surface magnetic field strength (Reiners et al. 2009;Vidotto et al. 2014).This is supported by the saturation in spot filling factors found by Cao & Pinsonneault (2022).Early models (Mac-Gregor & Brenner 1991) suggested that that  9 Ω  where  « 1 (linear dynamo) in the un-saturated regime and « 0 (saturated dynamo) in the saturated regime.In our dataset we find that the surface unsigned magnetic flux in the un-saturated regime Φ 9 Ω 1.25˘0.07 is matched by the   9 Ω 1.32˘0.14result of Vidotto et al. (2014).In comparison with the spin-down timescale and mass loss rate values, the surface unsigned magnetic flux appears to continue increasing past Ω sat , albeit at a slower rate as in Reiners et al. (2014Reiners et al. ( , 2022)).
We have applied a scaling on the wind torque values 9  w in order to match the gyrotrack torque 9  obs at solar maximum.In our work the scaling factor has the same value of 7 for all the stars modelled (and the Sun), and is thus independent of Ω.Similar torque scalings have a long history in the literature (Krishnamurthi et al. 1997).Solar wind torque estimates based on the solar surface magnetic field strength can be "20 times smaller than the torque required by the Skumanich relation and the Ω obs pq gyrotracks (Finley & Matt 2018;Finley et al. 2018).In See et al. (2019) the torque values of Finley & Matt (2018) were scaled by "25 to match the gyrotracks of Matt et al. (2015).Similarly, the Matt et al. (2012) wind torque formula, which was used in the Gallet & Bouvier (2013) gyrotracks and by Johnstone et al. (2015) was scaled by a factor of "2 and "11 respectively to match cluster rotation rate observations.In this work, we find that a single, constant 7ˆtorque scaling factor is sufficient for our Alfvén based wind model to agree with observations.
The need for the 7ˆtorque scaling may be rooted in underestimates of the stellar magnetic field strengths (Morin et al. 2010;Yadav et al. 2015;Lehmann et al. 2019) and solar magnetic field strengths, the latter in which a scaling of 2-4 is often applied to match the measured magnetic field at 1 au (Cohen et al. 2006;Oran et al. 2013;Linker et al. 2017;Pognan et al. 2018;Finley et al. 2018;Sachdeva et al. 2019;Riley et al. 2021).Zeeman-Doppler imaging reconstructs large-scale magnetic features on the stellar surface in an Ω-dependent way (Hussain et al. 2000;Morin et al. 2010), but only a fraction of the field strength is recovered (e.g.Lehmann et al. 2019).These two biases may be used to justify a scaling of 9  w including the 7ŝ caling applied in this work.On a similar note we find that scaling the surface magnetic field by a factor of "5 as in Evensberget et al. (2021Evensberget et al. ( , 2022Evensberget et al. ( , 2023) ) we also obtain gyrotracks with both exponential and Skumanich-type spin-down.We do not pursue this alternative approach here as the resulting wind mass loss rate values differ for the solar case (see Evensberget et al. 2023).We do however emphasise that the qualitative indications of a saturated and an un-saturated magnetic regime also appear in this model series.
Alternatively, the need for the 7ˆtorque scaling may be attributed to core-envelope decoupling (Stauffer & Hartmann 1986;Denissenkov et al. 2010;Gallet & Bouvier 2013, 2015), where the stellar convective envelope may spin at a different rate than the core while exchanging rotational momentum with it at a (configurable) coupling timescale .Core-envelope decoupling would change the direction of the coloured arrows in Fig. 1 and, depending on the length of , it could have the same effect as changing the torque scaling factor or scaling the ZDI magnetic field strength.A variable  was invoked by Gallet & Bouvier (2013, 2015) to create their observed-based gyrotracks extending from the stellar contraction era to the Sun's age.In the Baraffe et al. (2015) models "20 % of our stars' rotational inertia is in the envelope.This could produce an effect similar to required torque scaling for long , in that the envelope spins down with little angular momentum exchange with the core in parts of the stars' life.Core-envelope decoupling has theoretical support from hydrodynamic modelling (Endal & Sofia 1981;Rucinski 1988;Pinsonneault et al. 1989;Zahn 1992;Amard et al. 2016Amard et al. , 2019) ) but it has been suggested that the inclusion of magnetohydrodynamic instabilities such as the Tayler (1973) instability may effectively flatten the internal rotation curve (Fuller et al. 2019;Eggenberger et al. 2022).Some models of stellar spin-down (Collier Cameron & Li 1994;Matt et al. 2015;Johnstone et al. 2015) have applied a model of solid body rotation.The recently discovered 'stalled spin-down' (Agüeros et al. 2018;Curtis et al. 2019;Douglas et al. 2019;David et al. 2022) of K dwarfs is however attributed to core-envelope decoupling (e.g.Spada & Lanzafame 2020;Cao et al. 2022Cao et al. , 2023)), although a magnetic wind effect may also play a role (Curtis et al. 2020).

CONCLUSIONS
We have modelled the stellar spin-down rates of 27 young, solar-type stars based on wind models published in Evensberget et al. (2021Evensberget et al. ( , 2022Evensberget et al. ( , 2023) ) and whose surface magnetic fields were reconstructed by Folsom et al. (2016Folsom et al. ( , 2018) ) using Zeeman-Doppler imaging.We have compared our resulting wind torques and spin-down rates with the observation-based Ωpq gyrotracks of Gallet & Bouvier (2013).
Unlike many previous works, our wind torques are derived from fully three-dimensional wind models where the corona is heated by Alfvén waves.As our wind model (Sokolov et al. 2013;van der Holst et al. 2014) extends into the chromosphere, the coronal density and temperature emerge naturally as a consequence of the coronal heating by Alfvén wave turbulence.This means that we do not need to impose a desired mass loss rate by adjusting the coronal density, as is required in polytropic wind models (Vidotto et al. 2009;Matt et al. 2012;Réville et al. 2015;Finley & Matt 2017;Finley et al. 2018).Instead, the mass loss rate emerges naturally as a consequence of the coronal heating, and ultimately the surface magnetic field strength.
We find that the torque and the mass loss rate produced by our wind model (using the same solar configuration parameters across the entire age and rotation rate range) are compatible with features of the observationally derived gyrotracks Ω obs pq for both fast and slow rotators from the onset of the spin-down era until the Sun's age.We find signs of both the exponential ( 9Ω 9 Ω) spin-down rate for rapid rotators and the Skumanich-type ( 9Ω 9 Ω 3 ) spin-down rate for slow rotators.In our models the two regimes arise as a consequence of the ZDI observed magnetic field strength saturation at rapid rates of rotation.Our model is novel as it uses only a single free scaling parameter in order to agree with observations.By scaling all the wind torque based spin-down rates 9 Ω w by a factor of 7, so that the solar maximum value matches the 9 Ω d value predicted from the Skumanich (1972) relation Ω 9  ´1{2 , we find good agreement between the gyrotrackderived spin-down time-scales Ω{ 9 Ω obs and the wind torque-derived spin-down time-scales Ω{ 9 Ω w .Such scalings have long been used in evolutionary models (e.g.Krishnamurthi et al. 1997).The scaling may be needed because of changes in the effective rotational inertia due to core-envelope decoupling, or due to ZDI underestimates of the unsigned magnetic flux; both effects would lead to an increase in 9 Ω w .
Our work suggests that the wind torque saturation is caused by a saturation of the surface magnetic field strength, rather than a qualitative change in the magnetic field complexity for rapidly rotating stars.To further confirm or refute our findings it would be necessary to have more ZDI data of stars with very well constrained ages.ZDI data of stars with rotation periods in the range 0.2 d to 1 d and ą10 d would be of particular interest.3.

B: SPIN-DOWN TIMESCALES
To compute the spin-down timescales Ω{ 9 Ω w we use the scaled wind torque 9  w values from Evensberget et al. (2023) and the rotational inertia  from Baraffe et al. (2015).In this model the solar rotational inertia  d " 7.1 ˆ10 46 kg m 2 .The exact values used for each star is given in Table 1.We note that there is some spread in reported solar rotational inertia values, as  d " 5.8 ˆ10 46 kg m 2 in e.g.Cox (2000).
The rotational inertia  values are not part of the wind model input, but are used to compute the spin-down timescales Ω{ 9 Ω w as described in Section 2.1 with equation ( 1).Note that we do not consider the contraction torque term Ω 9 in this work so that in effect 9  w " 9 Ω w  in equation ( 1).In this work we also do not consider coreenvelope decoupling (e.g.Stauffer & Hartmann 1986).The spindown timescales Ω{ 9 Ω w for each star are computed from the values given in Table 1.These values are also plotted in Fig. 1 and the top panel of Fig. 2. Fig. B1 shows the residual ratio 9 Ω w { 9 Ω obs for each star.The residual ratio is computed from the gyrotrack spin-down timescale Ω{ 9 Ω obs and the wind spin-down timescale Ω{ 9 Ω w .The numerical values of Ω{ 9 Ω obs , Ω{ 9 Ω w , and 9 Ω w { 9 Ω obs are given in Table B1.The dark, coloured bars in Figure B1 are 95 % confidence intervals of the mean 9 Ω w { 9 Ω obs in each cluster.The confidence intervals are computed in the standard way (e.g.Draper 1998) using the Student's  distribution with  ´1 degrees of freedom, where  is the number of stars in each cluster.For the solar minimum and maximum models we do not compute or show a confidence interval as there is only one model for each of these cases.Similarly for Q BD-16 351 no confidence interval can be computed as we only have one model aged 42 ˘6 Myr.The confidence interval parameters are given in Table B2.
The confidence intervals can be seen to cover the 9 Ω w { 9 Ω obs " 1 line in Fig. B1 for each of the clusters except the Hyades where the confidence interval ends just above the line.Tentatively, it appears that the Hercules-Lyra cluster (green symbols) has a lower 9 Ω w { 9 Ω obs ratio than the other clusters.However, the confidence intervals are wide and the difference is not statistically significant.Similarly a trend of increasing 9 Ω w { 9 Ω obs with age is not statistically significant.We attribute the residual spread in 9 Ω w { 9 Ω obs to stellar magnetic cycles (see Evensberget et al. 2021) affecting the surface magnetic flux Φ (a source of scatter in the bottom panel of Fig. 2), but systematic errors in magnetogram detail with increasing Ω (Morin et al. 2010) may also play a role, as may stellar mass differences (see Table 1).This paper has been typeset from a T E X/L A T E X file prepared by the author.B2.
Table B1.Spin-down timescales computed in this work, and the residual ratio between the gyro-track spin-down timescale and the wind spin-down timescale.The residual ratios are plotted against stellar age in in Fig. B1.B2.Residual ratio statistics for each cluster.For each cluster of this study we give the age, followed by the (dimensionless) mean residual ratio 9 Ω w { 9 Ω obs and the 95 % confidence interval of the mean.These parameters can also be seen visualised in Fig. B1.In the solar case the range between our minimum and maximum models is 0.

Figure 1 .
Figure 1.Stellar rotational period and wind-derived spin-down rates 9Ω w plotted against stellar age for the stellar wind models inEvensberget et al. (2021Evensberget et al. ( , 2022Evensberget et al. ( , 2023)).The angular velocity of each star is plotted against the star's age with a coloured symbol.The spin-down rate 9 Ω w is indicated with an arrow extending from the symbol.Larger values of 9  w and thus of 9 Ω w produce more downwards-pointing arrows in the plot.In order for the symbols not to overlap each other we have added a small pseudorandom shift to the stellar ages.As the Hyades and Coma Berenices stars are close together in both Ω and age, we have added an inset in the top right corner which gives a zoomed in view of this region.The Sun symbol and the two black arrows extending from it represent the solar minimum and solar maximum models fromEvensberget et al. (2023).In the plot background we have added a shaded region corresponding to the region between the 25 th and 90 th percentile of stellar Ω obs pq gyrotracks as a function of age for a solar mass star fromGallet & Bouvier (2013).The plot background also contains a vector field (grey arrows) of 9 Ω obs values from the fitted model described by equation (2).The vertical dashed line denotes the onset of the Skumanich spin-down era, while the horizontal dash-dotted line separates the saturated and un-saturated spin-down regimes as noted in the top and right margins.

Figure 2 .
Figure 2.In the top panel we compare the spin-down timescale inferred from the wind models Ω{ 9 Ω w (coloured symbols) and the spin-down timescale Ω{ 9Ω obs (grey curve) from equation (2).The two Sun symbols connected by a black bar represents the solar maximum and minimum.The second panel shows the wind torque 9  w (1 erg " 10 ´7 N m) and the third panel shows the wind mass loss rate 9  w values (the commonly used solar mass loss rate value is 2 ˆ10 ´14  d yr ´1 " 1 ˆ10 9 kg s ´1).The bottom panels show the surface unsigned magnetic flux Φ which is determined directly from the magnetic maps driving the numerical models.The dashed black lines show the broken power law fit of equation (3).The vertical dash-dotted line shows Ω sat " 8.9 Ω d which separates the saturated (Ω ą Ω sat ) and un-saturated (Ω ă Ω sat ) spin-down regimes.

Figure A1 .
Figure A1.Mean Alfvén radius  plotted against the stellar rotation period as in Fig.2.The fit parameters for  a are given in Table3.

Figure B1 .
Figure B1.Residual ratio 9Ω w 9 Ω obs .The stellar symbol represent solar maximum and minimum.Values smaller than unity are associated with reduced rotational braking in Fig.1.The dark vertical lines show 95 % confidence intervals for the average 9 Ω w { 9 Ω obs in each cluster.Numerical values pertaining to this plot is given in TableB2.

Table 2 .
Wind modelling parameters used for all the stellar and solar models used in this work.The parameters that vary between wind models are the surface magnetic map, stellar rotation rate, mass, and radius.All other modelling parameters are fixed.Note that the stellar age is not a wind model input.

Table 3 .
Fitted coefficients for the broken power law { sat " pΩ{Ω sat q  .The first two numerical columns give the exponent  values on either side of the saturation point Ω sat " 8.9 Ω d .The final column gives the value of log 10  sat for  sat " pΩ sat q.The first three of these relationships are plotted in Fig.2.The average Alfvén radius  a is plotted in Fig.A1.The coefficient of determination  2 is given in the final column.