Rotational evolution of slowrotator sequence stars
Key Words.:
Stars: rotation – Stars: evolution – Stars: latetype – open clusters and associations: general – open clusters and associations: individual: Pleiades, M34, M35, M37, Hyades, Praesepe, Coma Berenices NGC6811 NGC6819Abstract
Context: The observed relationship between mass, age and rotation in open clusters shows the progressive development of a slowrotator sequence among stars possessing a radiative interior and a convective envelope during their premain sequence and mainsequence evolution. After 0.6 Gyr, most cluster members of this type have settled on this sequence.
Aims: The observed clustering on this sequence suggests that it corresponds to some equilibrium or asymptotic condition that still lacks a complete theoretical interpretation, and which is crucial to our understanding of the stellar angular momentum evolution.
Methods: We couple a rotational evolution model, which takes internal differential rotation into account, with classical and new proposals for the wind braking law, and fit models to the data using a Monte Carlo Markov Chain (MCMC) method tailored to the problem at hand. We explore to what extent these models are able to reproduce the mass and time dependence of the stellar rotational evolution on the slowrotator sequence.
Results: The description of the evolution of the slowrotator sequence requires taking the transfer of angular momentum from the radiative core to the convective envelope into account. We find that, in the mass range – , the coreenvelope coupling timescale for stars in the slowrotator sequence scales as . Quasisolid body rotation is achieved only after – Gyr, depending on stellar mass, which implies that observing small deviations from the Skumanich law () would require period data of older open clusters than is available to date. The observed evolution in the – Gyr age range and in the – mass range is best reproduced by assuming an empirical mass dependence of the wind angular momentum loss proportional to the convective turnover timescale and to the stellar moment of inertia. Period isochrones based on our MCMC fit provide a tool for inferring stellar ages of solarlike mainsequence stars from their mass and rotation period that is largely independent of the wind braking model adopted. These effectively represent gyrochronology relationships that take the physics of the twozone model for the stellar angular momentum evolution into account.
Conclusions:
1 Introduction
In the past decade, stellar rotation has received a great deal of attention, both because of its promising potential as a precise stellar age estimator (Barnes 2007, 2010; Meibom et al. 2015) and the large amount of good quality data increasingly available, either as a byproduct of planetary transit searches or from dedicated stellar rotation observational programs (see, e.g. Herbst & Mundt 2005; Hodgkin et al. 2006; von Braun et al. 2005; Messina 2007; Messina et al. 2008, 2010). Rotation periods of young stars can in principle be derived with a precision of 10 using photometric monitoring, with surface differential rotation being probably the most important limiting factor for individual old stars (see, e.g. Epstein & Pinsonneault 2014).
The realization by Barnes (2003) of the existence of a welldefined sequence of slowly rotating stars in the young open clusters colorperiod diagram has provided a criterion for establishing an empirical relationship between mass, age and rotation, which can be used to infer the age of stars that are sufficiently old to belong to that sequence. Barnes (2003) effectively disentangled the convergence of stellar periods towards a unique slowrotator sequence, seen as a progressive reduction of the scatter that is due to the presence of fast rotators (from very young clusters to the Hyades), from the rotational evolution of stars in this sequence, which dominates the colorperiod diagram from the Hyades onwards. Subsequently, Barnes (2007) provided a calibration of his rotationmassage relationship using the observed rotationmass distribution in open clusters and the solar age, rotation period and . More recently, Barnes & Kim (2010) provided evidence of a connection between the empirical rotational period functional dependence on and the convective turnover timescale . Based on this, Barnes (2010) derived a nonlinear model formulated in terms of the Rossby number () and two dimensionless empirical constants, which describes the rotational evolution towards the slowrotator sequence.
In this work we focus on the rotational evolution of stars already on the slowrotator sequence and study how well the current models are able to reproduce such evolution. We focus on the age interval between and Gyr, for which data with the richest details is available, providing strong constraints on the models. The observed clustering on this sequence in the massageperiod space suggests that it corresponds to some equilibrium or asymptotic condition that still lacks a complete theoretical interpretation and may be crucial to our understanding of the stellar angular momentum evolution.
Data that has become available in recent years has revealed an evolutionary behavior of the slowrotator sequence that seems difficult to reconcile with the original Barnes (2003) idea of a factorization of the time and mass dependence. Notably, in their study of rotation in M35 and M34, Meibom et al. (2009, 2011b) found that while Gtype stars on the slowrotator sequence spin down with a rate close to the Skumanich law, Ktype stars spin down more slowly.
On the theoretical side, recent work concentrated on improving the description of the magnetic braking due to coupling between the stellar wind and magnetic fields generated by an internal stellar dynamo. Matt et al. (2012) performed twodimensional axisymmetric magnetohydrodynamic simulations to compute steadystate solutions for solarlike stellar winds from rotating stars with dipolar magnetic fields, from which they derived a semianalytic formulation for the external torque on the star. Gallet & Bouvier (2013) derived the angular momentum loss rate using the Matt et al. (2012) equation for the Alfvén radius together with a dynamo prescription that relates the stellar magnetic field to stellar rotation, as well as a wind prescription that relates the massloss rate to the stellar angular velocity, obtained from current theory and calibrated to the presentday Sun. Given the current uncertainties on the stellar magnetic fields and on the wind massloss, however, Matt et al. (2015) derived a scaling for the dependence of the stellar wind braking on Rossby number and an empiricallyderived scaling with stellar mass and radius.
In this work, we include these various wind braking laws in the classical twozone description of the stellar rotational evolution and investigate to what extent they reproduce the slowrotator sequence evolution. The models are fitted to the data by means of a Monte Carlo Markov Chain (MCMC) method tailored to the case at hand. Using the parameters derived by the MCMC model fitting to the data, we derive period isochrones and tracks from 0.1 to 5 Gyr that are largely independent of the specific wind braking law adopted. These effectively represent gyrochronology relationships that take the effects of differential rotation in the period evolution into account.
The paper is organized as follows. In Sect. 2 we discuss the data used in the present analysis. In Sect. 3 we present an empirical description of the slowrotator sequence, which also illustrates the current stateoftheart, together with a novel nonparametric empirical fit of its evolution. In Sect. 4 we present the MCMC fitting method of twozone rotational evolution models with different wind braking prescriptions. In Sect. 5 the results of the MCMC fitting to the data are presented and discussed. The conclusions are in Sect. 6.
2 Data
Name  Age (Myr)  references 

