The breakdown of current gyrochronology as evidenced by old coeval stars

Empirically-determined timescales for star formation in giant molecular clouds are around 4-6 Myr (Turner et al., 2022). One can thus consider that stars in clusters older than, say, 10 times such timescales, should be well approximated as being coeval among themselves, as any age dispersion due to their formation epoch would not be larger than 10% of their current age. This of course is the reason why clusters are commonly used as calibrators for gyrochronology relations. For the purposes of the present work, by pairing stars in clusters older than 100 Myr we would be assembling pairs of coeval stars, which is a first step for our test. Moreover, if gyrochronology works, and current relations correctly capture that property, we would expect that pairs of stars belonging to the same cluster will show, by construction, a good agreement with gyrochronology.

As introduced in Section 2, our rotation period test takes pairs of stars as input. In order to apply our technique for a given star cluster, we randomly paired its members, so that these pairs retain the coevality inherited by belonging to the same (sufficiently old) cluster. Four clusters were used to form this sample: Praesepe (age $\sim$ 650 Myr) and NGC6811 (age $\sim$ 1 Gyr) from Godoy-Rivera et al. (2021a); Ruprecht 147 (age $\sim$ 2.5 Gyr) from Gruner & Barnes (2020) and Curtis et al. (2020); and NGC6819 (age $\sim$ 2.5 Gyr) from Meibom et al. (2015). We highlight the wide range of ages spanned by our cluster sample, as this is central to the gyrochronology examinations we perform below.

We compiled rotation period uncertainties from different sources. For Praesepe, as the original source (Rebull et al., 2017) did not provide rotation period errors, we adopted the median fractional error from the other sources ($\sigma_{P_{rot}}/P_{rot}\approx 7\%$) as a representative value. For NGC6811 we compiled the corresponding rotation period uncertainties from the original sources (Santos et al., 2019; Santos et al., 2021a; Curtis et al., 2019b; Meibom et al., 2011). For the Ruprecht 147 stars from Curtis et al. (2020), no period uncertainties were reported either, so we followed the same procedure as for Praesepe. Finally, the NGC6819 stars from Meibom et al. (2015) and the Ruprecht 147 stars from Gruner & Barnes (2020) had period errors already reported, and so we used them.

Note that for Praesepe and NGC6811 we took the stars that were astrometrically classified as probable cluster members by Godoy-Rivera et al. (2021a). In other words, these stars are fully consistent with being members of the clusters, but at the same time, rotational outliers are present in the sample (see figures 9 and 10 in Godoy-Rivera et al. 2021a). Similarly, for the Gruner & Barnes (2020) Ruprecht 147 sample, we discarded stars with low confidence in the rotation period measurement. These choices were made in order to perform meaningful comparisons with field stars, and the implications of the above on the gyrochronology tests are discussed in Section 4.

For photometric colors, we performed a cross-match between the cluster sample (as well as all the other samples described in the rest of this section), and the Gaia DR2 Catalog (Gaia Collaboration et al., 2018), to use the $G_{BP}-G_{RP}$ color as a proxy for mass. For consistency, we defined the primary for all our samples as the brighter star in the pair (i.e., the one with lower apparent $G$ magnitude). We de-reddened the Gaia photometry of the clusters by adopting a global $A_{V}$ extinction for each one, then computing $A_{G}$ and $E(\text{BP}-\text{RP})$ following the approach of Godoy-Rivera et al. (2021a) on a star-by-star basis (where the extinction coefficients depend on both $A_{V}$ extinction and $G_{BP}-G_{RP}$ color). A summary of the global V-band extinction values adopted and their sources, along with mean values for Gaia extinction $A_{G}$ and reddening $E(\text{BP}-\text{RP})$ for each cluster, is presented in Table 1.

Since all gyrochronology relations have a minimum color for which an age can be measured (colors bluer than this step onto the radiate-envelope regime), we only kept stars with²²2Being $(G_{BP}-G_{RP})_{0}=(G_{BP}-G_{RP})-E(\text{BP}-\text{RP})$ the color index corrected by reddening.$(G_{BP}-G_{RP})_{0}>0.562$. Additionally, we ignored fully convective stars – $(G_{BP}-G_{RP})_{0}>2.5$ –, as they tend to break apart from the usual $\sqrt{t}$ dependence, and in this regime the calibrations are considerably worse (see Angus et al. 2019; Pass et al. 2022; Lu et al. 2022). The color-cuts described above were applied directly to our data sample, i.e., on the Gaia DR2 colors. This is also the color-index in which the Angus et al. (2019) was calibrated.

