Magnetic activity evolution of solar-like stars:
I. $S_{\rm ph}$-Age relation derived from Kepler observationsastropy (Astropy Collaboration et al. 2013; Astropy Collaboration et al. 2018), kiauhoku (Claytor et al. 2020a; Claytor et al. 2020b), NumPy van der Walt et al. 2011, SciPy Jones et al. 2001–, STAREVOL (Siess et al. 2000; Amard et al. 2019), Diamonds (Corsaro & De Ridder 2014), matplotlib (Hunter 2007), pandas (pandas development team 2020; Wes McKinney 2010), KADACS (García et al. 2011; Pires et al. 2015)

From the kiauhoku analysis, ages were obtained for the full sample of stars. However there are some caveats to consider in order to select the most reliable ages. For this analysis, we wanted to remove as much as possible potential binary stars. Hence, we discarded stars for which Gaia Renormalized Unit Weight Error (RUWE) is larger than 1.2, as well as stars that are part of the Gaia DR3 non-single star catalog (Holl et al. 2022), which are either SB1 or eclipsing binaries. Finally we removed the stars where the CPCB1 flag (for close binary candidates based on the light curves, see S19 for more details) was set to 1. As aforementioned, for approximately 3000 stars, no luminosity from Gaia was available and we used surface gravity. This also means that the Renormalized Unit Weight Error (RUWE) was not defined for those stars. We thus discarded them as well. From the modeling results we removed stars with age older than 14 Gyr, stars whose posterior age uncertainty was greater than 60% (which is the optimal cut to remove stars whose age posteriors were poorly converged), and stars whose posterior luminosity mean and median were not mutually consistent (this occurred when the sampler favored two very different model fits equally, producing a bimodal posterior). If any other flag was set to 1 (see Section 3), we also discarded the stars from our analysis. This leads to a sample of 41,440 stars with reliable ages according to our criteria. The ages for these stars are provided in Table 1. In the remainder of this paper, we will refer to these ages as Age_kiauhoku.

In Figure 1, we plot the rotation period from S19 and S21 as a function of effective temperature and color-coded with Age_kiauhoku. We clearly see the boundaries of the gyrochrones. The range 4-5 Gyrs is shown in yellow. When going into more details, we find some fast old rotators with an age $>$ 5 Gyr and $P_{\rm rot}<$ 20 days mostly in the range of effective temperature 5,000-6,000 K. The stars with $P_{\rm rot}$ longer than 50 days with $T_{\rm eff}$ between 5,000 and 6,000 K correspond to subgiant stars that have already slowed down, due to radius expansion.

To check the reliability of our ages, a second set of ages was also computed with the same input parameters described above but with a different stellar evolution code. We used the stellar evolution grid of Amard et al. 2019, which was computed using the STAREVOL code (Siess et al. 2000; Lagarde et al. 2012; Amard et al. 2016a; Amard et al. 2019). The models include a self-consistent treatment of internal angular momentum transport, and angular momentum extraction by magnetised winds following Matt et al. 2015, calibrated to reproduce open clusters rotation period distributions (Amard et al. 2019, see the grids paper for more details). We use a MCMC maximum likelihood tool adapted from the SCePtER software by Valle et al. 2014 to interpolate in the grid of evolution models at intermediate rotation rates. The estimated values of age and mass are obtained by averaging the age and the mass of all the models with likelihood greater than $0.99\times\mathcal{L}_{\textrm{max}}$ similarly to what was done in Amard et al. 2020 to estimate only the stellar masses. The stellar parameters obtained with STAREVOL are given in Appendix A.

This allows us to quantitatively assess the impact of the different physics included in the models. The set of stars in common with reliable ages from both stellar evolution codes consists of 40,220 stars. We can see (Figure 2) that kiauhoku and STAREVOL ages are well correlated in general with a Spearman correlation coefficient of 0.89. However we discuss the differences in more detail below. In the lower panel, the difference between the two sets of ages is centered on 0 with some departure for younger stars. We also found that compared to the reported uncertainties from both stellar modeling approaches, the age differences are on average of -1.22 $\sigma$, where $\sigma$ is the sum of the quadratic uncertainties from each method. This suggests that the true uncertainty in the age is dominated by the systematic uncertainty inherent in choosing a specific model.

