Stellar Cruise Control: Weakened Magnetic Braking Leads to Sustained Rapid Rotation of Old Stars

Despite a growing sample of precisely measured stellar rotation periods and ages, the strength of magnetic braking and the degree of departure from standard (Skumanich-like) spindown have remained persistent questions, particularly for stars more evolved than the Sun. Rotation periods can be measured for stars older than the Sun by leveraging asteroseismology, enabling models to be tested against a larger sample of old field stars. Because asteroseismic measurements of rotation do not depend on starspot modulation, they avoid potential biases introduced by the need for a stellar dynamo to drive starspot production. Using a neural network trained on a grid of stellar evolution models and a hierarchical model-fitting approach, we constrain the onset of weakened magnetic braking. We find that a sample of stars with asteroseismically-measured rotation periods and ages is consistent with models that depart from standard spindown prior to reaching the evolutionary stage of the Sun. We test our approach using neural networks trained on model grids produced by separate stellar evolution codes with differing physical assumptions and find that the choices of grid physics can influence the inferred properties of the braking law. We identify the normalized critical Rossby number ${\rm Ro}_{\rm crit}/{\rm Ro}_\odot = 0.91\pm0.03$ as the threshold for the departure from standard rotational evolution. This suggests that weakened magnetic braking poses challenges to gyrochronology for roughly half of the main sequence lifetime of sun-like stars.

