Chemical Evolution in the Milky Way: Rotation-based Ages for APOGEE-Kepler Cool Dwarf Stars

Current theory of stellar spin-down or braking uses charged stellar winds in a magnetic field to carry angular momentum away from the star. The most basic braking laws have the form ${dJ}/{dt}\propto {\omega }^{3}$ (Kawaler 1988) for stars of a given mass, where J is the angular momentum and ω is the angular rotation velocity. In other words, faster-rotating stars spin down more quickly, while more slowly rotating stars spin down more slowly. This means that stars with similar masses born at nearly the same time will asymptotically approach the same spin rate, regardless of their initial angular speeds. Thus, we can map the rotation period of a star at late time robustly to its age, provided that the mapping is properly calibrated. This relationship is valid for stars with M ≲ 1.3 M_⊙ or ${T}_{\mathrm{eff}}\,\lesssim 6200$ K. Above these regimes the surface convection zone vanishes, affecting the magnetic field generation and thereby the spin-down process (this is known as the Kraft break, after Kraft 1967).

Describing the period distributions observed in clusters requires more complicated braking prescriptions than the Kawaler law. The discovery of rapidly rotating stars in the Pleiades (Soderblom et al. 1983; Stauffer et al. 1984) suggested that angular momentum loss must be weaker in some regimes. Observations of X-ray emission, a proxy for magnetic activity, also suggested that stars' magnetic fields saturate above some critical rotation speed ω_crit (e.g., Vilhu & Rucinski 1983; Wright et al. 2011). This saturation would allow for weaker magnetic braking for the fastest-rotating stars. Using magnetic saturation to modify the Kawaler law, Krishnamurthi et al. (1997) and Sills et al. (2000) were able to reproduce the period distributions in open clusters. This modified braking law has the form