While there is generally good agreement between the kiauhoku and STAREVOL ages, there is somewhat worse agreement in older ($>5$ Gyr) and in cooler ($<4500$ K) stars, with up to 5$\sigma$ difference (see Fig. 3). These discrepancies can be explained by the differences in angular momentum transport in the models underlying the two ages estimates. The most likely source of disagreement is the presence of internal differential rotation in the STAREVOL models, which the YREC-based kiauhoku models lack. Differential rotation implies that the torque is stronger at the age of the Sun (Amard et al. 2016b; Somers & Pinsonneault 2016, See), resulting in faster spin-down and a younger age given a rotation period. The effect is more pronounced at old ages, where the models have had more time to diverge in evolution, and in cool stars, which experience stronger differential rotation.

Other differences between the models include the weakening of magnetic braking (van Saders et al. 2016, e.g.,) in the YREC-based models at a given Rossby number, which the STAREVOL models do not include. The weakened braking allows for older ages at shorter periods, which may contribute to the age discrepancy, especially for the older kiauhoku ages. Additionally, both sets of models are calibrated to solar-mass stars, which could lead to different and discrepant behaviors in the coolest regime (Amard & Matt 2020, See e.g.).

We also note that STAREVOL predicts significantly younger stars below 1 Gyr compared to kiauhoku. The differences at these young ages could be explained by the different initial conditions used. Claytor et al. 2020a begin the evolution with a 0.28 Myr disc-locking time at a rotation period of 8.1 days, while Amard et al. 2019 start at 4.5 days after 5 Myr. By 5 Myr, the YREC-based models have spun-up under contraction to periods of 1.25 to 1.75 days. The initial conditions are only expected to affect the age estimates significantly for stars younger than 1 Gyr. Older than this, spin-down ensures that stars “forget” their initial conditions (Claytor et al. 2020a, e.g.,).

Finally, stellar metallicity may also contribute to the age discrepancies, as it has a non-negligible impact on the wind torque (Amard & Matt 2020). Recently Bonanno & Corsaro 2022 showed that metallicity has a significant impact on the rotational evolution of active stars in general.

Several catalogs of ages of the Kepler field stars have been compiled and we now compare them with our derived ages. The most classical method consists of isochrone fittings in a Color-Magnitude Diagram (CMD) or in a Kiel diagram where atmospheric parameters are fitted. This was done in a homogeneous way for the full or almost full Kepler sample in several Kepler star properties catalogs (Brown et al. 2011; Huber et al. 2014; Mathur et al. 2017; Berger et al. 2020a). In the most recent of these catalogs, which we refer to as B20, Gaia DR2 (Gaia Collaboration et al. 2018) luminosity was also used.

Another homogeneous and large-scale catalog was derived using gyro-kinematics relations from age-velocity dispersion (Lu et al. 2021, hereafter L21;) where they combined gyrochronology (using previous rotation periods measurements from Kepler) and vertical velocities from Gaia or LAMOST observations. It has been observed that vertical velocity increases with ages. However, by themselves, vertical velocities can only provide statistical ages. By combining them with gyrochronology, it is then possible to extract individual ages. This work led to a catalog of almost 30,000 stars with gyro-kinematic ages.

We also took advantage of the Gaia DR3 ages obtained with the Final Luminosity Age Mass Estimator (FLAME, Creevey et al. 2022, hereafter C22), that is available for 2945 of our targets. These ages are the most reliable for the Gaia solar-like stars as they are computed based on spectroscopic parameters as well as astrometry and photometry.

Finally, on a smaller scale but with a higher level of precision, we consider the catalogs of seismic ages for dwarfs and subgiant stars. These ages are obtained by finding the best-fit models to atmospheric as well as seismic parameters (global or with individual mode frequencies). A catalog based on the global seismic parameters was done for more than 400 stars (Serenelli et al. 2017) while the detailed analysis with individual frequencies was done for a smaller sample of 99 stars (Silva Aguirre et al. 2015; Silva Aguirre et al. 2017, hereafter SA17). As it was shown that the precision on stellar parameters is improved when using the additional information from the frequencies of the modes (Mathur et al. 2012), we make the comparison with the ages from the smaller sample of stars. SA17 provided ages from 6 different modeling pipelines and we compare our ages with the ones from the BAyesian STellar Algorithm pipeline for consistency between the two seismic samples. However we note that in terms of seismic ages, there can be discrepancies of 1 Gyr up to 4 Gyr between model pipelines.