In essentially all of these calibrators, rotation rates are measured by observing spot modulation due to dark starspots rotating in and out of view.The high photometric precision of the Kepler Space Telescope (Borucki et al. 2010), and the subsequent K2 mission (Howell et al. 2014), enabled predictions for magnetic braking to be tested on a wealth of open clusters and associations (see Cody et al. 2018) as well as a population of older field stars (McQuillan et al. 2014;Santos et al. 2021).
In addition to starspot modulation used to detect rotation, brightness modulations due to stellar oscillations are measurable in the high-precision, long-baseline Kepler time series photometry (Huber et al. 2011).Asteroseismology-the study of these oscillations-provides valuable information about the internal structure and evolution of stars.Specifically, stellar rotation rates can be measured from the mode frequencies (Nielsen et al. 2015;Davies et al. 2015;Hall et al. 2021) and ages can be inferred by comparisons with stellar models (Metcalfe et al. 2014(Metcalfe et al. , 2016;;Silva Aguirre et al. 2015;Creevey et al. 2017).
When the ages of older, sun-like field stars were asteroseismically measured with Kepler data, they were found to maintain surprisingly rapid rotation late into their main sequence lifetimes (Angus et al. 2015).To explain this sustained rapid rotation, it was proposed that stars diverge from the "standard spindown" model and enter a phase of "weakened magnetic braking" (WMB; van Saders et al. 2016van Saders et al. , 2019)).When stellar rotation was measured using asteroseismology rather than spot modulation, the observed rotation periods were consistently faster than predicted by the standard spindown model and evidence for WMB strengthened (Hall et al. 2021).Asteroseismology measures internal rotation rates in the stellar envelope, making it insensitive to surface differential rotation (Nielsen et al. 2015) and stellar inclination (Davies et al. 2015); additionally, asteroseismology can measure rotation rates for stars with weak surface magnetic activity and therefore undetectable spot modulation signals (Chaplin et al. 2011).These features allow asteroseismic rotation periods to avoid potential biases present in measurements from spot detection.
Careful analysis of pileups in the temperature-period distribution of sun-like stars also supported the WMB model.Studies of rotation rates in the Kepler field identified an upper envelope in stellar mass versus rotation period that matched a gyrochrone at ∼4 Gyr (Matt et al. 2015).An upper edge to the distribution could be caused by either a magnetic transition or detection bias in spot modulation (van Saders et al. 2019).Forward modeling of the Kepler field predicted a pileup of rotation periods in the weakened braking scenario that was not seen in the data, but van Saders et al. (2019) argued that errors in the measured effective temperatures were obscuring the feature.With refined measurements of stellar effective temperature, the predicted pileup in the temperatureperiod distribution was identified (David et al. 2022).
A study of sun-like stars with projected rotation periods measured from spectroscopic line broadening found them to be inconsistent with the Skumanich relation beyond ∼2 Gyr (dos Santos et al. 2016), supporting a departure from standard spindown.This sample was later revisited (Lorenzo-Oliveira et al. 2019), and the analysis suggested that the smooth rotational evolution scenario was favored, and if weakened braking takes place, it occurs at later times (≳ 5.3 Gyr).However, these measurements faced biases introduced by an uncertain distribution of inclinations, which can inflate rotation periods measured spectroscopically.
The physical mechanism that would lead to WMB remains uncertain, though some have proposed that a transition in the complexity of the magnetic field could reduce magnetic braking efficiency (Réville et al. 2015;Garraffo et al. 2016;van Saders et al. 2016;Metcalfe et al. 2016Metcalfe et al. , 2019)).Because the transition may to be rooted in the strength and morphology of the magnetic field, it is challenging to test with surface rotation rates measured through spot modulation, which require active stellar dynamos to drive starspot production (Matt et al. 2015;Reinhold et al. 2020).
To effectively use gyrochronology to estimate stellar ages, it is essential to understand when the transition to weakened braking occurs.Previous studies have provided estimates for the onset of WMB (van Saders et al. 2016(van Saders et al. , 2019;;David et al. 2022), but fully hierarchical modeling for the braking law has not been previously performed.As the departure from standard spindown depends on the dimensionless Rossby number and is predicted to be shared between all stars (van Saders et al. 2016), the problem is inherently hierarchical.Here, we provide new constraints on the evolutionary phase at which stars undergo weakened braking.We build on previous efforts (e.g.Hall et al. 2021) by modeling the rotational evolution of each star individually.
We apply a Hierarchical Bayesian Model (HBM) to constrain the population-level parameters for a WMB model.The use of an HBM has been shown to increase the precision of inferred stellar properties for high-dimensional models (Lyttle et al. 2021).Here, we model the weakened braking parameters as global properties shared by all stars, while simultaneously fitting individual stellar properties.We test the results of our fit using multiple model grids, and compare the performance of a WMB model to standard spindown.By comparing results between multiple model grids, we provide the first constraints on biases introduced by the choices of grid physics when modeling stellar rotational evolution.We find that weakened braking likely occurs before stars reach the evolutionary phase of the Sun.

DATA
We fit our rotational model to open clusters, the Sun, and Kepler field stars with asteroseismic measurements to ensure that we capture the early rotational evolution prior to the onset of weakened braking in addition to the behavior on the latter half of the main sequence.The seismic sample that best probes braking generally lies within 0.2 M ⊙ of the Sun and covers a wide range of ages.Stars hotter than 6250 K (∼1.2 M ⊙ ) lack deep convective envelopes on the main sequence, and do not undergo significant magnetic braking, and the seismic signals of stars cooler than 5000 K (∼0.8 M ⊙ ) have low pulsation amplitudes and are challenging to measure.We describe our calibrator sources in the following section.

Asteroseismic Sample
We also included a sample of Kepler field stars with asteroseismically-measured rotation rates and ages from Hall et al. (2021, hereafter Hall21).Rotation rates for main sequence stars can be challenging to measure with starspot modulation, particularly for older and less active stars, due to long rotation periods and diminished stellar activity.However, the rotational splitting of asteroseismic oscillation frequencies can be observed for stars in the end stages of the main sequence, and provides invaluable benchmarks for WMB.Hall et al. (2021) used asteroseismic mode splitting to measure rotation periods for 91 Kepler dwarfs.We augmented the Hall21 sample with two additional stars with asteroseismic rotation measurements in the wide binary system HD 176465 (KIC 10124866; White et al. 2017).The A and B components of this system are sometimes referred to by their nicknames Luke and Leia, respectively.The rotation periods reported in White et al. (2017) were derived by fitting asteroseismic mode splitting, following the same approach as Hall21.
We performed asteroseismic modeling for Luke & Leia and 47 stars from the Hall21 sample that fall within our desired mass range using version 2.0 of the Asteroseismic Modeling Portal1 (AMP; Metcalfe et al. 2009;Woitaszek et al. 2009;Metcalfe et al. 2023).This optimization method couples a parallel genetic algorithm (Metcalfe & Charbonneau 2003) with MESA stellar evolution models (Paxton et al. 2019) and the GYRE pulsation code (Townsend & Teitler 2013) to determine the stellar properties that most closely reproduce the observed oscillation frequencies and spectroscopic constraints for each star.The choices of input physics are nearly all the default choices in MESA release 12778, and the models include gravitational settling of helium and heavy elements (Thoul et al. 1994) as well as the two-term correction for surface effects proposed by Ball & Gizon (2014).The resulting asteroseismic sample is shown in panel (b) of Figure 1, while the stellar properties and rotation periods can be found in Table 1, which includes maximum-likelihood estimates of the age, mass, composition, and mixing-length from our AMP modeling.
With masses derived from asteroseismic modeling, we made mass cuts (0.8 M ⊙ ≤ M ≤ 1.2 M ⊙ ) to ensure our sample would fall within the bounds of our model grids.Previous studies have indicated that rotation periods in field stars < 7 days are likely due to non-eclipsing short-period binaries (Simonian et al. 2019(Simonian et al. , 2020)), and we therefore remove three stars (KIC 6603624, KIC 8760414, KIC 8938364) from the sample that showed rotation < 7 days at ages > 8 Gyr that we suspect are inconsistent with single star evolution.Panel (c) of Figure 1 shows the rotation periods and ages for our full sample of open clusters and asteroseismic field stars.

METHODS
We produced model grids for rotational evolution using two stellar evolution codes-Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2010Paxton et al. , 2013Paxton et al. , 2015Paxton et al. , 2018Paxton et al. , 2019) ) and Yale Rotating Stellar Evolution Code (YREC; Pinsonneault et al. 1989;Demarque et al. 2008).The ranges of stellar properties covered by our grid are detailed in Table 2, and we describe the model physics used to generate each grid in the following sections.

MESA Model Grid
We construct our MESA grid with identical input physics to the models used for asteroseismic inference (described in §2.2) in order to avoid biases introduced by the modeling (see Tayar et al. 2020).Our models used initial elemental abundances from Grevesse & Sauval (1998) and an atmospheric temperature structure following an Eddington T (τ) relation with fixed opacity.We smoothly ramp diffusion from fully modeled at M ≤ 1.1 M ⊙ to no diffusion at M ≥ 1.2 M ⊙ .We do not include core or envelope overshoot.We varied the mass M , metallicity [Fe/H], initial Helium abundance Y init , and mixing length parameter α MLT .
We calculated rotational evolution histories (as described in §3.3) for each combination of stellar properties and appended them to our grid.By default, MESA models do not output the necessary stellar parameters to perform rotational evolution, and it was necessary to adapt the outputs included in the grid.The additional parameters we include for each star were the total moment of inertia I tot , the moment of inertia of the convective envelope I env , the photospheric pressure P phot , and the  2017), all other rotation periods taken from Hall et al. (2021).T eff and [Fe/H] are adopted from the LEGACY (Lund et al. 2017;Silva Aguirre et al. 2015) and KAGES (Aguirre et al. 2015;Davies et al. 2016) catalogs (see Hall et al. 2021).Other stellar properties were derived for this work using asteroseismic mode fitting.
Table 2. Parameter boundaries of the MESA and YREC grids.
convective overturn timescale τ cz .We define τ cz as where H P is the pressure scale height at the convective zone boundary and v conv is the convective velocity one pressure scale height above the base of the convective zone.Stellar interiors in MESA models are divided into shells and the parameters are evaluated at a finite number of points.We identified the precise location of the base of the convective zone as a function of the star's mass fraction using the Schwarzschild criterion, and then interpolated between the values calculated at each shell boundary to more precisely identify the values of our desired parameters at each time step.

