Toward Precise Stellar Ages: Combining Isochrone Fitting with Empirical Gyrochronology

We present a new age-dating technique that combines gyrochronology with isochrone fitting to infer ages for FGKM main-sequence and subgiant field stars. Gyrochronology and isochrone fitting are each capable of providing relatively precise ages for field stars in certain areas of the Hertzsprung–Russell diagram (HRD): gyrochronology works optimally for cool main-sequence stars, and isochrone fitting can provide precise ages for stars near the main-sequence turnoff. Combined, these two age-dating techniques can provide precise and accurate ages for a broader range of stellar masses and evolutionary stages than either method used in isolation. We demonstrate that the position of a star on the HRD or color–magnitude diagram can be combined with its rotation period to infer a precise age via both isochrone fitting and gyrochronology simultaneously. We show that incorporating rotation periods with 5% uncertainties into stellar evolution models improves age precision for FGK stars on the main sequence and can, on average, provide age estimates up to three times more precise than isochrone fitting alone. In addition, we provide a new gyrochronology relation, calibrated to the Praesepe cluster and the Sun, that includes a variance model to capture the rotational behavior of stars whose rotation periods do not lengthen with the square root of time and parts of the HRD where gyrochronology has not been calibrated. This publication is accompanied by an open-source Python package (stardate) for inferring the ages of main-sequence and subgiant FGKM stars from rotation periods, spectroscopic parameters, and/or apparent magnitudes and parallaxes.


Introduction
Age is the most difficult stellar property to measure, and the difficulty of age dating is particularly acute for low-mass (GKM) stars on the main sequence (MS). Using conventional dating methods, uncertainties on the ages of these stars can be as large as the age of the universe. GKM dwarfs are difficult to age-date because most of their physical and observable properties do not change rapidly. This is represented in the spacing of isochrones on a Hertzsprung-Russell diagram (HRD) or color-magnitude diagram (CMD). On the MS, isochrones are tightly spaced, and even with very precise measurements of effective temperature and luminosity, the position of an MS star on the HRD may be consistent with the range of isochrones spanning several billion years. At the mainsequence turnoff, however, isochrones are spread farther apart, so that sufficiently precisely measured temperatures and luminosities can yield ages that are extremely precise (with minimum statistical uncertainties on the order of 5%-10%; e.g., Pont & Eyer 2004). The classical method for measuring stellar ages is isochrone placement, or isochrone fitting, where surface gravity changes resulting from fusion in the core (usually observed via luminosity, L, and effective temperature, T eff , or absolute magnitude and color) are compared with a set of models that trace stellar evolution across the HRD or CMD. CMD/HRD position has been thoroughly mapped with physical models and can be used to calculate relatively accurate (but not necessarily precise) ages, barring some systematic variations between different models, (e.g., Yi et al. 2001;Dotter et al. 2008;Dotter 2016). On the MS itself, there is little differentiation between stars of different ages in the L and T eff plane, so ages tend to be very imprecise. The method of inferring a star's age from its rotation period, called "gyrochronology," is much better suited for measuring ages on the MS because MS stars spin down relatively rapidly.
Magnetic braking in MS stars was first observed by Skumanich (1972), who, studying young clusters and the Sun, found that the rotation periods of solar-type stars decay with the square root of time. It has since been established that the rotation period of an MS star depends, to first order, only on its age and effective temperature or color (e.g., Barnes 2003). The convenient characteristic of stars that allows their ages to be inferred from their current rotation periods, and independently of their primordial ones, comes from the steep dependence of the spin-down rate on rotation period (Kawaler 1989). Stars spinning with high angular velocity will experience a much greater angular momentum loss rate than slowly spinning stars, and for this reason, no matter the initial rotation period, solar-type stars will have the same rotation period after around the age of the Hyades, 500-700 million years (Irwin & Bouvier 2009;Gallet & Bouvier 2015). After this time, the age of a star can be inferred, to first order, from its dust-corrected color (e.g., -B V or G BP -G RP ) and current rotation period alone (see Epstein & Pinsonneault 2014 for an analysis of how initial conditions affect gyrochronal ages).
The relation between age, rotation period, and mass has been studied in detail, and several different models have been developed to capture the rotational evolution of Sun-like stars. Some of these models are theoretical and based on physical processes, modeling angular momentum loss as a function of stellar properties as well as the properties of the magnetic field and stellar wind (e.g., Kawaler 1988Kawaler , 1989Pinsonneault et al. 1989;van Saders & Pinsonneault 2013;Matt et al. 2015;van Saders et al. 2016). Other models are empirical and capture the behavior of stars from a purely observational standpoint, using simple functional forms that can reproduce the data (e.g., Barnes 2003Barnes , 2007Mamajek & Hillenbrand 2008;Angus et al. 2015). Both types of model, theoretical and empirical, must be calibrated using observations. Old calibrators are especially important because new evidence suggests that rotational evolution goes through a transition at old age or, more specifically, at a large Rossby number, Ro (the ratio of rotation period to the convective overturn timescale). For example, stars shown to be old from Kepler asteroseismic data rotate more rapidly than expected given their age (Angus et al. 2015;van Saders et al. 2016). A new physically motivated gyrochronology model, capable of reproducing these data, was recently introduced (van Saders et al. 2016). It relaxes magnetic braking at a critical Rossby number of around 2, approximately the solar value. This model predicts that, after stellar rotation periods lengthen enough to move stars across this Ro threshold, stars conserve angular momentum and maintain a nearly constant rotation period from then until they evolve off the MS.
The gyrochronology models that capture post-Ro-threshold rotational evolution ) are the current state of the art in rotation dating. These models can be computed over a grid of stellar parameters and interpolated over to predict the age of a star. The process of measuring the age of a field star with these models is similar to inferring an age using any set of isochrones, with the difference being that rotation period is an additional observable dimension. Ages calculated using these models are therefore likely to be much more precise than using rotation-free isochrones because the rotation period provides an additional anchor point for the age of a star. We present here a complementary method that combines isochrones with an empirical gyrochronology model. The methodology is related to the van  model in that both use a combination of rotation periods and other observable properties that track stellar evolution on the HRD in concert. The main difference is that the gyrochronology model used here is an entirely empirically calibrated one, as opposed to a physically derived one.
One major advantage of using a physically motivated gyrochronology model is the ability to rely on physics to interpolate or extrapolate over parts of parameter space with sparse data coverage. However, rotational spin-down is a complex process that is not yet fully understood, and currently, no physical model can accurately reproduce all of the data available. For this reason, even physically motivated gyrochronology models cannot always be used to reliably extrapolate into unexplored parameter space. Physical models, when calibrated to data, can provide insight into the physics of stars; however, if accurate and precise prediction of stellar properties is desired, empirical models can have advantages over physical ones. For example, the data may reveal complex trends that cannot be reproduced with our current understanding of the physical processes involved, but may be captured by more flexible data-driven models. In addition, it is relatively straightforward to build an element of stochasticity into empirical models, i.e., to allow for and incorporate outliers or noisy trends. This is particularly important for stellar spindown because rotation periods can be affected by additional confounding variables which are not always observed (having a binary or planetary companion, for example). A further advantage of empirical models is that inference is more tractable: it can be extremely fast to fit them to data.
In this work, we calibrated a new empirical gyrochronology relation, fit to the Praesepe open cluster and the Sun, in Gaia G BP -G RP color. Gaia G, G BP , and G RP apparent magnitudes are now the most abundant photometric measurements, available for more than a billion stars (Gaia Collaboration et al. 2018). In fact, an important point of context for this work is the new availability of data relevant to gyrochronology and stellar ages. Gaia now provides broadband photometry and parallaxes for over a billion stars, and Kepler, K2, and TESS are providing rotation periods for hundreds of thousands of stars. Gyrochronology is becoming one of the most readily available age-dating methods, so continuing to improve gyrochronology relations and methods is important. Like any other gyrochronology relation, this new Praesepe-based model does not perfectly reproduce all of the observed data, and some simple modifications could make significant improvements, for example, by including a mixture model to account for outliers and binaries, and by removing the period-age, period-color separability to account for the different period-color shapes seen in clusters of different ages. We leave these improvements for a future project and, for now, test this new Praesepe-based model, which is built into an open-source Python package for gyrochronology called stardate (Angus et al. 2019). Stardate provides the framework for simultaneous gyrochronology and isochrone fitting, and because stardate is modular, it would be straightforward to update the gyrochronology relation in the future.
This paper is laid out as follows. In Section 2, we describe our new age-dating model and its implementation. In Section 3, we test this model on simulated stars and cluster stars, and in Section 4, we discuss the implications of these tests and future pathways for development. Throughout this paper, we use the term "observables" to refer to the following observed properties of a star: T eff , log g, observed bulk metallicity ([Fe/H]), parallax (p), photometric colors in different passbands phase (EEP), model metallicity (M/H), distance (D), and V-band extinction (A V ). These are the properties that generate the observables.