To implement these color-cuts for the Angus et al. (2015) and Spada & Lanzafame (2020) models, however, a few extra steps were required due to the specifics of each relation. The Angus et al. (2015) relation was calibrated with the $(B-V)_{0}$ color, so we transformed the stars’ Gaia DR2 $(G_{BP}-G_{RP})_{0}$ into $(B-V)_{0}$ using a 1 Gyr PARSEC isochrone (Bressan et al., 2012). We used the PARSEC models as they are empirically calibrated, and because the alternative analytical photometric relationships we tested (Jordi et al., 2010; van Leeuwen et al., 2018) yielded inconsistent results in the relevant color range. We also needed to impose a slightly higher minimum color for the pair components, as the Angus et al. (2015) relation has a lower color-boundary of $(B-V)_{0}=0.45$, which corresponds to $(G_{BP}-G_{RP})_{0}\approx 0.582$. On the other hand, the Spada & Lanzafame (2020) relation was calibrated with mass instead of color, so we used the same PARSEC isochrone to convert the stars $(G_{BP}-G_{RP})_{0}$ colors into mass. Moreover, we needed to impose stricter constraints to the colors and rotation periods of the stars. First, this model only reported gyrochrones for masses between $0.45M_{\odot}$ and $1.25M_{\odot}$, so we had to impose the equivalent filter for the Gaia color ($0.592<(G_{BP}-G_{RP})_{0}<2.346$). Second, the Spada & Lanzafame (2020) model only spans the ages between $0.1$ Gyr and $4.57$ Gyr, thus, we only kept stars for which their rotation periods (plus/minus their uncertainties) were between these two limiting gyrochrones.

To ensure that we are using well measured rotation periods, we required stars to have $P_{rot}>3\sigma_{P_{rot}}$. Additionally, in order to filter out potential close binaries and evolved stars from our sample, we imposed two further constraints to our stars. First, for each cluster we downloaded a PARSEC isochrone of its respective age, and eliminated all stars with absolute magnitudes more luminous than $0.38$ mag from the isochrone. This limit is calculated as $\approx 1/2\times 2.5\log_{10}(2)$, and corresponds to half of the magnitude difference between the isochrones of single stars and equal-mass photometric binaries (i.e., systems with twice the flux). We did this to identify and discard stars that may be part of close, unresolved binaries. Gruner & Barnes (2020) used a similar constraint for their Ruprecht 147 sample, thus, we also discard the stars they flagged as photometric binaries. Second, we dropped all stars more luminous than an absolute magnitude of $2.85$, as they probably include evolved stars.

Finally, in order to add robustness to our random pairing method, and increase the statistics of our sample, we repeated the process of randomly pairing the members of any given cluster a total of 20 times. This ensured a Cluster sample 20 times larger than the original, while simultaneously keeping the distributions of the stars’ properties (color, rotation period, etc.) identical, as each member was over-sampled the same number of times.

After our quality cuts, regarding the Angus et al. (2019) relation, we were able to run our test for 125 pairs (2500 pairs after the over-sampling process) for Praesepe, 75 (1500) pairs for NGC6811, 24 (480) pairs for Ruprecht 147, and 7 (140) pairs for NGC6819. Following Curtis et al. (2020), we combined the Ruprecht 147 and NGC6819 samples because of their practically identical ages, thus constructing a sample of 31 (620) pairs composed by $\sim 2.5$ Gyr-old cluster members. For the Angus et al. (2015) relation, the analogous numbers were: 124 (2480 after oversampling) pairs for Praesepe, 74 (1480) pairs for NGC6811, 24 (480) pairs for Ruprecht 147, and 7 (140) pairs for NGC6819. For the Spada & Lanzafame (2020) relation, the stricter constraints reduced the sample sizes, resulting in the final number of: 105 (2100) pairs for Praesepe, 65 (1300) pairs for NGC6811, 20 (400) pairs for Ruprecht 147, and 7 (140) pairs for NGC6819.

In Figure 2 we show the rotational sequences of the four clusters used (after applying our quality cuts). For each one, we also plot a gyrochrone (from Angus et al. 2019) of the corresponding cluster age, which are shown in the label for reference. As Angus et al. (2019) based their relation on Praesepe data, the 650 Myr line fits almost perfectly our Praesepe sample, by construction.