Pleiades  120  Hartman et al. (2010) 
M35  150  Meibom et al. (2009) 
M34  220  Meibom et al. (2011b) 
M37  550  Hartman et al. (2009) 
Hyades  600  Delorme et al. (2011) 
Praesepe  600  Delorme et al. (2011) 
Coma Berenices  600  Collier Cameron et al. (2009) 
NGC6811  1000  Meibom et al. (2011b); Janes et al. (2013) 
NGC6819  2500  Meibom et al. (2015) 
For this work we use literature data of stars in open clusters that show a welldefined slowrotator sequence, easily identifiable by eye as an overdensity in the periodcolor diagram, also known in the literature as the “Isequence” (Barnes 2003) or “ridge” (Kovacs 2015). The age at which this sequence begins to form is still rather undetermined, but it is evident that this is already well defined at the age of the Pleiades and present in all older open clusters. A summary is reported in Table 1. The data set includes clusters from 0.1 to 2.5 Gyr for which measurements in the relevant range are readily available. The age interval is sampled in an optimal way, in the sense that data for open clusters not considered here would not improve the age sampling in a significant way. Since our method relies on an empirical fitting of the slowrotator sequence to start with (see Sect. 3), data for field stars are not considered. Young loose stellar associations are also not considered in this work.
Rotational periods for the Pleiades are taken from Hartman et al. (2010). This data set includes 241 single stars in common with the compilation of probable Pleiades members by Stauffer et al. (2007), which provides near and midinfrared photometry together with a compilation of homogenised Johnson and magnitudes. For the Pleiades, following Stauffer et al. (1998) we adopt an age of 120 Myr.
Meibom et al. (2009) provide rotational periods and UBVRI photometry for 310 members of M35. The age of M35 has been estimated to 150 Myr (von Hippel et al. 2002), 175 Myr (Barrado y Navascués et al. 2001), and 180 Myr (Kalirai et al. 2003). Following Meibom et al. (2009), we adopt an age of 150 Myr.
Rotational periods and reddening corrected for M34 are taken from Meibom et al. (2011b). The sample contains periods for 118 stars. From the rotational period of slowrotator sequence stars, adopting the Skumanich law, Meibom et al. (2011b) estimated an age of 240 Myr, which we also adopt here.
Hartman et al. (2009) provide a clean sample of rotational periods for 372 stars in M37. Photometry for this cluster is from Kalirai et al. (2001). Following Hartman et al. (2008) we adopt an age of Myr.
The rotational period data set of Delorme et al. (2011) includes 56 members of the Hyades and 35 members of Praesepe. They also provide photometry collected from the literature. Following the discussion in Delorme et al. (2011), we adopt an age of 600 Myr for both clusters.
Collier Cameron et al. (2009) present a photometric survey of rotation rates in the Coma Berenices cluster containing data for 27 members with available photometry. Following their discussion, we adopt an age for Coma Berenices equal to the Hyades and Praesepe age (600 Myr) within the uncertainties.
Meibom et al. (2011a) provide rotational periods for 71 cluster members of NGC6811. Photometry for this cluster is provided by Janes et al. (2013). Following Meibom et al. (2011a, and references therein) we adopt an age of 1 Gyr.
Rotational periods and photometry for NGC6819 is taken from Meibom et al. (2015), who derived a gyro age of 2.5 Gyr in agreement with the estimate by Jeffries et al. (2013).
It is worth noting that the cluster ages adopted here do not differ significantly from the isochrone ages adopted by Kovacs (2015). In the present work, isochrone and rotational ages are assumed to coincide.
3 Empirical description of the slowrotator sequence evolution
3.1 Empirical models
The main assumption in Barnes (2003, 2007) is that rotation periods on the slowrotator sequence can be represented empirically by a relationship of the form
(1) 
where is the stellar age and the colour is taken as a proxy of the stellar mass. Other authors fitted the function to different data sets (e.g. Meibom et al. 2009; Mamajek & Hillenbrand 2008), maintaining the same form and factorisation as in Barnes (2007). Both Barnes (2007) and Mamajek & Hillenbrand (2008) derived their colorperiod relationship on a few well studied clusters (e.g. Persei, Pleiades, Hyades) while the fit by Meibom et al. (2009) is based on a single cluster (M35).
Barnes (2003) proposed that the mass dependence of the rotation periods in open clusters can be simply attributed to the moments of inertia of either the whole star or of the outer convection zone. This was later found inadequate by Barnes & Kim (2010), and led Barnes (2010) to the formulation of a different nonlinear empirical model, which provides the gyroage of a star as
(2) 
where and are two dimensionless constants, constrained by the observations, and is the initial rotation period, taken as the period on the zeroage main sequence. In the slowrotator sequence limit we are concerned with, this semiempirical formula gives
(3) 
which predicts that, asymptotically, with and . The Barnes (2010) model gives a description of the massdependent evolution of fastrotating stars to the slowrotator sequence. An alternative description is the metastable dynamo model (MDM) of Brown (2014), which assumes a stochastic transition between two different dynamo regimes and is rooted in an original idea of Barnes (2003).
The colorperiod diagrams for the clusters studied here are shown in Fig. 1. The stellar mass is estimated using the colormass relationship of Barnes & Kim (2010), based on the effective temperaturecolor transformation of Lejeune et al. (1997, 1998). Periods are scaled to the square root of age (in Myr), so that data points would overlap in diagrams at different ages if the rotational evolution were exactly Skumanichtype. The slowrotator sequence is easily identifiable in all clusters, with the only exception of the Pleiades and M35, for which some ambiguity may arise (see Sect. 3.2). From the comparison of the observed slowrotator sequence at different ages it is evident that:

Stars with slow down at a slower rate than predicted by the Skumanich law until Gyr;

From to Gyr, stars slow down at a faster rate than predicted by the Skumanich law, except for the NGC6811 stars in the mass range;

There is no evidence of stars of mass on the slowrotator sequence at ages younger than Gyr.
The current observed evolution of the slowrotator sequence from to Gyr shows that the factorization expressed by Eq. (1) does not hold, at least before Myr, because such a relationship would imply that the shape of the slowrotator sequence is maintained throughout its evolution.
The data available contains periods for stars on the slowrotator sequence with mass as low as for ages 0.6 Gyr, but no data for stars with is available at later ages.
In Fig. 1 we also report the empirical functions obtained by Barnes (2007), Meibom et al. (2009), and Mamajek & Hillenbrand (2008), as well as the inversion of Eq. (2) from Barnes (2010) and its asymptotic limit in Eq. (3). The comparison with observations shows that such relationships can describe the slowrotator sequence only on limited age and mass ranges.
With the exception of stars of in NGC6811, the slowrotator sequence lies above the Gyr sequence in the vs. diagram. In the framework of the twozone model (Sect. 4.1), this behavior could be attributed to the transport of angular momentum from the stellar core to the envelope in the earlier evolution after the ZAMS. We shall verify this in Sect. 4.
3.2 Nonparametric empirical fitting
cluster  age  