The comparison with the B20 isochrone ages (for 26,262 stars in common) shows a weak correlation, with a Spearman correlation coefficient of 0.27, as it can be seen in the top left panel of Figure 4. Some disagreement was expected given that isochrones are inaccurate or imprecise for main-sequence stars cooler than the Sun, which make up two-thirds of our sample. The luminosity of such stars does not evolve fast enough to constrain ages precisely, even with Gaia measurements. For cool stars, rotation period is a more precise tracer of age. The comparison of the gyro-kinematic ages with the isochrone ages also shows similar disagreement (see appendix of L21).

Compared to the Gaia DR3 FLAME ages (top right panel of Fig 4), the agreement is reasonably good with a Spearman correlation coefficient of 0.38, so a rather weak to moderate correlation.

In the bottom left panel of Figure 4, we can see that there is a better agreement between the gyro-kinematic ages from L21 and kiauhoku ages. The Spearman correlation coefficient of 0.53 suggests a moderate correlation. While gyro-kinematics ages were derived for individual stars, because they are computed for a given vertical velocity in bins of magnitude, temperature, rotation period, and Rossby number, the ages vary smoothly and a star with an unusual evolution or parameter would get averaged out and might have an over- or under-estimated age. We note that for stars older than $\sim$ 5 Gyr according to our models, the gyro-kinematics ages are smaller, similarly to what was seen in the comparison with STAREVOL ages (Sect. IV.2 and Figure 2). This could reflect the effect of smoothing of the parameter space used by L21. We note that some stars with gyro-kinematics ages are older than 14 Gyr (up to 25 Gyr) that have been removed.

The comparison with the seismic ages in the bottom right panel of Figure 4 shows a strong correlation with a Spearman correlation coefficient of 0.62. There are some stars with seismic ages between 2 and 4 Gyr that are older in our models and have small uncertainties (both from seismology and gyrochronology). These stars all have temperatures hotter than 6000 K and are near the Kraft break (Kraft 1967), where the outer convective envelope becomes thin and stars spin down substantially less over time. In this regime, stars’ rotation rates are more sensitive to their initial speeds, and the ages are more sensitive to the particular braking prescription used by the models.

We fit these stars with different spin-down models, varying the models’ starting speed and whether they experience weakened magnetic braking halfway through the main sequence. Our default kiauhoku/YREC models include weakened braking and have an initial rotation period of 8.1 days, consistent with observations of other seismic samples and young open clusters. Even using models with no weakened braking and slow launch condition of 13.8 days (van Saders & Pinsonneault 2013, consistent with the slow rotators of open clusters;), we recovered the seismic ages for only two of the discrepant stars.

In case the discrepancy was caused by the rotation period, we also tried relaxing the rotation constraint and fitting a YREC isochrone to the temperature, metallicity, and luminosity of the discrepant seismic stars. This reduced the disagreement for about half of the stars, although this was at least partly due to wider age posteriors because the models were less constrained.

Since the discrepancy with the seismic stars was not fully solved by varying the braking prescription or relaxing the rotation constraint, it is possible that the offsets arise because we use different models with different input physics compared to SA17. Tayar et al. 2022 showed that fitting different models to the same input parameters can result in large scatter in inferred ages, as high as 50% near the main sequence. We emphasize that while there may be scatter between different model-dependent ages, our homogeneous modeling procedure produces ages that are internally consistent and highly precise (15% median relative uncertainty).

Studies that are looking for magnetochronology relations usually search for a power-law relation, which translates into a simple linear relation in the log-log space of the magnetic activity proxy as a function of stellar age. Because of the way uncertainties are taken into account in the Bayesian framework, the coefficients found for a given relation do not correspond to the coefficients for the inverse relation. As we want to predict the ages from the magnetic activity proxy, we fitted a power law between the ages determined from the stellar evolution model and $S_{\mathrm{ph}}$, which in the log-log space can be expressed as:

where Age is in Gyr and $S_{\mathrm{ph}}$ in ppm.

We also fitted the inverse relation, i.e. $S_{\mathrm{ph}}$ as a function of age:

We will see that the coefficients can be different between the two ways of fitting the relations.

To have less scatter in our plot and to focus on the most reliable ages, we first selected stars with an effective temperature ranging between 5,700 K and 5,900 K (corresponding to 2,064 stars). We will refer to these stars as solar analogs. Figure 6 shows how the photometric magnetic activity proxy varies with the stellar age for these solar analogs. We clearly see a trend where the magnetic activity decreases as the star evolves. This behavior is expected and has been seen in other proxies of stellar magnetic activity such as Ca HK or magnetic field (Mamajek & Hillenbrand 2008; Vidotto et al. 2014, e.g. ). In Figure 6, we observe a plateau for ages above $\sim$ 5 Gyr, also seen by Masuda 2022 for a small sample of planet-host stars. We can wonder whether the plateau is real or rather a result of detection biases. From that same figure, we clearly see that in our sample of stars with measured rotation periods around 3-4 Gyr, we have $S_{\rm ph}$ values lower than the plateau, meaning that our detection limit could be lower than that observed plateau. In Section VII.4, we discuss more in details the possibility of a detection bias in our sample for low $S_{\mathrm{ph}}$ values.

In addition to the plateau after 5 Gyr, there is another one below 1 Gyr. So we performed a first Bayesian fit between 1 and 5 Gyr for Eq. (2) and found $a_{0}$ = $1.937^{+0.022}_{-0.014}$ and $a_{1}$ = -0.157 $\pm$ 0.009.

As mentioned above, we also performed another Bayesian fit, similar to the previous activity index-age relations, with Eq. (3). In that case we found $a^{\prime}_{0}$ = $8.513^{+0.004}_{-0.003}$ and $a^{\prime}_{1}$ = -2.170 $\pm$ 0.007 as shown by the dashed red line in Figure 6, which means that $S_{\rm ph}\propto({\rm Age})^{-2.17}$. This is quite different from what has been found for other proxies of magnetic activity, such as magnetic field or $\log R^{\prime}_{\rm HK}$, where the exponent varies between -0.5 and -0.6 (Soderblom et al. 1991; Vidotto et al. 2014; Lorenzo-Oliveira et al. 2018, e.g.). However given that $S_{\rm ph}$ measures a different type of magnetic activity feature in photometric data than the chromospheric one, we do not expect to find the same coefficients. We remind the reader that in the previous works, ages were obtained from simple isochrone fitting with different stellar evolution models. This can also impact the age determination and lead to different coefficients.

We note that the slope from Equation (3) does not correspond to the inverse relation of Equation (2), which is due to the different variable for which we optimize the inference process. Indeed in the first case of Equation (3), we use errors on $S_{\mathrm{ph}}$ while in the other case we use errors on age.

We then performed the fit by selecting different samples where we changed the range of the effective temperature for both G and K dwarfs and we found that the slope of the $\ln(S_{\rm ph})-\ln(Age)$ relation depends on $T_{\rm eff}$. When $T_{\rm eff}$ decreases, the absolute value of the slope decreases (see Figure 7), i.e. the relation becomes flatter as seen in Figure 5. This suggests that lower mass stars remain active for longer timescales compared to higher mass stars. The steep $S_{\mathrm{ph}}$-Age relation at young ages seen in Figure 5 is consistent with these results.

The fit is performed by adopting a specific Likelihood function that is formally equivalent to the standard Normal Likelihood given by Eq. (1), except for the treatment of the uncertainties. Following the work by Corsaro et al. 2013, the uncertainty on each data point is computed as a total uncertainty arising from that of each of the observables taken into account, including the one on Age. The uncertainty on a single data point $j$ is thus depending on the model parameters $\{\alpha_{i}\}$ and it is evaluated as:

where $\widetilde{\sigma}$ is the relative uncertainty because we are dealing with logarithmic quantities in the linearized model (e.g. $\widetilde{\sigma}_{\mathrm{Age}}\equiv\sigma_{\mathrm{Age}}/\mbox{Age}$).

We then took advantage of the new artificial intelligence tools that are commonly available and developed by the community, to estimate ages of stars from the same parameters used in the Bayesian fits earlier. This analysis was done for late-F, G, and K dwarfs together. We used Random Forest (RF) algorithms (Breiman 2001), as they have proved to be very useful in estimating stellar physical parameter from multi-parameters laws (Miller et al. 2014; Bugnet et al. 2018; Breton et al. 2021, e.g.). Random Forest regressors are based on the aggregation of a large number of random decision trees that are constructed from a training data set and internally validated to give a prediction based on the predictor for future observations. The RF method not only allows the use of a large number of parameters but also estimates their individual impact on the regression (see Section VI.2). Here we aimed at showing that it is possible to automatically estimate the age, based on the set of physical parameters [$S_{\rm ph}$, $P_{\rm rot}$, $T_{\rm eff}$, $\log(L/L_{\odot})$, [Fe/H]].

We split the sample of stars into two sets. First, a training set was used to train the algorithm to estimate the age from the physical parameters. Then, a test set composed of the remaining of stars was used to test and demonstrate the robustness of the method.

We tried different percentages of training sets from 20% to 80% and found very similar results for the accuracy (difference between the predicted age and the model one) in general with some improvements in the scatter. For the rest of the paper, we chose a training set of 50% of the stars that was a good trade-off to have a reduced scatter but still a large number of stars for the test set (more than 5,000 of them).

In order to take into account uncertainties on physical parameters [$\delta S_{\rm ph}$, $\delta P_{\rm rot}$, $\delta$[Fe/H], $\delta$log$L/L_{\odot}$], and to improve the training set, we randomly drew, for each of the stars ($\star$) in the training samples, $100$ artificial stars ($\star^{\mathcal{A}}$) by selecting their physical parameters as follow:

where $\textnormal{X}\in[S_{\rm ph}$, $P_{\rm rot}$, $T_{\rm eff}$, $\log L/L_{\odot}$, [Fe/H]], $\delta\textnormal{X}\in[\delta S_{\rm ph}$, $\delta P_{\rm rot}$, $\delta T_{\rm eff}$, $\delta\log L/L_{\odot}$, $\delta$[Fe/H] ] and $\mathcal{G}_{0\leq i\leq 100}$ random values following a standard normal distribution.

We also made two other runs with slightly different sets of parameters: one without $S_{\rm ph}$ that is similar to gyrochronology and one without $P_{\rm rot}$ that is similar to magnetochronology.

In Figure 8, we can see the importance of the different input parameters in the training of the RF algorithm. The so-called gini (Gini 1912) importance is a measure of how each parameter contributes to the homogeneity of the nodes and leaves in the trained random forest. The higher the importance value, the higher impact the parameter has during the training process. In the case of not including the rotation period in the RF (i.e. RF magnetochronology, top panel of Fig. 8), $S_{\rm ph}$ has the highest weight or importance, followed by the effective temperature, luminosity, and metallicity. This is in agreement with the results of the Bayesian fits done in Section V where $S_{\rm ph}$ has the largest coefficient.

In the case of RF gyrochonology (middle panel, Fig. 8), the highest importance is rotation period and then luminosity, while temperature and metallicity have very little impact.

Finally, for RF magneto-gyro-chronology, the picture changes drastically where the main feature that impacts the age estimate is still the rotation period followed by luminosity. The photometric magnetic activity proxy has very little impact, as much as [Fe/H]. This emphasizes that for our sample of late-F, G, and K dwarfs, rotation is the dominant factor for age determination and that the magnetic activity proxy is redundant compared to rotation period given the strong correlation between the 2 quantities (see Figure 7 of S21).

To obtain uncertainties on the ages predicted by the RF, we randomly draw, for each $\star$ in the test samples, $100$ $\star^{\mathcal{A}}$ by selecting their physical parameters following Eq. (10). The standard deviation of the distribution of the age predictions by the RF on the 100 $\star^{\mathcal{A}}$ is then measured for each $\star$ in the test sample.

