In This Day and Age: An Empirical Gyrochronology Relation for Partially and Fully Convective Single Field Stars

Gyrochronology, the field of age dating stars using mainly their rotation periods and masses, is ideal for inferring the ages of individual main-sequence stars. However, due to the lack of physical understanding of the complex magnetic fields in stars, gyrochronology relies heavily on empirical calibrations that require consistent and reliable stellar age measurements across a wide range of periods and masses. In this paper, we obtain a sample of consistent ages using the gyro-kinematic age-dating method, a technique to calculate the kinematics ages of stars. Using a Gaussian process model conditioned on ages from this sample (∼1–14 Gyr) and known clusters (0.67–3.8 Gyr), we calibrate the first empirical gyrochronology relation that is capable of inferring ages for single, main-sequence stars between 0.67 and 14 Gyr. Cross-validating and testing results suggest our model can infer cluster and asteroseismic ages with an average uncertainty of just over 1 Gyr, and the inferred ages for wide binaries agree within 0.83 Gyr. With this model, we obtain gyrochronology ages for ∼100,000 stars within 1.5 kpc of the Sun with period measurements from Kepler and Zwicky Transient Facility and 384 unique planet host stars. A simple code is provided to infer gyrochronology ages of stars with temperature and period measurements.


INTRODUCTION
Gyrochronology (Barnes 2003) is a method to agedate stars mainly using their rotation periods (P rot ) and mass/temperature (T eff ) measurements.It is based on the principle that stars lose angular momentum through magnetized winds and therefore, spin down with time (Kraft 1967).The simplest form of Gyrochronology relation is discovered by Skumanich (1972), stating that Unfortunately, This simple picture is heavily challenged by the emergence of large photometric surveys in the recent decade such as Kepler (Borucki et al. 2010), K2 (Howell et al. 2014), TESS (Ricker et al. 2015), MEarth (Berta et al. 2012a), andZTF (IRSA 2022a,b).These photometric surveys provided valuable data to measure stellar rotation in mass quantities (e.g.McQuillan et al. 2013McQuillan et al. , 2014;;García et al. 2014;Santos Corresponding author: Yuxi(Lucy) Lu lucylulu12311@gmail.com et al. 2019, 2021;Gordon et al. 2021;Lu et al. 2022;Holcomb et al. 2022;Claytor et al. 2023).These catalogs show sub-structures in the density distribution of stars in P rot -T eff space, suggesting not all stars spin down "Skumanich style".Some of the discoveries include: the upper boundary or pile-up of solar-like stars with intermediate ages (Angus et al. 2015;Hall et al. 2021;David et al. 2022) that could be caused by weakened magnetic braking (e.g.van Saders et al. 2016;Metcalfe et al. 2022) or perhaps the transition of latitudinal differential rotation (Tokuno et al. 2022); the intermediate period gap in partially convective GKM dwarfs (McQuillan et al. 2013;Gordon et al. 2021;Lu et al. 2022) most likely caused by stalled spin-down of low-mass stars (Curtis et al. 2020;Spada & Lanzafame 2020); the bi-modality of fast and slow-rotating M dwarfs that is difficult to explain with traditional models of angular-momentum loss (Irwin et al. 2011;Berta et al. 2012b;Newton et al. 2017;Pass et al. 2022;Garraffo et al. 2018); the abrupt change in stellar spin-down across the fully convective boundary (Lu et al. 2023, Chiti et al. in prep.).Therefore, modern-day gyrochronology heavily relies on empirical calibrations with benchmark stars such as those with asteroseismic ages (e.g.Angus et al. 2015;Hall et al. 2021), those in wide binaries (Pass et al. 2022;Silva-Beyer et al. 2022;Otani et al. 2022;Gruner et al. 2023, Chiti et al. in prep.), and open cluster members (e.g.Curtis et al. 2020;Agüeros et al. 2018;Gaidos et al. 2023;Dungee et al. 2022;Bouma et al. 2023).Asteroseismic ages can be accurate and precise to the 10% level with time series from Kepler.Unfortunately, asteroseismic signal strength/frequency decreases/increases dramatically as the mass of a star decreases, and no signals have been detected for low-mass M dwarfs.Open clusters are generally young as they typically dissolve in the Milky Way on a time-scale of ∼ 200 Myr.Much effort has been put into calibrating gyrochronology with wide binaries, however, no large catalog of consistent ages for wide binary stars currently exists.As a result, none of the above benchmark stars can provide a consistent sample of reliable ages for stars of vastly different masses and periods that can be used to calibrate empirical gyrochronology relations across a wide range of ages.
Recently, gyro-kinematic age-dating (Angus et al. 2020;Lu et al. 2021a), a method to obtain kinematic ages from stars with similar P rot -T eff -M G -Rossby Number (Ro; P rot divided by the convective turnover time), provide an opportunity to obtain a consistent benchmark sample for calibrating a fully empirical gyrochronology relation.One discovery using the ages obtained from this method is the fundamentally different spin-down law for fully and partially convective stars (Lu et al. 2023), as a result, it is important to obtain gyrochronology ages separately for partially and fully convective stars.By combining period measurements from Kepler and ZTF, we obtain gyro-kinematic ages for ∼ 50,000 stars and present the first fully empirical gyrochronology relation that is able to infer ages for single main-sequence stars of age 0.67-14 Gyr.In section 2, we describe the dataset, the method used to calibrate this gyrochronology relation, and the cross-validation test.
In section 3, we present the testing set and a catalog of ∼ 100,000 stars with gyrochronology ages.In section 4, we discuss the limitations, including the effect of metallicity, and future improvements.We de-reddened G BP − G RP , M G measurements from Gaia DR3 using dustmap (Green 2018;Green et al. 2018).The temperature is then calculated from G BP − G RP using a polynomial fit taken from Curtis et al. (2020).Ro is calculated as Ro=P rot /τ c , in which τ c is the convective turn-over time that depends only on the temperature of the star (See et al. in prep.).
We obtained rotation periods for ZTF stars with Gaia G band magnitude between 13 to 18 and G BP − G RP < 1 using the method described in Lu et al. (2022).For ZTF stars with G BP − G RP > 1, we adapted the rotation period measurements (before vetting) from Lu et al. (2022).From the full sample, we selected stars with agreeing periods from at least 2 seasons.By comparing 1,270 overlapping period measurements from ZTF and Kepler (Santos et al. 2021), we found an 81% agreement within 10% for stars with ZTF period measurements > 4 days (see Figure 1).As a result, we rejected stars with ZTF period measurements < 4 days.To roughly select dwarf stars, we also excluded stars with M G < 4.2 mag.This yielded ∼ 55,000 ZTF stars with RV measurements from Gaia DR3 (Gaia Collaboration et al. 2021).Combining ∼ 30,000 Kepler stars with M G > 4.2 mag from Lu et al. (2021a) with RV measurements from Gaia DR3, LAMOST (Cui et al. 2012), and inferred RV from Angus et al. (2022), we obtained a total of ∼ 85,000 stars with RV and relatively reliable period measurements (See Figure 2 top plot).(Santos et al. 2021) and ZTF (Lu et al. 2022, This work).We found a 81% agreement within 10% for stars with ZTF period measurements > 4 days.
We then excluded equal-mass binaries by fitting a 6 thorder polynomial (f 6 (T eff )) to the entire sample and only selecting stars with M G > f 6 (T eff ) − 0.4 (shifted by eye).We also excluded stars with Ro > 10.This left us with a final sample of 68,378 stars (ZTF: 49,928; Kepler:  (McQuillan et al. 2014;García et al. 2014;Santos et al. 2019;Lu et al. 2022, this work).The red dashed line shows the shifted 6 th -order polynomial (f6(T eff )) fitted to the entire sample that separates the equal-mass binaries from the rest of the sample.Middle: similar to the top plot but after excluding equal-mass binaries (a total of 68,378 dwarf stars).Bottom: period distribution of the 68,378 dwarf stars.
18,450).The period distribution for the final sample is shown in the bottom plot of Figure 2. The overall period distribution agrees with that of McQuillan et al. (2014); Santos et al. (2021), except for an over-density at ∼ 4,000 K with P rot < 10 days.Since we did not impose conservative vetting criteria, this over-density is most likely caused by systematic.We also see a systematic over-density at ∼ 30 days, this is a known systematic in ZTF, which is caused by the orbit of the moon (Lu et al. 2022).