1.10  1.05  1.00  0.95  0.90  0.85  0.80  0.75  0.70  0.65  0.60  
Pleiades  120  2.09  2.75  3.53  4.42  5.39  6.41  7.42  …  …  …  …  0.35  0.63 
M35  150  3.12  3.62  4.16  4.82  5.64  6.46  …  …  …  …  …  0.76  0.80 
M34  240  3.18  3.92  4.72  5.60  6.57  7.65  8.65  9.48  10.11  10.62  …  0.53  0.82 
M37  550  5.21  5.91  6.50  7.08  7.70  8.42  9.32  10.32  11.19  11.91  12.43  0.94  0.68 
Hyades  600  …  …  7.62  8.35  9.07  9.82  10.64  11.57  12.26  12.72  13.04  0.51  0.52 
Praesepe  600  …  6.53  7.36  8.16  8.90  9.55  10.02  10.48  11.09  11.77  …  0.63  0.89 
ComaBerenices  600  6.68  7.49  8.12  8.70  9.29  10.07  10.75  11.07  11.23  11.37  …  0.52  0.26 
HydesPraesepeComa  600  6.12  6.95  7.65  8.31  8.99  9.68  10.38  11.20  11.77  12.22  12.58  0.73  0.68 
NGC6811  1000  7.07  8.58  9.90  10.93  11.30  11.03  10.64  …  …  …  …  0.46  0.45 
NGC6819  2500  15.46  17.41  19.12  20.75  21.87  22.07  …  …  …  …  …  0.88  0.92 
In order to overcome the limitations of the empirical modeling discussed in the Sect. 3.1, we perform a nonparametric fit to the slowrotator sequence. To this end, we need a criterion for selecting the stars belonging to this sequence. For the older clusters, the selection is quite simple as only some outliers and a few remaining fastrotators need to be excluded. For the younger clusters, on the other hand, the slowrotator sequence overdensity is still easily identifiable in the colorperiod diagram, but the separation with stars approaching the slowrotator sequence is somewhat arbitrary and prone to subjective choices. Slowrotator sequence membership may become even more illdefined if, following Barnes (2010), we consider the sequence as an asymptotic limit.
We note, however, that nonparametric fits to the older slowrotator sequences produce residuals whose distribution can be approximated with a normal distribution. To trace back the sequence to the younger clusters with a criterion as homogeneous as possible, we set the width of the slowrotator sequence as the maximum width that produces a distribution of residuals with a normality probability (see below) of at least 30% (see Fig. 2 and Table 2). We use the local polynomial regression fit (LOESS, Cleveland et al. 1992) with smoothing parameter between and , and a Pearson normality test (Moore 1986; Thode Jr. 2002), as implemented in the statistical package R. In order to improve the fitting at the higher and lower end of the sequence, where the curvature is higher and a nonparametric fit would require a much denser data set, we use Eq. (3) as a normalization function, adjusting in order to mimic the slowrotator sequence at Gyr.
Our nonparametric fits are shown in Fig. 1, where the stars selected as members of the slowrotator sequence are also highlighted. The nonparametric fit is used to derive periods for a grid of stellar masses, reported in Table 2 together with the standard deviation and the normality test probability. The comparisons of the empirical distribution functions of the residuals with normal distributions are shown in Fig. 2.
Our selection of stars belonging to the slowrotator sequence in younger clusters is more restrictive than that adopted, for instance, by Barnes (2010). It should be considered, however, that the aim of the present work is to study in detail the rotational evolution of the slowrotator sequence, under the assumption that this represents some equilibrium state or asymptotic behavior, deliberately ignoring fast rotators or stars still approaching the slowrotator sequence. We also ignore periods much larger than the periods of the slowrotator sequence, assuming that they are either incorrect or that these stars are not cluster members. In other words, we assume that the peak of the overdensity of slow rotators in the colorperiod diagram corresponds to the maximum probability of having reached the equilibrium/asymptotic state and that stars in or near such a state have periods distributed normally around the peak. Such an approximation has the advantage of allowing a characterization of the width of the sequence using the standard deviation, which can then be used as input in our parametric fit described in Sect. 4. It turns out that such an approximation describes remarkably well the distribution in the older clusters, and that, by an appropriate selection, it can be extended to the younger clusters as well. For our purposes, it can safely be assumed that the standard deviation of the sequence includes both observational uncertainties and the intrinsic physical width due to the ongoing convergence onto the slowrotator sequence.
4 Twozone models Monte Carlo Markovchain fitting to the data
4.1 Twozone rotational evolution models
In our modeling of the rotational evolution of solarlike stars (whose structure is characterized by an inner radiative region, surrounded by a convective envelope), we adopt the theoretical framework of twozone models (TZMs; see, e.g. MacGregor & Brenner 1991; Keppens et al. 1995; Allain 1998; Spada et al. 2011). The main assumption of the TZM is that the radiative interior and the convective envelope of the star rotate as rigid bodies, with angular velocities and , respectively; the rotational state of the whole star at a given time is thus completely specified by these two quantities.
During the premain sequence, the rotational evolution is dominated by evolutionary changes to the stellar structure (i.e., overall contraction, appearance and growth of a radiative core), and the interaction with the circumstellar disc. From the zero age main sequence onwards, the interior structure changes very little, and the rotational evolution essentially depends on the interplay between the angular momentum loss at the surface, due to the magnetized stellar wind, and the angular momentum exchange between core and envelope. These main physical ingredients are included in our models as follows:

The radiative core grows at the expenses of the surrounding envelope, thus subtracting from it the angular momentum (Allain 1998): , where and are the mass and the radius of the core, respectively (in the following, we adopt the dot notation for derivatives with respect to time);

The stardisc interaction is assumed to be very efficient, capable of keeping the surface angular velocity constant during the entire disc lifetime, (the socalled disclocking approximation, see Koenigl 1991);

The rate of angular momentum loss from the surface due to the magnetized stellar wind is specified by the wind braking law (see Section 4.2 for details);

Over the coupling timescale , the amount of angular momentum is transferred from the radiative interior to the convective envelope (see MacGregor & Brenner 1991).
The TZM equations for and are therefore:

For :
(4) 
For :
(5)
In the equations above, we have introduced the moments of inertia, and , of the radiative zone and of the convective envelope, respectively. Note the appearance of terms involving the time derivatives of the moments of inertia (, etc.), as required by angular momentum conservation when significant structural changes take place.
We begin the integration of equations (4) and (5) when the star is still fully convective, when it is reasonable to assume that the star is rotating as a solid body. We thus set the initial conditions in terms of a single parameter, the initial period :
We extract the evolution of stellar structure quantities, , , etc., from stellar models of appropriate mass, calculated using the Yale Rotational stellar Evolution Code (YREC), in its nonrotational configuration (Demarque et al. 2008). The stellar models were calculated with initial composition and mixing length parameter calibrated on the Sun. Assuming the Grevesse & Sauval (1998) value of the solar metallicity, , our solar calibration gives , , , .
4.2 Wind braking laws
In our TZM framework, the wind braking law prescribes the dependence of on the stellar mass (either directly or indirectly, i.e., through a dependence on other stellar structure variables), and on the surface angular velocity . In this work, we focus on the following wind braking laws:

A hybrid braking law, which combines the classical dependence of Kawaler (1988) with the mass dependence suggested by Barnes & Kim (2010) and Barnes (2010) for the slowrotator sequence:
(7) where and are the moment of inertia of the whole star and of the Sun, respectively, and is an adjustable parameter with the same scaling as above.

The braking law proposed by Gallet & Bouvier (2013). Their starting point is the very general expression , where is the mass loss due to the stellar wind and is the Alfvèn radius, i.e., the radius of the (approximately) spherical region around the star where the wind is dominated by the magnetic field (Weber & Davis 1967). Gallet & Bouvier (2013) use an expression derived from the numerical simulation of Matt et al. (2012) to eliminate the Alfvèn radius from the braking law; their final result is (see Sect. 3.3 of Gallet & Bouvier 2013, for details):
(8) In the equation above, is the escape velocity at the stellar surface, is the surface magnetic field and is the stellar mass loss rate due to the wind. In the equation above, the constants , , and are fixed according to the numerical simulations of Matt et al. (2012); their values are: , , . The values of and are calculated using the BOREAS subroutine (Cranmer & Saar 2011) with the filling factor slightly modified in order to reproduce the average filling factor of the present Sun as in Gallet & Bouvier (2013).

The braking law proposed by Matt et al. (2015). They derived a physically motivated scaling for the dependence of the wind braking law on Rossby number and adopted an empirical scaling with stellar mass and radius. For our purposes, the Matt et al. (2015) braking law can be recast in the form:
(9) with
(10) and
(11) Following Matt et al. (2015), we adopt , which implies that this model also reproduces the empirical Skumanich (1972) law in the solidbody rotation limit. The saturation regime is equivalent to that introduced by Krishnamurthi et al. (1997); here we adopt as in Matt et al. (2015) (see also Sect. 5).
4.3 MCMC fitting of the TZM parameters
According to the TZM, the rotational evolution is completely determined by the stellar mass, which also selects the appropriate background stellar model, and by five parameters: the initial period , the disc lifetime , the coupling timescale , and the two parameters contained in the wind braking laws (6) or (7). We use a MCMC approach to constrain the values of these parameters.
The MCMC method is a powerful technique to obtain quantitative inferences on a set of model parameters from the observational data (see, e.g. Appendix A of Tegmark et al. 2004 for a concise introduction).
In the case at hand, the vector of model parameters for the TZM with, e.g, the braking law (6) is , while for each mass listed in Table 2 the vector of data points contains the periods from the nonparametric fits at the various cluster ages.
For the models with , we also consider the additional constraint based on the current solar rotation rate
Formally, Bayes’ Theorem relates the probability of the model, given the data, i.e., the posterior probability
(12) 
In the following, we ignore the normalization , which is a constant with respect to the models, and we use the notation for the likelihood. In essence, the MCMC method consists of sampling the posterior probability by means of a quasirandom walk in the space of the parameters , constructed as a chain of jumps guided by the values of the likelihood . Each step in the chain depends only on the previous one (Markov Chain), and is partly stochastic (Monte Carlo). The steps are performed according to the following twostage rule:

Propose stage: , with extracted from the jump function ;
The MetropolisHastings rule ensures that a trial step increasing the likelihood is always accepted (), while at the same time even steps which reduce the likelihood retain a (proportionally small) nonzero probability to be accepted. This procedure is completely specified by the choice of the functions and .
We use a multidimensional Gaussian form for the jump function:
(14) 
where are the jump sizes and the index (up to in our case) runs over the elements of the parameter vector . The jump function is critical to ensure that the steps in the chain are optimally correlated, providing a useful sampling of the posterior probability: the optimal regime is in between step acceptance (i.e., a purely random walk) and step rejection.
We define the likelihood in terms of the chisquare:
(15) 
The chisquare contains contributions from each cluster for which an estimate of the rotation period on the slowrotator sequence at the mass considered is available:
(16) 
where is the angular velocity of the convective envelope calculated from a run of the TZM with the set of parameters , while , are derived from Table 2 as follows:
In practical applications, running the MCMC procedure at the optimal acceptance regime ( of attempted steps accepted), is often hindered by the existence of correlations in the way the parameters influence the model. In our case, for example, we can anticipate the existence of degenerate pairs of values of and , i.e., a longer initial period and a shorter disc interaction, and vice versa, compensating each other to give a very similar evolution at sufficiently late times (say, a few hundred Myr onwards). Apart from computational efficiency, strong correlations can lead to oversampling of the degenerate regions of the parameter space, resulting in very biased, and practically useless, chains. We deal with this issue as recommended by Tegmark et al. (2004), by introducing a transformation from to an uncorrelated vector , and quantitatively checking the quality of the chain:

To transform to uncorrelated variables, we calculate the correlation matrix of the parameters
over a portion of the chain (typically, every steps), then diagonalize it to obtain the matrices and :
since is by definition real and symmetric, the eigenvalues are real and the matrix is orthogonal. The linear transformation
results in the very useful properties: , , i.e., the variables are (linearly) uncorrelated and normalized to one. Operationally, the vector is transformed into before substitution into the jump function (14); the trial step is performed in the variables , and then transformed back if the step is accepted. Since is normalized to one, the choice of the jump sizes is particularly simple: . The correlation matrix is recalculated, and the matrices and are updated, every steps.

To evaluate the quality of the chain, we calculate the autocorrelation of each parameter ,
where , etc., is the value of at the step of the chain. Typically, decreases with (its value at being unity by definition): we define the correlation length of the chain as the value for which drops below for the first time. The effective length of the chain is thus the number of steps divided by :
This is a measure of the number of independent points in the chain: if is sufficiently large, the average of over the chain, along with its standard deviation, can be considered representative of the bestfitting value of and its uncertainty, respectively. To ensure that we only use a portion of the chain sufficiently close to convergence, the first few hundred steps (burnin phase) are discarded.
Mass ()  Model  (Myr)  (Myr)  