In Figure 9 are reported such RF relative uncertainties in the case of magnetochronology (red crosses). For most stars in the test sample, the relative uncertainty is comparable to the original uncertainty of the age measurement (black crosses). The median uncertainties for the kiauhoku ages is 14.2% compared to 15% for the RF uncertainties. This result is also valid for gyro- and magnetochronology.

From the different relationships derived in Section V, we computed the predicted ages of the stars based on their stellar parameters. We then evaluated the accuracy of those ages compared to the modeled ones. This was done for the linear and multivariate fits for both magnetochronology and magneto-gyro-chronology relations. In the left panels of Figure 10, we can see the comparison between the ages predicted from relations presented in Section V and calibrated through a Bayesian approach and those coming from the stellar models in the case of magnetochronology (so no $P_{\rm rot}$ was used as input).

The top left panel shows the results for the linear fit (Eq. 2) for the solar analogs. We can see that there is a trend with an overestimation of the ages below 2 Gyr and an underestimation for stars older than 4 Gyr. In particular, this reflects that the analytical relation is not adequate enough to reproduce the age of the young stars. The median of the distribution of the differences between the predicted and modeled ages (bottom right panel) is of -6.3% with a Median Absolute Deviation (MAD) of 21.7%. We also notice that the distribution is skewed towards negative values. This means that for the solar analogs, a power law model (Eq. 2) between age and $S_{\mathrm{ph}}$, as studied in the past for other magnetic activity proxies, underestimates on average the ages with large uncertainties.

By adding the information on the rotation period to the linear fit (Eq. 7), the median difference between the model and predicted ages for the solar analogs is of 2.7%, with a smaller scatter of 19.1% (top right panel of Figure 10). Young ages below 2 Gyr are still not well retrieved. However, this analytical relation overestimates the ages on average with smaller differences in absolute value compared to the case without rotation. This clearly shows the importance of the knowledge of the rotation period to predict ages of solar analogs where no other stellar parameters are used in the analytical relation.

For the multivariate relation without rotation period ($\alpha_{5}$ = 0), we show the same comparison between the predicted and observed ages but for late-F and G dwarfs (middle left panel of Figure 10) and K dwarfs (bottom left panel of Figure 10). While the absolute values of the median differences with the stellar model ages are similar for both sets (-4.7% for the late-F and G dwarfs and 3.5% for the K dwarfs), the MAD is much larger for the K dwarfs (80.4% compared to 28.1%).

If we add the information on rotation period ($\alpha_{5}\neq$ 0), we clearly see an improvement on the ages prediction for both the late-F and G dwarfs (middle right panel of Figure 10) and K dwarfs (bottom right panel of Figure 10). For the former, the distribution of the differences is centered close to 0 (median of 0.1%) with a MAD of 14.1%. There is a slight bias of overestimated ages below 1.5 Gyr but on average the ages are well recovered. For the K dwarfs, the improvement is even more striking with a median difference of 0.2% and a MAD of 5.6% better than the late-F and G dwarfs. Such an improvement can be explained by the fact that rotation period is a stronger constraint of the ages of K dwarfs. Indeed rotation period is a stronger constraint on stellar age than the classical stellar parameters.

As mentioned in Section V, we also fitted the relations without $S_{\mathrm{ph}}$ to assess the importance of the magnetic activity proxy in the age prediction. The comparison between the predicted and kiauhoku ages is shown in Appendix D. While the results are similar to the case where we include $S_{\mathrm{ph}}$, the median differences and MAD values are slightly improved when the magnetic activity is taken into account.

We show in the top panel of Figure 11 the comparison of the predicted ages from the RF analysis without using the rotation period as input. This analysis was done on all late-F, G, and K dwarfs where the test set consisted of 5,452 stars. The predicted ages from the RF go close to 14 Gyr as the algorithm was trained on the full sample of late-F, G, and K dwarfs. We see a clear trend with increasing ages. Above 5 Gyr, the RF ages are biased towards smaller ages compared to kiauhoku ones whereas below 5 Gyr, the RF tends to overestimate the ages with a median difference of 15.9%. In particular for the old stars, the disagreement can be quite high due to the plateau that we observed above 5 Gyr, leading to some degeneracy. If we take into account the full range of model ages, the median difference goes down to 5.3% with a MAD of 48.8%.