YREC Model Grid
We construct our YREC grid following the settings laid out in van Saders & Pinsonneault (2013) and Metcalfe et al. (2020).We use the mixing length theory of convection (Vitense 1953;Cox & Giuli 1968) with the 2006 OPAL equation of state (Rogers et al. 1996;Rogers & Nayfonov 2002).Abundances were taken from Grevesse & Sauval (1998) and opacities from the Opacity Project (Mendoza et al. 2007).We define atmosphere and boundary conditions from Kurucz (1997).Nuclear reaction rates were drawn from Adelberger et al. (2011).Y init was fixed to a linear Helium-enrichment law anchored to the Sun with a slope of dY dZ ⊙ = 1.296 (see §5.4).We varied the same parameters as we did for the MESA grid, with the exception of Y init .
As with the MESA grid, we trace additional parameters to evaluate the angular momentum loss law.For each model at each timestep, we calculate the moment of inertia of both the star and its convective envelope, the photospheric pressure, and the convective overturn timescale.

Magnetic Braking Model
Prescriptions for magnetic braking often incorporate the dimensionless Rossby Number (Ro), defined as the ratio between the rotation period, P , and convective overturn timescale within the stellar envelope, τ cz , as a means to estimate magnetism across stars of different masses.We use the Rossby number in our rotation model due to its utility as a tracer for both the mass and composition dependence of spindown and magnetic field strength.We invoke a Rossby threshold, Ro crit , beyond which point stars depart from a simple power law spindown and conserve angular momentum (van Saders et al. 2016).We adopt the Matt et al. (2012) modification to the Kawaler (1988) braking law.We assume, as in van Saders & Pinsonneault (2013), that the magnetic field strength B scales as B ∝ P 1/2 phot Ro −1 , where P phot is the photospheric pressure, and that mass loss Ṁ scales as Ṁ ∝ L X ∝ L bol Ro −2 , where L X is the x-ray luminosity and L bol is the bolometric luminosity.
Our full model for rotational evolution is described by where Ro is defined as Ro = P τ cz , f K is the scaling factor for the strength of angular momentum loss during classical spindown, ω sat is the threshold at which angular momentum loss saturates for young stars, and with .
The term c(ω) is the centrifugal correction from Matt et al. (2012), and we assume c(ω) = 1, which is appropriate for slowly rotating stars.
To calculate the rotation histories for our grid, we take the outputs of non-rotating MESA and YREC models, and compute rotation periods with the rotevol code (van Saders & Pinsonneault 2013; Somers et al. 2017).We focus only on f K and Ro crit as they will be the most dominant parameters of a WMB law for the stars in our sample, which are old enough to have converged onto tight rotation sequences (Epstein & Pinsonneault 2013;Gallet & Bouvier 2015).We assume a disk locking period of 8.13 days and disk lifetime of 0.28 Myr, setting the initial rotation rates of our models (van Saders & Pinsonneault 2013).We fix ω sat to 3.863 ×10 −5 rad/s.Each of these parameters will be important at early (< 100 Myr) times, but will have negligible effects by the time stars reach the ages in our sample.We assume solid body rotation in our models, since the epoch of radial differential rotation in this mass range is again limited to young stars (Denissenkov et al. 2010;Gallet & Bouvier 2015;Spada & Lanzafame 2020).

Model Grid Emulator
With rotationally evolved model grids, we construct an emulator for rapid stellar evolution modeling.The general approach to this type of optimization problem is simple interpolation between tracks in a high-dimensional model grid (e.g.Berger et al. 2020).However, due to the size of the grid, number of parameters (4-5 per star and cluster, with 2 additional global braking law parameters), and large sample of potential targets, this approach becomes computationally expensive, particularly in the application of Bayesian inference through sampling the model.We therefore opt to train an artificial neural network (ANN) to map the stellar parameters of the grid to observable parameters of stars in our sample.
We define our MESA ANN with seven input parameters and four output parameters.Our inputs represent fundamental stellar properties: age, mass, metallicity, initial Helium abundance, mixing length parameter, braking law strength, and critical Rossby number.The ANN outputs are observable quantities: effective temperature, radius, surface metallicity, and rotation period.The YREC ANN has the above input parameters with Y init excluded, and identical output parameters.The remainder of this section describes the training and characterization of the MESA ANN.The process for training the YREC ANN is identical, and we compare the results when using different grids in §5.4.
Our model structure results in a neural network that acts as a stellar evolution emulator.Given some set of input stellar properties, the model will output the corresponding observable quantities.Because the emulation is rapid, the model can also be used to calculate likelihoods to infer input parameters-given some set of observed properties, we can sample prior distributions for the underlying stellar properties and retrieve posterior distributions, providing estimates for these values with uncertainties.
We construct an ANN with 6 hidden layers comprised of 128 neurons each (following the tuning process of Lyttle et al. 2021).Each hidden layer used an Exponential Linear Unit (ELU) activation function.Using TensorFlow (Abadi et al. 2016), we trained the model on an NVidia Tesla V100 graphics processing unit (GPU) for 10,000 epochs using an Adam optimizer (Kingma & Ba 2017) with a learning rate of 10 −5 .We trained the ANN in ∼8,000 batches of ∼16,000 points.The full model architecture is detailed in Appendix A.
Prior to training the ANN, we remove the pre-main sequence from the tracks in our grid, defined as the threshold at which the luminosity from nuclear burning exceeds 99% of the total stellar luminosity.We allow the tracks to begin evolving across the subgiant branch, as our sample includes stars at or approaching this evolutionary stage, but remove tracks that exceed a rotation period of 150 days.
In order to ensure that the mapping performed by the neural network does not introduce significant uncertainty to the inferred parameters, we divide the grid data into a training set and a validation set.The training set is composed of 80% of the models in the grid, drawn at random, and is used to generate the connections between the input model parameters and observed stellar properties.The remaining 20% of the grid is then used as a validation set to predict the observed parameters based on the provided input parameters, allowing us to characterize the neural network's ability to successfully predict well understood values.When compared to the measurement uncertainties associated with these parameters, the error introduced by the ANN is negligible, with typical fractional uncertainties of ∼10 −3 in the recovery of our validation set (see Figure 2).We also find negligible systematic offset for parameters in our validation set, indicating that the ANN is not introducing significant bias.