On the other hand, faster-than-expected rotational sequences can be seen for the 1 Gyr cluster NGC6811, and the 2.5 Gyr clusters NGC6819 and Ruprecht 147. More extensive studies of the NGC6811 cluster have revealed the apparent merging of its rotational sequence (for colors redder than $G_{BP}-G_{RP}\gtrsim 1$) with the Praesepe cluster (Curtis et al. 2019b; see also Meibom et al. 2011 and Agüeros et al. 2018). This is discussed further in Section 4.2.

We compiled a sample of Kepler-field wide binaries, given the availability of literature rotation periods for stars in that field (see Section 3.3). We do this by combining existing catalogs from dedicated searches. Our binary sample was composed of four sources: Deacon et al. (2016), Janes (2017), Godoy-Rivera & Chanamé (2018), and El-Badry & Rix (2018).

We de-reddened the Gaia colors following the same procedure as in Section 3.1. In this case, we queried the Bayestar19 map (Green et al., 2019) for extinction values for each star using the GALExtin website (Amôres et al., 2021). We then calculated the absolute magnitude for each star using the distances computed by Bailer-Jones et al. (2018).

The first row of plots in Figure 3 characterizes the Wide Binaries sample. The top-left plot is an HRD and the top-centre plot shows angular separation vs. proper motion difference. For guidance, in the background of these panels we respectively show the McQuillan et al. (2014) stars compared to the Kepler-field HRD, and the El-Badry & Rix (2018) binary pairs to illustrate their astrometric selection. The HRD shows that, by construction, most of the binaries components are MS stars, making them suitable for our gyrochronology tests. The centre plot shows that most pairs of the Wide Binaries sample fall between $\sim$ 5″ and 500″ of separation. These parameters are consistent with gravitationally bound systems, as indicated by the background pairs. On the centre plot, we draw a dashed line (with a $-1/2$ slope corresponding to a Keplerian orbit) to identify pairs that could be considered questionable due to their large proper motion differences (e.g., see Andrews et al., 2017; El-Badry & Rix, 2018; Ramírez et al., 2019). For our combined binary sample, we discard any pairs above this line, as they are probably chance alignments or one of their components is a close, unresolved binary. We do this to compensate for possible heterogeneities in the selection in the aforementioned source binary catalogs.

Stars in the Kepler field are arguably the best sample to use for our test, as the high precision photometry and observing baseline that encompasses many years translate into reliable rotation periods. In order to have a relative homogeneity of rotation period sources, we used primarily rotation periods from Santos et al. (2019), and adopted the values from McQuillan et al. (2014) in cases where stars were not found in the former. The wide binary lists from Janes (2017) and Godoy-Rivera & Chanamé (2018) compiled rotation periods too, so we used their values for stars where the two primary catalogs did not report rotation periods. ³³3Deacon et al. (2016) also compiled rotation periods, but they were exclusively from the McQuillan et al. (2014) catalog, so no extra rotation periods were provided.

While all these periods ultimately come from Kepler observations, the different sources calculate them using a variety of methods (e.g., Lomb-Scargle periodograms, wavelet analysis, auto-correlation function). Thus, we may suspect the rotation periods to have different biases and selection effects. To quantify this, we used the Wide Binaries and NGC6811 stars with more than one period source, and compared the rotation periods coming from the different references. We found that more than 98% of these stars had rotation periods that agree at the $95\%$ level or higher. This excellent agreement validates the reliability of the literature rotation periods, and we thus consider the prioritization of specific period sources to have a negligible effect on our work. ⁴⁴4The one caveat to this relates to the rotation period errors reported by McQuillan et al. (2014), but we further discuss this in Section 4.3.

After applying the same color, $P_{rot}/\sigma_{P_{rot}}$, and HRD filters used for the clusters (see Section 3.1; in this case using a 2 Gyr isochrone), we were able to run our test for a total of 40 Wide Binaries for Angus et al. (2019), 39 for Angus et al. (2015), and 19 for Spada & Lanzafame (2020). We show the former of these samples in the color vs. rotation period diagram on the top-right panel of Figure 3. Note that the smaller size of our Spada & Lanzafame (2020) Wide Binaries sample is due to the stricter constraints we had to apply due to the more restricted regime of applicability of their relation (see Section 3.1).

