A semiempirical model for magnetic braking of solartype stars
Abstract
We develop new angular momentum evolution models for stars with masses of to and from the premainsequence (PMS) through the end of their mainsequence (MS) lifetime. The parametric models include magnetic braking based on numerical simulations of magnetised stellar winds, mass loss rate prescription, coreenvelope decoupling as well as disk locking phenomena. We have also accounted for recent developments in modelling dramatically weakened magnetic braking in stars more evolved than the Sun. We fit the free parameters in our model by comparing model predictions to rotational distributions of a number of stellar clusters as well as individual field stars. Our model reasonably successfully reproduces the rotational behaviour of stars during the PMS phase to the zeroage mainsequence (ZAMS) spin up, sudden ZAMS spin down, and convergence of the rotation rates afterwards. We find that including coreenvelope decoupling improves our models especially for lowmass stars at younger ages. In addition, by accounting for the almost complete suppression of magnetic braking at slow spin periods, we provide better fits to observations of stellar rotations compared to previous models.
keywords:
Spinorbit Evolution  Angular Momentum Transfer  Magnetic Braking1 Introduction
Angular momentum evolution in lowmass stars is a result of a complex interplay between initial conditions during star formation, the stellar structure evolution, and the behaviour of stellar winds and magnetic fields. During the life of a lowmass star on the mainsequence (MS), angular momentum is lost to the magnetised wind (magnetic braking). By removing angular momentum, winds cause their host stars to spin down with time (e.g., Schatzman, 1962; Kawaler, 1988). Since rotation is the most important parameter that determines the strength of a star’s magnetic dynamo, this spin down leads to a decrease in the magnetic activity of lowmass stars as they age (e.g., Skumanich, 1972; Vidotto et al., 2014).
The study of angular momentum evolution has benefited a lot from recent measurements of rotation rates in starforming regions (e.g., Irwin et al., 2011), young open clusters (e.g., Hartman et al., 2010) and the field stars (e.g., McQuillan et al., 2014; Leão et al., 2015) covering an age range from to about . These results provide a detailed view of how surface rotational velocity changes as stars evolve from the PMS, through the ZAMS to the late MS. And a number of models have been developed to account for these observations (e.g., Irwin et al., 2007; Gallet & Bouvier, 2015; Johnstone et al., 2015). These models generally incorporate a wind braking law, disk locking and coreenvelope decoupling into their parametric models and in order to find an explanation for the observations.
In order to satisfy observational constraints, most of these models have to incorporate three major physical processes: stardisk interaction during the PMS, redistribution of angular momentum in the stellar interior, and angular momentum loss due to stellar winds. Stardisk interaction prevents the star spinning up during the PMS, although it is contracting at a fast rate. Recent theoretical advances have highlighted the impact of accretion/ejection phenomena on the angular momentum evolution of young suns (e.g., Matt et al., 2012; Zanni & Ferreira, 2013). The amount of angular momentum stored in the stellar interior throughout its evolution is usually unknown because observations have so far only revealed surface rotation, except for the Sun (e.g., Hiremath, 2013) and for a few evolved giants (e.g., Deheuvels et al., 2012).
It has recently been possible to gain a qualitative estimate of angular momentum loss rates due to magnetic braking with 2D and 3D numerical simulations of realistic magnetised stellar winds (e.g., Vidotto et al., 2011; Matt et al., 2012; Vidotto et al., 2014). Much of the difficulty to predict stellar wind torque arises from the uncertainty in our knowledge of the magnetic and stellar wind properties (both observational and theoretical). There is still much work needed in order to understand how wind properties depend upon stellar mass, rotation rate, and time.
The aim of the present study is to develop new angular momentum evolution models for a wide range of stars from F to Mtype and from a few to the end of their mainsequence life and incorporate some of the most recent advances.
The structure of this article is as follows: in Section 2, we rederive the parametrised model of Kawaler (1988), update it with available information and present our model for rotational evolution. In Section 3, we introduce our observational sample of over 1500 stars in 13 clusters and a number of individual field stars to derive observational constraints on the rotational evolution of stars. In Section 4, we calibrate the free parameters in our rotational evolution model using the observational constraints allowing us to find values that generate best fit models. Finally in Section 5, we compare the results of our rotational evolution model to observations.
2 Parameterisation of magnetic braking
The angular momentum evolution of isolated lowmass stars is controlled by the balance between three main physical mechanisms: angular momentum removal by magnetised stellar winds (here after “magnetic braking”), the stardisk interaction, and the angular momentum transfer within the stellar interior. In this section we discuss the corresponding model assumptions.
2.1 General formalism
As first proposed by Schatzman (1962), following the pioneering work of Parker (1958a), stars spin down because a coupling between their magnetic field and their ionised wind transfers angular momentum to the wind at the expense of that in the star. The stellar angular momentum lost by stellar wind in unit time can be estimated as a function of the stellar mass loss rate and rotation frequency :
(1) 
where is the stellar radius and is an equivalent radius that is to be determined. It can be thought as the distance after which wind material is no longer in bound to the star.
The radius depends mostly on the strength and geometry of the stellar magnetic field and the launching of the stellar wind. For example, Mestel (1986) proposes that for a radial magnetic field and , while for a dipolar magnetic field and .
In order to express as a function of , and we make two assumptions: (i) the wind velocity at is a fraction of the escape speed at that location,
(2) 
where , and (ii) the wind can be considered isotropic,
(3) 
Solving for it is easy for us to show that
(4) 
Given , this last equation is used to find .
With the above expressions for and , we obtain
(5) 
with for a radial field and for a dipolar field and is the surface escape velocity. The study of Matt & Pudritz (2008) shows that equation (5) has indeed the correct functional form for relatively slowly rotating stars, with the two parameters and which may be fitted to simulations. However, for more rapidly rotating stars the centrifugal acceleration changes the functional form slightly (Matt et al., 2012) to
(6) 
where Matt et al. (2012) find , , . And all quantities are in cgs units. is a coefficient that we will modify later. Equation 6 is equivalent to 5 in the limit that , and . We can thus rewrite equation (1) as
(7) 
2.2 Magnetic field parameterisation
Empirical evidence
The question of how to estimate the magnetic field strength is a difficult one. It has been known since Kraft (1967) that rapidly rotating stars have higher levels of magnetic activity than slowly rotating stars. At a given stellar mass, the magnetic field strength has approximately a power law dependence on rotation rate, such that for slow rotators and saturates at fast rotation. Following the work of Saar (1996) it is usually assumed that the surface average field strength is a strong function of the spin frequency. Wright et al. (2011) find . Furthermore, using the results of recent ZeemanDoppler Imaging (ZDI) studies of over 70 stars, Vidotto et al. (2014) showed that the largescale field strength averaged over the stellar surface, , scales with rotations as .
On the other hand, studies of the magnetic field in Mdwarfs indicate that the largescale magnetic field is inversely proportional to the Rossby number,
(8) 
where is the convective turnover timescale which appears to be a natural scaling factor for the stellar activity and magnetic field. For example, Noyes et al. (1984) showed that the mean level of chromospheric emission is inversely related to the rotation period. Also, Baliunas et al. (1996) found a significant correlation between the magnetic and rotational moments for the Mount Wilson sample. More recently, using a sample of stars in the range of to , for which the largescale surface magnetic field has been mapped using the Zeeman Doppler Imaging (ZDI) technique, See et al. (2017) have shown that the magnetic field strength and mass loss show less scatter when plotted against Rossby number rather than rotation period.
We therefore assume that magnetic field scales with the Rossby number as . See et al. (2017) and Folsom et al. (2016) provide predictions on the dependency of magnetic field strength on Rossby number using ZDI. In their figure 5, See et al. (2017) have plotted magnetic field strength against Rossby number. Using the bisector ordinary leastsquares method and after performing a cut for the critical Rossby number () at they found . Using linear regression we were also able to fit their sample with . In addition, using a different sample and the method, Folsom et al. (2016) have approximated . We consider as a free parameter to be determined.
The convective turnover timescale
The convective turnover timescale is defined as the ratio of the pressure scale height (or the mixing length) to the convective velocity near the base of the outer convective zone. Here, we calculate it using the CESAM stellar evolution code (Morel & Lebreton, 2008a), as described in Appendix A.
Fig. 1 shows how evolves as a function of time for stars of different masses. It is initially large (of order 20) and constant when the stars are fully convective. It then decreases on the premainsequence to become almost constant on the mainsequence. It is smaller for higher mass stars, with a significant drop above which coincides with the gradual suppression of the external convective zone for stars of these masses and higher. The fact that our differs slightly from those of (Landin et al., 2010) is due to a slightly different location for the calculation of . We use the pressure scale height and they use the mixing length of 1.5 times the scale height. The differential behaviour between stars of different masses is very similar, as it should.
In fig. 1, we also compare our results to those estimated empirically by Wright et al. (2011). By necessity, due to the large number of stars required, the latter estimates apply only to mainsequence stars. It is striking that although the is similar within an orderofmagnitude, the range, as a function of stellar mass, is much smaller in the empirical determination than in our calculations or those of Landin et al. (2010).
In this work, we choose to use based on our selfconsistent CESAM calculations. We note that the comparison with the Rossby number in the literature is not necessarily straightforward because of the different definitions and/or calculations used. However, with the exception of the premainsequence evolution phase, the different definitions and calculations lead to which is approximately proportional between different methods. Hence we can use proportionality relations between and other parameters but stress that the absolute values of are not directly comparable.
Magnetic field saturation at fast spin
Based on the slow rotational evolution of young rapidly rotating stars, there can be no serious doubt that both mass loss rate and the magnetic field saturate at fast rotation (e.g., Pizzolato et al., 2003). The reason is that, without saturation, the most rapidly rotating stars would spin down much faster than observed. Thus using a critical Rossby number () we can distinguish between the nonsaturated () and saturated () regimes as
(9) 
An estimate of can be derived by fitting rotational evolution models to the observational constraints. Donati et al. (2008) and Reiners et al. (2009) find . This determination of Romcrit is accurate only within a factor about and in addition, as discussed in Sect. 2.2.2, it is not defined in exactly the same way. Therefore, we treat as a free parameter and, by allowing a range , we look for a value which best fits observations in section 4.
Suppression of braking at slow spin
It has recently been shown by van Saders et al. (2016) that stars more evolved than the Sun rotate faster than previously expected. This indicates that magnetic braking is weaker in intermediateage and old stars. Therefore, they suggest the existence of another Rossby threshold, which we call , above which magnetic braking is suppressed. A possible explanation for this could be a change in field geometry from a simple dipole to a higher order field which could produce weakened braking (Réville et al., 2015). A realistic treatment should thus consider changes to equation (6) arising from this different field geometry, perhaps through a different . In the absence of any knowledge about the behaviour of these fields, we follow van Saders et al. (2016) and adopt a simple reduction of the angular momentum loss rate past this critical Rossby number through a change in the coefficient,
(10) 
Note that the above the critical Rossby number is arbitrary.
Following van Saders et al. (2016), we adopt ,because of our choice of convective turnover timescale , while this value is equal to 2.16 for van Saders et al..
2.3 Parameterisation of the mass loss
The question of how to consider mass loss is often left aside, probably because the change in angular momentum (equation 7) does not depend on mass loss when . This solution was preferred for that reason by Kawaler (1988) and Bouvier et al. (1997). For other an assumption must be made for . One can first note that, for , magnetic braking increases when mass loss decreases, a slightly problematic result at least asymptotically. For , braking increases with mass loss, an intuitively more satisfactory result. In the particular case of the numerical simulations by Matt et al. (2012) the spindown timescale . Thus the dependence on is not negligible.
As known since the pioneering work of Parker (1958b), mass loss is strongly linked to heating of the corona, which in itself is linked to the production of energetic photons and particles. Observations clearly indicate that the X and EUV irradiance of young stars (including the Sun) is orders of magnitude higher when young (e.g. Ribas et al., 2005). This activity is presumably mostly related to the fast rotation of the stars. Separately, Tout & Pringle (1992) show that for fully convective TTauri stars, very strong winds of the order of are expected and are proportional to . Fitting the spin rates of stars on the mainsequence, Johnstone et al. (2015) estimate that the wind mass loss rate scales with stellar parameters as .
See et al. (2017) estimate mass loss rates from magnetic field measurements of young stars and a potential field source surface model. Using their results, we consider mass loss to be correlated with the Rossby number, and like the magnetic field, saturate at fast rotation speed. Therefore, mass loss can be parameterised as
(11) 
See et al. (2017) find ranging between and with the bisector ordinary leastsquares method. Using linear regression, we were able to fit their data with . We treat as a parameter to be fitted within these extreme values.
2.4 Decoupling of the outer convective envelope
Up to now, we have implicitly assumed that the convective envelope is tied to the rest of the star which rotates as a solid body. In fact, as proposed by MacGregor & Brenner (1991) this is probably not the case and this has to be taken into account: the stellar wind acts on the star’s magnetic field lines which originate in the dynamo region, in the convective zone. The rotation state of the stellar interior then depends on the ability and speed at which angular momentum is transferred to the interior. MacGregor & Brenner (1991) proposed that a typical angular mixing timescale of the order of yr was required to explain the small amount of differential rotation observed with seismology between the solar convective zone and radiative interior. Allain (1998) applied this prescription to the case of Mdwarfs and showed the effect to be important to explain the efficient spinup of Mdwarfs in young clusters. We will see that interior angular momentum mixing also plays an essential role in Fdwarfs.
Following MacGregor & Brenner (1991), we assume that any difference in angular momentum between the outer convective zone (characterised by an angular momentum, moment of inertia and angular rotation linked by ) and the inner radiative layer (characterised by ) is mixed with a characteristic timescale . The equations for the evolution of the rotation frequencies in the outer convective zone (cz) and in the inner radiative zone (rad) are (Allain, 1998)
(12)  
and
(13)  
where is the mass of the convective zone and denotes the radius of the radiative layer. We have thus added another free parameter, (the characteristic timescale it takes for angular momentum to travel between radiative and convective zones), which from previous studies is expected to be of the order of yr (Allain, 1998). This timescale agrees both with the lack of differential rotation in our Sun and the fast rotation of some M and Gdwarfs in clusters younger than about .
2.5 Disk locking
The disklocking process arises from the observational evidence that stars that magnetically interact with their accretion disk during the first few Myr of PMS evolution are prevented from spinning up in spite of contracting towards the ZAMS. This behaviour is believed to result from the magnetic interaction between the young stellar object and its accretion disk, even though the details of the process are still to be understood (see Bouvier et al., 2014).
At the beginning, most stars have disks but lose them in a few Myr (see e.g. Haisch et al., 2001). From observations of the TaurusAuriga association, Bertout et al. (2007a) find that the average disk lifetime is a function of stellar mass (from about to ) and approximately equal to .
Following Bouvier et al. (1997) and Gallet & Bouvier (2015), we assume that fastest rotating stars at any age are those which lost their disk rapidly. Conversely, the stars rotating the slowest should have had their disks the longest. Furthermore, we assume that only the outer convective zone (source of the dynamo) is locked by the disk for a timescale which we also treat as a free parameter.
3 Observational Constraints and Models
3.1 Observational Constraints
The set of free parameters introduced in our rotational evolution model above needs to be constrained by observations in order for us to be able to study the distribution of stellar rotation rates at any age. For this reason, we use rotation periods of stars in a number of open clusters from the literature covering an age range of to . Determining the age of young stellar clusters can be very challenging, and in fact, for very young clusters about , age uncertainties may approximately be equal to the actual age estimate (Bell et al., 2013). Fortunately, the change in rotation rates at these very young ages are relatively mild so that modelling does not suffer much from these uncertainties. The ages of older clusters are better known, with uncertainties close to percent. Therefore, we neglect the age uncertainties in the rest of this paper.
Considering the fact that many of the stars in our sample are likely to have companions, one might be concerned that the fraction which are interacting would have a significant effect on angular momentum loss. Binaries may impact rotational evolution via gravitational or magnetic interactions. Stars in very close binaries can exert tidal forces on each other, spinning them up or down more rapidly than predicted for a single star. These systems are also close enough for one star to interact with the other’s largescale magnetic field. And at the earliest evolutionary stages, a companion may truncate the protoplanetary disk, minimizing the impact of magnetic braking and allowing the young star to spin faster than its single counterparts (e.g., Cieza et al., 2009). However, according to Duquennoy & Mayor (1991), only about percent of binaries have orbital periods short enough to be tidally interacting even on circular orbits. When considering eccentric planets, we estimate that another about percent with longer periods may be tidally interacting. In addition, the spin properties of spectroscopic binaries appear not to differ from other stars (see Meibom et al., 2011, concerning the M34 cluster). This indicates that while duplicity may be important in some cases, it does not play a major role in the process.
Cluster  age/Myr  /days  /days  Ref.  