Statistical Modelling
In order to efficiently optimize the braking law model parameters, we construct a hierarchical Bayesian model (HBM).The application of a similar HBM for constraining the distribution of Y init and α MLT has been demonstrated by Lyttle et al. (2021).We begin the construction of our model with Bayes' theorem-the posterior probability of our model parameters θ i given some set of observed data d i is where p(θ i ) is the prior on the model parameter θ i (for i parameters) and p(d i |θ i ) is the likelihood of the data given the model.We use our trained ANN to sample the prior distribution p(θ i ) for each parameter and evaluate an instance of the model µ i = λ i (θ i ), where λ i represents the ANN model.From this, we can represent the likelihood of each observation d i with uncertainty σ i given the model evaluation µ i as the normal distribution given N observed variables.
The hierarchical structure of our model allows us to prescribe various levels of pooling to different parameters.The WMB model parameters f K and Ro crit , for example, are assumed to be the same for all stars in our sample.For the ANNs trained on both the the MESA and YREC grids, we define the prior for Ro crit as Ro crit ∼ U(1.0, 4.5) and the prior for f K as f K ∼ U(4.0, 11.0) where θ ∼ X represents a parameter θ being randomly drawn from a distribution X, and U(a, b) is a uniform distribution bounded between a and b.The values of Ro crit and f K drawn from these uniform distributions are used to calculate the full set of model evaluations µ i for that step.The bounds for Ro crit and f K were centered near the solar Rossby number derived for our grids (for MESA: Ro ⊙ ≈ 2.05, f K,⊙ ≈ 5.89; for YREC: Ro ⊙ ≈ 2.33, f K,⊙ ≈ 7.52).
Other parameters are assumed to be unique to each star.For the YREC ANN, we constrain the mass, metallicity, mixing length parameter, and age.We constrain the same parameters for the MESA ANN with the addition of the initial Helium abundance.The prior distributions for these parameters are defined as truncated normal distributions, given by where N is the normal distribution, a and b are the lower and upper bounds, respectively, µ is the median and σ is the standard deviation.Here, µ and σ are taken from the observational constraints on the parameters and their uncertainties.For stars in clusters, we define a prior centered on the value reported in the corresponding reference (see §2.1) with a width set to the measurement uncertainty for age, metallicity, mixing length parameter, and rotation period (with the inclusion of Y init for the MESA grid).For the masses of cluster stars, we use a homology scaling relationship with T eff and set a broad prior (σ M =0.25 M ⊙ ), and for the mixing length parameter and initial Helium abundance we use uniform priors.For asteroseismic stars in our sample, all of the above properties are constrained by the asteroseismic fitting, and we use this asteroseismic value and its uncertainty as the center and width of the prior distributions, respectively.Our truncated distributions for all stars are bounded by the grid limits described in Table 2.
Finally, we include a third class of prior distributions in our model which are shared by some stars but not all.Each star within the same cluster is assumed to have the same age, metallicity, and initial Helium abundance, while these parameters should be fully independent for each target in the asteroseismic sample and for the Sun.These prior distributions share the same truncated normal form as the independent parameters, but can be selectively applied to specific subsets of the data.
With our priors and likelihoods defined, we sampled the model parameters.The ANN is compatible with automatic differentiation, allowing us to utilize No-U-Turn Sampling (NUTS; Hoffman & Gelman 2014).We constructed a probabilistic model with PyMC3 (Salvatier et al. 2016), then calculated the maximum a posteriori estimate as our starting point and sampled 4 chains for 5,000 draws with 1,000 tuning steps.We sampled chains long enough to ensure that the Gelman-Rubin R statistic (Gelman & Rubin 1992) was lower than 1.01 for all parameters indicating model convergence.The residuals from our fit, as well as an example of our model fit to the Sun, are shown in Appendix B.