Gyro-kinematic age data
We determined gyro-kinematic ages following the procedure described in Lu et al. (2021b), where the vertical velocity dispersion for each star is calculated from vertical velocities of stars that are similar in P rot , T eff , M G , and Ro to the targeted star.We then converted the velocity dispersion measurements into stellar ages using an age-velocity-dispersion relation in Yu & Liu (2018).The vertical velocities are calculated from Gaia DR3 proper motions (Gaia Collaboration et al. 2021) and RVs from various sources (data sample see Section 2.1.1)by transforming from the Solar system barycentric ICRS reference frame to Galactocentric Cartesian and cylindrical coordinates using astropy (Astropy Collaboration et al. 2013;Price-Whelan et al. 2018).The bin size to select similar stars to the targeted star in order to calculate gyro-kinematic ages was [T eff , log 10 (P rot ), R o , M G ] = [177.8K, 0.15 dex, 0.15 dex, 0.2 mag].This bin size is optimized by performing a grid search in the binning parameters (T eff , log 10 (P rot ), R o , M G ) and minimizing the total χ 2 in predicting individual cluster ages > 1.5 Gyr with M G > 4.2 and R o < 2 (data sample see Section 2.1.2).We did not use clusters with age < 1.5 Gyr in this process as gyro-kinematic ages for stars < 1.5 Gyr is heavily contaminated by binaries, and will overestimate cluster ages and produce unreliable results (See Figure 4 or A.1 in Lu et al. 2021b).Figure 4 shows the final optimization result.We excluded stars with gyro-kinematic age < 1.5 Gyr or > 14 Gyr as it is possible that a significant number of the youngest stars have not yet converged onto the slow-rotating sequence, and those that are very old are likely outliers.The sample of 46,362 stars with corrected gyro-kinematic ages between 1.5 and 14 Gyr and cluster ages between 0.67 Gyr and 4 Gyr are shown in Figure 5 top plot.