When adding the rotation period as input (middle panel of Figure 11) there is a very good agreement up to ages of $\sim$ 10 Gyr with less scatter compared to the previous RF run. For the full set, we find a median difference between the RF ages and the kiauhoku ones of 3.1% with a MAD of 8.9%. This is around half of the median uncertainties from the RF.

Finally, we also tested the RF algorithm for pure gyrochronology (bottom panel of Figure 11), where we used $P_{\rm rot}$ without $S_{\rm ph}$ along with the atmospheric parameters as inputs. The median difference between the predicted ages and the model ones is of 4.2% with a dispersion of 7.2%, smaller than the RF uncertainties. By taking into account the magnetic activity proxy, the median difference is improved but the dispersion increases. The use of both inputs $P_{\rm rot}$ and $S_{\rm ph}$ improves in average the age estimates with the RF.

We note that the ages computed in this work depend on models. They rely on a long list of choices of input physics and magnetic braking prescription, as well as calibrators needed to produce expected behavior. While we consider that all our choices are well-motivated, they are inevitably based on an incomplete knowledge of stellar physics and rotation.

One limitation of gyrochronological ages is the range of available calibrators. Spin-down relations are calibrated using stars for which precise ages and periods are known independently, which limits us to the Sun, a handful open clusters, and more recently a few tens of asteroseismic F, G, and K dwarfs (Barnes 2010; Angus et al. 2015a; van Saders et al. 2016; Curtis et al. 2020; Godoy-Rivera et al. 2021, e.g.,). The open clusters are all young with the oldest one of 4 Gyr, M67, being studied recently (Dungee et al. 2022); thus, until recently the Sun was our oldest anchor for gyrochronology at 4.57 Gyr. The few old, seismic dwarfs hint that braking weakens at ages older than the Sun (van Saders et al. 2016; Hall et al. 2021), but the physical processes that weaken the braking are still poorly understood. Furthermore, the benchmark clusters are all near solar composition. Any gyrochronology predictions for chemically different stars are then extrapolations into regimes where there are no calibrators. Spin-down is less constrained in these regimes.

Spin-down relations are also less reliable in regimes where stars lose their convective envelopes. With no surface convection, there is no dynamo to generate a magnetic field, so no braking occurs. In this regime (Kraft 1967, toward and above $\sim$6300 K, the Kraft break after), rotation is not meaningful as an age indicator. We note, however, that our use of luminosity in addition to rotation to constrain the age allows for luminosity to take over when rotation becomes less constraining.

Besides the evolution models themselves, the parameters we use as input to the MCMC fitting are also model dependent. Effective temperatures and metallicities are determined by fitting atmospheric models to spectra, which have different scales than those used in the evolution models. Furthermore, the luminosities are determined by fitting yet another set of models to Gaia parallaxes and 2MASS colors (Berger et al. 2020b). Some authors calibrate their stellar parameter estimates to a more standard scale (Holtzman et al. 2015, e.g., APOGEE;), but using different models with different parameter scales necessarily introduces uncertainty into the model estimates (Tayar et al. 2022).

When fitting evolution models to stars, we assume solar abundances of $\alpha$ elements. $\alpha$-abundances influence spin-down independently from metallicity alone, so this assumption adds uncertainty to the estimates as well. For example, Claytor et al. 2020a showed that assuming solar $\alpha$-abundance for a main-sequence star that actually has [$\alpha$/M] = 0.4 results in underestimating the age by as much as 20%.

Finally, our angular momentum evolution model in YREC assumes that stars spin down as rigid bodies, which may be a poor assumption at some stages of stellar evolution. Open cluster stars exhibit a temporary epoch of stalled spin-down in K dwarfs at 1 Gyr (Curtis et al. 2019) which could be explained by a two-zone angular momentum model such as that of Spada & Lanzafame 2020. The departure from expected behavior means that there are physical mechanisms our gyrochronology model does not capture.