RESULTS
We optimize the parameters of our model under two different assumptions-standard spindown and WMB.In the standard spindown framework, we assume stars follow a Skumanich-like angular momentum loss law, where J ∝ ω 3 at late times.Under the WMB assumption, stars lose angular momentum to magnetized stellar winds with the same relation as the standard spindown law until they reach a critical Rossby number Ro crit , at which point angular momentum is conserved.We use the MESA ANN as our primary emulator as its grid physics match the models used in the asteroseismic parameter estimates.In the standard spindown case, we only optimize for f K , and retrieve a constraint of f K = 6.11 ± 0.73.For the WMB model, we report f K = 5.46 ± 0.51 and Ro crit /Ro ⊙ = 0.91 ± 0.03.
Figure 3 shows the distribution of rotation periods predicted by our WMB model.We have divided the sample into equal-size bins in T eff because temperature captures the effects of both a star's mass and metallicity on its rotational evolution.The red shaded regions show the density of stars drawn from a simulated population of 1,000,000 stars under the best-fit WMB assumptions, generated with stellar properties drawn from uniform distributions for each parameter bounded by the edges of our sample using our MESA emulator.The width of the distribution is caused by the range of masses, metallicities, Helium abundances, and mixing length parameters within each T eff bin.Stars in clusters can be seen as groups with discrete, well-constrained ages below 2.5 Gyr, and are valuable calibrators for the early angular momentum loss J.In our model, this early J is captured by the braking law strength parameter, f K .Stars in our asteroseismic sample span a wide range of ages, particularly on the second half of the main sequence, and provide the constraint on Ro crit .In Figure 4, we show the comparison between the rotation periods predicted by both the standard spindown and WMB models (in blue and red, respectively).Each shaded region represents the density of points in a population of 100,000 simulated stars from our MESA emulator.The standard spindown model was fit to the full sample, without altering angular momentum loss beyond a Rossby threshold.The models produce similar constraints on f K , as the early rate of J is well-constrained by the clusters in both models.At older ages, the standard spindown model significantly overpredicts the rotation periods of stars in our asteroseismic sample.
The WMB model results in a smaller average deviation from the observed rotation periods.Figure 5 shows the the difference between predicted and observed rotation periods for our sample.The colored points show the uncertainty-weighted median within a 0.2 t/t MS bin.On average, the standard spindown model overpredicts rotation periods by 0.72 days for the full sample and 6.00 days for stars beyond the first half of the main sequence (t/t MS ≥ 0.5).Conversely, WMB underpredicts rotation periods by 0.31 days for the full sample and 3.18 days for stars past 0.5t/t MS .Isolating only the asteroseismic sample (at all ages), standard spindown overpredicts P rot by 4.66 days on average, and WMB underpredicts by 2.02 days.The corresponding fractional deviations for the asteroseismic sample are +17.73%for standard spindown and −9.09% for WMB.We perform a reduced chi-squared test to determine the goodness-of-fit for our models, and we find χ 2 ν,WMB = 1.07 and χ 2 ν,standard = 14.02.Because χ 2 ν,WMB ≪ χ 2 ν,standard , we conclude that the WMB model provides a better fit to the data.
Figure 5 shows the difference between predicted and observed rotation periods as a function of fraction of main sequence lifetime.For the first half of the main sequence, the standard spindown and WMB models both describe the observed rotation periods well.However, at roughly halfway though the main sequence (0.5t/t MS ), the standard spindown model deviates from the observed distribution and begins overpredicting rotation periods.Both models are consistent with the cluster data, which follow a tight spindown sequence that is nearly identical for the two models (see Figure 4).

DISCUSSION
We have provided refined probabilistic estimates for the onset of WMB, described by the parameter Ro crit .Our model indicates that stars enter a phase of weakened braking before reaching the Rossby number of the Sun (Ro crit = 0.91±0.03Ro ⊙ ).This result supports constraints by David et al. (2022), which found a sub-solar Ro crit when examining the pileup in the temperature-period distribution of Kepler stars.van Saders et al. (2016van Saders et al. ( , 2019) ) found that a critical Rossby number of Ro crit ≈ Ro ⊙ provided the best fit to the observed rotation periods, which agrees with our results within 2σ. .Difference between predicted and observed rotation periods for all stars in our sample (shown as gray points) as a function of fraction of main sequence lifetime t/t MS .The blue and red points represent the uncertainty-weighted median of δP within a 0.2 t/t MS bin for the standard and WMB models, respectively.Main sequence lifetime was estimated by identifying the age of core-H exhaustion in MESA models generated for each star.Roughly halfway through the main sequence lifetime, the standard spindown model begins significantly over-predicting rotation periods.The WMB model is consistent with the observed distribution until near the end of the main sequence, at which point it underestimates rotation periods.
The new constraints on weakened braking parameters provided here can be used as guidelines for where gyrochronology is likely to be accurate.Beyond Ro crit , rotation evolves only slowly with the changing moment of inertia, and stars can be observed with the same rotation period for Gyr timescales, challenging any gyrochronological estimate.We show that gyrochonological ages should be precise until ∼Ro ⊙ , corresponding to an age of ∼4 Gyr for sun-like stars.After the onset of WMB, age estimates should have significantly larger uncertainties due to the slowly evolving rotation on the second half of the main sequence.

WMB Model Performance
Towards the end of the main sequence, our model for weakened braking begins to underestimate rotation periods.This likely reflects our overly simple implementation of the transition from standard to weakened braking.The immediate shutdown of angular momentum loss beyond Ro crit is the simplest model which introduces the fewest new parameters.Given the limited sample of reliable calibrators spanning a wide range of T eff near the onset of WMB, any parameterization of a possible gradual transition, or a transition that does not completely shut down magnetic braking, is not well constrained.As more seismic constraints are placed on the ages and rotation periods near Ro crit , additional parameters that lead to a gradual transition, or J ̸ = 0 beyond the transition, can be tested.
The deviation between the WMB model and observed rotation periods could additionally be partially explained by small deviations in inferred model ages.At the end of a star's main sequence lifetime, even in the WMB framework when angular momentum is conserved, the rotation period increases steeply due to the changing stellar moment of inertia as the star's radius expands.Models for rotation increase on short time spans in parallel vertical tracks in rotation-age space as stars traverse the subgiant branch, with small separations between stars of different T eff .Improved asteroseismic modeling, or a larger sample of stars with asteroseismic parameter constraints, could better distinguish between these effects at the end of the main sequence.