Gaussian Process
Gaussian Processes (GP) is a generic supervised learning method designed to solve regression or classification problems.It has been applied frequently in timedomain astronomy (e.g.Foreman-Mackey et al. 2017;Angus et al. 2018;Gilbertson et al. 2020) as it can model the covariance between the noise in the data.Typically, a GP regressor is composed of a mean function (m; Equation 1), which is ideally physically motivated, and a covariance function (k; Equation 2) that captures the details that the mean function has missed.For a more detailed review of GP and its applications in astronomy, we direct the readers to Aigrain & Foreman-Mackey (2022).In this paper, we used the PYTHON package tinygp (Foreman-Mackey 2023) to construct our GP model.
Since there is an abrupt change in the spin-down law across the fully convective boundary (Lu et al. 2023), we fitted separate GP relations to the partially and fully convective stars.The division was made using the gap discovered in the color-magnitude-diagram (CMD).This gap is an under-density in the CMD near the fully convective boundary and can be approximated by a line connecting et al. 2018).It is thought to be caused by structural instabilities due to the non-equilibrium fusion of 3 He (van Saders & Pinsonneault 2012; Baraffe & Chabrier 2018;MacDonald & Gizis 2018;Feiden et al. 2021).
As fitting a multi-dimensional GP requires a large amount of computational resources, and it is not possible to fit to all ∼46,000 stars with gyro-kinematic ages within a reasonable amount of time, we constructed the final training sample by dividing the stars with gyrokinematic ages into bins with size [T eff , log 10 (P rot )] ∼ [50 K, log 10 (1 Days)] and calculating the median age in each bin if more than 10 stars were included.The fit was done separately for fully and partially convective stars as some of them overlap in T eff -P rot space.The temperature bin size is chosen based on the estimated uncertainty in temperature measurements of ∼ 50K, and the period bin size is chosen so that we can obtain enough training samples for the GP.The uncertainties associated with the training sample are measured with the standard deviation on the gyro-kinematic ages for stars in each bin.We then added all the individual cluster stars to the training sample and inflated their age uncertainty to be 0.5 Gyr to ensure a smooth GP fit.We found that using the true cluster age uncertainties reported in the literature, the GP over-fits the cluster data.The training sample for the partially (cir-  Classical gyrochronology relations assume the age of a star can be approximated with a separable function in temperature and period, and we constructed our mean function motivated by this relation.We formulated the mean function to be a double broken power law in P rot and T n for partially convective stars to capture the sudden increase of rotation periods of M dwarfs at ∼3,500 K and the plateauing of the rotation periods for G/K stars at ∼5,000 K.We define T n as the normalized temperature given by T n := (7000 − T eff )/(7000 − T break ) for the partially convective stars, and T n := (3500 − T eff )/500 for the fully convective stars, in which T break is the temperature at which the temperature power law changes.For fully convective stars, we used a single power law in T n and a double broken power law in rotation period, P or P rot , since the temperature range for the fully convective stars is small.
In equations, The mean function is defined to be: where the broken power law in rotation, f (P rot ), is defined to be: where P break is the rotation period at which the rotation power law changes.S P h (P rot , w P ) and S P l (P rot , w P ) are the smoothing functions in period space, defined to be: The broken power law in temperature, g(T n ), is defined to be: where S T h (T n , w T ) and S T l (T n , w T ) are the smoothing functions in temperature space, defined to be: The bold letters show the variables that will be fitted to the data, they are defined in Table 1.The smoothing functions (e.g. S P h and S T h ) can be viewed as switches for the broken power laws.w P and w T dictate how smooth the broken power laws are (e.g.w P =0 or w T =0 indicate a sharp transition between the power laws).Since the fully convective stars only span a small range in temperature, we used a single power law, so that g( Figure 6.Age predicted from the mean function calculated using the initial values shown in Table 1.The red points show the cluster star sample.The mean function is flexible enough to capture the cluster shapes.
For the covariance function of the GP model, we used a 2-D uncorrelated squared exponential kernel, meaning we assume no correlation between the temperature and period measurements.The function is defined to be: (2) where T eff , T ′ eff are two different data points in temperature space (same for log 10 (P rot ) and log 10 (P rot ) ′ ).l T and l P determine the length scale of the correlation between temperature measurements and period measurements, respectively.σ 2 determines the strength of the correlation.In other words, the covariance function determines how the response at one temperature and period point is affected by responses at other temperature and period points.
The initial values for the parameters used in the mean function and covariance function before optimizing are shown in Table 1. Figure 6 shows the mean function (background) calculated with the initial values and the cluster members overlayed on top (red points).
We built the GP model using tinygp (Foreman-Mackey 2023).tinygp is a PYTHON library for building GP models.It is built on jax (Bradbury et al. 2018), which supports automatic differentiation that enables efficient model training.We first optimized the parameters by maximizing the log-likelihood function, conditioned on the data described at the beginning of this section.The optimized parameters were then used as ini-tial inputs for the Markov chain Monte Carlo (MCMC) model to obtain the true distributions for the parameters.The priors are Gaussians centered around the optimized parameters with a width described in Table 1.We implemented the MCMC model in numpyro (Phan et al. 2019) for partially and fully convective stars separately.The best-fit parameters for partially and fully convective stars are shown in Table 1, and the corner (Foreman-Mackey 2016) plots are shown in Figure 7 and Figure 8 for partially and fully convective stars, respectively.

Cross-validation
To ensure our model did not over-fit the data, we performed the cross-validation test by first excluding a random 20% of the gyro-kinematic ages sample and optimized the model following the procedure described in the last section.The ages of these stars were then predicted using the trained model.We also carried out a leave-one-out cross-validation test for the cluster sample by excluding a single cluster at a time, retraining the model, and predicting the age of that cluster with the trained model.The cross-validation results are shown in Figure 9.The x-error bars for the cluster sample are taken from the literature, and the y-error bars are calculated by taking the standard deviation of the predicted ages of all the cluster members.The average standard deviation (y-error bars) for the cluster cross-validation result is 0.62 Gyr.The bias and variance for the cluster sample are -0.24Gyr and 0.43 Gyr, respectively, and those for the gyro-kinematic ages are 0.37 Gyr and 0.85 Gyr, respectively.The cross-validation results suggest that our model is able to predict ages within ∼ 1 Gyr for main-sequence stars with reliable P rot , G BP − G RP , and M G measurements.However, there exists a systematic at ∼1 Gyr in predicting gyro-kinematic ages, this systematic is most likely caused by the fact that the cluster sample between 0.67 Gyr to 1 Gyr occupies similar P rot -T eff space (see Figure 3), creating degeneracy in age predictions for stars younger than 1 Gyr.As a result, age predictions for stars < 1 Gyr might be biased.Stars around this age also occupies the P rot -T eff space where stars are expected to go through stalled spin-down (Curtis et al. 2020).Looking only at stars > 5,000 K greatly reduces the pile-up.

RESULT
Figure 10 shows the modeled isochrones for the cluster sample (left column) and the isochrones for 0.7 to 10 Gyr with 1.55 Gyr separations (right column) overlaid on the training sample of stars with gyro-kinematic ages that are partially convective (top row) and fully convective (bottom row).Since unlike most gyrochronol-   1.
because there are fewer stars in those bins.In addition, since fully convective stars could spin down faster than partially convective stars (e.g.Lu et al. 2023), the bin size used to calculate gyro-kinematic ages could induce blurring as it will include stars of different ages.Interestingly, there are stars with ages that match the M67 open-cluster age in the background gyro-kinematic age sample.This suggests some stars in this temperature and period range could be mis-classified as fully convective stars.However, looking at these stars in the CMD, they are far away from the gap that is typically used to distinguish partial and fully convective stars (Jao et al. 2018;van Saders & Pinsonneault 2012).One other possibility is stars in that temperature and period range can have multiple ages.Further study of the data and fully convective stars is needed to disentangle these scenarios.

Predicting ages for the LEGACY dwarfs
To test our model, we predicted ages for 51 LEGACY dwarf stars with asteroseismic ages derived from Kepler (Silva Aguirre et al. 2017), P rot , T eff , and M G data  2017).The ages and uncertainties for the gyrochronology ages were calculated by first calculating the ages for each star using 100 realizations of the model where the parameters are taken randomly from the MCMC fit.The 16th, 50th, and 84th percentile of the age predictions were then used to calculate the lower age limit, age, and higher age limit for each star.The crosses show the stars with M G < 4.2 mag, which are outside of our training set.The bias and median absolute deviation (MAD) for the entire testing sample are -0.07Gyr and 1.35 Gyr, respectively.This test suggests our model can estimate ages for single field dwarf stars with uncertainties just over 1 Gyr.
Since our model did not take into account the effects of metallicity, we investigated this by plotting the abso- lute difference between the LEGACY and gyrochronology age against the metallicity of the star (Figure 12 left plot).There is an obvious metallicity trend for stars with [Fe/H] < 0.0 dex, suggesting future work of incorporating metallicity into this model is necessary (also see Figure 9 in Claytor et al. 2020, for how metallicity can affect age determination using gyrochronology).However, metallicity measurements that currently exist for low-mass stars are either limited in sample size or inaccurate and imprecise due to the presence of star spots and molecular lines in the spectra (e.g.Allard et al. 2011;Cao & Pinsonneault 2022).As a result, we did not attempt to include training with metallicity in this work.
As mentioned in the introduction, stars likely stop spinning down due to weakened magnetic braking after reaching a critical Rossby number, Ro crit (van Saders et al. 2016).Recently, Saunders et al. (2023) fitted a magnetic braking model to asteroseismic and cluster data and concluded that Ro crit /Ro ⊙ = 0.91±0.03,which Ro crit ∼1.866 using MESA (Paxton et al. 2019).Indeed, the gyrochronology ages show large deviations from the asteroseismic ages for stars with Ro>1.866 (Figure 12 right plot).This suggests gyrochronology models should only be used to predict ages for stars with Ro<1.866.

Gyrochronology Ages for ∼ 100,000 stars
With this new gyrochronology relation, we predicted ages for ∼ 100,000 stars from Kepler (McQuillan et al. 2014;Santos et al. 2021) and ZTF (Lu et al. 2022, This work) with period measurements, in which the ZTF periods were vetted using a random forest (RF) regressor trained on Gaia bp-rp color, absolute G magnitude, ruwe, and parallax.We did this by first training the RF on the ZTF periods that are highly vetted (Lu et al. 2022).We then used the RF to predict the periods of the ZTF stars with measured periods described in section 2.1.Finally, we selected period measurements that agree within 10% of the predicted periods, which left us with 58,462 vetted ZTF periods with bp-rp color >∼ 1.3 mag and period >∼ 10 days.We excluded stars with Ro >1.866, this left us with a final sample of 94,064 stars with periods from Kepler and ZTF.We calculated the ages by using 100 realizations of the model with parameters taken randomly from the MCMC model after the chains had converged (same as what was done for the cluster isochrones in P rot -T eff space and stars with asteroseismic ages).We also tested how the measurement uncertainty in temperature and period could affect the ages by perturbing the measurements by 50 K and 10%*P rot , respectively, assuming Gaussian errors.We then recalculated the ages using these perturbed values.We performed this 50 times for each star and found that the age uncertainty caused by the measurement error was negligible compared to the uncertainty in the model parameters.Table 2 shows the column description for this final catalog.
Figure 13 shows the histograms for stars with inferred gyrochronology ages using the calibrated relation from this work.The black histogram shows the age distribution for the full sample of ∼ 100,000 stars, the red histogram shows those with T eff < 4000 K, and the blue histogram shows those with T eff ≥ 4000 K.The black dotted lines show the recent enhancement of star formation rate (SFR) shown in Ruiz-Lara et al. (2020) (5.7, 1.9 and 1.0 Gyr).The peaks in the histograms can correspond to the enhancements of the star formation rate in the Milky Way, changes in stellar spin-down, or system-Figure 10.Running median (solid lines) and the standard deviation (shaded region) of 100 realizations of the GP models from this work for partially convective (top row) and fully convective (bottom row) stars.The models are overlaid on the full sample with gyro-kinematic ages.The Jao's gap is used to distinguish between partially convective and fully convective stars.left column: modeled isochrones (solid lines; the shaded area representing the model uncertainty) for each cluster (points).right column: Isochrones between 0.7 Gyr and 10 Gyr with a 1.55 Gyr separation colored by age.atic bias.For example, a peak exist in the tail of the distribution at the time of SFR enhancement 5.7 Gyr ago.This is the first time SFR enhancement has been shown using gyrochornology.However, some peak also correspond to limitations in the gyrochronology model.For example, the peak around 2.5 Gyr exists only in stars ≥ 4000 K.This peak most likely corresponds to the stall in spin-down for partially convective stars (e.g.Curtis et al. 2020) that do not exist for fully convective stars (T eff < 3500 K; Lu et al. 2022).The stalling is thought to happen because the surface angular momen-tum loss is replenished by the core while the core and the envelope start re-coupling.Depending on the recoupling timescale, stars that span a range of ages will have very similar rotation period measurements, meaning they will have the same inferred age based on rotation and temperature alone.Future work should include other age indicators (e.g.stellar activity) to break this degeneracy.To infer gyrochronology ages for confirmed exoplanet host stars, we downloaded data from the NASA Exoplanet Archive 1 .We combined stars with period measurements publicly available from the NASA Exoplanet Archive and from this work and inferred ages with 100 model realizations as done in the rest of this paper.We excluded stars with age prediction < 0.67 Gyr and Ro > 1.866, which left us with 384 unique planet host stars.

Gyrochronology Ages for 384 Unique Planet Host Stars
1 https://exoplanetarchive.ipac.caltech.eduWithin these stars, 338 have new rotation period measurements from Lu et al. (2022) and this work.Figure 14 shows the age distribution of these stars, and the column description for this catalog is shown in Table 3 2020) (5.7, 1.9 and 1.0 Gyr).The peaks in the histograms can correspond to the enhancements of the star formation rate in the Milky Way (e.g., the peak in the tail around 5.7 Gyr), changes in stellar spin-down (e.g., the peak around 2.5 Gyr), or systematic bias (e.g., the peak around 1 Gyr).
(or stars with P rot and T eff measurements above those of the members of the Praesepe), Ro <1.866, and 3,000 K < T eff < 7,000 K. Inferring age for stars outside of this parameter space can lead to incorrect ages as the model is fully empirical, and stars with Ro >1.866 experienced weakened magnetic braking and stopped spinning down.However, Figure 11 suggests the model still has strong predicting power for stars with M G > 3.5 mag.
• A systematic exist at ∼ 1 Gyr for stars < 5,000 K.The cluster sample suggests the isochrones for stars between 0.67 Gyr and 1 Gyr (or even to 2.5 Gyr for low-mass stars due to stalling Curtis et al. 2020) in P rot -T eff space have significant overlaps (see Figure 4), as a result, stars with a range of ages but similar P rot and T eff measurements will have similar age inference of around 1 Gyr.
• Inferring ages ∼ 2.5 Gyr for partially convective stars could be inaccurate.Partially convective stars ∼ 2.5 Gyr start experiencing a stalled in their surface spin-down, most likely due to coreenvelope decoupling (Curtis et al. 2020).As a result, stars with a range of ages can overlap in P rot -T eff space and create prediction biases at ∼ 2.5 Gyr.Claytor et al. 2020).This means, all empirical gyrochornology relations available in the literature, calibrated on clusters or asteroseismic data, suffers from this bias.Figure 12 shows the absolute differences between the LEGACY ages and gyrochronology ages (∆Age) as a function of metallicity.An obvious trend is observed that gyrochronology ages for lower metallicity stars deviate more from the asteroseismic ages.