Some of these limitations have the effect of stretching or compressing our age scale, while others result in a systematic offset. Taken together, the effects mean that while our ages are internally consistent and precise, any individual age taken from this data set should be vetted cautiously before being used to support new claims.

We investigate whether the plateau above 5 Gyr seen in Fig. 6 is real. To evaluate what is the level of bias on the stars where we detected rotation periods in S19 and S21, we study here a sample of more than 100,000 solar-like stars observed by the Kepler mission for which no rotation periods could be reliably measured for different reasons (instrumental artifacts, low-amplitude signal, low inclination angle, or no stellar magnetic activity during the Kepler observations).

By definition $S_{\rm ph}$ is a proxy for magnetic activity measured from the rotational modulation in the light curves. This means that when it is not possible to retrieve a reliable $P_{\rm rot}$, we cannot derive $S_{\rm ph}$. In this section, in order to investigate the detection limit for these measurements, we define a variability index $V_{\rm ph}$, which can be obtained for stars regardless a $P_{\rm rot}$ was measured. We split the light curve in subseries of length 30 days to compute the standard deviation and we take the mean value of all subseries. Similarly to the $S_{\rm ph}$, we subtract the photon noise obtained from the formula by Jenkins et al. 2010b. We call this quantity the variability index, $V_{\rm ph}$, as the measured quantity may not be purely related to the magnetic activity of the star, in contrast to $S_{\mathrm{ph}}$.

In the comparison that follows, we also compute $V_{\rm ph}$ for stars with $P_{\rm rot}$ estimates. For these stars, $V_{\rm ph}$ is still related with magnetic activity as we find that $S_{\rm ph}$ and $V_{\rm ph}$ are very similar (see Appendix E). The comparison between $V_{\rm ph}$ for the rotation sample and the no rotation sample in Figure 12 shows that the stars with detected rotational modulation have larger $V_{\rm ph}$ in general than the stars without detection. This can also be seen in the histogram in the right panel of Figure 12. For active stars, indeed it is easier to observe rotational modulation in epochs of larger magnetic activity. Still, some stars without $P_{\rm rot}$ have large $V_{\rm ph}$ values similar to those of the stars with $P_{\rm rot}$. This could be related to instrumental artifacts, which might still be present in the light curves. A lack of modulation detection could also be due to a fully spotted star, a star observed pole-on with active latitudes close to the equator, or with high active latitudes and long-lived spots.

The rotation analysis done by S19 and S21 extended the rotation periods catalog towards slower rotators, which usually have low $S_{\mathrm{ph}}$ or $V_{\rm ph}$ values. Hence in this work, we derived ages older than the age of the Sun for stars with $P_{\mathrm{rot}}$ showing that we can still measure rotation periods for these old stars. This is slightly in contrast to the findings of Masuda 2022. Indeed the author studied a sample of 278 planet-host stars selected as having “robust” rotation periods from Mazeh et al. 2015 and computed ages with isochrone fittings based on the models of Angus et al. 2015b that also take into account surface rotation periods. He finds that there is a bias in the rotation detection and stars older than the Sun do not have a rotation periods measured. However, as pointed out by Masuda 2022 for the S21 sample used in our work, the situation might be different given the larger number of stars with longer rotation periods.

We looked at stars older than 5 Gyr in the $S_{\rm ph}$-age diagram to locate them in Figure 12. We found that most of the stars are above the lower edge, i.e. above the limit of the detection (see Figure 13). This points in the direction that the level of magnetic activity drops past a given age with a certain flat evolution afterwards. This does not discard the fact that old stars could be even less active.

From previous works, the magnetic activity of the Sun between the minimum and maximum of the solar cycle varies between 67.4 to 314.5 ppm (Mathur et al. 2019). We added those data points in the $S_{\mathrm{ph}}$-age diagram for solar analogs (see Figure 6). We can see that the magnetic activity level of the Sun agrees well with its counterpart for a similar age. Hence, it seems that the Sun behaves like Sun-like stars selected from spectroscopic parameters. This also agrees with the analysis of the variability of $S_{\mathrm{ph}}$ for solar analogs done by Santos et al. 2023.