ONC  
NGC2264  
CepOB3b  
NGC2362  
hPer  
Pleiades  
M35  
M34  
M37  
Praesepe  
Hyades  
NGC6811  
NGC6819 
It is generally common in the literature to report the rotation period of clusters as a function of their magnitudes or colours like , , , etc. This is while the rotational evolution models require the mass of the stars. Our sample includes clusters for which the statistical properties of the distributions of rotation rates are shown in Table 1. For four clusters (Praesepe, Hyades, NGC6811 and NGC6819) we have converted the commonly reported colours into masses. For the remaining clusters, the masses have been inferred directly from the references shown within Table 1.
In order to calculate the masses, we first extracted the age and metallicity of each cluster from the literature and using the CMD interface
Unlike clusters, masses are usually directly provided for field stars in the literature and that has been used for the objects in our sample. Table 4 shows the statistical properties of the distributions of rotation rates for the field stars used in our sample.
3.2 Numerical method
We solve our ensemble of equations (7) to (13) using a RungeKutta driver with adaptive step size control. All the quantities which depend on stellar evolution, namely radius, moment of inertia, convective timescale are calculated from a stellar evolution grid which extends from to in steps of and metallicities between to in steps of . This grid was calculated using CEPAM (Morel & Lebreton, 2008b) with a mixing length parameter that was fitted to the Sun. Practically this means that for our solar metallicity mass star reaches a radius of for an effective temperature of K in 4750 Myr. The small 200 Myr offset between the model age and the age of the Solar System is negligible given other sources of uncertainty (see e.g., Bonanno & Fröhlich, 2015).
The separation between the two zones defined as the deep radiative zone and the convective zone in equation (13) was obtained by looking at the first transition from a radiative to a convective zone, neglecting a possible central convective zone.
3.3 Comparison of existing models
At this stage, we are ready to compare the predictions of currently existing rotational evolution models to our observational sample. To perform this task, we implement the rotational evolution models of Bouvier et al. (1997) [B97] and Johnstone et al. (2015) [J15] into our numerical code. Next, we allow each model to evolve for Gyr and finally look at snapshots from their evolution at the same age as the clusters and field stars in our sample to compare how well models can predict the rotation rates. For the case of Gallet & Bouvier (2015) [GB15], because this work uses values of which are massdependent (see 3.3.2), we directly plot their published results. The assumptions of the different models are summarised in Table 2.
Model  Magnetic field  Mass loss  Saturation  Interior mixing  Ref  