0.85  KA  2.00  6.03  0.48  12250  85.5  11.1  8166 
KS  2.03  6.16  0.48  12250  86.8  10.7  8166  
KK  1.72  5.47  0.40  12250  79.7  10.6  8166  
KB  1.99  4.78  0.37  12250  84.6  10.8  8166  
MB  2.12  5.32  0.41  8166  83.1  11.5  8166  
GB  1.54  2.42  0.08  12250  85.0  9.1  8166  
0.90  KA  1.21  5.24  0.50  8166  51.3  9.4  8166 
KS  1.25  5.32  0.50  8166  52.6  9.4  8166  
KK  1.28  5.41  0.49  8166  53.5  9.1  8166  
KB  1.61  5.03  0.48  8166  59.6  9.3  8166  
MB  1.44  5.28  0.48  8166  52.3  9.9  8166  
GB  1.09  2.43  0.10  12250  67.4  7.1  12250  
0.95  KA  1.05  4.90  0.61  8166  34.9  9.4  8166 
KS  1.11  4.96  0.62  8166  35.5  9.1  8166  
KK  1.45  5.13  0.52  8166  38.3  7.9  8166  
KB  1.07  4.32  0.54  8166  35.1  9.1  8166  
MB  1.11  5.38  0.65  8166  34.9  9.4  8166  
GB  1.14  2.16  0.11  8166  49.0  6.4  8166  
1.00  KA  1.81  3.82  0.64  8166  19.1  8.5  8166 
KS  1.75  3.76  0.58  8166  19.0  8.1  8166  
KK  1.75  3.89  0.54  8166  20.0  7.4  8166  
KB  1.62  3.60  0.58  8166  18.7  8.1  8166  
MB  1.30  4.51  0.68  8166  17.3  8.3  8166  
GB  2.11  1.91  0.13  8166  34.1  5.7  8166  
1.05  KA  1.87  3.79  0.91  8166  15.6  7.9  8166 
KS  1.90  3.69  0.77  8166  15.6  7.5  8166  
KK  1.81  3.76  0.87  8166  15.4  7.8  8166  
KB  1.73  4.69  1.07  8166  15.3  7.4  8166  
MB  1.71  5.52  1.23  8166  15.0  8.3  8166  
GB  2.37  1.85  0.18  8166  23.0  5.2  8166  
1.10  KA  3.45  3.17  0.91  6125  12.0  6.6  6125 
KS  3.52  3.40  0.87  6125  12.8  6.4  6125  
KK  3.61  3.65  1.17  6125  12.9  7.3  6125  
KB  3.25  4.55  1.43  6125  12.5  6.9  6125  
MB  3.29  6.97  2.27  6125  12.4  7.9  6125  
GB  2.29  1.87  0.28  8166  14.1  4.3  8166 
5 Results and discussion
The empirical description of Sect. 3 implies that, once a star settles on the slowrotator sequence, it forgets most of its earlier rotational history. As suggested by Barnes (2003) or Brown (2014), fast rotators may have a different magnetic configuration than slow rotators and migrate to the slowrotator sequence after a change in their magnetic configuration. This scenario, which still lacks clear observational support (but see Garraffo et al. 2015), implies that stars may settle on the slowrotator sequence either after an initial rotational evolution like those described in Sect. 4.1 and 4.2, with parameters constant in time, or after an entirely different one, in which the parameters change with surface magnetic field configuration (e.g. Brown 2014). In the first case, convergence on the slowrotator sequence would be just a consequence of the dependence of the wind braking law, including the effects of the saturation regime, which is implicitly assumed in building our starting conditions. In the second case, both the overall scale factor and the exponent of in the wind braking law may change at some stage of the evolution, which would require a different approach than that adopted here. In both cases, however, once on the slowrotator sequence, the subsequent rotation period evolution is manifestly independent of the previous history, although the abundance of light elements may still depend on it (see, e.g. Bouvier 2008).
Our fitting, therefore, cannot constrain the stellar rotational history before the settling on the slowrotator sequence. In particular, it cannot constrain and even if the earlier rotational history follows one of the models listed in Sect. 4.1 and 4.2 with parameters constant in time. Nevertheless, such parameters are necessary in our modeling for advancing the rotational evolution from the birthline to the slowrotator sequence. We therefore treat and as nuisance parameters with reasonable priors and marginalize them out. The prior on is uniformly distributed over the observed period range in the ONC (Rebull et al. 2004). The prior on is assumed to be an exponential decay with an folding time of Myr (see, e.g. Fig. 11 in Hernández et al. 2008).
In our modeling, and are degenerate with respect to the slowrotator sequence evolution, in the sense that long and short can lead to the same period at the age of the Pleiades as short and long . Therefore, after checking compatibility by letting all parameters vary, we set Myr and marginalize out for each mass and model. It is worth stressing that this is just a convenient way of setting the initial conditions for the slowrotator sequence evolution and that no inference can be made on the actual distribution of and .
The duration of the phase of the early stellar evolution that is affected by the saturation regime (Eq. 6 or 9) depends on the assumptions made on the saturation threshold. As a consequence, the star may leave the saturation regime before or significantly after its settling on the slowrotator sequence. In fact, if we consider the relationship proposed by Krishnamurthi et al. (1997):
(17) 
with , or, equivalently, the saturation regime in the formulation of Matt et al. (2015) (see Sect. 4.2), stars of approximately solar mass will leave the saturation regime before settling on the slowrotator sequence (starting from the Pleiades age in our modeling), but stars of lower mass will stay in the saturation regime for a significant part of their earlier evolution on the slowrotator sequence. For this reasons, we adopt the strategy of testing different assumptions on , rather than trying to fit it as a free parameter. We therefore consider Eq. (6) with:

no saturation (model KA);

for all masses (model KS);

