MESA models with magnetic braking

Two magnetic braking models are implemented in MESA for use in the MIST stellar model grids. Stars less than about 1.3 $M_{\odot}$ are observed to spin down over time through interaction with their magnetized stellar winds (i.e., magnetic braking). This is the basis for gyrochronology, and fundamental to the evolution of lower mass stars. The detailed physics behind magnetic braking are uncertain, as are 1D stellar evolution models. Thus, we calibrate our models and compare to data from open clusters. Each braking model tested here is capable of reproducing the data, albeit with some important distinctions. The Matt et al. (2015) prescription matches the slowly rotating stars observed in open clusters, but tends to overestimate the presence of rapidly rotating stars. The Garraffo et al. (2018) prescription often produces too much angular momentum loss to accurately match the observed slow sequence for lower mass stars, but reproduces the bimodal nature of slow and rapidly rotating stars observed in open clusters fairly well. We find additional evidence that some level of mass dependency may be missing in these braking models to match the rotation periods observed in clusters older than 1 Gyr better.


INTRODUCTION
Stellar rotation rate is one of the fundamental properties of stars, determining all other properties of a star throughout their evolution, alongside stellar mass and metallicity. Rotation introduces centrifugal forces that strongly affect the structure of a star, altering temperatures and luminosities (an effect called gravity darkening), as well as inducing mechanical instabilities that induce mixing of stellar material throughout the star (von Zeipel 1924;Endal & Sofia 1976;Maeder 2009). Due to the immense challenges of full, 3D stellar modeling, we rely on simpler 1D models. The effects of stellar rotation have been studied for over a century at least. Models exploring stellar rotation in 1D have been calculated over the last several decades, e.g., Pinsonneault et al. (1989); ; ; Palacios et al. (2003); Charbonnel & Talon (2005); Denissenkov & Pinsonneault (2007). Numerous non-rotating grids (e.g., Baraffe et al. 1998;Pietrinferni et al. 2004;Demarque et al. 2004;Dotter et al. 2008;Bressan et al. 2012;Choi et al. 2016) with detailed physics and sometimes having their own specializations, e.g., as with the models of Somers & Pinsonneault (2015a); Somers et al. (2020) which models the structural effects magnetic fields and starspots, or e.g., binary evolution as in Eldridge et al. (2017). However, large 1D stellar model grids and isochrones that incorporate rotation effects over large ranges of mass, metallicity, and rotation rates, have only recently become available. This is not an exhaustive list, but see e.g., Ekström et al. 2012;Georgy et al. 2014;Gossage et al. 2018. Until recently, grids that fully model 1D stellar rotation have been limited by a lack of models that describe the spin evolution of low mass stars (i.e., with about M < 1.3M ) in this context. Implementing rotation in low mass 1D stellar models adds an additional layer of complexity. There are even fewer large grids of stellar models that incorporate the effects of rotation at low masses, but Amard et al. (2019) has produced some of the first. At masses below about 1.3M , stars develop a convective envelope that increases in depth from the surface of the star until masses ≤ 0.3 M , below which stars become fully convective. The convective envelopes of these low mass stars host surface magnetic fields that interact with stellar winds, extracting angular momentum, and causing them to slow down over time. The magnetic coupling between star and stellar wind, leading to a slowing of the stellar rotation rate is commonly called magnetic braking. This phenomenon was inferred early on by Kraft (1967) and Skumanich (1972) in nearby stars, but was theoretically anticipated beforehand as in e.g., Schatzman (1962); Brandt (1966); Weber & Davis (1967).
This characteristic of low mass stellar evolution may be exploited as an age determination technique known as gyrochronology (see e.g., Barnes 2003;Mamajek & Hillenbrand 2008). Thus, the accuracy of this age determination method is partially reliant on an accurate modeling of rotation in low mass stars, including the effects of magnetic braking. Accurately modeling magnetic braking is also crucial to studies of stellar populations 2 Gyr old, when the main sequence turn off (MSTO) is dominated by low mass stars that experience magnetic braking (Georgy et al. 2019;Gossage et al. 2019) and for studying their underlying rotation rate distributions. Furthermore, incorporating magnetic braking into 1D stellar models and applying them to data is a useful experiment to help constrain the detailed magnetohydrodynamic (MHD) simulations that inform them (e.g., Matt et al. 2012aMatt et al. , 2015Réville et al. 2015;Garraffo et al. 2015;Finley & Matt 2018;Garraffo et al. 2018;See et al. 2019). Note also that the models we present in this work self-consistently evolve models with rotation effects, e.g., internal angular momentum transport and rotational mixing processes calculated (like Amard et al. 2019). This is in contrast to other successful stellar models computed for gyrochronology that calculate rotation evolution from the structural properties of nonrotating stellar models (such as, van Saders & Pinsonneault 2013;Gallet & Bouvier 2013, 2015Lanzafame & Spada 2015), and thus do not self-consistently evolve the models with internal angular momentum transport, etc.
Furthermore, classical 1D stellar models can not reproduce the solar lithium abundance, and it has become clear that rotating stellar models, perhaps with additional physics (still under investigation) are required to explain it (Charbonnel & Talon 2005;Eggenberger et al. 2005;Somers & Pinsonneault 2016). The Sun has a greatly depleted surface abundance of lithium (measured at log( 7 Li/H) + 12 = 1.05 ± 0.1 now, in comparison to its assumed primordial value 3.26 ± 0.05 measured from meteorites, e.g., Asplund et al. 2009), and it is unclear why. Observations of solar twins (Carlos et al. 2019) and lithium abundances of stars in open clusters (Sestito & Randich 2005) have shown that stars near solar mass generally show a depletion of surface lithium abundance; although, there is an unexplained dispersion in the abundance at a given age, and the Sun appears to be amongst the most heavily depleted compared to similar stars (Carlos et al. 2019). It has been proposed that rotation-enhanced mixing could play a role in explaining this (e.g., Bouvier 2008;Eggenberger et al. 2010;Somers & Pinsonneault 2016, and see Bouvier 2020 for an overview), and comparison of v sin i and lithium abundance data may support this connection, e.g., Skumanich (1972); Bouvier (2008); Beck et al. (2017). Lithium is not the only light element observed to experience depletion. Beryllium and boron are observed to be depleted as well in F-and B-type stars, simultaneously and in a similar fashion (Charbonnel et al. 1994;Venn et al. 2002). Thus, reproducing light element abundances like lithium can serve as a useful constraint on low mass, rotating stellar models, and help us understand what its cause may be.
Our models are based on the framework of the MIST stellar model grid (Dotter 2016;Choi et al. 2016), which was built to cover a large range of masses (0.1 to 300 M ), evolution from pre-main sequence (PMS) to late evolutionary phases (such as white dwarfs), and a wide range of metallicities. The current iteration of MIST contains only non-rotating models, and models at a single higher rotation rate; although Gossage et al. (2018Gossage et al. ( , 2019 have done work towards expanding the grid to include more rotation rates. All the while, MIST has never included a proper treatment of rotation for stars with surface convective envelopes. In this work, we implement two magnetic braking models: Matt et al. (2015) and Garraffo et al. (2018). We examine their abilities to reproduce observed rotation periods of open cluster stars for a range of ages as implemented here. Our goal in doing this is to provide models available for studies, such as of the type listed above, but also to work towards converging our theories on the physics that drive magnetic braking. In so doing, hopefully we can make our 1D stellar models more accurate in this regime and uncover the physics behind rotation-driven phenomena. We have developed a preliminary grid of stellar models at solar metallicity to test implementations of magnetic braking in our models.
In Sec. 2, we provide the background of our models and some of the fundamental framework behind stellar rotation in MESA. We describe our implementation of the Matt et al. (2015) and Garraffo et al. (2018) braking models in Sec. 3. Our results are presented in Sec. 4, displaying the effects of each braking model in regards to reproducing observed open cluster rotation periods, effects on the Hertzsprung-Russell diagram, and surface lithium depletion. We discuss our results and discrepancies found in Sec. 5. Finally, we conclude our work in Sec. 6.

ROTATION IN MESA
Modeling the effects of stellar rotation on the evolution of a star is difficult. It is at least a 2D effect, but with current computing power, we rely on 1D stellar models. In this work we create our models with MESA r11701. The 1D stellar evolution code MESA is open source, and has been continually developed over the last decade (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019. The current implementation of rotation in MESA is described in Paxton et al. (2013) with updates in Paxton et al. (2019).
These models are an expansion of the MIST stellar model grid (Dotter 2016;Choi et al. 2016), and we adopt the physics mentioned therein. Until now, MIST has used an ad hoc treatment of rotation for stars with convective envelopes, scaling the rotation rate from its full value at stellar masses > 1.8 M towards zero at 1.2 M and below. In addition, models were previously made to rotate only once they reach the zero age main sequence (ZAMS); models are now initialized during the pre-main sequence (PMS) with rotation. The methodology described below has been developed to implement a more realistic treatment of stellar spin down for low mass stars that experience magnetic braking. For this work, we have calculated solar metallicity models, with masses ranging from 0.1 − 1.3 M .

Internal Angular Momentum Transport
MESA uses a diffusive equation (Endal & Sofia 1978;Pinsonneault et al. 1989; to handle angular momentum transport (Paxton et al. 2013). The equations for rotational chemical mixing (omitted here for brevity, but see  for the relevant formulation) are similar. The variation of angular velocity ω with time within the star takes the form with i, ρ, m, r as the specific moment of inertia, mass density, mass, radius, and ν as the turbulent viscosity. The turbulent viscosity is the sum of various diffusion coefficients, setting the strength of diffusive transport via e.g., convection and various hydrodynamical mixing processes from rotation induced instabilities. In addition, a constant, arbitrary diffusion coefficient may be set within MESA, representing some additional background mixing source(s), discussed below. There is some evidence that additional sources of diffusion may be necessary within this formalism, which could come from internal magnetic fields (Eggenberger et al. 2005), or possibly gravity waves (Charbonnel & Talon 2005) (both of which are physical processes that transport angular momentum in a non-diffusive manner, adding complexity, e.g., see Rogers 2015;Rogers & McElwaine 2017), as Denissenkov et al. (2010) and Somers & Pinsonneault (2016) mention. Additional sources of diffusion are required to reproduce the rotation profile of the Sun, which appears to be nearly solid body (Howe 2009). This additional source of diffusion is a free parameter (referred to as ν 0 going forward) that we tune to reproduce the solar rotation profile at the solar age for a 1 M model; this extra source operates throughout the entire star, with ν = ν hydro + ν 0 in Eq. 1 (where ν is the sum of all other diffusion coefficients). Due to the efficient transport of angular momentum in convection zones, solid body rotation develops there, while it generally does not within the radiative zone, allowing differential rotation.
We select a constant ν 0 for all of our models, regardless of mass for simplicity (although Lanzafame & Spada 2015;Somers & Pinsonneault 2016;Spada & Lanzafame 2020 find evidence of mass dependence). We discuss possible consequences of this choice in Sec. 5. We calibrate ν 0 to roughly reproduce the differential rotation observed in the Sun, where the core rotates 10% faster than the outer layers (e.g., Couvidat et al. 2003). We additionally use the spread in observed P rot for ∼ 1 M stars in open clusters ranging from 2 to 2500 Myr to calibrate ν 0 . Lower ν 0 produces weaker core-envelope coupling, and stronger differential rotation. The degree of differential rotation scales with rotation rate, and stronger differential rotation decreases angular momentum transport efficiency between the core and envelope, which hinders the efficiency of angular momentum loss. Thus, rapidly rotating stars spin down less efficiently than slow rotators, causing a spread in the predicted P rot of our models over time. Decreasing ν 0 causes stronger differential rotation, increasing this spread, which does not match open cluster data (e.g., see 2). In this way, cluster data at late ages (> 1 Gyr) provides a soft lower limit on the value of ν 0 in our models. We adopt a default, solar-calibrated value ν 0 = 2 × 10 4 cm 2 s −1 to satisfy both constraints.

BRAKING FORMALISMS
An early model for calculating the angular momentum lost by the Sun due to solar wind interactions came from Weber & Davis (1967). It was observed later in Kraft (1967) and Skumanich (1972) that stars spin down as they age. The equatorial surface velocity v eq appeared to be proportional to t −1/2 (with t as stellar age); this relation is now often referred to as the Skumanich law. Generally, more slowly rotating low mass stars in open clusters are observed to follow the Skumanich law, although more rapid rotators appear not to. Research has sought a model to describe the physical mechanisms that may cause stars to lose angular momentum as they age, with some early generalizations of the Weber & Davis (1967) model coming from Kawaler (1988) and Krishnamurthi et al. (1997).

Initial spin evolution
In understanding the physics of magnetic braking during PMS evolution, there are obfuscating, yet interlinked physical processes that must be addressed to properly model the rotation evolution of stars (Bouvier et al. 2014 gives an overview). One such process is disk locking. Disk locking is an observationally inferred process (Koenigl 1991;Shu et al. 1994) of the rotational evolution of these stars on the PMS, where stars appear to transfer angular momentum between their protostellar disk and themselves, such that they maintain a nearly constant rotation rate for some time (the disk locking time). Also see e.g., Matt & Pudritz (2005, 2008a; Matt et al. (2012b); Gallet et al. (2019) for more recent models of PMS star-disk interactions. To avoid initializing these models above their critical rotation rates (at which outward centrifugal force is equivalent to surface gravity), rather than starting at an initial angular velocity Ω i set by P rot,i , we initialize the models with a lower angular velocity, and allow them to spin up to their target P rot,i through contraction. Once the models have contracted enough, and have reached their target angular velocity (corresponding to P rot,i ), we force the stars to maintain a constant rotation period until the end of their disk locking time. This is similar to the procedures of e.g., Gallet & Bouvier (2013) and Amard et al. (2019). Incorporating a more physically motivated manner of initializing these models by following the developments of e.g., Matt et al. (2012b) orGallet et al. (2019 for disk locking would be an important next step, and is a goal of future work. The disk locking time influences the initial spin evolution of these stars, and is expected to last until about 2 − 5 Myr, or so. It is not well understood how this process works in detail, nor specifically how long the disk locking time should be, or what all of the conditions that could affect it might be. Some evidence points towards solar mass stars spinning up (through PMS contraction) between about 2 and 13 Myr, where solar mass stars with periods faster than 1 day do not appear in data for the Orion Nebular Cluster, NGC 6530, NGC 2264, nor NGC 2362 (1, 1.65, 2, and 5 Myr respectively; see e.g., Gallet & Bouvier 2013). Disks appear to dissipate largely within 13 Myr, the age of h Per (Currie et al. 2007;Fedele et al. 2010). As the physics of disk locking and PMS angular momentum evolution are highly uncertain, we simply select a representative disk locking time of 3 Myr for all models. We did test results with the variable disk locking time scheme used in Amard et al. (2019), but found it had minimal impact on our results.
With our disk locking time set to 3 Myr, we have used the observed rotation periods for stars in the range 0.95 − 1.05 (where colors were converted to masses via our evolution tracks) in the clusters NGC 6530 (Henderson & Stassun 2012) and NGC 2264 (Affer et al. 2013) at 2 and 3 Myr respectively, to inform our choice of initial rotation periods. Following the end of disk locking, at 3 Myr in our modelling, models contract during PMS evolution and freely spin up, decreasing their rotation period. We have selected initial rotation periods to match the rotation period range observed in h Per (Moraux et al. 2013) at 13 Myr, an age when all stars have largely lost their disks; this is similar to the proce-dure in Amard et al. (2019). Stars in h Per with masses estimated as ranging from 0.4 to 1.4 M generally have rotation periods ranging from about 0.3 to 10 days. We adopt an irregularly spaced grid of initial rotation periods (P rot,i , the initial period at which our disk locked models rotate for ages less than 3 Myr) for our models: 1.5, 3, 4.5, 6, 8, and 12 days, which through PMS spin up, roughly covers the distribution of rotation periods observed by Moraux et al. (2013) at 13 Myr. The most rapidly rotating stars in the Moraux et al. (2013) sample reach up to P rot,i ≈ 0.2 days; we have calculated models at P rot,i = 0.4, 0.6 and 0.8 days, but find that these models near solar mass reach super-critical rotation velocities. Therefore, we do not show our models with P rot,i < 1.5 days in our main results. We do find that these rapid rotators may be needed to reproduce the lithium abundance dispersion, including the low solar value (e.g., as observed in Carlos et al. 2019), as we discuss in Sec. 4.4; this issue is discussed further also in Sec. 5.5.
Although we adopt a flat distribution of rotation periods here for demonstrative purposes, the detailed distributions of stellar rotation periods seems to be dependent on mass dependent factors. Lower masses possess a more pronounced population of fast rotators than higher masses do (Moraux et al. 2013). Some of the fast rotators are likely the products of binary interaction (e.g. Douglas et al. 2017), but e.g., after removing binary candidates, Moraux et al. (2013) found that lower mass stars appear to possess a more significant population of fast rotators compared to their high mass counterparts. This may also be seen at ∼10 Myr in the Upper Scorpius (USco) association. Somers et al. (2017) and Rebull et al. (2018) describe a mass-rotation correlation wherein lower mass M dwarfs appear to rotate much faster than higher mass M dwarfs. This could mean, for example, that lower mass M dwarfs exhibit shorter disk locking times so that they spin up sooner than higher mass stars, or that some other mass dependent initial condition exists, which will likely need to be accounted for in models going forward.
In regard to this, rapidly rotating stars with masses 0.6 M have been difficult to understand. Recent models have worked to incorporate additional physics that utilizes data constrained by magnetic field activity, such as Wright et al. (2011); Wright & Drake (2016); Wright et al. (2018). Two braking models Garraffo et al. 2018) created to model the spin evolution of stars with surface magnetic fields are included in this study and described below.
3.2. Matt et al. (2015) braking Matt et al. (2015) produced a model that took advantage of more extensive rotation period data that had become available around the same time (Irwin & Bouvier 2009;Bouvier et al. 2014). Additionally, Matt et al. (2015) incorporated physically motivated stellar parameter scalings (from the MHD simulations of Matt & Pudritz 2008a;Matt et al. 2012a, as also included e.g., by Gallet & Bouvier 2013;van Saders & Pinsonneault 2013) as an update to previous stellar wind angular momentum loss laws. Following Matt et al. (2015) and the implementation described in Amard et al. (2019), we compute the torque due to magnetic braking aṡ , unsaturated (2) when the magnetic field is unsaturated, anḋ where and γ = 1 + (u/0.072) 2 from Eq. (8) of Matt et al. (2012a), where u is the ratio of rotation velocity to critical rotation velocity, v/v crit . The constants K, m, p, and χ are free parameters, calibrated to data. The values that we have adopted are collected in Table 1. The terms unsaturated and saturated in Eqs. 2 and 3 refer to two apparent regimes of stellar magnetic activity. These regimes have been found in e.g., (Wright et al. 2018, and references therein), where the X-ray luminosity of stars (a proxy for magnetic activity) appears to correlate strongly with the Rossby number R o = (P/τ c ), and becomes nearly constant below some critical value of R o (i.e., the saturated regime, where R o < R osat ). The parameter χ = R o /R osat relies on this critical value, R osat , which was measured in Wright et al. (2018) to be R osat = 0.14; the solar R o is believed to be around 2 (e.g., see See et al. 2016), and thus we take χ = 14. The constant m we set to 0.22, adopting the value from e.g., Matt et al. (2015); Amard et al. (2019). The constants K and p have been solar calibrated to reproduce the rotation period of the Sun (we adopt P = 28 days at the age of the Sun 4.6 Gyr). For K, we adopt a value 1.4 × 10 30 ; for p we adopt 2.6. The constant p was set to reproduce the spread of rotation rates observed in the Pleiades at 125 Myr for roughly 1 M stars according to data from Rebull et al. (2018). In calculating R o , we have adopted a different definition of the convective turnover time (τ c ) than Amard et al. (2019) did in their implementation of Matt et al. (2015).
Semi-empirical approaches have estimated expected values of τ c , such as in Wright et al. (2011), but it is not clear precisely where τ c should be calculated within the interior of a 1D stellar model. Cranmer & Saar (2011) derived an effective temperature dependent function to predict τ c , but this relationship may not be universally applicable to all stars. A number of stellar models base their convective turnover times on mixing length theory (Böhm-Vitense 1958), placing the calculation at some multiple (often 0.5 times) of the pressure scale height, as defined in Gilliland (1985). This roughly places the calculation within range of a convective eddy assumed to be related to the dynamo generating the magnetic field, tying the definition of τ c loosely to the magnetic dynamo. However, it is also unclear if fully convective stars generate surface magnetic fields with the same dynamo mechanisms that partially convective stars use (Mullan & MacDonald 2001;Reiners & Basri 2007;Irwin et al. 2011;Wright & Drake 2016;Wright et al. 2018). See Charbonnel et al. (2017) for a recent analysis of various definitions of the convective turnover time in relation to data.
We calculate τ c as where α MLT is the convective mixing length parameter (α MLT = 1.82 in our models; Choi et al. 2016), H P (r) is the scale height, and v c (r) is the convective velocity within the star at radius r. We calculate τ c one half a pressure scale height above the bottom of the outermost convection zone, defined as where r = r BCZ + 0.5H P (r), and where r BCZ is the position of the bottom of the outer convection zone. This is a simple approach that provides better agreement to the semi-empirical trend of τ c with mass in Wright et al. (2011) at an age of 1 Gyr (as assumed in that work) for our models than the definition of τ c at one half of a pressure scale height in Gilliland (1985) does.

Garraffo et al. (2018) braking
We have also included the braking model of Garraffo et al. (2018). This model is based on MHD simulations (Réville et al. 2015;Garraffo et al. 2015Garraffo et al. , 2016 that suggest the efficiency of magnetic braking depends on the complexity of the stellar magnetic field, not just its strength. The model combines an expression for angular momentum loss that is based on the Skumanich law, designatedJ dip , with a function Q J (n) to modulate this basic loss formalism by the magnetic field complexity, parameterized as n (Garraffo et al. 2016). The full angular momentum loss formalism iṡ representing Skumanich spin down (P ∝ t −1/2 ) that the model converges to when the magnetic field becomes dipolar. Here, c is a free parameter setting the overall strength ofJ dip , Ω is the angular velocity of the stellar surface, and τ c is the convective turnover time.
The function Q J (n) was derived in Garraffo et al. (2016) and serves to modulate the dipolar Skumanich angular momentum loss via magnetic field complexity. It takes the form where B is the magnetic field strength, and n parameterizes the magnetic field complexity. In practice, this second term only begins to matter in relatively complex magnetic fields where n > 7. Instead of following Garraffo et al. (2018), who imposed an upper limit at n = 7 so that this term may be ignored, and Q J (n) takes the form we simply ignore the term, and do not place an upper limit on n. We do this for simplicity because we currently do not calculate magnetic field strength in our models. The ignored term adds a constant to Q J (n), such that more complex magnetic fields are even more suppressed than they would be were the term included. We find that our models that reach n > 7 are often too suppressed regardless of the inclusion of this term, and so capping n has little consequence on our results. The magnetic complexity number n is itself a function of R o , with P being the rotation period of the star, P = 2π/Ω, representing the magnetic field complexity. In Eq. 10, n = 1 corresponds to a dipole field, while higher n represents more complexity and higher order multipoles. Also in Eq. 10, the free parameters a and b are free parameters. This equation was constructed to match expected magnetic field complexity trends based on Zeeman-Doppler-Imaging ( 2018), but our adopted saturation point lies at Ro,sat = 0.14, rather than their 0.11, as described in the text.
The parameters we adopt for this model were arrived at through calibration to cluster data, and are listed in Table 2. The parameter c controls the overall strength of angular momentum loss in the model, and was tuned to reproduce the solar rotation rate. We arrived at a larger value of a than in Garraffo et al. (2018). In Eq. 10, lower values of a mean that n reaches its minimum (becomes dipolar) at lower values of Ro; higher values of a mean that it takes higher R o before a dipolar field is achieved, effectively meaning that models take longer to reach a Skumanich type spin down, and remain spinning faster, longer. Lastly, we changed b to control magnetic suppression at higher R o (i.e., slower rotation and typically older ages). We adopt a smaller value of b so that our solar model reaches the solar rotation rate as well as matches the rotation periods of solar-like stars observed in open cluster data prior to 4.6 billion years. A larger b means that at later ages, stars regain magnetic complexity, and a smaller b means that stars stay mostly dipolar (and more closely follow Skumanich spin down) at later ages. Our adopted values are tuned mostly to reproduce solar rotation rates, but it is not clear that n should be universally assigned as in Eq. 10 for all stellar masses. As n (magnetic field complexity) increases, the portion of the stellar surface covered by open field lines decreases. Fewer open field lines leads to less outgoing stellar wind as complexity increases, reducing the capacity for magnetic braking to act in slowing the star down. Thus, more complex stars will tend to rotate faster, and less complex stars more slowly. The fast and slow branch rotators observed in open clusters may be explained through this relationship between field complexity and braking efficiency. Higher magnetic field complexity leads to greater suppression of the braking process. This has been shown previously in Garraffo et al. (2018), though models in that case were not evolved with these effects in place; rather, they were applied on pre-computed stellar models. Here, we self-consistently evolve the models with this braking scheme.

RESULTS
In this section we present the evolution of rotation period with time according to our models when using the   We reiterate that lower mass M dwarfs are observed to contain a larger fraction of rapid rotators, than higher mass stars do (e.g., Somers et al. 2017;Rebull et al. 2018). In this study, we have not attempted to match the observed distribution of rotation rates in such detail, rather we intend to display simply a representative range of observed rotation rates; matching distributions is a goal of future work. In Sec. 4.3 we discuss model behavior on the Hertzsprung-Russell diagram (HRD), and in Sec. 4.4 the surface lithium depletion seen in our solar mass models (panels (c) and (d), respectively in Figs. 2 and 4). In Figs. 2 and 4, the data that we compare to is near solar mass, in the range 0.95 to 1.05 M (based either directly on masses or on colors provided in the data tables; colors were converted to masses via our isochrones in the latter case). In interpreting our results, it is important to bear in mind that the angular momentum evolution of our models (and thus our results in general) is still linked to the diffusion parameter ν 0 (Sec. 2.1), representing the degree of core-envelope coupling. Core-envelope coupling affects the angular momentum transfer timescale, and can produce different rotation period evolution patterns. We do not perform an extensive parameter study of ν 0 in this work, but rather adopt a single value, calibrated to reproduce the differential rotation of the Sun (see e.g., Pinsonneault et al. 1989;Krishnamurthi et al. 1997;Denissenkov et al. 2010;Somers & Pinsonneault 2016 for studies more focused on exploring this). We further note issues related to core-envelope coupling, and their impact on our results in the following subsections.
Overall, our results are in line with those of similar studies, and the original papers of Matt et al. (2015) and Garraffo et al. (2018) in matching open cluster rotation periods (i.e., Figs. 2 and 2). Our results for Matt et al. (2015) models qualitatively match those shown in Fig. 7 of Amard et al. (2019), and our results from the Garraffo et al. (2018) model match the results in Fig. 4 of that respective paper (although Garraffo et al. (2018) performed population synthesis, so plot stellar densities, rather than simply positions on period-mass, or periodcolor diagrams). All results generically show too much angular momentum loss for sub-solar mass models at 1 Gyr and older ages, and carry similar morphologies to the results shown here in period-mass space, discussed further below.

Matt et al. (2015) comparisons
Results from our models employing the Matt et al. (2015) braking scheme show a plausible reproduction of the cluster data. Agreement is generally good between our 1 M model in panel (a) of Fig. 2 and cluster observations, given the presence of outliers (indicated by black crosses). In panel (b), we have included measured Rossby numbers for solar-like stars (sub-sampling the data from Vidotto et al. 2014 in the mass range 0.95 M ≤ M ≤ 1.05M ), and our calculated Rossby numbers agree fairly well over time.
For the evolution of rotation period regarding all masses, view Fig. 3. E.g., for the Pleiades (assuming an age of 125 Myr) by Rebull et al. (2016), the overall morphology in P rot -mass space agrees fairly well, but appears to overestimate the number of stars that rotate quickly, as the data shows a more collapsed sequence of slow rotators by this age. Many lower mass (≤ 0.6 M ) stars are just halting contraction and entering the ZAMS around this age and have spun up since the end of disk locking at 3 Myr. Higher mass stars have already completed contraction prior to this and begun to spin down under the influence of magnetic braking, converging towards the slow (Skumanich) sequence of stars. This model underestimates the slowly rotating stars around 0.5 M , and generally overestimates the number of rapid rotators with masses 0.8 M .
At the age of the Praesepe (assumed here to be 676 Myr), we compare to data from Douglas et al. (2017). Our models implementing the Matt et al. (2015) model allow for excellent reproduction of the slow rotators across the entire mass range. The lower mass, fast rotators are reproduced as well. The behavior is similar to what is shown in Matt et al. (2015), but Douglas et al. (2017) found that their Matt et al. (2015) models predicted too many quickly rotating stars with mass ≤ 0.8 M . Our models show similar behavior, with stars in this mass regime rotating fast. Essentially all of the stars 0.3 M ≤ M ≤ 0.8 M rotate too quickly, in comparison to the data which shows the majority of stars with these masses are slow rotators.
Faster rotation rates and lower masses mean shorter P rot and larger τ c (because of deeper convection zones at low mass), and thus smaller R o for these stars. The smaller R o keeps their magnetic fields in the saturated regime longer, delaying their spin down. This may be seen in Fig. 2, panel(b) for our 1 M model, and in Fig.  3, bottom row for all masses. So, while our implementation of Matt et al. (2015) is able to reproduce quickly rotating low mass stars, as observed in the data, it seems to overestimate the number of these stars, especially in the range 0.3 M ≤ M ≤ 0.8 M , or so.
Towards later ages, data from Curtis et al. (2019) of the open cluster NGC 6811 (estimated at 1 Gyr old) shows that stars have largely converged to a single sequence by this time. Meibom et al. (2015) also observed this in the 2.5 Gyr old NGC 6819 (not shown in Fig.  2). Using Matt et al. (2015), rotation periods tend to be too large for masses < 1 M . The models correctly predict a collapsed sequence of slow rotators by this age, but have generally experienced braking that is too effi-

Garraffo et al. (2018) comparisons
In Fig. 4, we show an overview of behavior for our solar mass model. The model does well matching the evolution of P rot up to 4.6 Gyr, as seen in panel (a). Likewise, it does a fair job matching the data of Vidotto et al. (2014) for measured Rossby numbers of solar-like stars (their data sub-sampled to be in the mass range 0.95 M ≤ M ≤ 1.05M here).
Refer to Fig. 5, top row, to see the evolution of P rot at selected ages for all masses. Using the Garraffo et al. (2018) braking scheme, models are able to reproduce the general morphology of P rot -mass space across time, but predict braking that is often too efficient at later ages for masses 0.3 M ≤ M ≤ 0.8 M (see Fig. 5). At 125 Myr, the slowly rotating stars near 0.5 M are captured better than when using Matt et al. (2015) braking. Slower models with roughly M > 0.5 M have large enough R o to enter a regime where n is closer to a dipole at this age (see Fig. 5 and Fig. 1), causing them to spin down more aggressively. Apparently the spin down is sufficient in the case for masses near 0.5 M , allowing them to reach periods closer to what is observed at this age.
The implementation of Garraffo et al. (2018) in our models qualitatively matches the morphology observed in the Praesepe, however it tends to predict too much spin down here as well with stars in the mass range 0.3 M ≤ M ≤ 0.8 M . Still, this model does well in reproducing the quickly rotating stars with masses < 0.3 M , and predicts a collapsed sequence of slow rotators by this age, (albeit rotating too slowly), as well as a sharp drop towards fast rotation with fully convective stars (M < 0.3 M ). Simultaneously capturing these low mass fast rotators, and a tight sequence of slowly rotating, higher mass stars allows this braking model to capture the bimodality observed e.g., in the Pleiades and the Praesepe, which is a relative strength compared to using our implementation of the Matt et al. (2015) formalism. In Fig. 5, notice the sharp transition from slow to fast rotation going towards lower masses, rather than the gradual, more continuous behavior seen with the Matt et al. (2015) braking model.
The reason models are able to reproduce the low mass fast rotators and tight sequence of slow rotators simultaneously under this braking model comes down to the parabolic dependency of the magnetic complexity (n, Eq. 10) on R o . In general our low mass models have   deeper convection zones and higher τ c , driving smaller R o (see Figs. 3 and 5, bottom rows), and faster rotation also decreases R o further. Due to their low R o (which may be seen in both Fig. 5 and Fig. 4), low mass models experience greater magnetic complexity (as per Eq. 10), suppressing their spin down until later ages. As seen in Fig. 10, there is a fairly narrow range of R o over which stars experience dipolar (strong) magnetic braking; higher mass models reach this threshold sooner than lower masses. The higher mass models (with generally higher R o ) spin down sooner, and re-enter a phase of increased magnetic field complexity (high n) as their R o increases, slowing their braking once again, and keeping them in a tight sequence of slow rotators. At the age of NGC 6811 (1 Gyr), models have spun down too much, similar to results from our implementation of Matt et al. (2015) braking, although the discrepancy is even more severe in this case. Again, models with mass ≥ 1 M are represented fairly well, and the models correctly predict a collapsed sequence of slow   Garraffo et al. (2018) braking model here. Note that the bimodal nature of rotation period data in e.g., the Pleiades and the Praesepe is reproduced; the model shows a tight sequence of slow rotators at higher masses (although often too slow), and jumps to faster rotation rates at lower masses, sparsely populating the region in between these regimes. Compare this to how the Matt et al. (2015) model produces many more rapid rotators at low masses that do not reflect the observed scarcity of such stars in Fig. 3, top row. rotators. A portion of masses below 0.6 M are predicted to remain quickly rotating, but the observations do not include stars in this mass range for comparison. This suggests that for masses < 1 M , there may be some additional mass dependence missing from this braking model, where they should be locked (by high n) at smaller R o than higher masses are (effectively, b in Eq. 1 would need to increase towards lower masses). However, and as previously mentioned in Sec. 4.1, this issue may be related to mass dependency in core-envelope coupling timescales, as suggested by Spada & Lanzafame (2020) (see Sec. 5.2 for more on this).

Effects on the Hertzsprung-Russell Diagram
In Figs. 2 and 4, panel (c) shows the evolution of our 1 M model on the HRD. There is a minimal effect going towards lower masses, and the effects are similar for masses up to at least 1.3 M included in this study. The inclusion of PMS rotation causes fast rotators to reach the ZAMS at a cooler temperature than slower rotators. This is slightly more pronounced in the case of Matt et al. (2015) braking, as those models take longer to spin down, and so are rotating faster near ZAMS than models do under our modeling with Garraffo et al. (2018) braking. Faster rotation rates cause stronger centrifugal force on the star, enhancing the effect of gravity darkening (e.g., von Zeipel 1924;Espinosa Lara & Rieutord 2011;Paxton et al. 2019), and causing stars that rotate fast enough to appear cooler upon reaching the ZAMS. This would suggest that low mass, rotating stars have some spread due to gravity darkening on the MS. Although, this effect is barely noticeable, except with the fastest rotation rates (P rot,i≤1.5 days).

Lithium burning
Additional mixing from stellar rotation has been proposed as a possible means of explaining the surface lithium abundances of ∼ 1 M stars. Observations of stars in open clusters (e.g., Sestito & Randich 2005), and our own Sun, reveal that solar mass stars appear to have a dispersion in their surface lithium abundances, hinting that something must cause more efficient lithium depletion in otherwise similar stars (see recent solar twin studies in, e.g., Thévenin et al. 2017;Carlos et al. 2019). Somers & Pinsonneault (2016) provides an overview of this, and also shows that including an additional angular momentum transport (parameterized by ν 0 here; refer to Sec. 2.1) may explain the dispersion. Rotating stellar models typically are not able to reproduce the solar lithium depletion via rotational mixing alone (e.g, Amard et al. 2016Amard et al. , 2019. As seen in Figs. 2 and 4, our default models do not reproduce the solar lithium abundance. We discuss this further in Sec. 5.3.

DISCUSSION
Our implementations of Matt et al. (2015) and Garraffo et al. (2018) allow us to compare a more traditional magnetic braking model with one that models the effect of magnetic complexity on stellar spin down. According to our results, magnetic complexity may be valuable in providing a mechanism for periods of reduced magnetic braking as stars evolve. This may aid in reproducing the solar lithium abundance, and in matching the behavior of rapidly rotating low mass stars at ages of the Praesepe and prior, and as a possible mechanism for the observed stall in angular momentum loss (Agüeros et al. 2018;Curtis et al. 2019) at later ages.

Discrepancies near 600 Myr
At the age of the Praesepe (assumed 676 Myr in this study), spin down via our implementation of Matt et al. (2015) appears to overestimate the presence of rapid rotators with 0.3 M < M < 0.8 M . Meanwhile, our models based on Garraffo et al. (2018) predict a more collapsed sequence of stars in this mass regime, but they tend to spin too slowly. In both cases, it is possible that re-calibrations of, or modifications to, the braking model may be necessary in this mass regime. Angular momentum loss appears to be too efficient in the case of Garraffo et al. (2018) braking, and too weak in the case of Matt et al. (2015).

Overestimated spin down after 1 Gyr
Our models generally experience too much angular momentum loss by 1 Gyr; compared to observations, the rotation periods of M < 1 M models are too long. Curtis et al. (2019) discuss that at 1 Gyr, there appears to be an epoch of stalled angular momentum loss in stars ≤ 1 M . Angular momentum loss seems to become less efficient around the age of the Prasepe for some time, as stars at the later age of NGC 6811 do not appear to have slowed down by very much more. Curtis et al. (2019) find that gyrochronology models tend to predict too much angular momentum loss, as our models do. The root of this issue is unclear, but may be due to metallicity effects (Angus et al. 2015), or perhaps a mass dependence on the angular momentum loss rate that is missing from current models (see also e.g., Agüeros et al. 2018).
One explanation for this discrepancy could lie in the parameter ν 0 (Sec. 2.1), which is related to the degree of core-envelope coupling within stars. Stronger coupling causes angular momentum transport to behave closer to a solid body (more efficient core-envelope transport). Somers & Pinsonneault (2016) studied this, and found evidence that the level of core-envelope coupling in stars may be mass dependent. Lanzafame & Spada (2015) studied this independently and found the same thing, and that this relationship could be crucial for modeling the slow rotator branch of period-mass diagrams. They discuss a trend suggesting that less massive stars experience stronger differential rotation (less core-envelope coupling). This would imply less efficient braking for lower mass stars as well. This is to say that the braking overall appears to be too strong in our current subsolar mass models, and reducing core-envelope coupling in these models could make them overall match the slow rotator branch better. This effect would not help up simultaneously match the fast rotator branch; that would still be due to e.g., saturated stellar winds (under the Matt et al. 2015 model) or increased magnetic complexity (under the Garraffo et al. 2018 model). As this is not currently accounted for in our models (we assume a single ν 0 , i.e., core-envelope coupling level, for all stars), it could be another source of the discrepancy and should be accounted for in the future. Recently, Spada & Lanzafame (2020) have presented a model with mass scaling where core-envelope coupling decreases towards lower masses, and have been able to model the slow rotator branch very well. A further physical explanation for this epoch of stalled spin down could be a period of greater magnetic field complexity, but ultimately any mechanism that dampens angular momentum loss rates could be at play, including, e.g., mass loss rates (as in See et al. 2019). However, if magnetic complexity were invoked, under the model of Garraffo et al. (2018), this would be reliant on the form of the function for n (Eq. 10), which could have its mass dependence modified to allow this. Eq. 10 was primarily formed to describe a general behavior of the models based on ZDI observations of stars, but is fairly unconstrained. If mass dependence were added, such that lower mass stars experience greater magnetic complexity after the age of the Praesepe, their spin down would be suppressed, achieving the observed behavior, similar to what is shown in our Fig. 6, for the case where a has been increased from our default 0.03 to a = 0.05. Recently See et al. (2019) have found that higher order modes in stellar magnetic fields may be more important at later ages as well, and van Saders et al. (2016Saders et al. ( , 2019 have found evidence for this too. If mass dependency in core-envelope coupling were introduced to our models, parameters would likely need to be re-calibrated, and it is unclear if a modified mass dependency in magnetic complexity would be required to reproduce the slow rotator branch; results from Spada & Lanzafame (2020) suggest that it may not be necessary for this purpose.

Reproducing the solar lithium abundance
As seen in Figs. 2 and 4, our solar mass models do not reproduce the solar lithium abundance, nor the observed dispersion in lithium abundances of solar-like stars (Carlos et al. 2019) under our default assumptions. This could be reconciled in at least two ways: either through prolonged rapid rotation through early ages, or strong differential rotation (as in Somers & Pinsonneault 2016). The historical arguments for rotationally induced mixing as a cause of lithium depletion involve several components. One is the observation that lithium depletion appears to follow a similar time dependence to stellar spin down (Skumanich 1972;Charbonnel et al. 1992). Additionally, the depletion has a dispersion that exists at fixed age, mass, and composition (as mentioned previously in Sec. 4.4). In absence of rotational mixing, standard stellar models predict lithium depletion in PMS stars, with deep convective envelopes, and not on the main sequence for solar-like stars. However, lithium depletion has been observed in open clusters (Soderblom et al. 1993;Sestito & Randich 2005;Thévenin et al. 2017;Carlos et al. 2019), in solar-like stars, taking place during the main sequence, calling for some mixing process (e.g., induced by rotation or other mechanisms) than standard stellar models operate with. Lithium is typically the most dramatically depleted light element in stellar surface abundances, but it is clear that some global mixing process is taking place within stars, as other light elements (boron and beryllium) can become depleted alongside lithium, especially in F-and B-class stars.
To test whether rapid rotation could facilitate greater lithium depletion in our models, we have increased the constant a in the Garraffo et al. (2018) braking model (Eq. 10) to a = 0.05 in Fig. 6, causing lower Rossby numbers to experience higher n (and thus have their braking suppressed more strongly at early ages). In this case, our solar model with P rot,i between 0.8 and 1.5 days appears to come close to reproducing the solar lithium abundance. In the case of Garraffo et al. (2018) braking, it would then be higher order magnetic field contributions that suppress magnetic braking at early ages, leading to an extended period of rapid rotation. Our models under the Matt et al. (2015) braking model could also achieve this through a lower p value, weakening the braking model's dependence on R o , and keeping them spinning faster for longer. This scenario would suggest that the Sun began life as a rapid rotator. This scenario is uncertain with our models however, because these rapid rotators reach critical rotation (Ω/Ω c ≈ 1), where rotating MESA models become unreliable (Paxton et al. 2013(Paxton et al. , 2019. More importantly, the scenario presented in Fig. 6 does not seem to be supported by data. However, for young stars of about solar mass and below, rapid rotators appear to be less depleted in lithium than slow rotators (Bouvier 2020). Physical interpretations of this trend vary, with some citing inhibition to convective mixing with greater rotation rates (Baraffe et al. 2017), or convective mixing inhibition from rotationally induced magnetic fields (e.g., Somers & Pinsonneault 2014, 2015bJeffries et al. 2017), or perhaps star-disk interactions on the PMS (Bouvier 2008;Eggenberger et al. 2012). Whatever the case, observations say that rapid rotation early on leads to weaker lithium depletion, not stronger, indicating our models are incomplete in explaining this phenomenon correctly. Furthermore, as shown by Delgado Mena et al. (2014); Carlos et al. (2019), the majority of solar twins at the age of the Sun are more depleted in surface lithium more than our slowly rotating models predict. I.e., our slowly rotating models have surface lithium abundances that are too high. It then seems that observations of solar twins tell us that (regardless of P rot,i ) all models need a boost in their mixing levels beyond what rotational mixing provides in order to match observations. These could come from various processes that alter the internal angular momentum transport (e.g., as explored by Somers & Pinsonneault 2016).
As found here and in Amard et al. (2016Amard et al. ( , 2019, Charbonnel & Talon (2005);Eggenberger et al. (2005) noted that other mixing processes are needed to supplement rotationally enhanced mixing. Mixing mechanisms like gravity waves and magnetic fields do not directly mix material, but rather facilitate greater hydrodynamical transport in conjunction with rotation induced mixing. Somers & Pinsonneault (2016) parameterized these effects through the degree to which they enhance coupling between the stellar core and envelope via enhanced angular momentum transport (similar to the concept of our ν 0 parameter). They found that solid body rotation (large ν 0 , and strong core-envelope coupling) paradoxically produces less efficient mixing than stars with differential rotation. They concluded that hybrid models that experience some degree of differential rotation (as observed in the Sun) are required to reproduce the solar lithium abundance. Too much differential rotation (small ν 0 ) produces too weak of a core-envelope coupling for solar mass stars to converge to the tight sequence of slow rotation observed at late ages (e.g., the small spread in periods at late cluster ages Meibom et al. 2015;Curtis et al. 2019, and see the final columns of Figs. 2, 4).
Our chosen default value of ν 0 = 2 × 10 4 cm 2 s −1 creates differential rotation reminiscent of what is observed in the Sun, and produces a spread of rotation rates that appears to match the spread observed in open clusters 10 1 10 2 10 3 10 4  Figure 6. Similar to Figs. 2 and 4, but just the evolution of (a) Prot and (b) surface lithium abundance (log( 7 Li/H) + 12) for our solar model using the Garraffo et al. (2018) braking model. We have changed a from our default 0.03 to 0.05, delaying braking until later ages, and allowing greater rotational mixing for demonstrative purposes. The Matt et al. (2015) could achieve this as well through a smaller p value. Colors for each solid line are the same as described in Fig. 2. Although effective in reproducing the solar lithium abundance, this scenario is not supported by observations. Note that the fastest rotators in panel (a) exhibit super-critical rotation, and so their behavior is erratic compared to our models with Prot,i > 1.5 days.
over time in our solar mass model (Figs. 2 and 4,panel (a)). If we were to decrease ν 0 somewhat from our default value, we would induce stronger differential rotation, driving stronger shears, and more hydrodynamic transport (e.g., Somers & Pinsonneault 2016) between the core and envelope, producing stronger mixing. Processes that may alter the internal angular momentum transport (such as gravity waves or internal magnetic fields) provide another route for reproducing the solar lithium abundance that could work regardless of the adopted braking model. We have tried various values of ν 0 , but they had little effect in enhancing lithium depletion with our current models, unless we re-calibrate the constants f c and f µ in the Pinsonneault et al. (1989) diffusion approximation (see, Choi et al. 2016 for additional details). In particular the parameter f c , which scales the efficiency of composition mixing to angular momentum transport, could be adjusted (as in Somers & Pinsonneault 2016) to produce qualitatively and quantitatively different sets of model outcomes, with more or less lithium depletion. We save this re-calibration for future work. While rotationally induced mixing may not be the direct driver of lithium depletion, rotation rate is still thought to be related to the depletion mechanism, and may account for the dispersion of abundances observed in solar twins, where the Sun appears to be especially low compared to most other stars (Carlos et al. 2019) at ∼4.6 Gyr. Observations by Sestito & Randich (2005) (also in Castro et al. 2016) provide a map of how lithium abundance may evolve with time for various stars in open clusters, including near solar masses. A range of initial rotation rates amongst stars could explain this by varying the efficiency of the lithium depletion process(es).
Assuming that clusters are born with some fraction of fast rotators (as in h Per and USco Moraux et al. 2013;Somers et al. 2017;Rebull et al. 2018), and that these open clusters truly show us an evolutionary sequence of rotation rate with time (e.g., Fritzewski et al. 2020, but c.f., Coker et al. 2016, observations suggest that clusters birth stars with a wide range of initial rotation periods. Observations e.g., from Hartman et al. (2010); Rebull et al. (2016) show that solar mass stars born with rapid rotation (P rot,i 0.6 days) are likely braked prior to ∼100 Myr, as very few such rapidly rotating stars near solar mass are observed beyond roughly this age in comparison to the distribution at 13 Myr (Moraux et al. 2013). So, although a tight sequence of slow rotators is observed for solar masses at later ages, a wide range of initial rotation periods for solar mass stars could still have existed within a cluster at early ages, leaving an imprint on the subsequent dispersion of lithium abundances at late ages.

The role of magnetic complexity
In our modeling, under Garraffo et al. (2018), magnetic complexity serves the practical role of stalling angular momentum loss. There is evidence that the efficiency of angular momentum loss changes over time, as in Agüeros et al. Presently, we have that n evolves through the evolution of a star's R o , which tends to increase as stars get older. As stars spin down, P rot increases, while along the main sequence, τ c hardly varies. Higher mass stars have smaller convective envelopes, and smaller τ c generally (see also the trend found in Wright et al. 2011), meaning they arrive on the ZAMS with higher R o , tending to appear more as magnetic dipoles than low mass stars. Thus, higher mass stars tend to spin down at earlier ages than lower masses, as may be seen in Fig.  2. This behavior helps our models using Garraffo et al. (2018) braking achieve a collapsed sequence of slowly rotating stars more readily than when employing Matt et al. (2015) braking, while also still predicting the presence of fast rotators with M < 0.3 M .
The behavior of n at higher R o according to Eq. 10 is meant to produce a period of high magnetic complexity again under the calibration used in Garraffo et al. (2018) (at what would tend to be later ages). The stalled angular momentum loss brought about by this could be an additional strength of including magnetic complexity in magnetic braking models. In the present case, we may need a modified mass dependence on the term a in Eq. 10, such that lower mass stars slow down, but re-enter a period of high magnetic complexity (suppressed angular momentum loss) again at lower R o than higher mass stars do. In other words, the function in Fig. 1 may need to shift leftwards, towards smaller R o as mass decreases. It is also worth noting that at low enough Rossby numbers the saturated regime (described in Sec. 3.2) may transition to a "super-saturated" regime (James et al. 2000;Wright et al. 2011;Jeffries et al. 2011;Argiroffi et al. 2016). Evidence for super-saturation is limited to a fairly low number of stars, but if the rapid rotation of super-saturated stars is due to inefficient angular momentum loss, the implicit link would be that supersaturation may be related to higher magnetic field complexity under the Garraffo et al. (2018) model, which is partially based on ZDI data that says magnetic field complexity increases with increasing rotation rate in active stars. Other dependencies of the angular momentum loss rate and magnetic field strength on stellar properties, such as atmospheric pressure and mass have been explored (e.g., by Pinsonneault 2013), and could provide additional physical updates to these models in the future. Whether higher order modes should play a significant role in modeling angular momentum loss is currently uncertain, as See et al. (2019See et al. ( , 2020 have shown that a majority of stars may not be significantly affected by high order magnetic field topologies, except perhaps at late ages. The results of Garraffo et al. (2018) were based on single mode MHD simulations, to conceptualize the effect of high order magnetic field topology on outgoing stellar winds (Réville et al. 2015;Garraffo et al. 2015Garraffo et al. , 2016. Finley & Matt (2017 conducted MHD simulations with mixed magnetic field modes (as would be found in nature) to study the extent to which higher order magnetic fields may affect outgoing stellar winds under more realistic conditions. They found that generally, the dipolar mode dominates the wind morphology for a wide range of mixed magnetic field topology configurations, but that certain configurations do show that high order magnetic fields dominate the wind morphology, so long as magnetic energy is concentrated moreso in the higher orders. Garraffo et al. (2016) compared their approximations to ZDI maps and MHD simulations for several real stars. They found very good agreement between the mass and angular momentum loss rates for the observed ZDI map MHD simulations and those of the pure modes that corresponded to the dominant magnetic field order of the observed ZDI maps. The decomposed observed ZDI maps showed fields dominated by higher order magnetic fields (n > 1, i.e., non-dipolar).
As analyzed by Lehmann et al. (2019), it appears to be difficult to accurately recover the magnetic fluxes that may be associated with the dipolar, quadrupolar, and octopolar magnetic field orders via ZDI due to the limits of spectral resolution. Given these difficulties, it is not easy to understand the partition of total magnetic flux among different field orders, and the true influence of higher order magnetic fields on actual stellar wind and angular momentum loss rates. If, as suggested by See et al. (2019See et al. ( , 2020, magnetic flux is primarily distributed in the dipole mode for a majority of stars, then high order magnetic complexity may play a reduced role on angular momentum loss rates. See et al. (2020) found that below a certain threshold of mass loss, non-dipolar magnetic fields negligibly influence angular momentum loss. They estimate mass loss rates for a sample of stars, and find many to be below this threshold, but note that estimating the mass loss rates is subject to significant uncertainty. Furthermore, Garraffo et al. (2015Garraffo et al. ( , 2016; Finley & Matt (2017); Garraffo et al. (2018); Finley & Matt (2018) point out that high order magnetic fields can reduce mass loss rates themselves (if they dominate the magnetic flux). Thus, it is not clear if estimated mass loss rates might be low because of high order magnetic fields in the first place. If high order magnetic fields comprise a dominant fraction of the magnetic field energy compared to the dipole mode, high order magnetic fields should have some effect on the angular momentum loss rate, and magnetic complexity should not be ignored in modeling magnetic braking.

Super-critical rotation
As mentioned in Sec. 3.1, we have calculated faster rotation rates (P rot,i = 0.4, 0.6 and 0.8 days) than what is shown in the bulk of our results. We have excluded these rapidly rotating models for the most part because near solar masses, they tend to reach super-critical rotation rates, at which point calculations of rotating models in MESA become uncertain. This does not necessarily mean that these stars would not exist in nature; indeed some solar mass stars rotating this fast are observed in the data of Irwin et al. (2008); Henderson & Stassun (2012); Affer et al. (2013), in clusters 2-8 Myr old or so. The number of rapidly rotating solar mass stars appears to diminish between about 13 and 100 Myr (Moraux et al. 2013). It would be straightforward to assume that the reason is because of magnetic braking, in lieu of a confident modeling of stellar evolution under critical rotation.
There is evidence that solar mass rapid rotators cease to rotate rapidly in open cluster rotation period data between ∼13 and ∼100 Myr(e.g., Hartman et al. 2010;Meibom et al. 2011Meibom et al. , 2015Núñez et al. 2015;Rebull et al. 2016;Douglas et al. 2017;Curtis et al. 2019, c.f., at early ages Moraux et al. 2013). It is likely that the rapid rotators have simply been braked, as rapidly rotating stars may be necessary to explain the observed dispersion in stellar lithium abundances (Sestito & Randich 2005;Delgado Mena et al. 2014;Castro et al. 2016;Carlos et al. 2019). As discussed above in Sect. 5.3, rotation is thought to be linked to the efficiency of the Li depletion mechanism (Bouvier 2020). Likely then, the braking of these rapid rotators takes place somewhere between 13 and 100 Myr; presumably, the more slowly rotating stars begin braking in a similar time frame. Our default implementations of Matt et al. (2015) and Garraffo et al. (2018) braking achieve this, rather than e.g., the scenario presented in our Fig. 6, though our fastest rotators remain unstable. We could tune the parameters in our implementations of the braking models to prevent our rapidly rotating solar mass stars from reaching critical rotation, but opt not to because they are not physically restricted from doing so, and they are not highly impactful on our results.

CONCLUSIONS
We have implemented the magnetic braking models of Matt et al. (2015) and Garraffo et al. (2018) in the MIST framework of the MESA stellar evolution code. Our goal was to test their predictions against observed stellar rotation period data and examine some of their implications for some aspects of stellar physics associated with rotation.
In terms of practicality, we find that the Matt et al. (2015) braking model provides a good overall description of stellar spin down with our models. However, we also see potential for the inclusion of magnetic complexity, as in the Garraffo et al. (2018) braking model to improve our modeling abilities. The Garraffo et al. (2018) formalism captures the bimodality of the observed rotation period distributions quite well, while the Matt et al. (2015) model does not. At old ages, both braking models predict rotation that is too slow, while prior to this, the Matt et al. (2015) model tends to predict rotation that is too fast for masses 0.3 M M 1 M , and the Garraffo et al. (2018) formalism predict rotation too slow in that mass regime at 600 Myr onwards. In spite of these shortcomings, it seems that magnetic complexity could play a unique role in characterizing the observed bimodal rotation period distributions in clusters younger than about 1 Gyr.
That both braking models tend to predict rotation that is too slow at ages > 1 Gyr, seems to suggest a missing mass dependant effect in the evolution of P rot in our models. Aside from explaining the observed bimodal distributions of low mass stars, increased magnetic complexity is a mechanism through which angular momentum loss may be halted by suppressing magnetic braking at late ages, keeping stars from spinning down too much. Recent findings suggest that a halt to angular momentum loss may be necessary at these late ages (van Saders et al. 2016;Agüeros et al. 2018;Curtis et al. 2019;van Saders et al. 2019). At least two mechanisms that could provide a mass dependent effect are a modified dependency in magnetic braking, or a mass dependency in core-envelope coupling time scales (Lanzafame & Spada 2015;Somers & Pinsonneault 2016;Spada & Lanzafame 2020). Of course, both effects could act on stars, and their interplay should be studied and constrained in future work on resolving spin down at ages greater than 1 Gyr. Thus, we plan to incorporate mass dependence in core-envelope re-coupling in future models to study these effects.
Studying the solar lithium abundance has provided a useful constraint on mixing efficiencies in these models, and in constraining the spin down rate of our models. Our default models show lithium depletion, but it is not sufficient to reproduce the solar lithium abundance. Extended periods of rapid rotation at early ages allow for a sufficient amount of lithium depletion in our solar model via enhanced rotational mixing. However, observations make this scenario appear unlikely. In regards to solar lithium depletion, there then appears to be a mixing process common to all solar mass stars (slow and rapid rotation) that is missing. E.g., Carlos et al. (2019) show that solar mass stars are generally more depleted than the majority of our models predict. Models may benefit from additional sources of mixing, such as gravity waves, or internal magnetic fields (Charbonnel & Talon 2005;Eggenberger et al. 2005;Denissenkov et al. 2010;Somers & Pinsonneault 2016), and/or other effects (e.g., variable convective mixing efficiencies Baraffe et al. 2017 linked to rotation rate).
We have laid groundwork for the extension of MIST (Dotter 2016;Choi et al. 2016;Gossage et al. 2018) model grids here with implementations of magnetic braking, important for allowing our rotating model grids to cover the complete stellar mass range. In future work, we plan to release an extended grid of models within the MIST framework utilizing the rotation methodology presented here.