B97  None  None  0.5  1  
GB15  to 500 Myr  0.22  1.7 to 8.5  2  
J15  None  0.2177  6.76  3  
This work  Myr  0.2177  6.34  4 

: This is an approximate relation for .
The B97 Model
The model of B97 for the stellar braking of lowmass stars ( ) with an outer convective zone stems from Kawaler (1988)’s formulation. They use a dynamo relation which saturates for fast rotators beyond the spin rate of and is described by the following equations:
(14) 
where is a calibration constant. The critical rotation rate in equation (14) is the only free parameter in B97. The values of the parameters used in this model are listed in table 2.
B97 can be made equivalent to our model by neglecting the term in equation (6), assuming and assuming that rather than . We thus obtain that and as a result . Assuming and G, and using cgs, B97 implies and .
The GB15 Model
Gallet & Bouvier (2013, 2015) develop an angular momentum evolution model for 0.5, 0.8 and 1.0 stars. They adopt Matt et al. (2012)’s prescriptions for the wind braking law and following their work choose and as the input parameters in equation 6. The mean magnetic field and mass loss relations used in their model are adopted from Cranmer & Saar (2011) with slight modifications. They are a function of stellar density, effective temperature, and angular velocity. As shown in Table 2, for angular velocities between 1.5 and 4 times that of the Sun, they find relatively steep dependences on Rossby number, i.e., for the magnetic field strength and for the mass loss rate. They also use a prescription that smoothly joins the nonsaturated and saturated magnetic regimes.
In addition, GB15 fit several other parameters. The value of is allowed to change as a function of mass. The value of the coreenvelope mixing timescale (our ) changes both as a function of stellar mass and initial rotation rate. Thus, for 3 mass bins and 2 rotation rates, 9 additional parameters must be found in addition to , , , , critical spin rate and disk lifetime . As a result, varies between for stars to for stars. The value of range for 10 Myr for fastrotating stars to 500 Myr for slowrotating stars.
The J15 Model
J15 construct a model to describe rotation and wind properties of lowmass mainsequence stars. Like GB15 they adopt Matt et al. (2012)’s prescriptions for the wind braking law and choose and . The mean magnetic field used in their study is adopted from Vidotto et al. (2014). They derive the mass loss rates by fitting their rotational evolution model to the observations by assuming that, in the unsaturated regime, the mass loss rate per unit surface area of a star has powerlaw dependences on its mass and rotation rate. Their model works with six parameters summarised in table 2. An important point to note about the J15 model is that it does not assume a decoupling of internal layers of stars but considers that they rotate as solid bodies. Since the model is limited to stars older than , the consequence of this assumption is that they do not consider the premainsequence phase.
Comparison of Models
Fig. 2 shows predictions of the evolution of stellar spin period made by the above models versus stellar mass. These have been plotted over observations of the rotation periods of the clusters and field stars of our sample at different ages (Tables 1 & 4). The two lines for each model represent the trend of the slow and fast rotators. All three models fitted percent of the stars to estimate their parameters and thus percent of the stars rotating the fastest are not included in their fit. The stars first rapidly spin up due to the contraction during the PMS stage. But after entering the mainsequence at about to , magnetic braking starts to slow the stars down.
Looking at fig. 2, the three models provide a fair fit to the ensemble of clusters that have been observed. Globally, the models do reproduce the observed spin down of stars with age and some of the observed trends in terms of mass. However, for the case of B97 we can see that, starting from two initial spin periods that cleanly encompass the observations at 1 Myr for the ONC, the situation deteriorates with more stars being outside this envelope, particularly for M35 at 150 Myr and for M37 at 550 Myr. B97 tends to explain the slow rotators in M37 and older clusters but not the less numerous fast rotators. Because of the early slow down of the fast rotators, this parameterisation creates poor fits to many of the stars in Pleiades at Myr as well as M35 at Myr.
GB15 obtains a quite good fit for the evolution of lowmass stars at the expense of using many free parameters. It is also important to stress that these fair fits are obtained by assuming a dependence of mixing timescale on rotation rate. However, this is a hypothesis that we view as unnecessary at this point and that we choose not to adopt. The J15 parameterisation (which assumes solidbody rotation) is more successful, including for the Pleiades. It generally explains the spin rates for about 90 percent of the observed stars. However, it also leaves out at least 10 percent of the fast rotators in M35 and M37.
Furthermore, these parameterisations are valid for stars with an outer convective zone but generally fail to reproduce the fastspinning Fdwarfs beyond the Kraft break (Kraft, 1967), as seen for M37, NGC6811 and NGC6819. This is not surprising because they were not specifically designed for that purpose, but it calls for a model that would work for a wider range of stars.
4 Parameter Calibrations
In this section we describe our method to calibrate our model parameters. We adopt the approach of trying to fit nearly all stars seen in clusters rather than focusing on a limited fraction which has usually been done in previous studies. In this way, we hope that the parameterisation can be used to identify more reliably outliers or stars which have been spun up by their exoplanet (e.g., Poppenhaeger & Wolk, 2014; Guillot et al., 2014). Also, using this approach, we are looking at the evolution of single clusters instead of individual mass bins.
Considering the collection of equations introduced in the previous sections and the assumptions made so far, we are left with free parameters which are listed in Table 3. Two of the parameters we do not vary. First, we do not calibrate because we are mainly focusing on clusters and stars at younger ages whereas mostly affects older stars. We adopt . Secondly, the parameter is constant because it has generally a small effect. Therefore, we are left with parameters to vary.
Parameter  Usage 