A New Empirical Gyrochronology Relation
We fit a broken power law to the 650 Myr Praesepe cluster in order to calibrate a gyrochronology relation that takes advantage of new data available from the K2 and Gaia surveys. This relation captures the detailed shape of Praesepe's rotation period-color relation and is calibrated to Gaia G BP -G RP color. Praesepe is the ideal calibration cluster because it has the largest number of members with precisely measured rotation periods, over a large range of colors (spanning spectral types A0 through M6; Rebull et al. 2017), of any open cluster. It is relatively tightly clustered on the sky, and many of its members were targeted in a single K2 campaign, from which rotation periods have been measured via light-curve frequency analysis (Douglas et al. 2017;Rebull et al. 2017). We compiled rotation periods, Gaia photometry, and Gaia parallaxes for members of the 650 Myr Praesepe cluster (Fossati et al. 2008) by identifying Praesepe members with measured rotation periods from Douglas et al. (2017) in the K2-Gaia crossmatch catalog provided athttps://gaia-kepler.fun/. This catalog crossmatched the EPIC catalog (Huber et al. 2016) with the Gaia DR2 catalog (Gaia Collaboration et al. 2018) using a 1″ search radius. The result was a sample of 757 stars with rotation periods, parallaxes, and Gaia G-, G BP -, and G RP -band photometry, shown in Figure 1. 12 Although this new model does not perfectly describe the rotation period of every star, it provides a better representation of rotational evolution than previous models such as the Angus et al. (2015) model (shown as the blue solid line in Figure 1).
In order to fit a relation to Praesepe, we removed rotational outliers bluer than G BP -G RP = 2.7 via sigma-clipping and fit a fifth-order polynomial to the remaining FGK and early M stars. We found that a fifth-order polynomial provided a substantially better fit than lower-order polynomials, which were not able to capture the sharp "elbow" in the rotation period-color relation. Additional orders provided either a worse fit, tending toward extreme values at the boundaries, or diminishing returns in goodness of fit. We also fit a straight line to the late M dwarfs (G BP -G RP >2.7) to capture the mass-dependent initial rotation periods of low-mass stars (Somers et al. 2017). We fit a separable straight line function to the period-age relation using the ages of Praesepe and the Sun. 13 This new Praesepecalibrated gyrochronology relation is for stars with G BP -G RP >2.7, where P rot is the rotation period in days, and t is age in years. Best-fit coefficient values are shown in Table 1. We calibrated this new relation using Praesepe and the Sun alone, without including other clusters, because different open clusters have slightly different period-color relationships Agüeros 2018;Curtis & Agüeros 2018;Curtis et al. 2019), and given that we used an age-color separable relation, adding in extra clusters was unlikely to improve the calibration. We did not include asteroseismic stars because most are slightly evolved and would require a gyrochronology relation that depends on logg. This new relation, fit to Praesepe and the Sun, does not perfectly predict the rotation periods of stars at all colors and ages, but it provides several improvements over previous empirical gyrochronology relations. First, it uses new K2 rotation period measurements to model the period-color relation of Praesepe in detail; second, it includes a model for the rotational behavior of M dwarfs; and third, it is calibrated to GaiaG BP -G RP color: a directly observable quantity and the most widely available photometric color index.
Equation (1), describing the rotational evolution of FGK and early M stars, is most closely analogous to previously calibrated empirical gyrochronology relations (e.g., Barnes 2003Barnes , 2007Mamajek & Hillenbrand 2008;Barnes 2010;Angus et al. 2015). It describes stars with Sun-like magnetic  (3)). The rotation periods of stars bluer than around 0.56 dex and redder than around 2.7 dex in G BP -G RP are modeled as a broad log-normal distribution with a standard deviation of 0.5 dex, added to observational uncertainties. This figure was generated in a Jupyter notebook available athttps://github.com/RuthAngus/ stardate/blob/master/paper/code/Fitting_Praesepe.ipynb.  (1) and (2) Coefficient Value dynamos that follow a "Skumanich-like" magnetic braking law, i.e., their rotation period increases with the square root of their age. It does not describe stars hotter than around 6250 K, which have a thin convective layer and a weak magnetic dynamo, nor does it describe fully convective stars, which take a long time to converge onto the Skumanich (1972) braking law (Krishnamurthi et al. 1997). In addition, stars with Rossby numbers larger than around 2 do not show Skumanich-like magnetic braking (van Saders et al. , 2019, so Equation (1) does not describe these stars. It also does not describe the rotation periods of subgiants or giants, whose rotation periods are influenced by their expanding radii, changing winds, and core-to-surface differential rotation (e.g., van Saders & Pinsonneault 2013;Tayar & Pinsonneault 2018). Finally, this relation does not describe the rotation periods of dynamically or magnetically interacting binaries, which often rotate more rapidly than isolated stars at the same age and color (Douglas et al. 2016). In order to include non-Skumanich type stars in our combined gyrochronal and isochronal model, we designed a composite gyrochronology relation, which describes a mean and variance model for the rotation period distributions across the HRD. Rotation periods are modeled differently for stars of different photometric color, EEP, Rossby number, age, and metallicity. The rotation period model for each of these groups is described below.
1. The rotation periods of late F, GK, and early M dwarfs (0.56<G BP -G RP <2.7), with Rossby numbers less than 2, are modeled as a log-normal distribution, with mean given by Equation (1) and variance given by the squared inverse of their observational uncertainties. 2. The rotation periods of fully convective stars (G BP -G RP >2.7) are modeled with a log-normal distribution with mean given by Equation (2), and an extra standard deviation of 0.5 dex, added to observational uncertainties. This distribution reflects the observed rotation periods of late M dwarfs in the Praesepe cluster, which span a broad range at every color. Late M dwarfs with masses 0.3 M e , temperatures 3500 K, and G BP -G RP 2.7 exhibit weak magnetic braking until at least after the age of Praesepe (∼650 million years). 3. The rotation periods of F-type and hotter stars, with G BP -G RP <0.56 dex, are modeled as a log-normal distribution with mean log 10 P rot =0.56 and an additional standard deviation, added to observational uncertainties, of 0.5 dex. Stars more massive than around 1.25 M e , with temperature6250 K and G BP -G RP  0.56, do not spin down appreciably over their mainsequence lifetimes because they do not have the deep convective envelope needed to generate strong magnetic fields. This model for the mean and variance of hot star rotation periods is based on stars hotter than 6250 K in the McQuillan et al. (2014) sample, which have a mean log 10 P rot of 0.56 dex and a standard deviation of around 0.5 dex. Both hot and cool stars retain rotation periods that are similar to their primordial distribution (see e.g., Matt et al. 2012;Somers et al. 2017). 4. The rotation periods of stars with large Rossby numbers are modeled as follows. The age at which a star's Rossby number would exceed 2 is calculated by inverting Equation (1). If a star's age is greater than this, its mean rotation period is given by P max =2/τ, where τ is the convective turnover timescale, calculated via stellar mass using Equation (11) from Wright et al. (2011). 5. The rotation periods of subgiants are described with a log-normal distribution with mean given by Equations (1) or (2), or 0.56, depending on its color, and an additional standard deviation of 5 dex. This is not an accurate model of subgiant rotation periods (see, e.g., van Saders & Pinsonneault 2013), but the highly inflated variance makes it a weakly constraining one. Isochrone fitting provides precise ages for subgiants, so by inflating the variance of the gyrochronology relation at large EEP, we allow a star's position on the HRD/CMD to dominate the age information over its rotation period. This is useful because we have not yet built an accurate rotation-age relation for subgiants into our model. 6. We model stars younger than around 250 Myr with a lognormal distribution with mean function given by Equations (1) or (2) (1) or (2), or 0.56, depending on its color, and a standard deviation of 0.5 dex. Gyrochronology is not calibrated at these extreme metallicities, due to a lack of suitable metalpoor and rich calibration stars. Rather than assume that the same gyrochronology model can be naïvely applied to these stars, we take a more conservative approach and model them with a broad Gaussian distribution.
Inflating the variance of the rotation period distribution for stars with non-Skumanich magnetic braking behavior has two purposes: (1) in the cases of hot stars and fully convective stars, it allows the broad distributions of rotation periods observed in clusters and the field to be matched, and (2) it downweights the age information provided by rotation periods in regions of the HRD/CMD where rotation periods are not information-rich or the gyrochronology model is inaccurate or uncalibrated. If the observational uncertainties on rotation periods are 5% on average, which corresponds to 0.05/ln 10∼0.02 dex, adding 0.5 amounts to a 25× increase in standard deviation, or a ∼600× increase in variance. So, practically speaking, when the standard deviation is inflated to 0.5 dex or more, ages are almost entirely inferred via isochrone fitting.
We used the following composite gyrochronology model to infer ages from rotation periods: where σ is the relative period uncertainty, divided by ln 10, on individual rotation period measurements, and σ P is an additional scatter that is a function of EEP, age, metallicity, and color. It takes a maximum value of 0.5 for hot stars and fully convective stars, 5 for subgiants and giants, and a minimum value of zero for late F, GK, and early M dwarfs. The variance model is shown in Figure 2. Sigmoid functions were used to provide smooth transitions between regions of low and high variance. Sharp changes in variance would produce sharp changes in likelihood, which would cause the posterior distributions over stellar parameters to be more difficult to sample. The sigmoid functions shown in Figure 2 reach half their maximum values at G BP -G RP =0.56 and 2.5 for hot and cool stars, respectively; EEP=454 for subgiants; age = 250 Myr for young stars; [M/H]=−0.2 for metal-poor stars; and 0.2 for metal-rich stars. The logistic growth rate, or steepness, of the sigmoid functions is 100 dex −1 for both color transitions, 0.2 EEP −1 for the EEP transition, 20 dex −1 for the age transition, and 5 dex −1 for the both metallicity transitions. The additional standard deviation is additive, so if a star is, e.g., hot, evolved, and metal poor, the additional standard deviation of its rotation period rises to 6.