as in Eq. (17) (model KK).
Note that an implicit mass dependence (through ) enters both model KK, as a dependent saturation threshold, and Eq. (7) (hereafter model KB), as a factor in the wind angular momentum loss. On the other hand, Eq. (9) (hereafter model MB) contains a complex mass dependence that is both explicit and implicit, through the dependence on and , and the dependence of the saturation threshold, respectively.
The wind braking law (8) (model GB hereafter) contains, in principle, no adjustable parameters. However, in order to compare it with models derived from Eqs. (6), (7), and (9), where is treated as a free parameter, we assume variable and check whether the fitted shows some residual correlation with mass.
By fitting or , we both test how well the models describe the mass dependence of the rotational evolution, and correct the mass dependence using the observational data.
After setting the initial conditions, we are left with two free parameters, for the models KA, KS, KK, KB, and MB and for model GB.
We fit these parameters using the MCMC method described in Sect. 4.3 and the data of Table 2.
From this data set we exclude data for stellar mass , since, for our purposes, the sampling in age in this mass range is still insufficient
The results of the fitting of the MCMC parameters are reported in Table 3, Fig. 3 and 5. Fig. 3 shows that the vs. relation obtained for the KB models produces the lowest linear Pearson correlation coefficient, compatible, within the uncertainties, with no correlation. It is, indeed, remarkable that at Gyr the slowrotator sequence follows the shape predicted by Barnes (2010) for the asymptotic slowrotator sequence (Eq. 3) with an appropriate rescaling of the free parameter . In a solidbody rotation regime, this shape comparison would manifestly imply the mass scaling proportional to , implemented in model KB. Although stars at 0.55 Gyr have not yet reached this regime, our fit supports the validity of such mass scaling, at least in the mass range 0.85–1.1 . The confirmation of the general validity of such an empirical mass scaling awaits the acquisition of period data for at ages well beyond the current 0.6 Gyr limit.
Model MB also produces a significantly lower vs. correlation than the KA, KS, KK, and GB model. However, the rise at indicates that the mass dependence of the MB wind braking law is less accurate than model KB at larger mass.
Model KA, i.e. the original Kawaler (1988) model with no saturation, tends to underestimate the wind braking at lower mass with respect to higher mass (or vice versa). Our fit compensates for this behavior with a systematic variation of with mass. A saturation threshold constant in mass, as in model KS, helps to correct this behavior, but only by a negligible amount for . A massdependent saturation threshold, as in Eq. (17), is somewhat more effective, but a clear vs. correlation remains.
The fitted values of for model GB are very close to the one given in Gallet & Bouvier (2013), which derives from the original simulation of Matt et al. (2012), but show a clear correlation with . A strong correlation of the fitted with mass (using percentiles of the whole period sample) was also found by Gallet & Bouvier (2015), who interpreted it as due to a change in magnetic configuration with mass.
A comparison of the wind braking mass dependence of models KA, KB, and MB in the range is shown in Fig. 4. Comparing the latter with Fig. 3, we deduce that the slope of the mass dependence of both models KB and MB is compatible with the observations in the range . For , however, the slope of model MB is too steep and the observations favor a slope more similar to that of model KB. Note, furthermore, that the effective mass dependence of model MB depends on the saturation regime, although this has only a minor impact on the fit, as in model KK.
Age  