Calibration to the Sun  
Accounting for rapid rotators  
a  
d  
Timescale of angular momentum transfer  
between radiative and convective zones  
Timescale of keeping stellar rotation  
constant due to the presence of the disk  
Discerning between critical and  
subcritical regimes  
Suppression of magnetic field at old ages 
Because of the large number of parameters and the inhomogeneous data to be used, we proceed as follows. We start from the parameters that, from initial tests, are shown to have higher effects on the evolution of angular momentum. We choose to fit individual clusters by eye and find one optimal solution per parameter.
In the next subsections, we present the results of each parameter search around our preferred values for the other parameters. For simplicity we present these results only for the cluster(s) which led to our choice.
4.1 Calibration to the Sun: and
We use the ensemble of relations defined by equations (7), (9), (11) & (13) to calculate the evolution of our Sun and constrain as a function of the other parameters of the model. is adjusted so that the model correctly reproduces the Sun’s rotation period at its current age.
We adopt an initial rotation period of 8 d as typical of the rotation of TTauri stars (e.g. Irwin & Bouvier, 2009b). The corresponding initial CESAM model has been evolved from an initially extended state for Myr. At that point, it is a fully convective premainsequence star with a radius of .
Fig. 3 shows the resulting spin evolution of the Sun from various sources in the literature and for our preferred model. The details of the parameters of each model are presented in Table 2. For , following the results obtained by Matt et al. (2012), we adopt a value of and, as described earlier because of the small effect it has on the overall evolution of angular momentum, we choose not to vary this parameter.
By construction (through an optimization of the parameter ), all parameterisations reproduce the solar rotation but the evolution histories differ. In the premainsequence evolution phase, the spin period is independent of the parameterisation: the spinup is controlled by the stellar contraction. The evolutionary tracks begin to deviate from each other after 30 Myr around the beginning of the mainsequence. The J15 parameterisation produces the strongest braking for solarmass stars younger than the Sun. In this phase, the GB15 parameterisation yields spin periods which can be shorter by up to a factor . The parameterisation that we present here results in a spin period which can be almost one order of magnitude smaller than J15 around an age of 200 Myr. For older stars, the trends are inverted, with GB15 producing the slowest spinning stars and J15 predicting slightly faster rotation speeds. However, our parameterisation yields still faster rotation speeds because of the assumption of the suppression of magnetic braking when the Rossby number exceeds that of the Sun. This yields a nonsmooth curve for our model in fig.3.
4.2 Disk lock parameter:
It was shown in section 2.5 that the disk lock effect and the value of the parameter controlling it in the model is very important in terms of specifying the initial content of angular momentum in stars. Because magnetic braking has a timescale of at least 10 Myr, is fully decoupled from the other parameters and may be determined independently.
In addition, it was also highlighted in 2.5 that the fastest rotating stars at, say, are those which lose their disks rapidly. As a result we should expect the fastest rotators to have a to . Therefore, following Haisch et al. (2001), we started by adopting a disk lifetime of , corresponding to the mean overall disk lifetime observed in young clusters and then tested different values in its proximity to the slow rotator branch. In addition, for the fast rotators we tested values close to . Bertout et al. (2007b) suggested that there is a dependence between and stellar mass. However, our calculations of spin rates in young clusters do not enable us to distinguish between the Bertout et al. relation and one with a value of that is independent of stellar mass. As a result we simply assume a constant timescale for all masses but acknowledge that for small mass fast rotators, there may be a slight tendency for an even smaller .
Based on fig. 4 we have adopted a value of and for best fit to the fast and slow rotators, respectively.
4.3 Magnetic field and mass loss parameters: , and
According to equations (9) and (11), the magnetic field of the star and its mass loss are controlled by three parameters in our model, which describes the dependence of the magnetic field to the star’s Rossby number, which plays the same role for mass loss and which determines when the magnetic field saturates and magnetic braking is weakened.
In the nonsaturated regime our model predicts that the stars have an angular momentum decay rate of . This shows and are degenerate. In addition, because of the calibration of our models to the solar spin rate, constant depends on which therefore also makes degenerate with and .
In order to explore the extent of this degeneracy we used the results of See et al. (2017) and Folsom et al. (2016) who provide predictions on the dependency of mass loss and magnetic field strength on Rossby number using Zeeman Doppler Imaging. It should be stressed that the estimates obtained for & parameters by See et al. (2017) and Folsom et al. (2016) have been found by adopting a value of . This shows that the parameter is also degenerate. Hence, we vary this parameter alongside & .
On this basis, we choose to test values of , and to explore the degeneracy between these parameters and look for cases that could provide fair fits to our clusters. In appendix B we show how , & affect the models against the stars in the M35 and M37 clusters.
In fig. 5 we explore the degeneracy between parameters , & . It should be noted that since none of the cases in the limit were successful, we do not include them in fig. 5. In each panel we have tested different values of the parameter versus and for different values. The plots show whether or not a certain set of parameters fits our sample at two different ages; 150 Myr and 550 Myr. These two ages have been chosen because they contained the most stars among our clusters and in all cases if a model provided a good fit for these two, it was also successful at the other ages.The red and black colours show the results of the tests at 150 Myr and 550 Myr respectively.
The three symbols represent how well a model could explain the rotation profile of a cluster. Circles represent fits that were fair, i.e. covered most of the stars in a cluster. For example, in fig. 11 the curve represents a fair fit for the 150 Myr cluster. Diamonds show a good fit, i.e. The majority of the stars in a cluster are situated between the two fastest and slowest rotating envelopes while the curves are not far away from the stars either. For instance, in fig. 11 the is considered as a good fit while although the curve covers more stars, it is located far away from the envelope and thus not regarded as a good fit. Finally, crosses show tests which were ruled out.
Fig. 11 to 13 show that the slope at 550 Myr is not steep enough to explain both lowmass fast rotators and the lack of massive fast rotators. Hence, none of the models were accepted as good fits at 550 Myr. The models which created fair fits at both (and thus all) ages are shown with two concentric black and red circles. We find that only in 6 out of 576 cases we can find a successful fit. This also gives an acceptable limit for the range of each parameter.
Fig. 6 compares these successful models at 150 Myr and 550 Myr. The degeneracy between the three parameters , and is apparent. We point out that, for the 150 Myrold cluster, our envelope for fast rotating stars is slightly below the cloud of points. On the other hand, at 550 Myr, we cannot explain a few fast rotating stars below 0.6 and our braking model seems to be slightly too weak for fast rotating stars above 0.8 . On the other hand, our model is quite successful at explaining most of the slow rotating stars at both ages (and, as we will see afterwards, for older clusters). Overall, while our solutions appear to be capturing the dominant features of the spin–mass–age relation, we still have problems for fast rotating stars. This indicates that either the observations are incomplete or more probably that magnetic braking is more complex than assumed here.
Since all the successful models are extremely close to one another, we cannot select one over the others based on these spin–mass–age relations. We therefore choose our preferred model based on the proximity of the parameters to what has been found in the literature. These parameters are , & . We adopt them for the rest of the paper.
4.4 Decoupling of internal layers:
The timescale over which angular momentum is transferred between the inner core and the outer convective zone, the mixing timescale , is a free parameter that we adjust to fit observations at this stage. Several previous studies have included a decoupling timescale in their models and found it to be of order of a few . J15 have focused on ages above and found that they were able to fit their sample without assuming coreenvelope decoupling. We also performed a series of tests to check whether we could indeed describe observations without . In order to do this, we repeated our tests in section 4.3. However, our results showed that a major defect of such models is that they are generally unsuccessful at reproducing a good fit to both young and old clusters with the same parameters. We could not find a set of parameters that described all ages successfully and thus infer that coreenvelope decoupling is a necessary component of our model.
GB15 on the other hand, assume to depend on mass and initial rotation simultaneously. They argue that the reason for adopting a different mixing timescale for the slow and fast rotators, is that angular momentum mixing may depend on the width of the tachocline (between the radiative core and the convective zone) which is influenced by shear. We acknowledge that this is a possibility, but it introduces yet another source of freedom to the model.
Fig. 7 shows the result that we get after adopting GB15’s . It also shows that we were able to fit our sample using a constant for all masses and rotation regimes. We did not attempt to fit all possible combinations of but instead tried fixed and found that 150 Myr improved the fit in the sense that it could explain more stars in the lowmass fast rotators branch at young and intermediate ages. However, we are also getting a larger range of spin periods at old ages. This indicates that the solution is more complex and beyond our model.
Finally, we tested how different orders of magnitude of our adopted would fit observations (fig. 7). An interesting point in the behaviour of the mixing timescale parameter is how either lower or higher values generally predict slower rotation rates.
5 Results: Confronting The Model With Observations
Based on the investigations performed in the previous sections we are now at a stage where we can provide a model with our preferred set of parameters (see Table 2) for predicting the rotational evolution of stars vs. their mass.
5.1 Comparison to Cluster Observations
Fig. 8 summarises the spin period evolution of fast and slow rotators as a function of mass and at different ages. The two solid black lines represent the slow (upper curve) and fast (lower curve) rotator population envelopes. The dotted red lines represent models with intermediate initial spin periods of , & days. Our models provide a good overall fit to the observational constraints.
As previously discussed in 4.2, the spin rates of the outer convective zones of the stars are assumed to be held constant during the first for solarmass fast rotators and for solarmass slow rotators. We have also assumed a linear variation of between these two for the intermediate models in fig. 8.
After the disk disappears, stars spin up due to their PMS contraction until they approach the mainsequence and magnetic braking begins to be important. As a result, rotation periods which are uniformly distributed between about 1 and 16 days in the ONC at 1 Myr start decreasing at 5 Myr in NGC2362, in particular for fast rotators which supposedly lost their disk early and began their spinup. At 13 Myr the highmass stars have spun up significantly due to premainsequence contraction. There are a few outliers with days spin period and a mass of , but these could easily be explained by a or so error on the mass determination. After entering the mainsequence, at the age of about 30 Myr, stars spin down through magnetic braking. Fig. 8 shows that after about 100 Myr slow rotators appear to converge towards a relatively welldefined sequence whereas fast rotators remain scattered relatively significantly. The envelopes for the fast and slow rotators appear to globally capture the evolution of these populations. We see that for the M37 cluster, we slightly underestimate the spread in rotation rates. The difference, however, is small. Johnstone et al. (2015) also report a similar issue in their predictions for this cluster.
We acknowledge that between 600 Myr and 940 Myr, our model seems to predict a slightly too fast envelope for fast rotators. At older ages, up to 2.5 Gyr, the slow rotators branch of our model reproduces the observations well. However, there are very few fast rotators. This raises the question of whether our envelope for fast rotators is correct. In fact, as shown in fig. 2, the B97 and GB15 parameterisations tend to have a much faster convergence so that, after about 500 Myr, the initial conditions have been lost. It is not clear, however, if the observations require that. The clusters which have been observed at these relatively old ages contain somewhat fewer stars and it seems plausible that the absence of these fast rotating stars is due to the fact that they are rarer than the slowrotating ones. In addition, we can see that the Hyades does possess a number of lowmass, fast rotating stars which are not explained by the B97 and G15 parameterisations. We hence believe that our parameterisation provides a plausible conservative estimate for rapidly rotating stars.
Furthermore, our model also captures the essence of the Kraft break the fact that Ftype stars, above about , rotate much more rapidly than Gtype stars. This was not the case of the B97, G15 and J15 parameterisations (see fig. 2). The agreement between model and observations is not as good for these massive stars as for the lowermass stars however.
Finally, we point out that the increase in spin periods seen in the NGC6819 cluster is due to stars leaving the mainsequence. Because they expand, conservation of angular momentum requires that their spin rate decreases.
5.2 Bimodality in Rotation Rates
It has long been known that in most of the intermediateage clusters such as M35 & M37, many of the stars lie on one of two principal sequences, identified as C and I by Barnes (2003). These sequences are not necessarily clearly defined with a significant number of stars found in between them.
The Isequence consists of relatively slow rotators. With increasing age, the Isequence becomes more tightly defined, and moves to longer periods. Stars on the Csequence are rapid rotators. In young clusters ,, the Csequence has members that span the mass range of solarlike stars. In older clusters the Csequence fades away, starting with the highermass stars. By the age of the Hyades () it seems as if the only remaining fast rotators are lowmass stars. However, as discussed previously, this may be an artefact from the relatively limited number of stars associated with these old clusters.
Our model reproduces this feature qualitatively. This is because slow rotators with have an angular momentum decay rate whereas for fast rotators with , (assuming ). Extremely fast rotators, characterised by , are found on the PMS and for spin periods shorter than about a day. For these, the magnetic braking is even weaker, .
Thus, the stronger dependence of magnetic braking on for slow rotators leads to a faster convergence towards a relatively welldefined relation. On the other hand, the weaker dependence on for fast rotators with already maximal field strength implies that the scatter in spin rates is conserved for longer. We can see, however, by focusing on the intermediate initial spin rates (red curves) in fig. 8 that although our model captures the broad behaviour of these populations, it fails to reproduce the fact that the slow rotators appear to rapidly outnumber the fast rotators. This probably requires a different physics from that assumed here (see e.g., Brown, 2014).
5.3 Application to individual field stars
Our model parameters have been calculated based on cluster observations and a calibration to the Sun. We now turn to the comparison of our model predictions to those observed on individual field stars, listed in Table 4. Most of those were obtained by an analysis of Kepler observations combining photometrically determined spin periods to precise determination of age, mass and metallicity from astroseismic modelling (van Saders et al., 2016). The ensemble spans a range of age between about 1 and more than 10 Gyr, masses between and and [Fe/H] between and . It is thus well suited for our purposes.
Star  age [Ma]  [days]  Ref.  