CONCLUSION
Gyrochronology is one of the few promising methods to age-date single main-sequence field stars.However, gyrochronology relies strongly on empirical calibrations as the theories behind magnetic braking are complex and still unclear.The lack of a relatively complete sample of consistent and reliable ages for old, low-mass mainsequence stars with period measurements has prevented the use of gyrochronology for relatively old low-mass stars beyond ∼ 4 Gyr (the age of the oldest cluster with significant period measurements Dungee et al. 2022).
By combining period measurements from Kepler and ZTF, using the gyro-kinematic age-dating method, we constructed a large sample of reliable kinematic ages expanding the P rot -T eff space that is most suitable for gyrochronology (4 days < P rot < 200 days; 3,000 K < T eff < 7,000 K).By using a Gaussian Process model, we constructed the first calibrated gyrochronology relation that extends to the fully convective limit and is suitable for stars with ages between 0.67 Gyr and 14 Gyr.Cross-validation tests and predicting ages for dwarf stars with asteroseismic signals suggest our model can provide reliable ages with uncertainties on the order of 1 Gyr, similar to that of isochrone ages (e.g.Berger et al. 2023, Figure 9).In this paper, we provide ages for ∼ 100,000 stars with period measurements from Kepler and ZTF, of which 763 are exoplanet host stars with a total of 1,060 planets.
Systematic exist at stellar age ∼ 1 (for T eff < 5,000 K) and 2.5 Gyr (for partially convective stars) due to the fact that stars with a range of ages overlap in T eff -P rot space, most likely due to stalling caused by coreenvelope decoupling.This causes the model to infer similar ages for stars of a range of ages.Adding other age indicators such as stellar activity in the future could potentially break the degeneracy in T eff -P rot space for stars of certain range of ages.Obvious metallicity bias exists for this model (see Figure 12 left plot; deviation of ∼2 Gyr from the asteroseismic ages as the metallicity of the star reaches -0.5 dex), as a result, future work should incorporate metallicity measurements.
6. ACKNOWLEDGMENTS Y.L would like to thank Joel Ong for suggesting the title.R.A. acknowledges support by NASA under award #80NSSC21K0636 and NSF under award #2108251.This work has made use of data from the European Space Agency (ESA) mission Gaia,2 processed by the Gaia Data Processing and Analysis Consortium (DPAC).3Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.This research also made use of public auxiliary data provided by ESA/Gaia/DPAC/CU5 and prepared by Carine Babusiaux.
This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.
This research was done using services provided by the OSG Consortium (Pordes et al. 2007;Sfiligoi et al. 2009), which is supported by the National Science Foundation awards #2030508 and #1836650.
This research has also made use of NASA's Astrophysics Data System, and the VizieR (Ochsenbein et al. 2000) and SIMBAD (Wenger et al. 2000) databases, operated at CDS, Strasbourg, France.