Simultaneously Fitting Gyrochronology and Isochrones
The previous part of this section describes the model for the mean and variance of rotation periods as a function of their ages (and colors, EEPs, masses, and metallicities). In what follows, we describe how this model was combined with a stellar evolution model to infer stellar ages (and other parameters) via gyrochronology and isochrone fitting simultaneously. Our goal was to infer the age of a star from its observable properties by estimating the posterior probability density function (PDF) over age, where t is age, m x is a vector of apparent magnitudes in various bandpasses, P rot is the rotation period, andp is parallax. Spectroscopic properties (T eff , log g, and [Fe/H]) and/or asteroseismic parameters (Δν and ν max ) may also be available for a star, in which case they would appear to the right of the "|" in the above equation because they are observables. In order to calculate a posterior PDF over age, other stellar parameters must be marginalized over. These parameters are distance (D),  where the likelihood of the data given the model is , , , , M H , 7 x V rot and the prior PDF over parameters is The priors we used are described in the Appendix. We assumed that the process of magnetic braking is independent of hydrogen burning in the core, outside of the dependencies that are captured in the model. This assumption allowed us to multiply two separate likelihood functions Figure 2. The additional rotation period scatter, σ P , added to the observational period uncertainties in the model (see Equation (3)). The standard deviation was increased for early F and hotter stars (G BP -G RP <0.56), late M dwarfs (G BP -G RP >2.7), and evolved stars (EEP420) in order to downweight the age information supplied by rotation periods and reproduce the observed rotation period distributions. We also increased the variance for stars younger than around 250 Myr, because the rotation periods of these stars typically have not yet converged onto a tight gyrochronology sequence, and for very high-and low-metallicity stars (−0.2[Fe/H]0.2), because the gyrochronology relations have not been calibrated at these extreme values. Downweighting the gyrochronal likelihood by the inverse variance (1/σ 2 ) allowed the ages of these stars to be mostly inferred via isochrone fitting. Sigmoid functions were used to provide smooth transitions between regions of low and high variance.
14 The inferred metallicity, [M/H], is a model parameter that is different from the observed metallicity, [Fe/H], which would appear on the right side of the "|" in Equation (4). together: one computed using an isochronal model and one computed using a gyrochronal model. We assumed that the probability of observing the measured observables given the model parameters was a Gaussian and that the observables were identically and independently distributed. The isochronal likelihood function was where O I is the vector of n observables:p, m x plus spectroscopic and/or asteroseismic observables if available, and Σ is the covariance matrix of that set of observables. I is the vector of model observables that correspond to a set of parameters: t, E, [M/H], D, and A V , calculated using an isochrone model. We assumed there is no covariance between these observables, and so this covariance matrix consists of individual parameter variances along the diagonal with zeros everywhere else. The gyrochronal likelihood function was where P O is a 1D vector of observed logarithmic rotation periods and P P is the vector of corresponding logarithmic rotation periods, predicted by the model. Σ P comprises individual rotation period measurement uncertainties, plus an additional variance that is a function of EEP and G BP -G RP color, added in quadrature. This variance accounts for the stochastic nature of the rotation periods of very hot and very cool stars and allowed us to predominantly use isochrone fitting to measure the ages of subgiants. The full likelihood used in our model was the product of these two likelihood functions, The inference processes proceeded as follows. First, a set of parameters-age, EEP, metallicity, distance, and extinctionand observables for a single star were passed to the isochronal likelihood function in Equation (9). Then, a set of observables corresponding to those parameters were generated from the MIST model grid using isochrones.py (Morton 2015) and compared to the measured observables via the isochronal likelihood,  iso (also computed using isochrones.py). The parameters were also passed to the gyrochronology model (Equation (3)), where t, E, [M/H], D, and A V were used to calculate the G BP -G RP color and mass from the MIST model grid and, in turn, the rotation period via the gyrochronology model. EEP and G BP -G RP were also used to calculate the additional rotation period variance, added to the individual period uncertainties. This model rotation period was compared to the measured rotation period using the gyrochronal likelihood function of Equation (10). The gyrochronal loglikelihood was added to the isochronal log-likelihood to give the full likelihood, which was then added to the log prior to produce a single sample from the posterior PDF.
Ages were inferred with stardate using Markov Chain Monte Carlo. The joint posterior PDF over age, mass, metallicity, distance, and extinction was sampled using the affine invariant ensemble sampler, emcee (Foreman-Mackey et al. 2013), with 50 walkers. Samples were drawn from the posterior PDF until 100 independent samples were obtained. We actively estimated the autocorrelation length, which indicates how many steps were taken per independent sample, after every 100 steps using the autocorrelation tool built into emcee. The MCMC concluded when either 100 times the autocorrelation length was reached and the change in autocorrelation length over 100 samples was less than 0.01, or the maximum of 500,000 samples was obtained. This method is trivially parallelizable, as the inference process for each star can be performed on a separate core. The age of a single star can be inferred in around 1 hr on a laptop computer.