(Myr)  1.10  1.05  1.00  0.95  0.90  0.85 
100.  2.65  3.04  4.04  4.62  5.22  6.38 
125.  2.78  3.18  4.18  4.76  5.36  6.51 
160.  2.95  3.37  4.41  4.98  5.55  6.68 
200.  3.17  3.59  4.67  5.23  5.77  6.87 
250.  3.44  3.90  5.01  5.55  6.06  7.13 
320.  3.85  4.33  5.49  6.02  6.49  7.51 
400.  4.33  4.85  6.08  6.59  6.99  7.95 
500.  4.97  5.52  6.83  7.32  7.64  8.54 
630.  5.80  6.41  7.82  8.30  8.54  9.33 
800.  6.91  7.57  9.15  9.60  9.77  10.42 
1000.  8.17  8.90  10.67  11.15  11.24  11.75 
1250.  9.67  10.48  12.49  13.03  13.12  13.47 
1600.  11.62  12.52  14.86  15.54  15.66  15.90 
2000.  13.64  14.63  17.34  18.19  18.45  18.64 
2500.  15.93  17.03  20.14  21.23  21.68  21.92 
3200.  18.85  20.06  23.68  25.05  25.78  26.20 
4000.  21.98  23.21  27.32  28.96  29.98  30.66 
5000.  26.17  26.96  31.50  33.40  34.72  35.70 
The resulting is largely independent of the model adopted (Fig. 5). The largest difference is between model GB, which is the only one with a wind braking law not strictly proportional to , and all the others. Note that our derived from model GB is very close to that derived by Gallet & Bouvier (2015) for the 90th percentile of the whole period distribution. As discussed below, a different dependence results in a different angular momentum storage in the radiative interior, and therefore a different coupling timescale, but remains a steep function of mass in all cases. Adopting model KB with an average , we obtain a vs. relationship which is wellrepresented by the powerlaw scaling:
(18) 
In Fig. 5 this relationships is compared with the obtained from each of the models tested in the MCMC fit.
The evolution of the envelope and core angular velocities is plotted in Fig. 6. For clarity, only the results of models KB, MB, and GB are shown. Since the mass dependence has been corrected by the fit, all models with a wind braking law proportional to produce essentially the same and evolution.
In the – Gyr age range, the evolution deduced from model GB also overlaps with that deduced from all the other models. Differences with the other models become significant after Gyr, but remain small even at the solar age (see Fig. 6). Model GB model predicts, however, a larger angular momentum storage in the core at early ages with respect to the other models.
The comparison of Figs. 1 and 6 explains the behavior of the slowrotator sequence in the vs. diagram. For at 0.1 Gyr, is well below the
(19) 
relationship with Gyr, and the corresponding value from the non parametric fit (Sect. 3.2). The evolution proceeds by approaching such relationship until Gyr, then decreases more steeply than Eq. (19). In the timeframes of Fig. 1, this corresponds to a slowrotator sequence approaching the M37 sequence between and Gyr, then departing from it between and Gyr towards increasing . The behaviour for is qualitatively similar, but, in this case, remains quite close to Eq. (19) until Gyr, then decreases significantly faster than Eq. (19) afterwards. This explains why the slowrotator sequence for remains close to the M37 sequence in the vs. diagram until Gyr, and then departs from it again from 0.55 to 2.5 Gyr toward higher values.
The comparison between and shows that a quasisolid body rotation regime, and therefore a regime in which the wind braking dominates the rotational evolution, is achieved well after 1 Gyr. Furthermore, the comparison between and Eq. (19) with Gyr shows that a quasiSkumanich evolution is achieved well after 1 Gyr. This implies that in order to verify a departure from the proportionality of the wind braking law, more data at ages well above 1 Gyr are needed.
Period isochrones in the periodcolor diagrams are calculated with model KB adopting and correcting for likely biases according to the powerlaw fitting shown in Fig. 5. These isochrones are essentially equivalent to those computed using any of the wind braking models tested after correcting the mass dependence through the fitted and correcting in the same way. The comparison with observations, shown in Fig. 7, outlines the observational origin of the small discrepancies found at the ages of the Pleiades and M35 for , where photometry and reddening uncertainties, particularly crucial for a cluster like M35, affect more significantly the results. Fig. 7 also shows that our model reproduces satisfactorily the period distribution in NGC6811 as well, even if it had been originally excluded from the MCMC fit. Even the reconstructed period at , where the distribution apparently bends towards lower periods, lies on the upper envelope of the observed distribution and it is therefore compatible, in a statistical sense, with the observed distribution. Overall, however, the evolution of the slowrotator sequence is reproduced with unprecedented accuracy.
Period isochrones recomputed for a logarithmically spaced age grid are reported in Table 4 and Fig. 8. Interpolation of Table 4 provides a tool for inferring stellar ages of solarlike mainsequence stars from their mass and rotational period, effectively representing an alternative gyrochronology relationship that takes the physics of the twozone model for the stellar angular momentum evolution into account.
6 Conclusions
Using a twozone model MCMC fitting to the rotational periods vs color data from selected open clusters, we have compared new proposals for the wind braking law with the classical Kawaler (1988) one modified by Chaboyer et al. (1995). In our tests, we coupled the Kawaler (1988) braking law with different hypotheses regarding the saturation regime, including that of Krishnamurthi et al. (1997). The wind braking laws tested also included that of Gallet & Bouvier (2013), one derived from the work of Barnes & Kim (2010) and Barnes (2010), and the new proposal by Matt et al. (2015).
The MCMC fitting has been performed on the slowrotator sequence between and Gyr. We proposed a new quantitative criterion for identifying the slowrotator sequence based on the symmetry of the data distribution around the density peak in the periodcolor diagram. Following this criterion, it is possible to assign a width of the period distribution around the density peak using the standard deviation. The value of the peak density as a function of colors (mass) is estimated using a nonparametric fit and, together with the standard deviation, constitute the data used in the MCMC fitting and in the definition of the likelihood function.
The main idea behind this work was to see how accurately a twozone model with parameters constant in time can describe the rotational evolution on the slowrotator sequence. In fact, it is evident that the slowrotator sequence evolves in a regular and predictable way and it is then natural to explore the extent to which current models can reproduce this sort evolution.
We find that the difficulties in reconciling the observed behavior of the slowrotator sequence evolution between 0.1 and 0.6 Gyr with the original Barnes (2003) idea of a factorization of time and mass dependence (Meibom et al. 2009, 2011b) can be simply attributed to a transfer of angular momentum from the radiative core to the convective envelope in the early slowrotator sequence evolution. The twozone model MCMC fitting leads to a robust (largely independent of the wind braking model adopted) estimate of the coreenvelope coupling timescale as a function of mass in the range . In this mass range we find that scales as on the slowrotator sequence.
The importance of the core angular momentum storage and transfer to the convective envelope, in which enhanced viscosity that is due to (magneto) hydrodynamical instabilities is expected to play an important role, has been outlined by many authors (see, e.g. references in Sect. 4.1). In this work we derive, for the first time, the empirical mass dependence of the coreenvelope coupling timescale in the slowrotator regime. Our empirical mass dependence of represents a useful term of comparison for the theoretical description of angular momentum transport in stellar interiors. By focusing on the slowrotator sequence, in which the TZM with constant parameters, as described in Sect. 4.1, is expected to capture the essential physics of the angular momentum evolution, we avoid the complications that arise in dealing with the fastrotator regime, in which the TZM with constant parameters may break down (e.g. Brown 2014). The validity of the TZM for the slowrotator regime is confirmed by the successful fit to the cluster data currently available. The wind braking law, for which the empirical dependence on stellar mass and angular velocity still gives better results than purely theoretical ones (see e.g. Matt et al. 2015), remains one of the critical aspects of this modeling.
Our results outline and explain the invalidity of a factorization like the one expressed by Eq. (1), on which most of the proposed gyrochronology relationships are based, at least for the early (massdependent) evolution of the slowrotator sequence. On the other hand, our results mean that the search for gyrochronology relationships like the one recently proposed by Kovacs (2015), which is not based on known physics, are not justified. Kovacs (2015) proposal was motivated in fact by the disagreement between gyrochronology relationships based on Eq. (1) and stellar evolution isochrone fitting, but the ages assumed for the clusters included in our analysis, which do not differ significantly from the isochrone ages assumed in Kovacs (2015), are fully compatible with our modeling of the rotational evolution of the slowrotator sequence.
With a suitable choice of the free parameters, we tested the mass dependence of the models under scrutiny. We found that the mass scaling of the wind braking law implied by the models of Barnes (2010) and Barnes & Kim (2010) is the most consistent with the data, supporting the proposal that the convective turnover timescale is a key parameter in such a mass dependence. A definitive confirmation of the general validity of such mass dependence, however, would require period data of older cluster and lower stellar mass than is available to date. The mass scaling of the recent model by Matt et al. (2015) is found consistent with the data for , but it decreases too steeply with mass for . All the other models considered show a clear residual correlation with mass, indicating an incorrect mass scaling.
Small deviations from the dependence of the wind braking law, such as those predicted by the Gallet & Bouvier (2013) model, are also compatible with observations in the 0.1–2.5 Gyr range. In fact, quasisolid body rotation, the regime in which the wind braking dominates the angular momentum evolution, is achieved after 1–2 Gyr, depending on mass, so that small deviations from the proportionality can be compensated for by a different angular momentum storage in the stellar core to produce, essentially, the same evolution. After correcting the mass scaling, differences between the Gallet & Bouvier (2013) model and models based on the proportionality becomes not negligible after 2.5 Gyr, although they remain small even at the age of the Sun. A corollary of this is that verifying small deviations from the proportionality would require data for clusters much older than 2.5 Gyr and that, in the 2–5 Gyr age range, deviations from the Skumanich law for are expected to be small.
From our MCMC model fitting to the data, we derive period isochrones and tracks from 0.1 to 5.0 Gyr valid for stars with . These are essentially independent of any assumption regarding the mass scaling of the wind braking law and compatible with small deviations from the proportionality. By interpolating between the grid point, these tracks can be used to infer stellar ages from their mass and rotational period, providing alternative gyrochronology relationships that take the physics of the twozone model into account.
Acknowledgements.
We are grateful to an anonymous referee, whose comments have contributed to improving the paper. ACL thanks the Cosmic Magnetic Fields Branch of the LeibnizInstitut für Astrophysik Potsdam (AIP) for their kind hospitality. ACL acknowledges the support received from the European Science Foundation (ESF) for the activity entitled ’Gaia Research for European Astronomy Training’. FS acknowledges support from the Leibniz Institute for Astrophysics Potsdam (AIP) through the Karl Schwarzschild Postdoctoral Fellowship.Footnotes
 offprints: A. C. Lanzafame
 We adopt the values Gyr, s.
 The standard notation denotes the conditional probability of event , given that event is observed.
 Note that the strong mass dependence discussed in the following implies that, in order to constrain the TZM parameters at , data at ages Gyr would be required, which is not available yet.