Figure 2 .
Figure 2. Top: MG-T eff for the full sample of ∼ 85,000 dwarf stars with period measurements from ZTF and Kepler(McQuillan et al. 2014;García et al. 2014;Santos et al. 2019;  Lu et al. 2022, this work).The red dashed line shows the shifted 6 th -order polynomial (f6(T eff )) fitted to the entire sample that separates the equal-mass binaries from the rest of the sample.Middle: similar to the top plot but after excluding equal-mass binaries (a total of 68,378 dwarf stars).Bottom: period distribution of the 68,378 dwarf stars.

Figure 3 .
Figure 3. Full cluster sample (small points) and the final 660 cluster stars (large points) used in training.The cluster ages range from 670 Myr to 4 Gyr, and the final training set is selected with 3-sigma clipping to exclude fast-rotating stars that have not converged onto the slow-rotating sequence.
cles; 1,109 data points) and fully convective (crosses; 96 data points) stars colored by the median gyro-kinematic or the cluster ages are shown in Figure5bottom plot.

Figure 5 .
Figure 5. Top: 46,362 stars with corrected gyro-kinematic ages between 1.5 and 14 Gyr (background histogram) and individual cluster stars (circles) with literature ages between 0.67 Gyr and 4 Gyr.Bottom: GP training set for the partially (circles; 1,109 data points) and fully convective (crosses; 96 data points) stars colored by the median gyrokinematic ages in each [T eff , log10(P rot )] ∼ [50 K, log10(1 Days)] bin or the cluster ages.