Results
In order to demonstrate the performance of our method, we conducted three sets of tests. In the first, we simulated observables from a set of stellar parameters for a few hundred stars using the MIST stellar evolution models and the gyrochronology model of Equation (3). The ages predicted with our model were compared to the true parameters used to generate the data. In the second, we tested our model by measuring the ages of individual stars in the NGC 6819 open cluster, and in the third, we tested our model on Keplerasteroseismic stars.

Test 1: Simulated Stars
For the first test, we drew masses, ages, bulk metallicities, distances, and extinctions at random for 1000 stars from the following uniform distributions: T eff ; log g;F; parallax; and apparent magnitudes B, V, J, H, and K and Gaia G, G BP , and G RP were generated from these stellar parameters using the MIST stellar evolution models. We added a small amount of noise to the "observed" stellar properties in order to reflect optimistic observational uncertainties for isochrone dating. We added Gaussian noise with a standard deviation of 25 K to T eff , 0.01 dex to [Fe/H] and log g, and 10 mmag to B, V, J, H, and K magnitudes. These are just one choice of uncertainties that we could have adopted and are extremely optimistic. The uncertainties on predicted ages will depend strongly on all observational uncertainties; however, because this analysis is designed to show the relative improvement in stellar age precision when rotation periods are included, we chose to use best-case spectroscopic parameters. The noise added to Gaia G-band photometry ranged from 0.3 mmag for stars brighter than 13th magnitude to 10 mmag for stars around 20th magnitude (Evans et al. 2017;Gaia Collaboration et al. 2018). Noise added to the Gaia G BP and G RP bands ranged from 2 mmag for stars brighter than 13th magnitude to 200 mmag for stars fainter than 17th. Unphysical combinations of stellar parameters were discarded, resulting in a final sample size of 841 simulated stars. Figure 3 shows the position of these stars on an HRD (with log g on the y-axis instead of luminosity to improve the visibility of the MS), colored by their age. Rotation periods for these stars were generated using the gyrochronology relation described in Equation (3). We added a 5% Gaussian noise to all stellar rotation periods to represent realistic measurement uncertainties of 5%. The median uncertainty on rotation periods calculated from Kepler light curves, provided in the McQuillan et al.
(2014) catalog, is 1%. However, the Aigrain et al. (2015) injection and recovery study showed that true rotation period uncertainties are often slightly larger than this, and the noise distribution of rotation periods can be highly non-Gaussian (e.g., Aigrain et al. 2015;Angus et al. 2018). Figure 4 shows the rotation periods of 841 stars generated from the gyrochronology model.
We took two approaches to inferring the ages of these simulated stars: first using isochrone fitting only, and second using isochrone fitting combined with a gyrochronology model (stardate). Because the posterior PDFs of stars are often multimodal, we found that the choice of initial positions of the emcee walkers influenced the final outcome because walkers occasionally got stuck in local minima. We found that the following set of initial parameters worked well, though not perfectly: EEP=330, t=9.56 Gyr, [M/H]=−0.05, D= 269 pc, and A V =0.0 15 . Figure 5 shows the results of combining gyrochronology with isochrone fitting for the simulated sample. The stars' true ages are plotted against their predicted ages, with ages inferred with gyrochronology and isochrone fitting in color, and ages inferred using isochrone fitting only plotted in light gray. The different panels show the results for different types of stars: FGK dwarfs that are still undergoing magnetic braking, FGK dwarfs that have ceased magnetic braking (their Rossby number is around 2), M dwarfs, and evolved stars. The selection criteria for these groups are in the panel headings. The FGK dwarfs with low Rossby numbers showed the largest improvement: the median age precision for this group (defined as the standard deviation of the posterior as a percentage of the median age) was 8% when using a combination of isochrones and gyrochronology, and 22% using isochrone fitting alone. This equates to an almost 3× improvement in age precision. The age rms of this group was 0.8 Gyr using isochrones and gyrochronology, and 2 Gyr using isochrones only. Despite the fact that stars with Rossby numbers of 2 have stopped spinning down and their rotation periods no longer evolve with age, the ages of these stars can still be relatively precisely constrained with isochrone fitting. The median age precision for stars with large Rossby numbers was 11% with gyrochronology and isochrone fitting, and 13% with isochrone fitting only. The precision of M dwarf ages did improve overall when their rotation periods were included, but this improvement was entirely driven by the early M dwarfs. The precision of this group improved overall from 33% to 22% and rms from 5.4 to 4.6 Gyr. The precision of ages inferred for evolved stars changed very little when rotation periods were included in the inference process. This is because the variance of the rotation-age relation was inflated by a large amount in the gyrochronology model (Equation (3)), making rotation periods almost entirely uninformative for this group. The median age precision of subgiants from both gyrochronology and isochrone fitting and isochrone fitting alone was 7%.
An important caveat associated with these results is that they strongly depend on the uncertainties adopted for all observables Figure 3. The simulated star sample plotted on an HRD, colored by age (top panel) and rotation period (bottom panel). HRD positions were calculated using MIST isochrones via the isochrones.py Python package, and rotation periods were generated using Equation (3). This figure was generated in a Jupyter notebook available athttps://github.com/RuthAngus/stardate/blob/ master/paper/code/Simulate_data.ipynb. Figure 4. Data simulated from the rotation period model. Late F, GK, and early M dwarfs (stars with 0.56<G BP -G RP <2.7 follow the Praesepe-calibrated gyrochronology relation (dashed gray lines), with the exception of old, slowly rotating stars with large Rossby numbers whose rotation periods are fixed at 2× their convective overturn time. The rotation periods of early F (G BP -G RP <0.56) and late M dwarfs (G BP -G RP <2.7) and subgiants (EEP420) were generated from a log-normal distribution with standard deviation given by Equation (3). The top panel shows the rotation periods versus G BP -G RP colors of simulated stars, colored by their age, and the bottom panel shows the same stars colored by their equivalent evolutionary phase (EEP). The gray lines describe the mean gyrochronology model at ages 1, 3, 5, 7, 9, 11, and 13 (rotation periods rise with age). This figure was generated in a Jupyter Notebook, available athttps://github.com/RuthAngus/stardate/blob/master/ paper/code/Simulate_data.ipynb.
used in the analysis: T eff , [Fe/H], log g, G, G BP , G RP , J, H, K, B, V, and rotation period. Changing the uncertainties on these observables will affect the uncertainties on inferred ages in different ways. This demonstration is not intended to reflect the typical age uncertainties that will result for all stars in realityit merely exemplifies the increase in stellar age precision that results from one specific choice of uncertainties. Needless to say, estimating the uncertainties on observables accurately can be as important as accurately measuring the observables themselves. This simulation experiment was designed to show the theoretical improvement in age precision when gyrochronology is incorporated into isochrone fitting. However, it does not demonstrate the accuracy of this method because the test data were simulated from the same model used to infer ages. The results of this experiment are therefore extremely accurate by design. When applying this method to real data, the results will only be accurate if the model is accurate. In other words, stardate, like any age-dating method, provides modeldependent ages. Stellar ages calculated with stardate depend on both the accuracy of the MIST models and the accuracy of the gyrochronology model (Equation (3)). In order to test the accuracy of stardate, we applied it to real data, as described in the following section.