Constructing a random pairs sample is important, as it gives us a statistical control group composed of stars that are not coeval. These samples are typically used in binary studies (e.g., Lépine & Bongiorno, 2007; Andrews et al., 2019; Espinoza-Rojas et al., 2021). For this purpose, two different random pair samples were constructed.

• •

  Random Pairs [1]: random pairs generated from the components of the wide binaries sample. We took all the components of the wide binaries selected for the rotation period test (for their respective gyrochronology model) and paired them randomly, ensuring that all real pairs from the original sample were destroyed. This guaranteed us a sample of pure random pairs with no real binaries in it. It was also generated so the stars in this sample of random pairs have, by construction, identical distributions of magnitudes, colors, $P_{rot}$, and $\sigma_{P_{rot}}$ as the wide binaries (because they are the same stars). This avoids potential systematic biases that could appear if we compare two samples with different distributions of stellar properties. As with the Clusters sample, we repeated the pairing process 20 times, over-sampling each star an equal amount of times, so the underlying distributions are kept identical.

• •

  Random Pairs [2]: random pairs from McQuillan et al. (2014) stars. The McQuillan et al. (2014) catalog contains rotation periods for over 34,000 stars in the Kepler field. We paired randomly every star in the catalog, constructing a sample of over 17,000 random pairs, with no repeated stars. We also de-reddened their colors through the same procedure as for the Wide Binaries. The result, after applying the previously mentioned color, $\sigma_{P_{rot}}$ and HRD constraints, was a sample of 11,011 random pairs. This approach has the advantage of generating a very large number of pairs, making the sample significantly less sensitive to statistical variations, and avoiding the need for over-sampling a small population of stars. A drawback of this sample, however, is that it does not follow the same distribution of stellar properties as the Wide Binaries and Random Pairs [1]. This is reflected in the color and period distributions shown in the bottom-left and bottom-right panels of Figure 3. Another aspect to consider for this sample, is that by being constructed exclusively from the McQuillan et al. (2014) catalog, it inherits its $\sigma_{P_{rot}}$ values. We find that these tend to be significantly smaller than the $\sigma_{P_{rot}}$ values of the other sources, which has important effects on the x-parameter we have previously defined. We discuss this effect in Section 4.

With this, we have constructed two samples of random pairs with pros and cons that complement each other. The middle and bottom rows of Figure 3 show plots for the Random Pairs [1] (red points) and [2] (green points) samples, respectively. The centre column plots show that the random pairs have much wider separations and larger proper motion differences than the Wide Binaries, both of which are to be expected. Note that, as Wide Binaries and Random Pairs [1] are composed of the same stars (only the pairings are different), the HRD and color vs. rotation period plots for the upper and middle rows are virtually identical, except for the connecting lines in the latter panel.

Since $\Delta P_{rot,gyro}$ is the difference of two variables with relatively similar distributions (i.e., $\Delta P_{rot,gyro}=P_{rot,measured}-P_{rot,expected}$ for the secondary of a given pair), we set the context of our expectations by exploring the distribution that results from subtracting two random variables that follow the same distribution. In the left panel of Figure 4, we plot (in black) a uniform distribution of rotation periods, and (in purple) the distribution of the rotation periods from McQuillan et al. (2014).⁵⁵5In order to keep consistency with the data, the rotation periods used are those of the components of the Random Pairs [2] sample. In the right panel, we plot (in black) the resulting distribution obtained by subtracting two random variables that followed the uniform distribution, and (in purple) the one obtained by subtracting two random variables that followed the McQuillan et al. (2014) distribution.

The black distribution on the right-hand panel has a triangular shape, it is symmetric and centered at zero, and it represents the floor for a $\Delta P_{rot}$ distribution, i.e., the case where all pairs of stars have unrelated components, and thus their rotation period values should be completely independent (the components of some random pairs may have similar ages, and thus small $\Delta P_{rot}$, but this happens by chance, not systematically). It is not a uniform distribution, despite the distributions of the quantities in the difference being uniform themselves, and shows a peak at zero. This means that, even among pairs of unrelated, non-associated stars, a preference for seemingly correlated rotation periods (and therefore ages) is to be expected.

The purple distribution is very similar to the black one, although slightly taller and narrower (somewhat resembling a more Gaussian shape). This demonstrates that we should expect the $\Delta P_{rot,gyro}$ distribution of each of our samples (even the random pairs) to be centered on zero and fairly symmetric. This also justifies our assumption of taking the absolute value in the definition of the x-parameter (Equation 2).