Figure 7 .
Figure 7. Parameter posterior distributions for the mean function of the Gaussian Process model for the partially convective stars after the MCMC has converged.The parameter descriptions are shown in Table1.

Figure 8 .
Figure 8. Same as Figure 7 but for the fully convective stars.available from Santos et al. (2021).Figure 11 shows the 1-to-1 comparison between the LEGACY asteroseismic ages and the gyrochronology ages from our model colored by T eff (left; Curtis et al. 2020) and [Fe/H] (right; Silva Aguirre et al. 2017).The uncertainties for the asteroseismic ages were calculated by taking the standard deviation of the age predictions from various pipelines from Silva Aguirre et al. (2017).The ages and uncertainties for the gyrochronology ages were calculated by first calculating the ages for each star using 100 realizations of the model where the parameters are taken

Figure 9 .
Figure 9. Cross-validation results for the 20% gyrokinematic age sample (black histogram) and individual clusters (red points).The systematic at 1 Gyr indicates existing bias in predicting stars younger than 1 Gyr old.

Figure 11 .
Figure 11.Testing result for 51 LEGACY stars with asteroseismic ages (not included in our training sample) colored by T eff (left; this work) and [Fe/H] (right; Silva Aguirre et al. 2017).Stars with MG < 4.2 mag (outside of the training sample) are shown in crosses.This result suggests our model can estimate ages for single dwarf field stars with uncertainties just over 1 Gyr.