Test 2: Open Clusters
In order to test our model on real stars with known ages, we selected a sample of stars in the 2.5 Gyr NGC 6819 cluster. We compiled Kepler-based rotation periods (Meibom et al. 2015), Gaiaphotometry, and Gaia parallaxes for members of the NGC 6819 cluster. Figure 6 shows the period-color relation of this cluster, and Figure 7 shows the results of inferring the ages of individual cluster members using a combination of gyrochronology and isochrone fitting (via stardate) and isochrone fitting alone. The ages of F stars in this cluster (G BP -G RP ∼5.5-6.5) were relatively precisely constrained by isochrone fitting alone because, at 2.5 Gyr, they are approaching the MS turnoff. For these hot stars, ages inferred with gyrochronology and isochrones were similar to ages inferred with isochrones and similarly precise, showing that isochrones provide a lot of age information for these stars, and rotation periods do not add significantly more information. The G and early K stars in this cluster (G BP -G RP .65) were not precisely recovered from isochrone fitting alone-the isochrone-only age posteriors tend toward the prior, which is a uniform distribution between 0 and 13.8 Gyr. The median age of stars in the cluster was 4.27±0.48 Gyr when only isochrone fitting was used. In contrast, including gyrochronology when inferring the ages of G and K stars in this cluster Figure 5. The true vs. predicted ages of simulated stars. Ages calculated by combining gyrochronology and isochrone fitting with stardate are shown in color, and ages calculated with isochrone fitting only are shown in gray. The different panels show the results for stars with G BP -G RP <2.2 (FGK dwarfs) that are still braking magnetically (Ro < 2), stars with G BP -G RP <2.2 that have stopped braking magnetically (Ro 2), stars with 2.2<G BP -G RP (M dwarfs), and evolved stars (EEP>420). Gyrochronology is highly effective for FGK stars and ages inferred with both gyrochronology and isochrone fitting are more accurate and precise than ages inferred via isochrone fitting only for this group. Neither gyrochronology nor isochrone fitting can provide precise ages for M dwarfs, so the ages of these stars are imprecise regardless of the age-dating method. This figure was generated in a Jupyter Notebook available athttps://github.com/RuthAngus/stardate/blob/ master/paper/code/Results_plots.ipynb.
significantly improved age precision. The median age of stars in this cluster was 2.63±0.16 Gyr when ages were inferred with a combination of isochrone fitting and gyrochronology, using the newly calibrated Praesepe-based gyrochronology model. The previously calibrated Angus et al. (2015) model resulted in a median stellar age of 2.66±0.21 Gyr, which is still consistent with the established cluster age of 2.5 Gyr. The median age of stars in the cluster using uncorrected photometry was slightly underestimated at 1.86±0.22 Gyr. This suggests that, despite the fact that V-band extinction is marginalized over during the inference process, correcting for extinction before ages are estimated will reduce bias introduced by dust.
We found that 5% rotation period uncertainties resulted in the most accurate ages for NGC 6819. The uncertainties on the measured rotation periods, provided in Meibom et al. (2015) and shown in Figure 6, were likely underestimated for some stars. Underestimated rotation period uncertainties can result in inaccurate age estimates. This raises the question: how should uncertainties on rotation periods be estimated? The likelihood is weighted by the inverse variance, so uncertainties on the rotation period control the relative information provided by gyrochronology, isochrones, and the prior. If rotation period uncertainties are either too large or too small, the resulting age estimate will be imprecise and/or inaccurate. It is difficult to measure uncertainties on rotation periods directly: standard techniques such as Lomb-Scargle periodograms and autocorrelation functions do not provide them. Ideally, rotation period uncertainties should capture both the measurement precision, and the physical uncertainty introduced by the latitudinal movement of star spots on the surface of a differentially rotating star. For example, Donahue et al. (1996) demonstrated that the seasonal variation in measurements of G and K star rotation periods is a function of period, D µ  P P rot rot 1.3 0.1 . This variation is presumably caused by a latitudinal drift in the dominant active regions, which traces the stellar cycle over several years, in combination with latitudinal differential rotation. This suggests that latitudinal spot drifting does not significantly affect stellar rotation periods when stars are young-for example, the scatter of rotation periods about the mean gyrochronology model in Praesepe is only around 5%. However, it is likely that this effect will become more important at older ages. A thorough exploration of how rotation period uncertainties, from both measurement uncertainty and physical variation, affect stellar ages via gyrochronology is key to understanding the power of gyrochronology as an age-dating method. For now, we leave this exploration for a future study.