In Figure 5, we show the $\Delta P_{rot,gyro}$ distributions resulting from applying our test with the Angus et al. (2019) relation to the four samples constructed in Section 3.

As it can be seen, the Clusters sample produces the best results among our four samples, with its $\Delta P_{rot,gyro}$ distribution displaying the highest peak and the narrowest width of all distributions, i.e., a better statistical agreement with gyrochronology. Of course, this is the expected outcome, given that the gyrochronology relations have been calibrated using the color vs. rotation period distributions measured for stars in clusters. Indeed, Angus et al. (2019) used the Sun and the Praesepe cluster to constraint their relation.

In Table 2 we present the fraction of pairs that satisfy $x\leq 1,2,3$. As can be inferred from Equation (2), this means the fraction of pairs for which $\Delta P_{rot,gyro}$ falls inside $1,2,3$ times its uncertainty. Here, the fraction of pairs with $x$ smaller than a given value is always higher for the Clusters sample than for the other samples, further demonstrating their better agreement with gyrochronology.

In Figure 6, we present the comparison between the cumulative distributions of the x-parameter for the samples. Thus, the cumulative probability where $x=1,2,3$, corresponds to the values reported in Table 2. The distribution of the Clusters climbs much more rapidly in comparison with the other distributions, reflecting its better concordance (i.e., a narrower peak in Figure 5) with the gyrochronology relations by construction.

In the hypothetical case that the agreement of a given sample of pairs with a gyrochronology relation were “perfect” considering the measurement uncertainties, then we would expect the distribution of $\Delta P_{rot,gyro}$ to be Gaussian, with a width defined by propagating the typical rotation period uncertainties. Thus, in Figure 6 we also show a half-normal distribution (to account for the absolute value factor in the x-parameter) as a black line, and include the corresponding fractions at a given $\sigma$ in Table 2.

Both in Figure 6 and Table 2, the Clusters sample exhibits cumulative distributions below the half-normal distribution (e.g., the probability of $x\leq 1$ in Table 2 for the half-normal distribution is $0.683$, slightly higher than the values for the Clusters). This shows that, while the Clusters present the best agreement with gyrochronology compared to the other samples, it is below of what we could expect under the “best-case scenario” of the Gaussian distribution.

We now examine our results in further detail by comparing the clusters among each other, with respect to the Angus et al. (2019) relation. There is a clear trend in Table 2 for the younger Praesepe cluster to show the best agreement with the relation among the clusters, then the older NGC6811, and finally the oldest $\sim 2.5$ Gyr clusters Ruprecht 147 and NGC6819, with the worst agreement. This suggests a possible degradation of the gyrochronology relations beyond the $\sim 1$ Gyr of age, and we further discuss this in Section 5.

Both Random Pair [1] and [2] samples show very similar $\Delta P_{rot,gyro}$ distributions in Figure 5, although the latter one shows a slightly higher peak. Their overall agreement validates both methods of random pair generation.

However, the differences arise when we compare them in the distributions of the x-parameter in Table 2 and Figure 6. Although the CDF of $x$ for the Random Pairs [1] raises more rapidly with increasing $x$ than the one for Random Pairs [2], this does not necessarily mean a worse agreement with the gyrochronology relations for the latter sample. This is because the rotation period uncertainties $\sigma_{P_{rot}}$ of McQuillan et al. (2014) (the source of the Random Pairs [2] sample) tend to be much lower than the ones from the Wide Binaries sample, and thus for the Random Pairs [1] sample. The median rotation period fractional uncertainty for the Random Pairs [2] is $\sigma_{P_{rot}}/P_{rot}\sim 0.008$, while the median for the Random Pairs [1] sample is $\sim 0.076$, almost and order of magnitude larger. Measurement errors of $\sigma_{P_{rot}}/P_{rot}<0.05$ are likely underestimated and do not account for the uncertainties inherited from phenomena such as differential rotation, spot evolution, and the lifetime of active regions (e.g., Giles et al., 2017; Santos et al., 2021b; Basri et al., 2022). All of these result in much higher values of $x$ for Random Pairs [2] (as the denominator in Equation 2 is smaller), making the comparison of the x-parameter between the Wide Binaries and Random Pairs [1] samples (which have identical $\sigma_{P_{rot}}$ distributions) much more appropriate than with Random Pairs [2].