Assessing the Asteroseismic Constraint
To illustrate the impact of the asteroseismic sample on our ability to constrain Ro crit , we fit our model to two subsets of the data: one comprised of only clusters and the Sun, capturing the early rotational evolution, and one that adds the asteroseismic stars.Figure 6 shows a Kernel Density Estimate (KDE) of the sampled marginal posterior distributions for Ro crit when fit to each of these samples.When fit to only clusters and the Sun, Ro crit has little to no likelihood below the solar value, and is unconstrained beyond the solar value.This aligns with our expectations, as the young cluster sample has repeatedly been shown to follow standard braking (Barnes 2007(Barnes , 2010;;Mamajek & Hillenbrand 2008;Gallet & Bouvier 2015;Meibom et al. 2011Meibom et al. , 2015)).When the asteroseismic sample is included, the posterior becomes tightly constrained near the solar value.This exercise clearly demonstrates why the effects of WMB were not identified until a large enough sample of stars with precise rotation periods and ages spanning the main sequence were available.

Consistency with Solar Twins
A recent study by Lorenzo-Oliveira et al. (2019) proposed tension between the weakened magnetic braking model and an observed population of "solar twins."The stars in this sample have typical masses within ±0.05 M ⊙ of solar and metallicities with ±0.04 dex of solar.Rotation periods were not directly measured for the majority of stars in this sample, instead the projected rotational velocity v sin i of each star was estimated from spectral line broadening.This was converted to a projected rotation period, P rot / sin i, using stellar properties derived from from Gaia DR2 (Gaia Collaboration et al. 2018) and ground-based spectroscopic data.
If a system is observed directly edge on (i = 90 • ), the projected rotation period will match that measured from photometric spot modulation or asteroseismic mode splitting.The primary effect of rotation axis inclination away from 90 • is to shift the projected rotation period to a higher value (see panel (a) of Figure 7).Lorenzo-Oliveira et al. (2019) undergo a selection process of simulating projected rotation periods given some random orientation between 0 and 90 • , comparing their measured population against these simulations, and reducing their sample to stars they found most likely to be seen edge on based on the agreement (see  description of their approach).As only a fraction of the observed sample is likely to be observed directly edge on, the fastest rotation periods in the solar twins sample represent a lower envelope to the true distribution of rotation periods of the sample.
We test the standard spindown and WMB models the solar twins sample, seen in panels (b) and (c) of Figure 7.We calculate P rot / sin i for our MESA emulator model tracks, drawing inclinations randomly from a uniform distribution between 0 and 1 in cos i.The stellar properties of our model grid were drawn from uniform distributions bounded by the parameter cuts described in Lorenzo-Oliveira et al. ( 2019)-mass and metallicity were bounded by 0.8 M ⊙ ≤ M ≤ 1.2 M ⊙ and −0.04 ≤ [Fe/H] ≤ +0.04, and unconstrained parameters were given broad uniform priors (0.22 ≤ Y init ≤ 0.28, 1.4 ≤ α MLT ≤ 2.0).We note that fixing Y init and α MLT to solar-calibrated values has negligible impact on the model fit.We find that the standard spindown model overpredicts projected rotation periods beyond the age of the Sun.The WMB model predicts the observed population with minor deviations from entirely edge-on inclinations.We find that the WMB model reasonably reproduces the behavior observed in the solar twins, and does so better than the standard spindown model.