Figure 12 .
Figure 12.Absolute differences between the LEGACY ages and gyrochronology ages as a function of metallicity (left) and Rossby number (right).The red dotted lines show where the difference is 0. The uncertainties are calculated assuming Gaussian uncertainty (σ 2 = σ 2 LEGACY + σ 2 gyro ).There exists an obvious metallicity trend for stars with [Fe/H] < 0.0 dex, in which the ages can deviate up to ∼ 2 Gyr as metallicity goes down to ∼ -0.5 dex.Age prediction significantly worsens for stars with Rossby number > Rocrit, which is ∼ 1.866 according to Saunders et al. (2023). .

4.
LIMITATIONS & FUTURE WORKSome possible limitations and biases of this model include:• This model should only be applied to stars with M G > 4.2 mag, P rot < 200 Days, ages > 0.67 Gyr

Figure 13 .
Figure 13.Histograms of stars with inferred gyrochronology ages from this work.The black dotted lines show the recent enhancement of star formation rate (SFR) shown in Ruiz-Lara et al. (2020) (5.7, 1.9 and 1.0 Gyr).The peaks in the histograms can correspond to the enhancements of the star formation rate in the Milky Way (e.g., the peak in the tail around 5.7 Gyr), changes in stellar spin-down (e.g., the peak around 2.5 Gyr), or systematic bias (e.g., the peak around 1 Gyr).