All things considered, both Random Pairs samples show worse agreement with gyrochronology than the Clusters and Wide Binaries. This implies that the standard gyrochronology relations do have some predictive power in identifying coeval populations of field stars.

Ultimately, the results from the Random Pairs samples show a baseline for the non-zero agreement with gyrochronology that can be expected for unassociated pairs of field stars. This could be considered as analogous to the existence of chemical doppelgangers found in studies of Galactic archaeology (e.g., Ness et al., 2018), although in our case the dimensionality of the agreement is significantly smaller.

The Wide Binaries sample shows a better agreement with the tested relations in Figure 5, in comparison with the samples of random pairs. The distribution of $\Delta P_{rot,gyro}$ is both narrower and taller than the ones from the Random Pairs. However, it also shows a worse agreement with gyrochronology than the combined Clusters sample.

As already mentioned in Section 4.3, the best comparison between wide binaries and non-coeval pairs would be that made against the sample of Random Pairs [1], as their $\sigma_{P_{rot}}$ distributions are, by construction, identical. Here, it is shown that the cumulative distribution of $x$ raises substantially more rapidly for the Wide Binaries than for the Random Pairs [1]. Moreover, in Table 2, we can compare quantitatively the $x$ cumulative distributions for $x=1,2,3$. For instance, for the Angus et al. (2019) relation, we have that the ratio of Wide Binaries to Random Pairs [1] that satisfy $x\leq 1$ is $\approx 2.0$ ($0.375/0.186$). At face value, our results support empirically the initial assumption that wide binary components are, indeed, coeval.

Although the Wide Binaries sample shows better results than the Random Pairs, its agreement with gyrochronology is clearly below the Clusters sample. As we would expect many of our field binaries to be of several Gyr of age, the effect of weakened magnetic braking in older stars (van Saders et al., 2016; Hall et al., 2021; Masuda et al., 2022) is likely to be an important factor that affects our sample. By the same line of reasoning, our method should be sensitive to such departures from standard gyrochronology, and therefore it can be used to examine future calibrations that will account for this weakened magnetic braking. Today, given the lack of suitable clusters, only wide binaries (Chanamé & Ramírez, 2012), stars suitable for seismic analyses (García et al., 2014; Hall et al., 2021), and possibly other seemingly coeval structures (e.g., Andrews et al., 2022; Kounkel et al., 2022) offer hope to better calibrate gyrochronology at several Gyr.

Gyrochronology relations usually depend directly on color (used as a proxy for stellar mass), such that for a given age, the $P_{rot}$ will increase toward redder colors (lower masses). Given this intrinsic dependence, for a given set of gyrochrones that span a range of ages, the gyrochrones will tend to be more clustered at bluer colors (higher masses), and they will diverge and be more separated at redder colors (e.g., see Figure 2). Thus, bluer stars with very different ages may appear with similar rotation periods. This could introduce a bias in our test, as small uncertainties in the rotation period of a blue primary star can imply a large $\Delta P_{rot,gyro}$ for a redder secondary.

To examine the effect of this in our tests, we divided each of our samples into two sub-samples of equal size, by small and large color differences between the primary and secondary. We found that, in general, pairs with similar colors showed smaller $\Delta P_{rot,gyro}$ and x-parameter values than pairs with different colors. This happens with both coeval and non-coeval pairs.

To weigh the dependence of our test on the color of the primary, we flagged pairs with primary colors that were bluer than a certain threshold (we tested several limits between $(G_{BP}-G_{RP})_{0}=0.7$ and 1.0 mag, and the results were very similar), and defined the reference gyrochrone for $\Delta P_{rot,gyro}$ by fitting it to the redder secondary instead. This resulted in $\Delta P_{rot,gyro}$ values that tended to be closer to zero than the original sample, for samples with both coeval and non-coeval star pairs. However, the x-parameter values were almost identical. This is because, by definition (see Section 2), the x-parameter value is basically independent of which star of the pair is treated as the primary. Thus, at least partially, our x-parameter does take this dependence on the primary’s color into account.

Ultimately, our samples are too small to refine our test as to fully take color into account. In order to mitigate the apparent bias in pairs of stars with similar colors, future efforts with larger samples will be needed.