Test 3: Kepler Asteroseismic Stars
In order to test our method in the regime where both isochrone fitting and gyrochronology become important, we recovered the ages of the 21 asteroseismic stars analyzed in van . These 21 stars were observed in Kepler's short cadence mode and are a mixture of dwarfs and subgiants. Their asteroseismic ages were calculated from the analysis of the frequencies of individual oscillation modes Metcalfe et al. 2014;Silva Aguirre et al. 2015;Ceillier et al. 2016) and their rotation periods from their Kepler light curves (García et al. 2014). We crossmatched these stars with the Gaia catalog to obtain parallaxes and apparent magnitudes in the Gaia G, G BP , and G RP bandpasses. We also added J, H, and K 2MASS magnitudes from the Kepler input catalog (Brown et al. 2011) and used spectroscopic effective temperatures, spectroscopic metallicities, and rotation periods reported in Table1 of van Saders et al. (2016). The Gaia photometry is extremely precise, and we found that artificially inflating the uncertainties on Gaia apparent magnitudes by a factor of 10 substantially improved the quality of fit, both in terms of MCMC convergence and agreement with asteroseismic age measurements. In Figure 8, we show the ages of these 21 stars inferred using isochrone fitting and gyrochronology, against their asteroseismic ages, calculated using the Asteroseismic Modeling Portal (AMP; Metcalfe et al. 2009Metcalfe et al. , 2012Metcalfe et al. , 2014. In this figure, the colored symbols show ages inferred using a combination of gyrochronology and isochrone fitting, Figure 6. The Kepler-based rotation periods of members of the 2.5 Gyr NGC 6819 open cluster. The raw G BP -G RP colors are shown in red, and the dustcorrected colors are shown in black. The dashed line shows a gyrochronology model that was fit to the Praesepe cluster and the Sun in this work, interpolated to 2.5 Gyr. The solid blue line shows a previously calibrated gyrochronology model (Angus et al. 2015). This figure was generated in a Jupyter Notebook available athttps://github.com/RuthAngus/stardate/blob/master/paper/code/NGC6819. ipynb. Figure 7. The inferred ages of members of the NGC 6819 open cluster as a function of their G BP -G RP color. Ages of stars inferred using a combination of isochrone fitting and gyrochronology (Praesepe and Sun calibration) with dereddened Gaia G, G BP , and G RP photometry (black circles) and uncorrected, raw, photometry (red squares). Even though V-band extinction is marginalized over in the inference process, reddening can still bias ages. Blue triangles, pointing up, show ages inferred using isochrone fitting and gyrochronology, with the Angus et al. (2015) gyrochronology model. Orange triangles, pointing down, show ages inferred using isochrone fitting only. The ages of F stars (stars bluer than 0.7) were precisely constrained by isochrones and including gyrochronology makes little difference in their inferred ages. The age precision of G and K dwarfs (stars redder than 0.7) was improved by including gyrochronology. The median age of stars inferred using the gyrochronology model calibrated to Praesepe and the Sun (black circles) was 2.65±0.13, which is consistent with the established cluster age (2.5 Gyr). This figure was generated in a Jupyter Notebook available athttps://github.com/RuthAngus/ stardate/blob/master/paper/code/NGC6819.ipynb.
implemented with the stardate Python package. The black and gray triangles show ages inferred from gyrochronology only, where the mass and color of the stars were not inferred, but fixed to be the asteroseismic mass and the Gaia G BP -G RP color. These gyrochronal ages were calculated with age as the only free parameter, without marginalizing over stellar mass, G BP -G RP color, or extinction. As a result, these gyrochronal ages are not the same as the ages given by the maximum of the gyrochronal likelihood function used in the combined age model; they simply represent an approximation to the gyrochronal age. Gray triangles are shown for stars where gyrochronology is not applicable because the stars are either too metal poor, too metal rich, or too evolved. Black triangles are shown for stars where gyrochronology is applicable. The white circles show ages inferred from isochrones only. Dashed lines connect the three different age measurements for the same stars.
Although the ages of all 21 stars shown in Figure 8 were inferred with a joint isochronal and gyrochronal model, most (all but eight) were either too evolved, too metal poor, or too metal rich for gyrochronology to contribute any information to the ages. These metal-poor/-rich or evolved stars lie in a regime where the variance on their rotation period was artificially inflated because the gyrochronology relations are not thoroughly understood or well calibrated. The rotation periods of the remaining eight stars did contribute to their inferred ages to some degree; however, isochrones still dominated the age information for some of them because these stars are relatively old and/or relatively hot.
In general, there is relatively poor agreement between the asteroseismic ages and the ages inferred using stardate. Much of this discrepancy is driven by differences in the isochronal ages, which is likely attributable to differences between the MIST stellar evolution models and those used in the AMP analysis: a combination of the Aarhus stellar evolution code (ASTEC Christensen-Dalsgaard 2008a) and the adiabatic pulsation code (ADIPLS Christensen-Dalsgaard 2008b). We compared nonrotating, solar-metallicity MIST isochrones for middle-aged stars with solar-metallicity BaSTI isochrones (Pietrinferni et al. 2004;Hidalgo et al. 2018) and found that, for stars between 4 and 8 Gyr, in the same effective temperature range as the asteroseismic stars, the age discrepancy between the two sets of models can be as large as 1-2 billion years. The MIST isochrones lie above the BaSTI isochrones on the HRD, leading to a systematic underprediction of ages.
The gyrochronal ages, where gyrochronology is applicable, do not show excellent agreement with the asteroseismic ages either. The four hot stars to the left in Figure 8 are rotating more slowly than predicted by the Praesepe-based gyrochronology models, and as a result, their gyro-ages are older than their asteroseismic ones. For these hot stars, two out of four have ages that are still consistent, or close to consistent, with their asteroseismic ages. The third star from the left is an anomalously slow rotator for its age and mass, and van  also found this star to be surprisingly slowly rotating. In contrast, the star with a black triangle symbol (indicating that gyrochronology is applicable) farthest to the right is rapidly rotating for its mass and age, even when weakened braking is taken into account. This star is around solar mass (1.0 ± 0.03 M e ) and has an asteroseismic age older than the Sun (7.28 ± 0.51 Gyr), yet rotates with a period of only 19.8±1.3 days. This is KIC 9098294, a single-lined spectroscopic binary with an orbital period of around 20 days