observed  
KIC10644253  \textcolorgreen✳  \textcolorgreen✳  1  
KIC9139151  \textcolorgreen✳  \textcolorgreen✳  1  
KIC10454113  \textcolorred✷  \textcolorred✷  1  
KIC3427720  \textcolorgreen✳  \textcolorgreen✳  1  
KIC5184732  \textcolorgreen✳  \textcolorgreen✳  1  
KIC10963065  \textcolorgreen✳  \textcolorgreen✳  1  
AlphaCenA  \textcolorgreen✳  \textcolorgreen✳  3  
Sun  \textcolorgreen✳  \textcolorgreen✳  2  
KIC8006161  \textcolorgreen✳  \textcolorgreen✳  1  
KIC6196457  \textcolorgreen✳  \textcolorgreen✳  1  
KIC11401755  \textcolorgreen✳  \textcolorgreen✳  1  
70OphA  \textcolorgreen✳  \textcolorgreen✳  3  
KIC6116048  \textcolorgreen✳  \textcolorgreen✳  1  
KIC6521045  \textcolorgreen✳  \textcolorgreen✳  1  
KIC10586004  \textcolorgreen✳  \textcolorgreen✳  1  
KIC9955598  \textcolorgreen✳  \textcolorgreen✳  1  
KIC11244118  \textcolorred✷  \textcolorred✷  1  
AlphaCenB  \textcolorgreen✳  \textcolorgreen✳  3  
16CygB  \textcolorred✷  \textcolorgreen✳  1  
16CygA  \textcolorred✷  \textcolorgreen✳  1  
KIC9098294  \textcolorred✷  \textcolorgreen✳  1  
KIC8349582  \textcolorgreen✳  \textcolorgreen✳  1  
KIC3656476  \textcolorred✷  \textcolorgreen✳  1  
KIC7871531  \textcolorred✷  \textcolorgreen✳  1  
55Cnc  \textcolorred✷  \textcolorred✷  4 
In order to compare observations and models, we proceed as follows. For each star, we calculate a nominal model for the stellar age, mass and metallicity and an initial spin period of d (taken as the geometrical mean of the values observed in the ONC). We also calculate models for ages, masses and metallicities which differ by around the observational values, and models for the nominal stellar parameters but initial spin periods of and d, respectively. We note the difference between the spin periods thus obtained and those of the nominal model , , and respectively. We then calculate our model uncertainties from a quadratic mean, . The model generates asymmetries so negative and positive uncertainties are calculated separately. We then compare the observed and modelled values within uncertainties. A model is successful for one star if the observed and modelled values overlap. It is unsuccessful otherwise.
Following the approach taken by van Saders et al. (2016), two kinds of model are calculated. A standard one with a magnetic braking that continues unhindered at large Rossby numbers (slow spin rates) and one that is stalled for . As shown in Table 4, the former model is successful for stars (excluding the Sun, which is used for the calibration). It fails mostly for old, slowly spinning stars. The latter model that we adopted is successful for stars, which thus points to a tapering of magnetic braking for slowly spinning stars, in line with the results of van Saders et al. (2016).
The three cases that still present a mismatch with our preferred model are:

KIC10454113: This midlife, solarmetallicity, F star (with an estimated K based on our model calculations) is rotating much more slowly ( days) than expected from our models which predict spin periods between about 2 and 8 d. The difference is significant and to be compared to our results for the NGC6819 cluster which also showed that two stars were rotating more slowly than our predictions. At this point, we do not know how to explain, on one hand, the fast rotators in NGC6811 (0.94 Gyr) and, on the other hand, KIC10454113 and the slow rotators in NGC6819 (2.5 Gyr).

KIC11244118: Another star that our model has difficulty explaining is this 6.4 Gyr F star, with a metallicity of (+0.35) which is rotating faster ( days) than our model predictions (30 to 41 days). We notice that this star is characterised by .

55Cnc: This is a binary system consisting of a Gtype star (55 Cancri A), known to host at least five planets, and a smaller red dwarf (55 Cancri B). The main star spins in 30 to 48 days, faster than the model prediction of 56 to 65 days, even when accounting for the suppression of magnetic braking after .
The engulfment of a hot Jupiter (e.g. Poppenhaeger & Wolk, 2014; Guillot et al., 2014) could be invoked to explain the anomalously fast rotation of KIC11244118 and 55Cnc. However, given that the list is limited to only 24 stars and given that hot Jupiters are quite rare, this is highly unlikely. For KIC10454113, the fact that the inferred spin period is two to three times longer than the model predictions is even more puzzling. We are therefore not able to explain the spin rates of these stars and recommend further studies, in particular to confirm their characteristics.
The special case of Cen A & B illustrates the importance of a precise characterization of the targets. These two stars are known for their proximity to the Sun and in the case of Cen B for the possible presence of an Earthmass planet in a short orbit that has recently been invalidated (Rajpaul et al., 2016). Delorme et al. (2011) indicate, from astroseismic constraints, an age of Gyr which yields a spin period clearly in excess of the observations for both stars. However, the detailed study of Bazot et al. (2012) on Cen A yields a younger age of Gyr. With this age, our preferred model reproduces the observed spin periods of both components.
6 Conclusion
We have presented a parametric semiempirical model for the magnetic braking of stars with masses between 0.5 and 1.6 and from premainsequence to the end of their mainsequence lifetime. The parameters prescribe the stellar mass loss rate and magnetic field amplitude as a function of the physical characteristics of the star obtained from a stellar evolution model and of its spin rate. In the spirit of previous studies (B97, GB15), we also assume that the spin of the outer convective layer of stars is maintained constant by their circumstellar disks and that the slowest (fastest) spinning stars are those which had a disk for the longest (shortest). Two critical Rossby numbers define when the magnetic field saturates and stops increasing with increasing spin rates (for ) and when the magnetic field changes form, leading to a decrease in angular momentum loss efficiency (for ). Contrary to previous models which parameterised the magnetic field intensity and mass loss as a function of the rotation rate, we tie them to the star’s Rossby number calculated at the base of the outer convective zone. This feature enables applying the same model to a wide variety of stellar types. All models are calibrated to the Sun, in order to reproduce a d spin period after Gyr of evolution of a solar metallicity star.
The model parameters were determined from an analysis of the observation of young clusters with ages between 1 Myr and 2.5 Gyr. For nonsaturating magnetic fields , we determined that and led to a fair fit of the observations for a wide range of stellar types (M to F). We found that the suppression of magnetic braking in F stars is naturally explained by the shortening of the convective turnover timescale at the base of their (small) outer convective zone. We also investigated the degeneracy between the parameters controlling magnetic field strength and mass loss rate (parameters , & ). We found that in only 6 out of 576 cases we could find a successful fit. For instance, the following values could all provide fair fits: , , ; , , ; & , , . Therefore, our parametric model provides a way to analyse the spin rate of stars and to identify outliers which may have been spun up or down by external processes.
One of the consequences of our parameterisation is that the convergence in spin rates with age is slower than generally assumed. At the age of the Sun, we predict a spread of rotation rates of about 10 percent that depends on the angular momentum acquired by stars at the end of the circumstellar disk phase. This spread increases at lower masses.
Following van Saders et al. (2016), we could confirm that stars of old age and slow rotation rates see their braking suppressed. This is plausibly a consequence of a change in the magnetic field geometry but should be investigated further.
Although our model provides a fair fit to the spin evolution of stars between and , we acknowledge that it does not provide the final answer to the problem. For instance a remaining issue that needs to be addressed is that our model predicts too many fast rotators at older ages which is not in line with observations.
We are providing our model both for fast and slowrotating stars as a function of mass and age. Our model is thus a useful basis to analyse and predict stellar spin periods. A more exhaustive exploration of the ensemble of solutions would be desirable in order to understand and constrain the physical mechanisms at play. Beyond that, progress on the evolution of the spin rates of stars will require a combination of dedicated studies and observations of spin rates and magnetic field properties of stars as a function of age, mass and metallicity.
Finally, although many of the stars in our sample are likely to have companions, as we showed in section 3.1 binarity does not play a major role for the purposes of this study.
Acknowledgements
We thank the anonymous referee for providing useful comments on the paper that have helped us to improve the original manuscript. We also thank Colin Johnstone, Habib G. Khosroshahi and Raphaël Raynaud for their comments and help on the paper. We acknowledge support from the Centre for International Scientific Studies and Collaboration (CISSC) and from the Cultural Section of the French Embassy through a Gundishapur fellowship, as well as the OCA BQR and the IPM School of Astronomy.
Appendix A On the convective turnover timescale
a.1 Calculation method
In the present work, we used , the ratio of the pressure scale height to the convective velocity calculated with the mixing length approach. We used the CESAM stellar evolution code (Morel & Lebreton, 2008a), and these quantities were calculated at half a pressure scale height over the base of the convective zone.
The definition of is based on the idea that dynamo action should take place at the base of the convective zone (e.g., Kim & Demarque, 1996), but where exactly this should be measured is vague. It is thus rather the differential behaviour, from one star to the next, that is important to capture. In that sense, the parameter which affects most the value of is the depth of the convective zone itself. This is obvious in Figure 9, which shows that is almost uniquely defined by the relative mass of the outer convective zone and is almost insensitive to stellar mass, metallicity and age (once the dependence on has been removed).
In fact, can be well fitted for a wide range of premainsequence and mainsequence evolutions by the following relation:
(15)  
where the value of is in seconds and denotes mass of the external convective zone divided by the total stellar mass.
a.2 Mainsequence approximation
Stars on the premainsequence see their physical properties evolve rapidly and as a consequence, becomes a strong function of age. However, when they get on to the mainsequence, the time dependence becomes much weaker, and, as shown by Figure 10, a simple relation then links , stellar mass and stellar radius:
This relation holds at least up to 1.0 M. For higher stellar masses, however, the retreat of the external convective zone implies that drops, which is at the heart of the existence of the Kraft break, i.e., a much weaker magnetic field and magnetic braking for Ftype stars than Gtype stars.
Appendix B Exploring the magnetic field and mass loss parameter space
In this appendix we discuss how the three parameters , & affect model predictions. Figures 11 through 13 show how various values for , and affect the spin period evolution models against the stars in the M35 (150 Myr) and M37 (550 Myr) clusters.
In figure 11 we have varied from to . It should be noted that since values above were completely unsuccessful at fitting the observations they have been removed from plots. Small values of result in a too fast aggregation of the two rotating regimes, while bigger values are slightly more successful at reproducing the slow rotators branch. A value of successfully covers both the fast and slow rotators and seems to create the best fit to the observations. This is in close agreement with the usually adopted value of in the literature.
In addition, the results of our tests for the parameter are shown in figure 12. The results differ quite significantly on the choice of : Small values below 1.2 result in a lower envelope that is detached from the cloud of fastspinning stars. Large values above 1.6 yield results in contradiction with the existence of lowmass fast spinning stars and also predicts a slow rotation envelope that is detached from the cloud of points for smallmass stars. Based on these results, we choose which is between the two values predicted by Folsom et al. (2016) and See et al. (2017) who find .
In figure 12, we also show our two critical Rossby numbers and . We can see in the case of these two young clusters (and also for other clusters) that for a majority of stars, the relatively slow rotators, are located between these two lines. The fast rotators are, however, characterised by and are thus in a regime in which the magnetic field saturates (i.e., it is not proportional to the spin rate). This feature is crucial to account for the presence of fast rotating solartype stars even at relatively old ages.
Finally, figure 13 shows the results of testing parameter . It can be seen that this parameter affects all stellar masses but has a higher effect on fast rotators. Based on this figure we choose which is close to the predicted value by See et al. (2017) (d=1.49) and also in agreement with the theoretical discussion of Tout & Pringle (1992) who predicted mass loss to be a function of .
Appendix C Models
We provide our rotational evolution models for stellar masses to and with five initial rotational periods; , , , and days in electronic format. The models start from zero to 10 Gyr. Each model contains the following parameters:
AGE: age of the star [],
H: the total angular momentum [],
PSPIN1: rotation period of the primary [],
M1: mass [],
R1: radius [],
MTI1: moment of inertia [