Following our discussion in Section 4.2.1, we now assess the impact that close binaries may have in our results. This arises from the possibility that the rotation periods of the stars we study, and therefore their $\Delta P_{rot,gyro}$ and x-parameter values, may be affected by the presence of unresolved companions at close separations that have survived our HRD cuts. We do this by taking advantage of the recently published Gaia EDR3 and DR3 catalogs (Gaia Collaboration et al., 2021, 2022).

First, we studied the subset of stars with high Renormalized Unit-Weight Error (RUWE), which can be used to identify unresolved companions due to the astrometric noise they imprint in the Gaia solutions (Belokurov et al., 2020). We searched the Gaia DR3 RUWE information for our stars, and found values for virtually 100% of them. We classified those with RUWE$>1.4$ as having high-RUWE (e.g., see El-Badry et al. 2021, Berger et al. 2020a). Across all samples, $\sim$ 20% of the pairs were classified as having at least one component with high-RUWE. When comparing the pairs’ x-parameter distributions, however, we noticed that the ones with high-RUWE components have nearly identical distributions as the full samples. Therefore, we conclude that while these high-RUWE stars may host close companions, these are not having an impact on the rotation periods of the stars we study. This is in agreement with the recent results by Avallone et al. (2022) (e.g., see their figure 13).

Second, we studied the subset of stars that can be classified as radial velocity (RV) variables. We again searched in Gaia DR3 and found that $\sim$ 45% of the pairs in our samples have at least one component with RV data. To identify the RV-variable stars, we followed the prescription of Katz et al. (2022) (see their section 3.7). Out of the stars with RV data, we found the fraction of pairs with at least one RV-variable star to be $\sim$ 5–15% across all samples. When comparing the pairs’ x-parameter distributions with their respective full samples, these pairs with RV-variable stars showed noticeably worse agreements with the gyrochronology models. Further inspections showed that these RV-variable stars tend to be concentrated at bluer colors (presumably inherited from the Gaia RV selection function), and more importantly, at shorter rotation periods. For instance, for the most numerous Random Pairs [2] sample, the RV-variable stars predominantly have periods shorter than 5 days. This is again in agreement with Avallone et al. (2022), where spectroscopic binaries are heavily concentrated at short periods. These stars are likely in the regime where tidal interactions are affecting their rotation rates, hence directly affecting the parameters we study (e.g., see also Simonian et al. 2019, 2020; Angus et al. 2020).

For simplicity, in the analysis that follows we kept all the targets in our full samples. This is motivated by the unnoticeable effect that the high-RUWE stars have in our results, and the small fraction that the RV-variable stars represent compared to the full samples ($\sim$ 3%). In other words, our conclusions are not impacted by this choice. Nonetheless, we highlight that when applying gyrochronology to field MS stars, or performing analyses similar to ours in larger samples, the existing RV-variability data should be checked as to flag potential rotational outliers (see also Daher et al. 2022).

A final, alternative way of quantifying the agreement of our pairs with gyrochronology, is to directly compare the gyrochronology-derived ages of their members, computed from the relations tested. We followed the method of Otani et al. (2022), and defined an age ratio “tolerance” $\Delta$, so we can test if the ages of the components of a given pair are consistent within a certain $\Delta$. For simplicity, instead of defining distributions for the age ratio given by the gyrochronology equations, we compared directly the computed ages. In summary, we calculated the fraction of pairs in a certain sample that satisfied

For instance, when adopting a value of $\Delta=30\%$, we found the following fractions of the samples to be “coeval”: $\sim 73\%$ for Praesepe, $\sim 46\%$ for NGC6811, $\sim 36\%$ for our 2.5 Gyr clusters, $\sim 50\%$ for our Wide Binaries, $\sim 23\%$ for our Random Pairs [1], and $\sim 25\%$ for our Random Pairs [2]. For a tolerance of $\Delta=10\%$, we found “coeval” fractions of: $\sim 38\%$ for Praesepe, $\sim 20\%$ for NGC6811, $\sim 17\%$ for our 2.5 Gyr clusters, $\sim 20\%$ for our Wide Binaries, $\sim 8\%$ for our Random Pairs [1], and $\sim 8\%$ for our Random Pairs [2]. Note that even in the case of the youngest Praesepe cluster, the agreement is not perfect, as discussed in 4.2.1.

These results follow a similar trend to the ones shown by the other results in this section, validating our $\Delta P_{rot,gyro}$ and x-parameter tests, and further suggesting a decline in the reliability of gyrochronology for old stars (Section 5).