Accounting for Grid Bias
We test our model fit using neural networks trained on grids of models generated by two stellar evolution codes, MESA and YREC.This provides an opportunity to independently validate our results as well as test for any bias introduced by the choice of grid.To date, most investigations of WMB have used ages and rotational evolution that were inferred using reasonable, but different, underlying stellar evolution models.Our MESA grid was constructed with input physics matching the asteroseismic modeling, avoiding the cross-grid bias when fitting the MESA-trained neural network to the asteroseismic observations.While we have matched the physics in the seismic and rotational models, we have not performed the fits simultaneously, which we reserve for future work.The primary difference between the construction of the grids was to vary Y init as an additional dimension of the MESA grid, while calculating it with a fixed He-enrichment law in the YREC grid.We used a relation to compute the value of Y init for a model in the YREC grid given its metallicity [Fe/H] given by Fe/H] + 1 where Y P is the primordial Helium abundance, the slope of the Helium enrichment law that matches the solar value is dY dZ ⊙ = 1.296, and the solar metal fraction is Z X ⊙ = 0.02289 (Grevesse & Sauval 1998).
The ANN for the YREC grid was trained identically to the process for the MESA grid described in §3.4, and we constructed the probabilistic model following the process described in §3.5.For the YREC ANN, the value of Y init fit by our asteroseismic modeling with MESA was not used as a constraint on the model likelihood, while it was for the MESA ANN.The choice to include Y init as a free parameter, as well as the differences between how different stellar evolution codes calculate quantities used in our modeling, have the potential to introduce systematic biases in the resulting model fits.Here, we compare between the results inferred by emulators trained on different model grids.
Most braking laws include a strong Ro dependence, and thus a dependence on the convective overturn timescale τ cz , and there is no single agreed upon means of calculating this value (see Kim & Demarque 1996).Furthermore, changes in grid physics can result in different values of τ cz , even in solar-calibrated models.To account for this, we normalize Ro crit by a grid-dependent solar Rossby number Ro ⊙ .To calculate Ro ⊙ for each grid, we produced solar-calibrated stellar evolution tracks and compute the Rossby number at the age of the Sun.For each model grid, we also compute the value of f K that reproduced solar rotation at solar age under the standard spindown assumption, and apply this as a normalization factor when comparing the inferred values of f K in our WMB models.We notate this solar-normalized braking law strength as f ′ K .These normalization factors allow us to compare directly between the braking law parameters inferred from the ANN trained on each model grid.
The left panel of Figure 8 shows the marginal and joint posterior distributions for the braking law parameters when fit with the MESA and YREC ANNs.The black dashed line shows the solar Rossby number, Ro ⊙ .Both MESA and YREC return values of Ro crit below Ro ⊙ , indicating that the onset of WMB occurs before the age of the Sun for a solar analog.The inferred braking law parameters have slight offsets, but agree within 1σ.To assess the impact of leaving the Y init parameter free, we also performed probabilistic modeling with the MESA ANN with Y init set to the He-enrichment law described above.We show the updated posterior distributions for this fit compared to the YREC ANN in the right panel of Figure 8.Using the YREC emulator model, we retrieve constraints on the braking law parameters of f ′ K = 0.86 ± 0.07 and Ro crit = 0.94 ± 0.04.When Y init is left as a free parameter, the MESA emulator model returns f ′ K = 0.77 ± 0.07 and Ro crit = 0.91 ± 0.03.When we fix the Helium enrichment law to that used in the YREC grid, the MESA emulator model reports f ′ K = 0.80 ± 0.07 and Ro crit = 0.94 ± 0.03.We note that all models consistently return a value of Ro crit below the solar Rossby number.
When holding Y init fixed to the YREC He-enrichment law, we find closer agreement between the braking law parameters inferred by our model fitting, with Ro crit in near-perfect agreement.This implies that Y init provides additional constraints on the braking law parameters, and its inclusion as a grid dimension can influence the result.Y init is a challenging property to measure for sun-like stars, and yet affects our inferred value of Ro crit at the ∼1σ-level.We conclude that uncertainty in the Helium enrichment law should be treated as a systematic uncertainty in the inference of Ro crit .

Future Applications
In this study, we focus only on f K and Ro crit due to the age distribution of our sample.In the future, the same approach described here could be applied to a sample of targets which span earlier phases of evolution (i.e.young open clusters), at which time braking law assumptions, such as the disk-locking timescale, disk lifetime, ω sat , and internal angular momentum transport must be treated more carefully.
We limited the range of our input model grid to cover the parameters of our sample in order to reduce the computational time required for model generation and neural network training.The framework for the ANN emulator could easily be applied to a grid spanning a wider range of stellar properties, and would provide a useful tool for quickly evaluating stellar evolution tracks or simulating stellar populations.To reduce training time, the grid resolution could be selectively increased to reach a precision threshold.Scutt et al. (2023) suggested that parameter spacing can be modified in different regions of the grid to improve ANN precision.
Asteroseismic pulsation frequencies are often generated alongside stellar models using tools such as GYRE (Townsend & Teitler 2013).These pulsation frequencies, particularly the large frequency spacing (∆ν), can be included in the grid dimensions (e.g.Lyttle et al. 2021) and applied as further likelihood constraints for models.Ideally, some combination of the above additions could be implemented to produce a broadly applicable stellar evolution emulator that does not require generating or interpolating large model grids.

CONCLUSIONS
In summary, our primary conclusions are: 1. We present evidence for weakened magnetic braking in old stars.Using a neural network as a stellar evolution emulator, we perform probabilistic modeling to produce posterior distributions for the parameters of the weakened braking model.We find that the weakened braking model provides the best fit to the observed distribution of rotation periods.
2. We show that the most likely weakened braking scenario diverges from standard spindown at a slightly earlier evolutionary phase than the Sun (Ro crit /Ro ⊙ = 0.91 ± 0.03).We caution that our WMB model is a simplified case in which angular momentum loss is fully switched off at a critical Rossby number, and likely does not fully capture the time evolution of the stellar dynamo.The relatively sparse calibrator sample near Ro crit means that it remains challenging to infer the precise onset of WMB relative to the Sun's evolution.
3. Our method for emulating stellar evolution with a neural network enables rapid evaluation of stellar models, making it possible to fit braking law parameters while properly accounting for the uncertainties in the stellar parameters of our calibrator sample.By modifying the braking law used to generate our training set, we could test other effects at early times, such as the impact of internal angular momentum transport or disk-locking.
4. We report mild disagreement between the constraints on WMB parameters when using different underlying model grids.This indicates that the choice of grid physics and which parameters are varied in the model can impact the inferred model parameters.For our choices, the impact is at the 1σ level.
5. The WMB model appears compatible with the solar twins sample.The standard spindown model predicts slower rotation than observed in the solar twins stars during the second half of the main sequence, while their rotation periods can be described by the WMB model with modest deviations from a fully edge-on population.
6. Our constraint on the Ro crit at which stars enter a phase of weakened braking suggests that gyrochronology faces challenges when estimating stellar ages for much of the main sequence lifetime.For sun-like stars, gyrochronological age estimates are likely unreliable beyond an age of ∼4 Gyr.For more massive stars (≳ 1.1M ⊙ ), gyrochronology relations appear break down even earlier, at an age of ∼2.5 Gyr.Even after a star has entered the weakened braking phase, a reasonable range for its age can be estimated from its rotation period, and our constraint on Ro crit enables gyrochronological modeling that will provide a realistic uncertainty on the stellar age.
The growing population of stars with precisely measured ages and rotation periods from asteroseismology is shedding essential light on the evolution of stellar rotation.Improved direct observations of magnetic field strength can add additional constraints on the braking law parameters.As more stars are added to this sample, the transition to WMB can be constrained to higher precision.