$\begin{eqnarray*}&&\displaystyle \frac{{dJ}}{{dt}}=\left\{\begin{array}{l}-K{\omega }_{\mathrm{crit}}^{2}\omega {\left(\displaystyle \frac{R}{{R}_{\odot }}\right)}^{\displaystyle \frac{1}{2}}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{-\tfrac{1}{2}},\omega \gt {\omega }_{\mathrm{crit}}\\ -K{\omega }^{3}{\left(\displaystyle \frac{R}{{R}_{\odot }}\right)}^{\displaystyle \frac{1}{2}}{\left(\displaystyle \frac{M}{{M}_{\odot }}\right)}^{-\tfrac{1}{2}},\,\omega \leqslant {\omega }_{\mathrm{crit}}\end{array}\right.,\end{eqnarray*}$

where

$\begin{eqnarray*}&&{\omega }_{\mathrm{crit}}={\omega }_{\mathrm{crit},\odot }\displaystyle \frac{{\tau }_{\mathrm{cz}}}{{\tau }_{\mathrm{cz},\odot }},\end{eqnarray*}$

τ_cz is the convective overturn timescale, and K is a constant calibrated so that the law reproduces the Sun's rotation at the solar age. R and M are the stellar radius and mass, respectively, and the subscript ⊙ denotes solar values.

This form of the braking law ties the magnetic braking to the properties of the convection zone (e.g., Cranmer & Saar 2011; Wright et al. 2011), introducing an additional mass dependence. However, to second order, metallicity also affects the properties of the convection zone: metallicity increases opacity, which increases the depth of the surface convection zone and convective overturn timescale. Furthermore, metal content impacts the main-sequence lifetimes of stars at fixed mass, meaning that stars of different metallicities have different amounts of time over which to lose angular momentum on the main sequence. Two stars with the same mass and age, but different compositions, will therefore have different rotation periods.

Empirical braking laws generally neglect metallicity because there are too few calibrators to span that dimension sufficiently. Most open clusters used to calibrate spin-down relations have solar metallicity, so using these relations to predict angular momentum loss for stars with nonsolar abundances requires extrapolation. However, with a full evolutionary model, one can account for the impact of metallicity on stellar structure and lifetime naturally and then apply a braking law to the model. This will result in a metallicity-dependent spin-down relation, and while there is still a lack of calibrators beyond solar metallicity, the extrapolation to other metallicities is physically motivated.

Using the braking law in Matt & Pudritz (2008) and Matt et al. (2012), van Saders & Pinsonneault (2013) produced such a model-based spin-down relation. They presented a torque of the form

$\begin{eqnarray}&&\displaystyle \frac{{dJ}}{{dt}}=\left\{\begin{array}{l}{f}_{K}{K}_{M}\omega {\left(\displaystyle \frac{{\omega }_{\mathrm{crit}}}{{\omega }_{\odot }}\right)}^{2},\,{\omega }_{\mathrm{crit}}\leqslant \omega \displaystyle \frac{{\tau }_{\mathrm{cz}}}{{\tau }_{\mathrm{cz},\odot }}\\ {f}_{K}{K}_{M}\omega {\left(\displaystyle \frac{\omega {\tau }_{\mathrm{cz}}}{{\omega }_{\odot }{\tau }_{\mathrm{cz},\odot }}\right)}^{2},\,{\omega }_{\mathrm{crit}}\gt \omega \displaystyle \frac{{\tau }_{\mathrm{cz}}}{{\tau }_{\mathrm{cz},\odot }},\end{array}\right.\end{eqnarray} \tag{ 1 }$

with

$\begin{eqnarray*}&&\displaystyle \frac{{K}_{M}}{{K}_{M,\odot }}={R}^{3.1}{M}^{-0.22}{L}^{0.56}{P}_{\mathrm{phot}}^{0.44}.\end{eqnarray*}$

Here f_K is the solar calibration constant and R, M, L, and P_phot are the radius, mass, luminosity, and photospheric pressure, respectively, in solar units. To summarize, van Saders & Pinsonneault (2013) took the Matt & Pudritz (2008) formulation, which uses the magnetic field strength and mass-loss rate, and parameterized them in terms of photospheric pressure, rotation velocity, and convective overturn timescale based on empirical relations (Pizzolato et al. 2003; Wood et al. 2005; Hartman et al. 2009). Notice that, for a star of known mass, composition, and age, the other parameters in K_M can be determined from a stellar model, allowing for extrapolation of the spin-down law to stars with nonsolar composition. While this partially accounts for metallicity in the braking law, the exponents of the braking law parameters are anchored in empirical calibrations that are largely based on stars with Sun-like composition. We do not treat the composition dependence of these exponents, nor any metallicity dependence of the initial conditions, both of which may be important. Instead, we focus only on the underlying stellar evolution models. The remaining values of f_K and ω_crit, as well as initial conditions (i.e., initial disk period and disk-locking timescale), are typically calibrated using solar and cluster data (van Saders & Pinsonneault 2013).

The braking law in Equation (1) has been tested in several regimes against other spin-down models. For example, van Saders & Pinsonneault (2013) showed that this law reproduces the rapid rotation of stars near the low-mass side of the Kraft break better than other Kawaler law variants. Similarly, Tayar & Pinsonneault (2018) demonstrated the strength of Equation (1) in reproducing the rotation periods of intermediate-mass core helium burning stars. For stars more like the Sun, most braking laws perform equally well since they are all anchored to the same set of Sun-like calibration stars. Because of its track record in reproducing observed periods of stars of various masses and evolutionary stages, we adopt the braking law in Equation (1) for our study.

For our analysis we make use of the grid of stellar rotation models created using the Yale Rotating stellar Evolution Code (YREC; Demarque et al. 2008) from van Saders & Pinsonneault (2013) and van Saders et al. (2016). These models do not intrinsically include rotation. Instead, rotational evolution is determined by integrating Equation (1), while the models are used to provide the time and metallicity dependence of the stellar parameters appearing in the braking law. Thus, our rotational evolution does not include the effects of rotational mixing or evolutionary feedback from starspots, both of which may be important for the most active and rapidly rotating stars. We used the same input physics as van Saders et al. (2016) but recomputed the model grid using the Castelli & Kurucz (2004) atmosphere tables, which include nonsolar abundances of α-elements (O, Ne, Mg, Si, S, Ar, Ca, and Ti).

Our model grid spans the range of mass, metallicity, and α-abundance detailed in Table 1, while each evolution track has been run to either the helium flash or 30 Gyr, whichever occurs first. The models utilize the nuclear reaction rates of Adelberger et al. (2011) and chemical evolution of the form $Y={Y}_{\mathrm{primordial}}\,+\tfrac{{\rm{\Delta }}Y}{{\rm{\Delta }}Z}Z$, with $\tfrac{{\rm{\Delta }}Y}{{\rm{\Delta }}Z}=1.4$. To model convection, YREC employs mixing-length theory (MLT; Vitense 1953; Cox & Giuli 1968) with no convective overshoot. We use opacities from the Opacity Project (Mendoza et al. 2007), supplemented with low-temperature opacities of Ferguson et al. (2005), for a Grevesse & Sauval (1998) solar mixture of elements, and the OPAL equation of state (Rogers & Nayfonov 2002). A solar calibration sets the values of the mixing-length parameter α = 1.86, solar hydrogen mass fraction X_⊙ = 0.7089, and solar metallicity Z_⊙ = 0.0183. The rotational evolution uses solid-body rotation and the angular momentum loss of Equation (1), with f_K = 6.575, ω_crit = 3.394 × 10⁻⁵ s⁻¹, disk-locking timescale τ_disk = 0. 281 Myr, and initial disk period P_disk = 8.134 days (the "fast-launch" conditions of van Saders & Pinsonneault 2013). We calculate the convective overturn timescale using the local MLT convective velocity at one pressure scale height above the convective boundary. We note that the effective temperatures for our sample (see Section 2.2) are generally in a regime where the weakened braking observed by van Saders et al. (2016) should not be important; most of our stars are not massive enough to have evolved to half of their main-sequence lifetimes within the age of the Galaxy. Our tracks have been condensed in the time domain to a more practical set of equivalent evolutionary phases (EEPs) according to the method of Dotter (2016).

Table 1.  YREC Model Grid Nodes

	Min	Max	Step
Mass (${M}_{\odot }$)	0.3	2.0	0.01
[M/H] (dex)	−1.0	0.5	0.5
[α/M] (dex)	0.0	0.4	0.4

Download table as:  ASCIITypeset image

We study the sample of Kepler Cool Dwarfs from APOGEE-1 (Majewski et al. 2017) as described in the target selection of Zasowski et al. (2013). APOGEE-1, one of the programs in phase three of the Sloan Digital Sky Survey (SDSS-III; Eisenstein et al. 2011), spectroscopically surveyed the entire Milky Way with resolution R ∼ 22,500 in the near-infrared. While APOGEE was designed primarily to constrain dynamical and evolution models of our Galaxy, in collecting data across all Galactic regions, the survey also targeted the Kepler field. The APOGEE-1 data set therefore contains stars for which we also have rotation periods from Kepler. The stars in our sample were specifically targeted as an ancillary program to the core APOGEE survey. They were selected using Kepler Input Catalog (KIC; Brown et al. 2011) parameters to have ${T}_{\mathrm{eff}}\leqslant 5500$ K, $\mathrm{log}g\geqslant 4.0$ dex, and 7 < H magnitude < 11. With this intentionally simple selection function, the target sample contained 1241 stars. However, three plates were not observed as planned, so only 930 were observed.

The three panels of Figure 1 show where the stars in our sample are located on the H-R diagram, rotation–temperature space, and composition space. Following the method of García et al. (2014a), Ceillier et al. (2016, 2017), and Santos et al. (2019), we derived rotation periods from Kepler light curves obtained with KADACS (Kepler Asteroseismic Data Analysis and Calibration Software; García et al. 2011) using three high-pass filters at 20, 55, and 80 days and interpolating gaps in the data using the inpainting techniques of García et al. (2014b) and Pires et al. (2015). The rotational analysis employs wavelet decomposition based on Torrence & Compo (1998) and as implemented in the A2Z pipeline (Mathur et al. 2010), autocorrelation function (ACF, e.g., McQuillan et al. 2014), and the product of the two (composite spectrum, CS; e.g., Ceillier et al. 2017), obtaining three period estimates for each target, one for each KADACS filter. Aigrain et al. (2015) used injection/recovery methods to test the accuracy of this and other period-detection methods; they found that between 88% (noisy simulated stars) and 92% (noise-free stars) of the periods recovered using this technique were accurate to within 10%. Targets with reliable rotation period estimates are automatically selected if the different estimates agree and the heights of the ACF and CS peaks are larger than a given threshold (for details see Santos et al. 2019). For the remainder of the sample, we proceed with visual inspection. The rotation estimate that we provide and use for the subsequent analysis is that retrieved from the wavelet decomposition. Of the 930 stars in our sample, we recovered rotation period estimates for 532 stars, 65% being automatically selected. The remaining stars were discarded. We detected rotation periods for more stars using this method than using autocorrelation alone: based on the detection fractions in McQuillan et al. (2014), we would expect ∼450 detections. Our period estimates have a median uncertainty of ∼2 days and a median relative uncertainty of ∼8%.

Zoom In Zoom Out Reset image size

Periods were not reliably detected in the remaining stars for a number of reasons, likely due to noise, unfavorable orientations of the rotation axes, long rotation periods (many tens to hundreds of days), and/or small amplitudes of spot modulation. The amplitude of the spot modulation depends on the stellar inclination angle and spot latitudes, but it also depends on spot area coverage, which is expected to depend on stellar age. Thus, both period and amplitude of the spot modulation are markers of age, so we suspect that we preferentially detect rotation for young stars. This is evident in the bottom panel of Figure 1, where most of the metal-poor, α-rich stars (which tend to be old) are preferentially undetected in rotation.

We adopted temperatures, metallicities ([M/H]), and α-element abundances ([α/M]; α includes O, Mg, Si, S, Ca, and Ti but is dominated by Mg, O, and Si in the dwarfs) from APOGEE Data Release 14 (Holtzman et al. 2018). These parameters were derived from spectra by the APOGEE Stellar Parameters and Chemical Abundances Pipeline (ASPCAP; García Pérez et al. 2016). ASPCAP performs a temperature-dependent calibration to the abundance measurements under the assumption that stars in a given open cluster should have homogeneous abundances. We used these calibrated values for 461 of our stars with detected rotation periods; the remaining 71 stars had no calibrated spectroscopic parameters available, usually due to the proximity of the stars' temperatures or metallicities to ASPCAP atmosphere grid edges. For these stars we used the uncalibrated values. The temperatures have a median uncertainty of 100 K, while the metallicities and α-abundances have median formal uncertainties of 0.06 and 0.02 dex, respectively. The spectroscopic temperatures are, on average, about 20 K cooler than their photometric counterparts reported by Berger et al. (2018), with an rms scatter of 147 K. These offsets are consistent with those reported by Holtzman et al. (2018) for dwarfs in APOGEE DR14. Serenelli et al. (2017) noticed a dispersion in ASPCAP metallicity when compared to values derived from optical spectroscopy; following their practice, we added a 0.1 dex uncertainty in quadrature to the ASPCAP formal uncertainties in both [M/H] and [α/M]. For the uncalibrated parameters where uncertainties were not reported, we adopted 100 K uncertainty in temperature and 0.1 dex uncertainty in both metallicity and α-abundance.

As an extra parameter to constrain the evolutionary states of stars in our sample, we adopted luminosities derived by Berger et al. (2018) from Gaia parallaxes. Our sample contains 31 stars that Berger et al. (2018) are missing, and for these we used ASPCAP-derived surface gravities to constrain evolutionary phase.

Single red giants are predicted to have very long rotation periods, and detected rotation signals in them are typically related to binary mass transfer (Tayar et al. 2015). We therefore focus on main-sequence period–age relations, discarding stars with luminosities above 20 L_⊙ from our analysis. We retain subgiants to gauge the usefulness of our braking law for non-main-sequence stars. Of the stars with no luminosity information, we removed those with $\mathrm{log}g\lt 3$. Using the luminosity and surface gravity cuts, 15 stars were removed from our sample. Finally, 40 of the remaining stars were common to the giant candidate catalog of R. A. García et al. (2019, in preparation). These stars' light curves and periodograms were visually inspected to ensure that there was no contaminating power excess at the frequencies from which we derived rotation periods. Upon visual inspection, another 15 stars were cut, leaving 502 for further analysis.

Because our original sample was selected using KIC surface gravities, the giant contamination represents cases in which the original KIC was ineffective at separating evolved and unevolved stars. The fraction of contaminating giants and subgiants in our original sample of 930 stars is 25%. This increases to 45% when we consider stars redder than K_p − J = 2, in contrast to 74% found in Mann et al. (2012) for stars with KIC-determined $\mathrm{log}g\gt 4$, Kepler magnitude K_p < 14, and K_p − J > 2. The difference is likely due to the magnitude cut in our sample (7 < H magnitude <11): our faint-end limiting magnitude results in a smaller observed volume and less giant contamination than Mann et al. (2012) observed.

To obtain ages, we used MCMC sampling to map the stellar observables of metallicity, temperature, luminosity, and period onto the fundamental parameters of mass, bulk composition, and age. We ran chains for 502 stars using the Python emcee package developed by Foreman-Mackey et al. (2013, 2019). emcee samples a posterior probability distribution using an ensemble of walkers, each running its own Markov chain.

For priors, we used uniform constraints on mass, metallicity, and α-abundance only to ensure that the sampler remained within our model grid edges. For the log-likelihood, we used a χ² in the form of

$\begin{eqnarray*}\begin{array}{rcl}\mathrm{ln}{P}_{\mathrm{like}} & = & -\displaystyle \frac{1}{2}({\chi }^{2}+{\chi }_{\mathrm{evo}}^{2})\\ {\chi }^{2} & = & \displaystyle \sum _{i}^{N}{\left(\displaystyle \frac{{x}_{\mathrm{obs},i}-{x}_{\mathrm{mod},i}}{{\sigma }_{{x}_{\mathrm{obs},i}}}\right)}^{2},\end{array}\end{eqnarray*}$

where the ${x}_{\mathrm{obs},i}$ are the Kepler rotation period and ASPCAP-derived ${T}_{\mathrm{eff}}$, [M/H], and [α/M], and the ${x}_{\mathrm{mod},i}$ are the sampled YREC model parameters. The second term ${\chi }_{\mathrm{evo}}^{2}$ was used as a relatively weak constraint on the evolutionary state, and only when the sampled model luminosity was different from the observed luminosity by more than a factor of two. It took the form

$\begin{eqnarray*}{\chi }_{\mathrm{evo}}^{2}=\left\{\begin{array}{ll}{\left(\displaystyle \frac{\tfrac{1}{2}{L}_{\mathrm{obs}}-{L}_{\mathrm{mod}}}{{\sigma }_{{L}_{\mathrm{obs}}}}\right)}^{2}, & {L}_{\mathrm{mod}}\lt \displaystyle \frac{1}{2}{L}_{\mathrm{obs}}\\ 0, & \displaystyle \frac{1}{2}{L}_{\mathrm{obs}}\leqslant {L}_{\mathrm{mod}}\leqslant 2{L}_{\mathrm{obs}}\\ {\left(\displaystyle \frac{2{L}_{\mathrm{obs}}-{L}_{\mathrm{mod}}}{{\sigma }_{{L}_{\mathrm{obs}}}}\right)}^{2}, & {L}_{\mathrm{mod}}\gt 2{L}_{\mathrm{obs}}.\end{array}\right.\end{eqnarray*}$

This form encouraged the walkers to remain within a factor of two of the observed luminosity without allowing luminosity the same constraining power as rotation period and temperature. Thus, if the period and luminosity provided contradictory constraints (e.g., in the case of an unresolved wide binary system, the observed luminosity would be greater than that of a single star, but the observed period may be uncontaminated), the period information would dominate the likelihood value.

An approximation was necessary to make the composition data compatible with the model grid we used. Namely, about half of our stars had α-abundances that were less than solar (with the lowest being [α/M] = −0.44), but the Castelli & Kurucz (2004) atmosphere tables are available only for [α/M] = 0 or +0.4. Rather than extrapolate to subsolar values, we modeled any star with an α-abundance less than solar using a solar value. Since α-enhancement has an effect on the angular momentum loss, this will have induced a slight bias in our determined ages for these stars. We discuss this bias in Section 4.1.

In sampling effective temperature, rotation period, metallicity, and α-abundance, we obtained posterior distributions of ages for our sample of 483 cool Kepler dwarfs. We report for the first time gyrochronological ages inferred using full evolutionary models. The ages and other relevant stellar parameters for our sample are listed in Table 3. We consider the ages in this table to be reliable given our model assumptions. The distribution of ages in our reliable-period sample and the distribution of relative uncertainties are shown in Figure 3.

Zoom In Zoom Out Reset image size

Table 3.  Input Data and Derived Rotation-based Ages for Stars in the APOGEE-Kepler Cool Dwarfs Sample

KIC ID	${T}_{\mathrm{eff}}$	σT	[M/H]	σ[M/H]	[α/M]	σ[α/M]	Prot	σP	L	${\sigma }_{L}^{-}$	${\sigma }_{L}^{+}$	Age	${\sigma }_{\mathrm{age}}^{-}$	${\sigma }_{\mathrm{age}}^{+}$
	(K)	(K)	(dex)	(dex)	(dex)	(dex)	(days)	(days)	(L⊙)	(L⊙)	(L⊙)	(Gyr)	(Gyr)	(Gyr)
1432745	4601	99	0.171	0.118	−0.024	0.102	21.81	1.70	0.168	0.003	0.003	2.88	0.35	0.38
1724975	5259	108	−0.011	0.117	−0.042	0.101	10.61	0.75	0.854	0.018	0.019	1.03	0.10	0.11
1996721	5020	101	−0.005	0.115	0.012	0.101	31.93	3.42	0.295	0.007	0.007	6.04	0.97	1.12
2018047	4049	103	−0.375	0.116	0.103	0.103	43.02	2.45	0.075	0.001	0.001	8.79	0.89	1.02
2156061	4284	96	0.146	0.117	−0.047	0.103	29.38	3.06	0.450	0.116	0.192	4.45	0.72	0.81

Note. We list the Kepler IDs, APOGEE effective temperatures, bulk metallicities, and α-enhancements with formal uncertainties, rotation periods with uncertainties estimated from Kepler light curves, luminosities with uncertainties inferred from Gaia parallaxes by Berger et al. (2018), and the median of the age posteriors with 68% credibility limits.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Simonian et al. (2019) found that more than half of stars with observed periods between 1.5 and 7 days are actually synchronized binaries, noting that 10 days is the boundary longer than which no synchronization should take place (but see also Fleming et al. 2019). For this reason we plot stars with periods shorter than 10 days separately from those with longer periods.

Unlike the case of solar analogs, we clearly predict a large number of old ages for our targets, with no obvious break in the distribution. The age distribution peaks around 1–2 Gyr, which is younger than that found by Haywood et al. (2013) and Buder et al. (2019), but consistent with the result of Silva Aguirre et al. (2018), who also studied stars in the Kepler field. Even if the fast rotators are not binned separately, the peak remains unchanged. The median relative age error is 14%, which is better than the 20%–30% range usually seen for isochrone ages of similar stellar populations (e.g., Buder et al. 2019, who studied over 7000 dwarf and main-sequence turnoff stars in the Galactic Archaeology with HERMES survey, GALAH). For a discussion of systematic bias in our age estimates, see Section 4.1.

We do not expect our sample to be representative of the distribution underlying the Galactic stellar population for several reasons. First, our selection criteria remove stars that have KIC surface gravity $\mathrm{log}g\lt 4$ dex. Since these stars are evolved and therefore preferentially older, this cut biases our sample toward younger stars. Moreover, requiring that our stars' rotation be detectable biases the sample toward more active stars, which also tend to be younger. Furthermore, although the original sample was selected with simple cuts in magnitude, temperature, and surface gravity in the KIC, the KIC itself is somewhat biased in its determination of stellar parameters.

Our sample of stellar ages is uniquely biased in comparison to other such Galactic samples: our ability to infer ages is tied directly to our ability to detect rotational modulation. This modulation is stronger in younger, more rapidly rotating stars. At fixed rotation period, it is stronger in stars with cooler temperatures and deeper convective envelopes (see McQuillan et al. 2014; van Saders et al. 2019). The α-enhancement also tends to deepen convective envelopes at fixed metallicity, although it is a second-order effect in comparison to the bulk metallicity (Karoff et al. 2018). Taken together, it means that we expect our sample to be biased against old, metal-poor objects, which is clearly seen in the bottom panel of Figure 1. This is in contrast to asteroseismology, which is not particularly biased against old stars, and with isochrone methods, which tend to struggle with young and intermediate-age stars on the main sequence in this temperature range.

Van Saders et al. (2019) suggested that stars with Rossby number $\mathrm{Ro}\gt 2$ might become undetectable in spot modulation based on the solar-temperature stars in the McQuillan et al. (2014) field sample, creating a sharp "detection edge" at long periods. For the metallicities and α-enhancements of stars undetected in spot modulation in our sample, we would have expected all stars cooler than 4500 K to fall at Ro < 2 for t < 12 Gyr, and thus to fall below this detection edge.

For those stars with measured rotation periods, we calculated the Rossby number (Ro = P_rot/ τ_cz) using the observed period and the model convective overturn timescale from our MCMC stellar parameter estimation. To estimate a distribution of Rossby numbers for the stars not detected in rotation, we used the following recipe:

• 1.  
  We split the stars into an α-rich and α-poor sequence using the same [M/H]–[α/M] selection as Silva Aguirre et al. (2008, see their Figure 8).
• 2.  
  We randomly sampled an age for each star from the appropriate α-rich or α-poor age distribution of Silva Aguirre et al. (2008, see their Figure 9).
• 3.  
  Using the star's ASPCAP abundances, sampled age, and mass from Mathur et al. (2017), we interpolated a Rossby number for each star from our stellar model grid.

Figure 4 shows a histogram of Rossby numbers for the rotationally undetected stars, compared with the MCMC-estimated Rossby numbers for our "good" sample. While we can say little about the robustness of a Rossby number estimate for an individual undetected star, statistically the set of Rossby numbers should provide some insight into the sample of rotationally undetected stars. Strikingly, the undetected stars are not overwhelmingly predicted to have high Rossby numbers: the undetected stars, like the detected stars, mostly have Ro < 2, save a small but notable tail of undetected stars with predicted Rossby numbers Ro > 2. This suggests that the nondetections are only mildly biased against the oldest stars; most are likely undetected owing to unfavorable spot patterns or inclinations.

Zoom In Zoom Out Reset image size

Another view of detectability is encapsulated in Figure 5, which shows a 2D histogram of the fraction of stars detected in rotation as a function of effective temperature and α-enhancement. We include all stars in our sample with ASPCAP $\mathrm{log}g\gt 4$ so that the fractions are unaffected by contamination by giant stars. Like McQuillan et al. (2014), we find higher detection fractions at cooler temperatures. The decreasing detectability with increasing α-enhancement represents the difficulty of detecting rotation in progressively older and more slowly rotating stars.

Zoom In Zoom Out Reset image size

Here we examine whether we recover the various chemical enrichment trends observed in studies that use ages determined from isochrones or asteroseismology.

3.3.1. Trends in α-abundance with Age

Using our rotation-based ages, we investigated the abundance of α-elements as a function of age to look for trends in Galactic chemical evolution. Our results are shown in Figure 6 and are plotted over results of the similar studies of Haywood et al. (2013), Silva Aguirre et al. (2018), and Buder et al. (2019) to compare age trends. Each comparison study uses a different combination of α-elements; to ensure accurate comparisons, we plot our ages against the appropriate combination of elemental abundances from APOGEE. We consider trends of individual α-elements with age in the Appendix.

Zoom In Zoom Out Reset image size

We recover the general trend of increasing α-abundance with age. The rate of α-enhancement is consistent with the rates observed by Haywood et al. (2013) and Silva Aguirre et al. (2018) for stars younger than ∼8 Gyr, but the sharp upturn in slope at about 8 Gyr seen by Haywood et al. (2013) is absent from our sample. Instead, we observe a spread of α-abundances for old stars more consistent with the results of Silva Aguirre et al. (2018) and Buder et al. (2019). The difference is likely due to the sample selection used by the different studies: Haywood et al. (2013) selected stars in the solar vicinity, while Silva Aguirre et al. (2018) and Buder et al. (2019) study stars in APOGEE and GALAH, which probe beyond the solar neighborhood and detect a wider variety of stars.

Perhaps the most drastic difference between our sample and the comparison samples is that we find very few stars in the old, α-enriched sequence. Figure 7 shows this same result in composition space, this time colored by age. The lower, young, α-poor sequence is well populated, but the upper α-rich sequence is nearly empty. This is likely due to the fact that old, α-enhanced stars are relatively inactive and are therefore difficult to detect in rotation (see again the bottom of Figure 1). Metallicity may play a role too: old stars tend to be relatively metal-poor and thus have thinner convective envelopes. Presumably, this would result in weaker magnetic fields, less starspot activity, and fewer detections at fixed age. In addition, almost all the stars in our sample are located within 0.5 kpc of Earth. The α-rich "thick"-disk stars tend to be farther away, mostly at distances greater than 1 kpc (see Figure 4 of Hayden et al. 2015), well outside of Kepler's peak detectability.

Zoom In Zoom Out Reset image size

3.3.2. Young, α-rich Stars

Our sample includes three young (<6 Gyr), α-rich ([α/M] > 0.13) stars, selected to have periods P > 10 days. These stars are present in numerous other surveys (Chiappini et al. 2015; Martig et al. 2015; Silva Aguirre et al. 2018; Buder et al. 2019) where ages are derived by a variety of other means including isochrone fitting, asteroseismology, and the [C/N]–mass relationship in giants. Chiappini et al. (2015) suggested that this population is formed in the region near corotation with the Galactic bar and is then subject to strong radial migration. They may also represent a blue straggler population (Martig et al. 2015). We caution that if these stars are indeed merger products, our rotation-based ages are unreliable: if mass transfer spins up the stars, they will appear young when ages are inferred through rotation–age relations. Interestingly, these stars are also unusually kinematically hot (see Figure 8).

Zoom In Zoom Out Reset image size

3.3.3. Old, α-poor Stars

3.3.4. Metallicity Dispersion

3.3.5. Metallicity Distribution

We examine the kinematics of our sample as a function of age using the combination of positions, proper motions, parallaxes, and radial velocities from the second Gaia data release (Gaia Collaboration et al. 2018). We draw associate Kepler targets with their Gaia counterparts following the method of Berger et al. (2018), which compared the magnitudes of the Kepler and Gaia source matches, removed stars with uncertain (σ_π/π > 0.2) parallaxes, and implemented the astrometric quality cuts described in detail in Berger et al. (2018), Lindegren et al. (2018), and Arenou et al. (2018). In total, 381 stars in our sample survive the Berger et al. (2018) cuts, have reported Gaia RV values (Gaia Collaboration et al. 2018), and have a "good" label. We convert the Gaia five-parameter astrometric solution and line-of-sight velocities into U (radial), V (azimuthal), and W (vertical) velocities following Johnson & Soderblom (1987). We adopt a peculiar velocity for the Sun of $({U}_{\odot },{V}_{\odot },{W}_{\odot })=(11.0,12.24,7.25)$ km s⁻¹ from Schönrich et al. (2010). We estimate velocity errors by propagating the full Gaia 5D $(\alpha ,\delta ,\varpi ,{\mu }_{\mathrm{ra}},{\mu }_{\mathrm{dec}})$ covariance matrix through the coordinate transformations, assuming an independent Gaussian error distribution for the line-of-sight velocities.

We show the velocities with respect to the local standard of rest in Figure 8 as a function of rotation-based age. The velocity dispersion in the sample is higher in the older stars, consistent with a picture in which stars are dynamically heated over time. Our sample shows much the same behavior as the M-dwarf sample of Newton et al. (2016, their Figure 13), who also observed an increased velocity dispersion in slow rotators, corresponding to old stars.

The uncertainties associated with model-based gyrochronological ages include effects of composition and initial conditions. For composition, changing the bulk metallicity or detailed abundances changes the depth of the stellar convection zone, which in turn affects the rate of spin-down. Meanwhile, it has been demonstrated that a star "forgets" its initial conditions after some sufficient amount of time has passed, usually less than a gigayear for G and K dwarfs, but up to several gigayears for the coolest M dwarfs (e.g., Epstein & Pinsonneault 2014; Gallet & Bouvier 2015). Any lingering dependence on initial conditions remains to be quantified. Using our model grid, we performed three tests to quantify independently the effects of metallicity, α-abundance, and initial conditions on our derived ages.

First, we investigated the effects on age if metallicity is ignored, i.e., if one assumes that a star has solar composition. We attempted to mimic the scenario one might encounter when attempting to use gyrochronology relations on a poorly studied field star whose temperature and period are known, but whose metallicity, detailed abundances, or initial rotation period are not. We employed the following procedure:

• 1.  
  Interpolate ages for two stars from the model grid using the same randomly sampled mass, α-abundance, and rotation period, but with different metallicities: one with [M/H] = 0.5, and one with a solar value (i.e., [M/H] = 0).
• 2.  
  Compare the ages for the two stars, treating the metal-rich star as the "real" star and the star with solar abundance as an "incorrect" model.
• 3.  
  Repeat until all regions of the model grid have been sufficiently sampled.

We repeated this test for a metal-poor case as well, using a "real" star with [M/H] = −0.5.

The second test was like the first, except that the α-enhancement was used as the manipulated variable. In other words, we sampled two stars of the same mass, metallicity, and rotation period, but differing α-abundances of [α/M] = 0 (the solar case) and 0.4 (the enhanced case).

Finally, in our third test we sampled stars of the same mass, metallicity, and α-abundance, but with different initial conditions. The control sample employed the standard "fast-launch" conditions of disk-locking timescale ${\tau }_{\mathrm{disk}}=0.281\,\mathrm{Myr}$ and disk period P_disk = 8.134 days. The test sample, on the other hand, used "slow-launch" conditions of τ_disk = 5. 425 Myr and P_disk = 13.809 days (van Saders & Pinsonneault 2013). These values were calibrated using stars in the Pleiades (125 Myr old) and M37 (550 Myr).

Each test consisted of 50,000 random samples, and the results are shown in Figures 9 and 10. In all plots, the value (t_inferred − t_true)/t_true indicates by what fraction the age is overestimated (or underestimated) under the given false assumption. The panels in a given row show the results for the same sample in two different ways. The left panels of Figure 9 show the age bias as a function of mass, while the right panels show a ${T}_{\mathrm{eff}}$–period diagram with bins colored by the age bias. For the α-abundance and launch condition tests, the bias depended on metallicity, so the left panels of Figure 10 have metallicity on the vertical axis and are colored by age bias.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If one assumes solar metallicity for a metal-rich, 0.75 M_⊙ star (top of Figure 9), the age will be overestimated by ∼15%. The bias is much larger for stars more massive than the Sun, i.e., near the Kraft break. Here a subtle change in metallicity changes whether the star has a convective envelope and spins down over time, so the inferred rotation-based age is drastically affected. These effects are slightly worse, but in the opposite direction, for the metal-poor case: for a metal-poor, 0.75 M_⊙ star, the age will be underestimated by ∼40%. In addition, the Kraft break occurs at lower masses in metal-poor stars, so the bias is also worse than the metal-rich case for stars slightly less massive than the Sun. On the low-mass end, the bias is reversed: for a metal-rich M dwarf, assuming solar metallicity results in underestimating the age by ∼30%. However, this is also the regime in which a spread of initial rotation periods persists. Both of these cases are beyond regimes where gyrochronology has been validated: above the Kraft break, the lack of spin-down makes any kind of rotation dating useless, while M dwarfs suffer from a lack of calibrators. If we consider the somewhat better-studied regime of 0.5–1.0 M_⊙, the bias ranges from −30% (for metal-poor stars) to 15% (for metal-rich stars) when comparing the inferred ages of two stars of identical mass and period, ignoring metallicity. For the observer's case where temperature and period are known (say, for example, they are solar) and metallicity is assumed to be solar, the inferred age for a truly metal-rich ([M/H] = 0.5) star will be overestimated by ∼15%, while the inferred age for a truly metal-poor ([M/H] = −0.5) star will be underestimated by ∼30%. Similarly, two stars with identical ${T}_{\mathrm{eff}}=4500\,{\rm{K}}$ and P_rot = 40 days will have ages different by ∼5% in the metal-rich case and by ∼10% in the metal-poor case. Metallicity is therefore an important consideration in full evolutionary spin-down models; more calibrators at a range of metallicities should be pursued.

The top row of Figure 10 shows the results of the α-abundance test. If one assumes solar α-abundance for an α-enhanced, solar-metallicity star of 0.75 M_⊙, detailed abundances enter the systematic error budget at the 10%–20% level in a model grid. This bias remains fairly constant in the mass ranges where gyrochronology is useful, but the bias is more drastic in stars near the Kraft break, similarly to the metallicity tests. This suggests that even if the metallicity is known perfectly, rotation-based ages may be biased by as much as 20% if detailed abundances are not known to some level. Since the atmosphere tables used in our model grids are defined only for [α/M] = 0 and 0.4, we cannot say without extrapolation how this bias behaves in α-poor stars.

Finally, the bottom row of Figure 10 shows the bias in age if fast-launch models are used but the system is truly slow-launch. For late G and K dwarfs, assuming the incorrect launch condition generally does not result in more than a 10% bias in age because sufficient time has passed to allow them to "forget" their initial rotation conditions. For late K and early M dwarfs, incorrectly assuming initial conditions comes into the systematic error budget at the ∼10%–20% level, with the bias growing worse toward M dwarfs, which take longer to forget their initial conditions. The bias is also larger for stars near the Kraft break; since stars in this regime do not spin down with age, two identical stars with differing launch conditions never have the same period on the main sequence. This results in an undefined age bias, as shown by the blank bins in the bottom left panel of Figure 10.

While we have tested the sensitivity of our assumed braking law to initial conditions and metallicity, we have neglected a number of confounding factors. We assume only one form of the period–age relation. While all period–age relations reproduce the behavior of solar-temperature stars by construction, they tend to make diverging predictions at cooler temperatures (see Figure 1 in Matt et al. 2015). This is largely a consequence of the scarcity of empirical calibrators for cool stars. The asteroseismic calibrator sample extends only to $\sim 0.8\ {M}_{\odot }$ (van Saders et al. 2016), and open cluster sequences are also truncated at roughly the same mass in systems older than ∼1 Gyr (Meibom et al. 2015; Agüeros et al. 2018; Curtis et al. 2019). A detailed assessment of systematic age offsets produced by different literature braking laws is beyond the scope of this article. Here instead we have focused on whether, under the assumption of a particular braking law, we can recover the galactochemical trends inferred using other age metrics.

Furthermore, we have neglected the impact of internal angular momentum transport in our models. To maintain solid-body rotation, the transport in the interior must be fast in comparison to the angular momentum loss from the surface, which is not necessarily true during all phases of evolution. In so-called "two-zone" models (MacGregor & Brenner 1991), the core and envelope are treated as separate zones, with angular momentum transport occurring between them on a characteristic timescale τ_c–e. Two-zone models can display slower rotation than the solid-body case (at fixed age) while the core has had insufficient time to couple (the effective angular momentum reservoir being depleted by magnetized winds is smaller), but faster than solid-body rotation once the core and envelope recouple. Recoupling timescales in Sun-like stars are τ_c−e ≈ 10–55 Myr (Denissenkov et al. 2010; Gallet & Bouvier 2015; Lanzafame & Spada 2015; Somers & Pinsonneault 2016), while for 0.5 M_⊙ stars τ_c−e ≈ 10²–10³ Myr (Gallet & Bouvier 2015; Lanzafame & Spada 2015). During the preparation of this manuscript, Curtis et al. (2019) published rotation periods for K and M dwarfs in the 1 Gyr old cluster NGC 6811 and showed that spin-down appears to stall between the ages of Praesepe (∼700 Myr) and NGC 6811 (1 Gyr), which they attribute to potential core–envelope decoupling. However, van Saders et al. (2019) showed that the longest-period stars from McQuillan et al. (2014) with ${T}_{\mathrm{eff}}\lt 5000\,{\rm{K}}$ have periods fully consistent with a braking law of the solid-body form that we have adopted here for stars roughly the age of the Galactic disk, suggesting that the impact of this stalled braking does not persist for old ages. The impact of core–envelope decoupling on the observed rotation periods in K and M stars is uncertain, due again to a lack of calibrators with known ages. While a detailed treatment of core–envelope decoupling is also beyond the scope of this article, we caution that if it is important, the age scale we present here is likely stretched or compressed with respect to the true underlying age scale, particularly at ages of hundreds of megayears (Sun-like stars) to several gigayears (M dwarfs).