Conclusions
We present a statistical framework for measuring precise ages of MS stars and subgiants by combining observables that relate, via different evolutionary processes, to stellar age. Specifically, we combined HRD/CMD placement with rotation periods, in a hierarchical Bayesian model, to age-date stars based on both their hydrogen burning and magnetic braking history. The two methods of isochrone fitting and gyrochronology were combined by taking the product of two likelihoods: one that contains an isochronal model and the other a gyrochronal one. We used the MIST stellar evolution models and computed isochronal ages and likelihoods using the isochrones Python package. We fit a new broken power-law gyrochronology model to the Praesepe cluster and included a modification recommended by van Saders et al. (2016) that accounts for weakened magnetic braking at Rossby numbers larger than 2. The rotation periods of hot stars, cool stars, evolved stars, young stars, and metal-poor and -rich stars were modeled with a broad log-normal distribution. We tested stardate on simulated data and cluster stars and demonstrated that combining gyrochronology with isochrone fitting improves the precision of age estimates for FGK dwarfs by a factor of 3 over isochrone fitting alone, assuming 5% measurement uncertainties on rotation periods. Incorporating rotation periods into stellar evolution models also improves the precision of the EEP parameter and because EEP, combined with age and metallicity, determines the mass, radius, and log g of a star, this means rotation periods can improve the precision of all stellar parameters. Although V-band extinction is marginalized over during inference, correcting photometry for dustextinction before analysis or including it as a prior can improve the accuracy of stellar ages measured with stardate.
We also tested stardate on a set of 21 Kepler asteroseismic stars . We found that Figure 8. A comparison of stellar ages inferred using asteroseismic modeling with ages inferred using a combination of isochrone fitting and gyrochronology. Colored circles show ages inferred using isochrone fitting and gyrochronology combined via the stardate software package. Black triangles show the ages of all stars inferred via gyrochronology only, and white circles show ages inferred via isochrone fitting only. discrepancies between ages measured with stardate and ages measured with asteroseismology are likely produced by differences between the MIST and BaSTI stellar evolution models. Asteroseismic and cluster stars provide an opportunity for calibration, but given the high dimensionality of the gyrochronology relations (i.e., rotation period depends on age, mass, metallicity, surface gravity, etc.), many stars with precise ages, spanning a range of properties, are still needed to reliably calibrate them.
In cases where gyrochronology predicts inaccurate stellar ages, it is either because models are not correctly calibrated, because the rotation periods or rotation period uncertainties are themselves inaccurate, or because of rotational outliers. For example, stardate may predict inaccurate ages for stars in close binaries whose interactions influence their rotation period evolution. Rotational outliers are often seen in clusters (see, e.g., Douglas et al. 2016;Rebull et al. 2016Rebull et al. , 2017Douglas et al. 2017), and many of these fall above the main sequence on a CMD, indicating that they are binaries. In addition, measured rotation periods may not always be accurate and can, in many cases, be a harmonic of the true rotation period. For example, a common rotation period measurement failure mode is to measure half of the true rotation period. The best way to prevent an erroneous or outlying rotation period from resulting in an erroneous age measurement is to allow for outlying rotation periods using a mixture model, a feature that could be built into stardate in the future.
The optimal way to age-date stars is by combining all their available age-related observables. This could ultimately include activity dating via flare rates and chromospheric activity indices, kinematic dating, and chemical dating. Of all the established age-dating methods, gyrochronology and isochrone fitting are two of the most complementary. The two methods are optimal in different parts of the HRD: gyrochronology works well for FGK dwarfs and isochrone fitting works well for subgiants and hot stars, so combining the two methods results in consistently precise ages across a range of masses, ages, and evolutionary stages. In addition, using both methods at once circumvents the need to decide which method to use a priori. It eliminates the circular process of classifying a star based on its CMD position (M dwarf, subgiant, etc.), then deciding which age-dating method to use, then inferring an age which itself depends on the classification that was made. It is important to infer all stellar properties at once because they all depend on each other. Stardate is applicable to a large number of stars: FGKM dwarfs and subgiants with a rotation period and broadband photometry. This already includes tens of thousands of Kepler and K2 stars and could include millions more from TESS, LSST, WFIRST, PLATO, Gaia, and others in the future.
The code used in this project is available as a Python package called stardate, doi:https://doi.org/10.5281/ zenodo.2712419. It is available for download via Github 16 or through PyPI 17 Documentation, available athttps://stardate. readthedocs.io/en/latest/. All codes used to produce the figures in this paper are available athttps://github.com/RuthAngus/ stardate. This paper is based on code with the following Git hash: f739562c1546e117e9bb217e1732c62b41be8061.
Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is Priors We used the default priors in the isochrones.py Python package. The prior over age was where A is log 10 (Age [yr]). The prior over EEP was uniform with an upper limit of 800. We found that adding this upper limit reduced some multimodality caused by the giant branch and resulted in better performance. The prior over true bulk metallicity was based on the galactic metallicity distribution, as inferred using data from the Sloan Digital Sky Survey (Casagrande et al. 2011). It is the product of a Gaussian that describes the metallicity distribution over halo stars and two Gaussians that describe the metallicity distribution in the thin and thick disks: where H F =0.001 is the halo fraction, μ halo and σ halo are the mean and standard deviation of a Gaussian that describes a probability distribution over metallicity in the halo and take values −1.5 and 0.4, respectively. The two Gaussians inside the square brackets describe probability distributions over metallicity in the thin and thick disks. The values of the means and standard deviations in these Gaussians are from Casagrande et al. (2011). ξ is the integral of everything in the square brackets from -¥ to ¥ and takes the value ∼2.507. The prior