APPENDIX
A. NEURAL NETWORK STRUCTURE As described in §3.4,our artificial neural network was generated with six hidden layers of 128 neurons.A model summary can be found in Table 3

B. MODEL VALIDATION
To validate the performance of our model, we calculated the difference between the observed value and the median of the posterior sampled distribution for each parameter in our grid.Figure 9 shows this value, where δX = X predicted − X observed for a parameter X.We find good overall agreement between predicted and observed values, with no significant systematic offsets.
For each star in our sample, we retrieve full posterior distributions for each parameter.In Figure 10, we show the sampled marginal and joint posterior distributions for the Sun.

Figure 1 .
Figure 1.(a) Hertzsprung-Russell diagram showing our sample of calibrators in open clusters.Model tracks generated by our emulator are shown as gray lines.(b) Hertzsprung-Russell diagram of our asteroseismic sample from Hall21.We derived the stellar properties shown here with asteroseismic modeling.(c) Observed rotation period plotted as a function of stellar age.We color points by their effective temperature.Asteroseismic stars are shown as circles and open cluster members are marked by triangles.

Figure 2 .
Figure 2. Uncertainty introduced by the MESA ANN emulator.The histograms for P, R, and T eff show the (predicted−truth)/truth value for our training set, and the bottom right panel shows predicted−truth for the surface metallicity to account for points where [Fe/H] surface,truth ≈ 0. The median µ and standard deviation σ of these distributions are shown in the top right corner for each parameter, and µ is marked by the solid vertical line.The error incurred by the ANN is negligible compared to the uncertainty on the observed values.

Figure 3 .
Figure3.Stellar rotation period versus age, shown in three bins each spanning 300 K in T eff .Asteroseismic measurements and cluster stars are shown by points-black points represent rotation periods with fractional uncertainties σ P /P ≤ 25% and gray points show σ P /P > 25%.The Sun is marked by the ⊙ symbol.Red contours represent the distribution of rotation periods within a given T eff bin predicted by our MESA emulator model, produced from a sample of one million emulated stars with stellar properties randomly drawn from uniform distributions bounded by our sample, and f K and Ro crit fixed to the median values of the posterior distributions.

Figure 4 .
Figure 4. Same as Figure 3, with the additional comparison to the standard spindown model.The contours represent the distribution of predicted rotation periods within a given T eff bin, with red showing our WMB model and blue showing a standard Skumanich-like spindown model, both generated with our MESA emulator.The value of f K for the standard spindown model is the median value from the posterior of a fit to our full sample with no Ro crit constraint.
Figure5.Difference between predicted and observed rotation periods for all stars in our sample (shown as gray points) as a function of fraction of main sequence lifetime t/t MS .The blue and red points represent the uncertainty-weighted median of δP within a 0.2 t/t MS bin for the standard and WMB models, respectively.Main sequence lifetime was estimated by identifying the age of core-H exhaustion in MESA models generated for each star.Roughly halfway through the main sequence lifetime, the standard spindown model begins significantly over-predicting rotation periods.The WMB model is consistent with the observed distribution until near the end of the main sequence, at which point it underestimates rotation periods.

Figure 6 .
Figure 6.Comparison between posterior distributions for Ro crit from models fit to different subsets of the data.When fit to only clusters and the Sun (shown in red), Ro crit is unconstrained beyond the solar Rossby number.With the inclusion of the asteroseismic sample (shown in blue), Ro crit is tightly constrained just below the solar Rossby number.The y-axis has been arbitrarily scaled for clarity.

Figure 7 .
Figure 7. (a) Projected rotation period, P rot / sin i, of the solar twins sample versus age.The colored lines show tracks from our MESA emulator for a solar-calibrated model with a range of stellar inclinations, evolved under WMB assumptions.Models that are not observed edge on have their projected periods shifted to higher values.(b) The solar twins sample compared to a standard spindown model with a range of stellar inclinations.We generated a population of 1,000,000 stars with parameters drawn from uniform distributions within ±0.05 M ⊙ of solar for M, ±0.04 dex of solar for [Fe/H], and inclinations, i, drawn from a uniform distribution in cos i. Y init and α MLT were drawn from uniform distributions covering our model grid.(c) Same as panel (b), but with the WMB model.f K and Ro crit were fixed to values fit to our full sample.The standard model overpredicts rotation periods of the solar twins sample beyond the age of the Sun, while they are consistent with WMB when accounting for inclinations.

Figure 8 .
Figure 8.(a) Corner plot showing the marginal and joint posterior distributions for the global parameters of our WMB model.Blue shows the samples from the fit using a neural net trained on a grid of MESA models, and red shows the samples from a fit using the YREC-trained neural net.The solar Rossby number Ro ⊙ is shown as a dashed black line.The median values of each distribution are shown as dashed lines in their respective colors in the top and right panels.(b) The same posterior distributions, now with the He-enrichment law in the MESA probabilistic model fixed to the relation used when generating the YREC grid.The primary difference between the grids used to train the emulator models is the varied Helium abundance Y init in the MESA grid.When fixed to the YREC enrichment law, the constraints on WMB global parameters are in closer agreement.

Figure 9 .
Figure 9. Difference between observed values of stellar properties and predictions from our probabilistic model (predicted−observed), plotted against stellar age.The red line shows the median of the difference, with the standard deviation shown by the shaded region.

Figure 10 .
Figure 10.Corner plot showing marginal and joint posterior distributions for the Sun from our sampled probabilistic model.
§2 of Lorenzo-Oliveira et al. 2019 for a full .

Table 3 .
Model summary for our ANN.