References
 Allain, S. 1998, A&A, 333, 629
 Barnes, S. A. 2003, ApJ, 586, 464
 Barnes, S. A. 2007, ApJ, 669, 1167
 Barnes, S. A. 2010, ApJ, 722, 222
 Barnes, S. A. & Kim, Y.C. 2010, ApJ, 721, 675
 Barrado y Navascués, D., Deliyannis, C. P., & Stauffer, J. R. 2001, ApJ, 549, 452
 Bouvier, J. 2008, A&A, 489, L53
 Brown, T. M. 2014, ApJ, 789, 101
 Chaboyer, B., Demarque, P., & Pinsonneault, M. H. 1995, ApJ, 441, 876
 Cleveland, W. S., Grosse, E., & Shyu, W. M. 1992, Local regression models, Statistical Models in S (Wadsworth & Brooks/Cole)
 Collier Cameron, A., Davidson, V. A., Hebb, L., et al. 2009, MNRAS, 400, 451
 Cranmer, S. R. & Saar, S. H. 2011, ApJ, 741, 54
 Delorme, P., Collier Cameron, A., Hebb, L., et al. 2011, MNRAS, 413, 2218
 Demarque, P., Guenther, D. B., Li, L. H., Mazumdar, A., & Straka, C. W. 2008, Ap&SS, 316, 31
 Epstein, C. R. & Pinsonneault, M. H. 2014, ApJ, 780, 159
 Gallet, F. & Bouvier, J. 2013, A&A, 556, A36
 Gallet, F. & Bouvier, J. 2015, A&A, 577, A98
 Garraffo, C., Drake, J. J., & Cohen, O. 2015, ApJ, 807, L6
 Grevesse, N. & Sauval, A. J. 1998, Space Sci. Rev., 85, 161
 Hartman, J. D., Bakos, G. Á., Kovács, G., & Noyes, R. W. 2010, MNRAS, 408, 475
 Hartman, J. D., Gaudi, B. S., Holman, M. J., et al. 2008, ApJ, 675, 1233
 Hartman, J. D., Gaudi, B. S., Pinsonneault, M. H., et al. 2009, ApJ, 691, 342
 Hastings, W. K. 1970, Biometrika, 57, 97
 Herbst, W. & Mundt, R. 2005, ApJ, 633, 967
 Hernández, J., Hartmann, L., Calvet, N., et al. 2008, ApJ, 686, 1195
 Hodgkin, S. T., Irwin, J. M., Aigrain, S., et al. 2006, Astronomische Nachrichten, 327, 9
 Janes, K., Barnes, S. A., Meibom, S., & Hoq, S. 2013, AJ, 145, 7
 Jeffries, Jr., M. W., Sandquist, E. L., Mathieu, R. D., et al. 2013, AJ, 146, 58
 Kalirai, J. S., Fahlman, G. G., Richer, H. B., & Ventura, P. 2003, AJ, 126, 1402
 Kalirai, J. S., Ventura, P., Richer, H. B., et al. 2001, AJ, 122, 3239
 Kawaler, S. D. 1988, ApJ, 333, 236
 Keppens, R., MacGregor, K. B., & Charbonneau, P. 1995, A&A, 294, 469
 Koenigl, A. 1991, ApJ, 370, L39
 Kovacs, G. 2015, ArXiv eprints
 Krishnamurthi, A., Pinsonneault, M. H., Barnes, S., & Sofia, S. 1997, ApJ, 480, 303
 Lejeune, T., Cuisinier, F., & Buser, R. 1997, A&AS, 125, 229
 Lejeune, T., Cuisinier, F., & Buser, R. 1998, A&AS, 130, 65
 MacGregor, K. B. & Brenner, M. 1991, ApJ, 376, 204
 Mamajek, E. E. & Hillenbrand, L. A. 2008, ApJ, 687, 1264
 Matt, S. P., Brun, A. S., Baraffe, I., Bouvier, J., & Chabrier, G. 2015, ApJ, 799, L23
 Matt, S. P., MacGregor, K. B., Pinsonneault, M. H., & Greene, T. P. 2012, ApJ, 754, L26
 Meibom, S., Barnes, S. A., Latham, D. W., et al. 2011a, ApJ, 733, L9
 Meibom, S., Barnes, S. A., Platais, I., et al. 2015, Nature, 517, 589
 Meibom, S., Mathieu, R. D., & Stassun, K. G. 2009, ApJ, 695, 679
 Meibom, S., Mathieu, R. D., Stassun, K. G., Liebesny, P., & Saar, S. H. 2011b, ApJ, 733, 115
 Messina, S. 2007, Mem. Soc. Astron. Italiana, 78, 628
 Messina, S., Distefano, E., Parihar, P., et al. 2008, A&A, 483, 253
 Messina, S., Parihar, P., Koo, J.R., et al. 2010, A&A, 513, A29
 Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087
 Moore, D. 1986, Tests of the chisquared type. GoodnessofFit Techniques. (Marcel Dekker, New York.)
 Rebull, L. M., Wolff, S. C., & Strom, S. E. 2004, AJ, 127, 1029
 Skumanich, A. 1972, ApJ, 171, 565
 Spada, F., Lanzafame, A. C., Lanza, A. F., Messina, S., & Collier Cameron, A. 2011, MNRAS, 416, 447
 Stauffer, J. R., Hartmann, L. W., Fazio, G. G., et al. 2007, ApJS, 172, 663
 Stauffer, J. R., Schultz, G., & Kirkpatrick, J. D. 1998, ApJ, 499, L199
 Tegmark, M., Strauss, M. A., Blanton, M. R., et al. 2004, Phys. Rev. D, 69, 103501
 Thode Jr., H. 2002, Testing for Normality. (Marcel Dekker, New York.)
 von Braun, K., Lee, B. L., Seager, S., et al. 2005, PASP, 117, 141
 von Hippel, T., Steinhauer, A., Sarajedini, A., & Deliyannis, C. P. 2002, AJ, 124, 1555
 Weber, E. J. & Davis, Jr., L. 1967, ApJ, 148, 217