Figure 14 .
Figure 14.Left: P rot -T eff diagram of exoplanet host stars colored by their gyrochronology ages inferred from this work.Right: histogram of the gyrochronology ages inferred from this work.

Facilities:
Gaia, Kepler, TESS, PO:1.2m (ZTF), Exoplanet Archive Software: Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), dustmaps (Green 2018), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), Pandas (McKinney et al. 2010), tinygp (Foreman-Mackey 2023), numpyro (Phan et al. 2019) APPENDIX A. VISUALIZING THE GYROCHRONOLOGY MODEL We can visualize the gyrochronology model by looking at the best-fit mean and covariance functions separately.Figure15shows our full model in the training parameter space (left column), the mean function prediction (second column), and the covariance function correction (third column) for partially convective (top row) and fully convective (bottom row) stars.The last column shows the age uncertainty associated with the model parameters.The uncertainty is calculated based on 100 realization of the model with parameter drown from the MCMC posterior distribution.The large fractional uncertainty for partically convective stars around 6,000 K is both caused by the young age and the overlapping isochrones in the cluster training data (see Figure10).

Figure 15 .
Figure 15.Full model in the training parameter space (left column), the mean function prediction (second column), and the covariance function correction (third column) for partially convective (top row) and fully convective (bottom row) stars.The last column shows the age uncertainty associated with the model parameters.

Table 1 .
Initial values for maximizing the log-likelihood function, Gaussian prior width used in the MCMC fits, and and final values for the mean function (Equation1) and GP squared exponential kernel (Equation2) parameters after the MCMC fitting.
eff -log 10 (P rot ) space, with the size of the grids to be [T eff , log 10 (P rot )] = [52 K, log 10 (1.1 Days)].We then selected all (P rot , T eff ) points that had predicted ages within 5% of the desired age.The running median (solid lines) and standard deviation (shaded area) of values that are close to the mean function outside of the range of the training data.Moreover, since obtaining gyro-kinematic ages requires binning stars in similar P rot , T eff , and M G , they are less reliable at the edges

Table 2 .
Catalog description of the gyrochronology ages for ∼ 100,000 stars derived from this work.This table is published in its entirety in a machine-readable format in the online journal.Gyr gyrochronology age upper uncertainty Age err m Gyr gyrochronology age lower uncertainty

Table 3 .
Catalog description of the gyrochronology ages for 384 exoplanet host stars derived from this work.This table is published in its entirety in a machine-readable format in the online journal.