© The Authors 2024

1 Introduction

Determining the ages of individual stars is among the most difficult problems in astrophysics because they cannot be measured directly. The most direct method to date is nucleochronometry: radioactive dating of the oldest meteorites provides a solid estimate of the Solar System’s age, and therefore of the age of the Sun (e.g. Chaussidon 2007), while long-half-life isotope abundances inferred from the spectra of a few very old metal-poor halo stars give a rather straightforward access to their age (e.g. Christlieb 2016). For other stars, several methods are used in the literature to estimate ages (see Soderblom 2010, for a review), even though no single method is valid for all age ranges or spectral types.

One can classify the different techniques into model-dependent ones, which make use of stellar evolutionary models, and empirical relations between the age and a given stellar observable. In the latter type, the underlying physics of the empirical relations is usually not fully understood and needs to be calibrated on samples of stars with high-quality ages (obtained by model-dependent methods). One example is gyrochronology, which uses the empirical relation between a star’s rotation period and its age (Barnes 2007), or the correlation between the activity or lithium abundance of F, G, K stars and age. Another possibility, which has recently been given a boost with the advent of large spectroscopic surveys, is to use abundance-age relations, which display linear dependences on age partly explained by stellar and/or Galactic evolution. This could be the case of the [α/Fe]-age relation (Haywood et al. 2013; Bensby et al. 2014; Haywood et al. 2015; Ciucă et al. 2021; Katz et al. 2021), the CNO abundances (which are mostly driven by stellar evolution; Masseron & Gilmore 2015; Lagarde et al. 2017), or the so-called chemical clocks (e.g. [Y/Mg] or [Sr/Mg], which are driven more by a smooth Galactic evolution of s-process elements; Nissen 2015, 2016; Tucci Maia et al. 2016; Delgado Mena et al. 2019; Jofré et al. 2020; Casamiquela et al. 2021). As a generalisation, spectroscopic ages can be derived using supervised machine-learning techniques that take advantage of the full chemical information of stars trained on high-precision datasets, such as those based on asteroseismology (Ciucă et al. 2021; Hayden et al. 2022; Anders et al. 2023; Boulet 2024). This approach, however, has a significant dependence on which stellar population the training is performed on, and the age estimations are intrinsically entangled with Galactic chemical evolution; thus, there is a potential redundancy when using them to study the chemical evolution of the Milky Way.

The methods that are considered most reliable are those that use stellar evolution models; that is, asteroseismology or isochrone placing, for which the underlying physics is rather well understood (but see Lebreton et al. 2014a,b) and which rely on the fewest assumptions. Asteroseismology provides a precise way to constrain ages through the measurement of solarlike oscillations. Until the advent of the PLATO mission (Rauer et al. 2014), there are only limited samples of stars (≲10 000) with asteroseismic constraints and only in selected fields of the CoRoT (Baglin et al. 2006), Kepler (Borucki et al. 2010), K2 (Howell et al. 2014), and TESS (Ricker et al. 2015) missions (see e.g. Miglio et al. 2017, for a review). On the other hand, placing isochrones becomes ‘easy’ in open and globular clusters, since they have many coeval mono-metallic members. This allows a distribution fit over the whole range of mass being fit, leading to a generally well-constrained situation¹. That is why stellar clusters provide the primary benchmarks to study age-related properties. In fitting a cluster age, one has many stars distributed in mass in an HR diagram, and it is the full behaviour of the models over that mass range that is being fit, leading to a highly precise result overall. For individual field stars, the distribution fit is not possible, so isochrone placement becomes challenging because multiple isochrones can pass through a given point in the HR diagram given the degeneracies in the different observables and the underlying photometric uncertainties. Indeed, this method can only be applied with good precision (≲15%) for individual stars near the main-sequence turn-off (MSTO) or in the subgiant branch (SGB), where evolutionary models of different masses separate in the HR diagram, providing less degeneracy (Jørgensen & Lindegren 2005). Therefore, both methods, asteroseismology and isochrone placement, are complementary in terms of the spectral types of stars for which precise ages can be obtained.

The advent of Gaia, with its precise photometry and parallaxes for a vast number of stars, has enabled massive stellar parametrisation. With Gaia DR3 (Gaia Collaboration 2023), stellar parameters for 471 million sources were estimated from low-resolution blue photometer (BP) and red photometer (RP) spectra (Andrae et al. 2023), together with the chemo-physical parametrisation for 5.6 million stars estimated from the Radial Velocity Spectrometer (RVS) data (Recio-Blanco et al. 2023). The FLAME module then produced luminosities, radii, masses, and ages, among other parameters, for 284 million stars in Gaia DR3. Ages, computed mainly using atmospheric parameters (T_eff, log g, and [M/H]) and absolute magnitude as inputs, were produced for half of the sample (Creevey et al. 2023).

Many other studies in the literature have used Gaia in combination with other external data to produce more precise ages (e.g. McMillan et al. 2018; Sanders & Das 2018; Wang et al. 2023; Stone-Martinez et al. 2024). For instance, Xiang & Rix (2022) recently computed ages for subgiant and MSTO stars that were observed as part of the seventh data release (DR7) of the largescale spectroscopic survey LAMOST (Cui et al. 2012; Zhao et al. 2012). They made use of absolute magnitudes in the K band, and spectroscopic estimates of T_eff, [Fe/H], and [α/Fe] abundances from the DD-payne pipeline in LAMOST DR7. Queiroz et al. (2023) recently used the StarHorse code (Queiroz et al. 2018; Anders et al. 2022) to derive stellar ages for millions of MSTO and SGB stars observed by several spectroscopic surveys, including APOGEE (Majewski et al. 2017), GALAH (De Silva et al. 2015), and Gaia RVS. In this case, in addition to Gaia DR3 data and spectroscopically determined atmospheric parameters (effective temperature, surface gravity, and metallicity), it uses infrared photometry from 2MASS and AllWISE to infer stellar astrophysical parameters, including ages. Kordopatis et al. (2023) determined Bayesian isochrone age estimates for stars observed by Gaia’s RVS. They used GSP-Spec (Recio-Blanco et al. 2023) calibrated atmospheric parameters, 2MASS and Gaia-EDR3 photometry, and parallax-based distances to compute ages, initial stellar masses, and reddenings for 5 million stars with spectroscopic parameters in Gaia-DR3.

This study aims to determine ages using only Gaia DR3 photometry (G,G_BP,G_RP) and parallaxes, in combination with recent interstellar extinction 3D maps, and chemical abundances from the data release (DR8) of LAMOST. In contrast with previous age estimates, here we do not make use of a particular stellar parametrisation (e.g. spectroscopic T_eff or log g) because the information of the HR diagram is provided only through the absolute Gaia magnitudes. Indeed, one of the strong points of the Gaia mission is to deliver very precise and accurate photometry and parallaxes, which makes it possible to date stars in an HR diagram with the lowest possible uncertainties. In turn, spectroscopic information on the metallicity and chemical composition (the α abundance in particular) allows us to break some degeneracies of isochrone fitting. We used the recent code SPInS (Lebreton & Reese 2020), which uses a Bayesian framework with an MCMC sampler coupled with an interpolation scheme, to parameterise stars using stellar evolutionary models.

The paper is structured as follows. Section 2 explains the basics of our method (for a detailed paper, we refer the reader to Lebreton & Reese 2020). Section 3 describes the selection of the LAMOST and Gaia DR3 data used in this paper. In Sect. 4, we present a detailed comparison with open and globular clusters to validate our method, and in Sects. 5 and 6 we show and discuss the results of applying our method to the field star sample described in Sect. 3. Finally, we conclude the paper in Sect. 7.

2 The method

We aim to obtain stellar parameters, in particular ages, of a large sample of stars (selected as described in Sect. 3). To this end, we used the code Stellar Parameters INferred Systematically (SPInS Lebreton & Reese 2020), a public Python pipeline² that takes different types of inputs (e.g. photometric, spectroscopic, interferometric, and/or averaged asteroseismic) to provide the age, mass, and radius (among others) of a star, relying on a grid of evolutionary tracks. In brief, the code works in a Bayesian framework to provide the posterior probability distribution function (PDF) of the inferred stellar parameters from a set of observational constraints, a grid of stellar models (see Sect. 2.1), and a set of priors. The PDF is sampled using an MCMC solver based on the emcee Python package (Foreman-Mackey et al. 2013), coupled with a versatile interpolation scheme for the stellar models (see Sect. 2.2).

2.1 Stellar grids and priors

We used the BaSTI³ grid of stellar models (Hidalgo et al. 2018), calculated for a solar-scaled heavy elements distribution and updated input physics, including atomic diffusion of helium and metals, overshooting of convective cores, and mass loss. The corresponding solar mixture was that of Caffau et al. (2011) complemented by Lodders (2010). Evolutionary tracks are provided for a set of masses and metallicities (see below), with a helium abundance derived assuming a helium-to-metal enrichment ratio of ΔY/ΔZ = 1.31 (see Hidalgo et al. 2018, for more details). For each stellar model of a given age, mass, and chemical composition, the luminosity and effective temperature are provided, as well as Gaia EDR3 magnitudes, M_G, M_GBP, and M_GRP. On the SPInS website, the previous BaSTI stellar evolution tracks (Pietrinferni et al. 2004) are provided in a format directly readable by SPInS. However, any stellar evolution grid, calculated with any stellar evolution code, can be used if it is written in SPInS input format.

For the sake of homogeneity, in order to always use the same models for all stars, we have not used the α-enhanced model grid, but rather the solar-scale grid, scaling the input metallicities to mimic the α-element enrichment. This was done via the commonly used relation derived by Salaris et al. (1993)⁴.

The whole BaSTI grid contains a total of 1120 evolutionary tracks with initial masses in the range of M₀ ∈ [0.1 M_⊙, 15 M_⊙] and metallicities of [M/H] ∈ [−3.197, +0.3]. The grid covers all evolutionary stages from pre-main-sequence up to either the first thermal pulses on the AGB or C-ignition, or to the age of the Universe, depending on the mass. In the grid prepared for SPInS, we have excluded high masses of M > 10 M_⊙, the pre-main-sequence parts of each track (t_adim < 0.05), and the phases beyond the RGB tip, given the low probability of these objects being found in our sample (described in Sect. 3). t_adim is an adimensional age parameter that goes from 0 to 1, that has a value dependent on the evolutionary stage, and that is thus homogeneous for all masses. The stellar evolutionary model grid used has a total of 978 732 points.

Interpolation on the grid was done in the (log(M₀/M_⊙), [M/H], t_adim) space, at each MCMC iteration. On the one hand, this involves the interpolation between evolutionary tracks, which is coded as a linear barycentric interpolation on a simplex defined by a Delaunay tesselation on the grid of models. On the other hand, the interpolation along the evolutionary track is linear between adjacent points, performed on t_adim, with the purpose of combining models at the same evolutionary stage. Then, a transformation to physical age (t) was implemented in SPInS for each MCMC step to provide a sample on physical age. We refer the reader to Lebreton & Reese (2020) for more technical details about the interpolation scheme.

SPInS allows the use of priors on the grid parameters, such as the initial mass function (IMF), the metallicity distribution function, or the star formation rate. In the case of this work, we have not imposed any strict prior on metallicity, mass, or age to avoid biases in the statistical interpretation of the resulting trends. This is motivated by the fact that overall systematic offsets can be present in the derived ages, depending for instance on the set of evolutionary tracks or the photometric transformations. This means that ages larger than the age of the Universe are allowed, as much as this is allowed by the evolutionary tracks.

2.2 Configuration of the MCMC sampler and evaluation of the convergence

There are several options concerning the MCMC sampler integrated into SPInS that need to be fixed in the code on a case-by-case basis. Firstly, we chose to do the initialisation of walkers based on a Gaussian distribution centred around the best-fitting model. This is a convenient option to reach quick convergence, since it already places the starting points near a preferred position, although later the walkers can move and explore more distant points.

Other free parameters such as the number of walkers, burnin, and production steps have a large dependence on the complexity of the posterior and need to be properly set to allow a correct exploration of the parameter space. This is a crucial issue because, if the PDF is multimodal (as is often the case in an HR diagram), it can be tough to sample it with a standard MCMC, and this can result in biased solutions in some regions of the HR diagram. A possible solution to this problem is to use a parallel-tempering ensemble MCMC that runs in a modified posterior given by a transformation of the likelihood as $${\mathcal{L}}^{\mathrm{\prime }}(x)=\mathcal{L}(x{)}^{(1/T)},$$(1)

where T is a parameter usually called ‘temperature’ (Vousden et al. 2015). Using higher T values allows MCMC chains to explore the parameter space more easily because the likelihood is flatter and broader. The parallel-tempered MCMC is implemented to be used in SPInS via the package ptemcee⁵. In our case, we have seen that this strategy dramatically improve the convergence of the chains. However, the usage of parallel tempering implies that, for a fixed number of iterations, the computation times increase with the number of temperatures, reaching several minutes per star. Therefore, we have found a need to evaluate the convergence on the fly for each star, in order to accelerate the computational run.

Evaluating the convergence of the walker chains is an essential step in an MCMC analysis. Even though it is formally impossible to guarantee the convergence of an MCMC sampler, there are some diagnostic tools to evaluate if we are obtaining an accurate approximation of the PDF. We chose to use a criterion based on the integrated autocorrelation time (τ), whose basic idea is that the chain has sampled long enough when the walker has traversed the high-probability parts of the parameter space many times in the length of the chain. Following the discussion of Hogg & Foreman-Mackey (2018), a small value of τ compared with the length of the chain can be used as a sign of convergence.

In this work, we have implemented in SPInS a robust on-the-fly computation of the autocorrelation time⁶ done while sampling the MCMC every n_steps. We considered the sample to have reached convergence when the mean of τ in the three main dimensions of the grid (log (M₀/M_⊙), [M/H], t) reaches τ < N/100, where N is the number of steps that the chain has sampled. With the achievement of the convergence criteria, we stopped the MCMC sampler and drew the statistical analysis of the resulting PDF. This strategy is particularly useful for running SPInS in a massive way because it allows us to save computation time once convergence is reached, and provides an automatic evaluation of the goodness of the PDF when the MCMC has finished.

A new version of SPInS with improved features, including the on-the-fly computation of the autocorrelation time and the automatic convergence evaluation used in this paper, is publicly available. Our tests show that with 10 temperatures, 10 walkers, and 1000 burn-in steps, most of the stars in our sample (MSTO and subgiants, see Sect. 3) reach convergence in around 2000 production steps. With this configuration on an Apple M1 processor, SPInS takes around 20 seconds per star using the native SPInS parallelisation on two processes. On top of this, in this study we have parallelised among 20 nodes the sample of MSTO and SGB stars detailed in Sect. 3.3.

3 Data selection

SPInS is very flexible in terms of input observational constraints; any parameter included in the grid of evolutionary models can essentially be used as a constraint. In this paper, except for spectroscopic metallicity estimates, we have exclusively used photometric and astrometric constraints: M_G and (G_BP − G_RP)₀, obtained from Gaia DR3 photometry and parallaxes, coupled with a reddening estimate from the recent three-dimensional extinction maps provided by Vergely et al. (2022). These maps are based on the inversion of large spectroscopic and photometric catalogues including Gaia DR3, 2MASS, and AllWISE.

To build our final sample, we performed an initial query in the Gaia DR3 archive, selecting stars brighter than magnitude G = 18 and in a 3 kpc region around the Sun. Being able to retrieve accurate absolute magnitudes and intrinsic colours is essential for the isochrone placement method. We performed cuts on the relative uncertainties in the parallax and in the three magnitudes. We performed cuts on Galactic latitude to avoid high-extinction regions in the Galactic plane, where the extinction maps are less accurate. Additional cuts were applied on the normalised unit weighted error (RUWE), the percentage of successful Image Parameter Determination (IPD) windows with more than one peak (ipd_frac_multi_peak), and the amplitude of the IPD goodness of fit, to avoid binary or multiple stars. The complete query is:

SELECT * FROM gaiaedr3.gaia_source WHERE
parallax_over_error > 10 AND
phot_g_mean_flux_over_error>50 AND
phot_rp_mean_flux_over_error>20 AND
phot_bp_mean_flux_over_error>20 AND
phot_g_mean_mag < 18 AND
parallax > 0.33 AND
ipd_frac_multi_peak < 2 AND
ipd_gof_harmonic_amplitude < 0.1 AND
ruwe<1.1 AND (b<-10 OR b>10) AND dec>-10.

We have additionally performed a cut in declination because we later crossmatched the sample with LAMOST (see Sect. 3.2), which is entirely contained above −10°. This query gives a total of 23 million stars.

3.1 Building the intrinsic colour-magnitude diagram

For the computation of interstellar absorption, we chose the map from Vergely et al. (2022) that covers a volume of 3 kpc × 3 kpc × 800 pc around the Sun at a resolution of 10 pc. The individual A_V values per star were computed with a linear 3D interpolation of the absorption density towards the line of sight, followed by an integration along the line of sight. Some of the stars are located outside the map at |Z| > 800 pc. Accordingly, we consider the absorption to be zero at higher |Z| values. Absolute magnitudes were computed using an inversion of the parallax.

We applied zero-point corrections on Gaia DR3 parallaxes, as is recommended by Lindegren et al. (2021). In all cases, we obtained the apparent magnitudes G, BP, RP from the photometric fluxes and corrected for Gaia colour zero points following Riello et al. (2021). To correct for interstellar extinction, we used the empirical relations provided by Danielski et al. (2018) with coefficients updated for Gaia DR3⁷. For instance, the extinction coefficient for the G band is k_G((G_BP − G_RP)₀, A₀) = A_G/A₀, where (G_BP − G_RP)₀ is the star’s intrinsic colour, and A_G and A₀ are the extinction in the G band and at λ₀ = 550 nm, respectively. For stars that have (G_BP − G_RP) ≤ −0.06, which is the limit of validity of Danielski et al. (2018)’s formulae, we used Wang & Chen (2019)’s law. For A₀, we used the value of A_V. The extinction maps used do not have detailed uncertainty values; thus, we warn the reader that the derived uncertainties on the absolute magnitudes are a lower limit.

3.2 Metallicities and α abundances

The previous sample was finally crossmatched with the DR8 catalogue for A, F, G, and K stars of the LAMOST spectroscopic survey to use the [Fe/H] and [α/Fe] abundances as additional constraints, which significantly breaks degeneracies to get the physical parameters from evolutionary tracks.

LAMOST DR8 contains:

• a spectroscopic parameter catalogue of 6.6 million stars from the low-resolution survey (LRS), which provides atmospheric parameters, radial velocities, and iron and α-element abundances obtained with a spectral resolution of 1800 in the wavelength range of 3700–9000 Å.

• a spectroscopic parameter catalogue of 1.2 million stars from the medium-resolution (7500) survey (MRS), which provides, additionally, overall α-element abundances and individual abundances for certain elements. In this case, two sets of parameters and abundances are listed, one from the LASP pipeline (Wu et al. 2014), and another from a convolutional neural network (CNN) that up to DR8 was called Data-Driven Payne (DD-Payne). Soubiran et al. (2022) did a thorough comparison of the [Fe/H] results among spectroscopic surveys, finding that DD-Payne for LAMOST DR5 tends to give larger biases than LASP in the metal-poor regime ([Fe/H] < −1) when compared to higher-resolution studies. In this work, we have used both estimations of [Fe/H] and [α/Fe] (LASP and CNN) to compute two sets of MRS ages.

Cuts on the [Fe/H] uncertainty of <0.05 and the radial velocity uncertainty of <5 kms⁻¹ were applied to both samples to filter low-quality values and problematic spectra. Once we had added the Gaia selection (see previous subsection), we ended up with a total of 490 233 stars in the LRS sample and 71 473 in the MRS sample.

Figure 1 shows the heliocentric X–Y and X–Z distributions of the LAMOST DR8 LRS sample, together with the corresponding colour-magnitude diagram (CMD) corrected for extinction. The intrinsic CMD in the bottom left panel, coloured by density, exhibits a thin main sequence, which also makes visible an equal-mass binary sequence in the faint magnitudes. Clear turn-off, subgiant, and red giant branches can be seen, together with a prominent red clump of helium-burning stars. Below the red clump, there is a sign of the red giant branch bump, caused by a discontinuity in the luminosity of hydrogen shell-burning stars (see e.g. King et al. 1985). The bottom right plot of Fig. 1 shows the same intrinsic CMD, but now binned and coloured depending on mean [Fe/H], which highlights a gradient with metallicity in the horizontal direction, where metal-rich stars tend to be redder and metal-poor ones bluer. This is particularly remarkable for the main sequence, red giant branch, and red clump, as was expected from stellar evolution models. The secondary red clump is visible as a small compact group towards the bluer and fainter part of the main red clump. This corresponds to more massive red clump stars, which are probably younger (Girardi 1999), and this is consistent with the plot, since one would expect them to be overall more metal-rich.

Fig. 1 Distribution of the initial sample of 352 k stars in the LAMOST DR8 LRS sample obtained from the procedure indicated in Sect. 3. Top: galactic X–Y and X–Z distributions coloured by stellar density. The Sun is at (0, 0). Bottom: intrinsic CMD coloured by density (left), and binned with colours representing mean [Fe/H] (right).

3.3 Selection of main-sequence turn-off and subgiant stars

It is widely known (e.g. Pont & Eyer 2004; Takeda et al. 2007) that ages for individual stars derived by isochrone fitting are only reliable for the MSTO and SGB regimes because these positions in the Hertzsprung-Russell diagram have fewer degeneracies in the evolutionary models for different masses. For stars in the red giant branch and the main sequence, SPInS typically obtains very flat or ill-defined PDFs that are uninformative of the age of the star (see tests done with synthetic stars in Lebreton & Reese 2020).

Thus, we selected the MSTO and SGB regions in bins of [Fe/H], as is shown in Fig. 2, similar to the MSTO and SGB selection done in Queiroz et al. (2023). The faint end of the MSTO was selected with a horizontal cut at magnitude M_G = 4.5, and the redder limit of the SGB was selected with a linear function of M_G = 20 · (G_BP − G_RP)₀ + b_i, where b_i = [18, 17, 16, 15] for the metallicity bins [Fe/H] = (0.5, 0.0),(0.0, −0.5),(−0.5, −1.0),(−1.0, −2.5), respectively.

The selection gives 243768 stars in the LRS and 35 096 stars in the MRS. Both samples consist of stars with exquisite photometric and astrometric quality with parallax errors of 1% in mean. This provides very good uncertainties for the absolute photometry that are in general smaller than 0.04 in M_G, and of the order of 0.01 in G_BP − G_RP (see Fig. 3). Metallicity errors have values of 0.03 dex on average, for both samples.

4 Validation sample: Open and globular clusters

We validated the method of determining the ages of individual field stars using open and globular clusters. Even though they are also model-dependent, star cluster ages are among the best anchors for validating age estimates of field stars, since they represent, in general, mono-age, mono-metallicity populations, and thus it is possible to do a distributed isochrone fit across the entire mass range. Indeed, the observed CMDs of star clusters serve as important calibrators for stellar evolutionary models.

In our case, we selected a set of validation stars in stellar clusters over a wide range of ages (see Table 1), similarly selected as our sample of field stars. We then computed the ages of the individual cluster member stars, independently of the cluster they belong to, then compared them with cluster ages obtained in the literature. Most of the literature cluster ages were obtained with a distributed isochrone fit on Gaia photometry, using different stellar isochrones and fitting methods. We also included ages determined using eclipsing binaries, asteroseismic ages for giants, and using the white dwarf cooling sequence, as is described in the following subsection.

4.1 Cluster and star selections

We first selected a suitable sample of clusters, starting from the open cluster catalogue of Cantat-Gaudin et al. (2020), which we crossmatched with Gaia DR3. We restricted the sample of stars using the filters on photometric and astrometric quality from the Gaia catalogue specified in Sect. 3, as was done for the field stars. To allow for better statistics in the comparison, we only kept the clusters that have a significant number of selected stars (>500), and that have an age and extinction determination in the Cantat-Gaudin et al. (2020) catalogue. This results in a sample of 9 k stars in 13 open clusters. We also searched for [Fe/H] determinations of the selected clusters in the literature (Casamiquela et al. 2021; Netopil et al. 2016), prioritising those coming from LAMOST (Zhong et al. 2020) for consistency with Sect. 3.

We noticed that for some open clusters (NGC 1039, NGC 2287, NGC2516, NGC 3532, NGC2447, NGC 2632, NGC 2682, and NGC 188), the extinction values derived by the automated procedure from Cantat-Gaudin et al. (2020) are underestimated by 0.05 up to 0.1 mag with respect to other literature values. In particular, Gaia Collaboration (2018) and the recent determinations from Tsantaki et al. (2023) provide very coherent values among them for the clusters in common, and also compared to previous literature studies. The assumed extinction has an impact in the derived ages, so we decided to use the A_V values from Tsantaki et al. (2023) and Gaia Collaboration (2018) for the eight clusters mentioned.

We list in Table 1 the clusters’ physical characteristics from the literature used in this work, including mean distances, absorption (A_V), and [Fe/H] abundances. We also include a non-exhaustive list of literature ages for these clusters, taken mainly from the studies of Cantat-Gaudin et al. (2020), Gaia Collaboration (2018) and Tsantaki et al. (2023). We find it relevant to include in this list ages coming from eclipsing binaries for NGC 6819 (Brewer et al. 2016) and NGC 188 (Meibom et al. 2009), as well as asteroseismic ages for giants obtained by Rodrigues et al. (2017) for NGC 6819.

As a test for old ages, we used globular cluster members from the catalogue from Gaia DR3 Vasiliev & Baumgardt (2021). We applied a strategy analogous to the case of open clusters, filtering stars according to their photometric and astrometric quality. We additionally required the cluster to be closer than 3 kpc and to have a well-populated turn-off with high-probability members (> 0.7) at magnitude G < 18. This selection yields only two clusters, NGC 6397 and NGC 6121 (M 4). We discarded NGC 6121 for this validation, since it is highly extincted due to its location behind the Upper Scorpius star-forming region. Its CMD shows a very broad main sequence and MSTO, pointing to a significant differential reddening. Therefore, we kept as the only suitable globular cluster NGC 6397, which is located at 2.488 ± 0.019 kpc with an overall reddening of E(B − V) = 0.18 (Baumgardt & Vasiliev 2021). We adopted [Fe/H] = −1.99 ± 0.01 and [α/Fe] = 0.46 ± 0.04 (from [Mg/Fe]) based on the mean values from a high-spectroscopic-resolution analysis of 13 member stars (Carretta et al. 2009). We also list in Table 1 its properties, including several age determinations from the literature obtained with isochrone fitting (Correnti et al. 2018; VandenBerg et al. 2013), as well as an age derived from its white dwarf cooling sequence (Torres et al. 2015).

We obtained an intrinsic CMD for the cluster sample in the same way as for field stars, described in Sect. 3.1. We applied the zero-point corrections on Gaia DR3 parallaxes to obtain the distances for the individual stars of the open clusters. For the globular cluster, we did not apply this correction but instead adopted for all stars the mean distance carefully derived by Baumgardt & Vasiliev (2021), to avoid enlarging the scatter in the SGB of the cluster. We used the mean absorption values per cluster from Table 1, and we applied the zero-point corrections and empirical calibrations to derive an extinction-corrected CMD.

We then selected MSTO and SGB stars of the clusters in the same way as was done for the field stars in Sect. 3.3 and Fig. 2. For NGC 6397, we additionally set a limit on bright magnitudes (M_G > 2) to exclude horizontal branch stars. This gives a selection of 4374 stars in 14 clusters.

Fig. 2 Intrinsic CMD of the initial sample of LRS (top) and MRS (bottom) stars, divided in four [Fe/H] bins, coloured by the density of stars (logarithmic scale). The number of stars is plotted in each panel. The grey lines represent the cuts done to select MSTO and SGB, as is explained in Sect. 3.3. To guide the eye, we additionally plot in blue two isochrones representative of each [Fe/H] bin at two different ages.

Fig. 3 Distribution of uncertainties in absolute magnitude and colour for the selection of MSTO and SG stars in the MRS (bottom) and LRS (top) samples.

4.2 Results per cluster and literature comparison

We ran SPInS star by star with the observational constraints (M_G,(G_BP − G_RP)₀) and metallicities from the literature, as was indicated in the previous subsection. Since we used solar-scaled mixture evolutionary models, in the case of the globular cluster NGC 6397 we took into account its [α/Fe] determination to rescale the metallicity to mimic the α-element enrichment, as is explained in Sect. 2.1.

The results from SPInS give a PDF sampled by the MCMC of the three fitted parameters (age, metallicity, and mass), for which we can take the median value as the best estimate, and the 16^th and 84^th percentile to compute the uncertainties. We show in Fig. 4 the results of the age determination star by star for four of the analysed clusters. We plot the intrinsic CMDs per cluster, where each analysed star is coloured according to its median age. We can extract several conclusions on the performance of our age determination method, which we describe below:

• We notice that stars in the lower main sequence (far from the turn-off) tend to have overestimated ages. This is no surprise since, for a given mass, an evolutionary track stays in nearly the same place of the main sequence for a long time. These stars will therefore present broader PDFs in the age dimension, which overall will result in larger median ages, particularly for the less massive members of young clusters. These lower main-sequence stars will also have larger uncertainties derived from the 16^th/84^th percentiles of the PDF, which means that we can filter them out using a cut on their relative uncertainties.

• We can also see the biases due to the presence of blue stragglers (BSs) in old clusters such as NGC 6819 and NGC 6397, which will tend to provide younger ages than the real cluster age.

• Finally, we also notice that unresolved binaries (particularly those that belong to the equal-mass binary sequence) tend to provide older ages than a single star at a given magnitude because its locus in the CMD corresponds well with an old turn-off star. We have not done any hard cut on the RUWE in this sample of stars because it penalised most of the stars belonging to the turn-off of a few clusters (particularly NGC 2287 and NGC 2447), and thus biases the final results.

The discrepancies in the cases of BSs and equal-mass binaries are expected simply because we cannot use standard evolutionary models to describe multiple stellar systems and exotic objects. In this test case using clusters, we take advantage of the fact that we can easily identify these stars by eye in a CMD and remove them manually from the sample. From the total of 4374 stars in our sample of clusters, we identified 58 BSs and 464 binaries, which represents 11% of the sample. We have not found a way to filter them out in a general way for the entire range of possible ages in the scientific case of field stars. However, the fraction of BSs and unresolved binaries found in the clusters allows us to set an order of magnitude of the ‘contamination’ rate that one can expect from a sample of field stars selected in the same way.

In the bottom plots of Fig. 4, we show the histogram of the median age obtained per star in grey, and in blue we only represent stars filtered according to the following criteria: (i) We remove stars that give large uncertainties on age, restricting the set of stars to those with relative uncertainties lower than 20% and uncertainties better than 500 Myr (a cut on absolute age uncertainty is needed to filter old stars with large uncertainties), (ii) we exclude the identified BSs and binaries. We notice a general improvement in the consistency with literature ages by applying these filters. In Fig. 5, we plot the median age values per cluster obtained from the blue histograms, with error bars representing the 16^th/84^th percentiles, coloured by the mean literature age.

To have a better description of the involved uncertainties in a Bayesian manner we can describe the age PDF of the cluster as the sum of the individual PDFs of the member stars. This allows us to keep the full information given by SPInS in the final age distribution, including stars with multiple solutions. We show the age distribution for each cluster in Fig. 5 representing in the form of violin plots the 95% of the cumulative distribution functions. As was done for the previous case, we have filtered the stars previously identified as BSS, unresolved equal-mass binaries or stars with uncertainties larger than 20% or 500 Myr. We also show in diamond symbols the literature values for each cluster listed in Table 1. In this way we are able to see the details of the posterior distributions, which are highly skewed, sometimes presenting multiple peaks.

In general, we see a good consistency between our age estimates and the literature ages. For young clusters (<1 Gyr), even though the peaks of the distributions are very close to the literature values, we have a tendency to slightly overestimate the ages of these clusters due to a tail towards old ages. These tails are the consequence of the uniform cut in absolute magnitude performed for our sample selection for all clusters. A cut at M_G < 4.5 makes us include a larger proportion of stars in the main sequence for young clusters compared to old clusters (see for instance NGC 3532 in Fig. 4). As explained above, ages for young main-sequence stars are in general overestimated because of their broad PDFs. This overestimate is larger when taking a simple median of the stars per cluster, instead of combining up the individual PDFs (violin plots), where there is a clear peak of the PDF distribution very close to the literature values. For clusters older than ~500 Myr these tails are smaller because there is a larger proportion of stars in the turn-off, which dominate the cluster’s PDF. For the oldest case, the globular cluster NGC 6397, we find that the peak of our distribution is in nice agreement with literature estimates, which are very consistent among them, though a clear gradient of age across the SGB is seen for this cluster in Fig. 4.

For intermediate-age clusters we tend to obtain relatively broad distributions, particularly for NGC 188 (for which literature ages are also quite diverse), but also for NGC 2682 and NGC 6819. These three clusters are also more distant than the sample of field stars analysed in Sect. 5, which is essentially limited to 1 kpc. Thus, we expect larger uncertainties in photometry, manifested in the wider main sequences, which unavoidably give broader age distributions.

Particularly for these three clusters, we find a large proportion of turn-off stars with double peaks in the PDFs (see for instance Fig. A.1). We have developed an algorithm to automatically detect the cases where the PDF shows multi-peaks (see Sect. A.1 for a detailed explanation), which we have run on the full sample of cluster stars. We have found that multi-peak solutions are in general found for stars placed near the turn-off loop. This feature causes the stellar tracks to overlap particularly for tracks of mass >1M_⊙ at solar metallicities. In general, a filter in the overall age uncertainty removes a large fraction of multi-peak stars, but those for which the two peaks are relatively close are not removed.

Overall, we find a good agreement of the cluster ages results’ from SPInS when compared to literature results. When we look at the comparison over the full age range in Fig. 6, we see a small tendency from our study to overestimate the ages with respect to literature, particularly for clusters younger than 1 Gyr, with a mean deviation of 200–500 Myr, depending on whether we use the mode or the median of the distribution for each cluster.

Fig. 4 SPInS results for three of the analysed clusters: NGC 2516 (0.2 Gyr), NGC 3532 (0.4 Gyr), NGC 6819 (2.2 Gyr) and NGC 6397 (12.6 Gyr). Top: intrinsic CMDs of the member stars (grey) and the selection of stars analysed by SPInS coloured by the obtained age (from the median value of the PDFs). Bottom: age histograms of all analysed stars in grey; and a filtered subsample is shown in blue, selected as: relative uncertainties in age better than 20% (and better than 500 Myr), excluding BSs, and excluding stars in the equal-mass binary sequence. Literature ages are marked with vertical dashed orange lines.

Fig. 5 Comparison between the age quoted in the literature from studies indicated in Table 1 (diamond symbols), and the age distribution per cluster obtained with SPInS in the form of violin plots. The horizontal black line depicts the mode of the age distribution. We also plot with filled circles and error bars the age estimate coming from the median age per star instead of the full PDF.

5 Results: Field stars

We use SPInS to determine the masses and ages of the selected MSTO and SGB stars with the configuration as explained in Sect. 2. We use as observational constraints: absolute magnitudes from Gaia DR3 (M_G,(G_BP − G_RP)₀), and metallicity from LAMOST. Uncertainties in these parameters are assumed to follow Gaussian distributions.

As for the case of the globular cluster in Sect. 4.1, we re-scale all the metallicities using the [α/Fe] provided by LAMOST to mimic the α element enrichment, to be able to use solar-scaled evolutionary tracks (see Sect. 4.2). For the MRS we compute two sets of ages, one using [α/Fe] and [Fe/H] from the main LAMOST pipeline (LASP), and another one using the abundances from LAMOST DR8 obtained using the label-transfer method based on a CNN, trained on APOGEE. The uncertainty propagation of the scaled metallicities using the Salaris et al. (1993) formula is not straightforward because it is likely that the uncertainties in the [α/Fe] and [Fe/H] of LAMOST are correlated. Since LAMOST does not provide a specific correlation matrix between the two values, here we have propagated the uncertainties analytically as if the two quantities were uncorrelated, i.e. the quadratic sum of the [α/Fe] and [Fe/H] multiplied by the partial derivatives of the Salaris et al. (1993) formula. The partial derivative with respect to [Fe/H] is always 1, and the one with respect to [α/Fe] is always smaller than 1 taking into account the values of [α/Fe] in our sample. We have decided to assume a value of exactly 1 for the latest derivative in order to slightly inflate the final errors, and possibly mitigate the fact that we are assuming independent measurements for [α/Fe] and [Fe/H].

We provide age estimates for 35 096 and 243 768 stars based on the metallicities and α element abundances of the LAMOST DR8 MRS and LRS samples, respectively. For the case of the CNN values of MRS, the sample includes 34 779 stars with valid abundances from CNN.

In Fig. 7, we show the comparison of the ages derived for the MRS-LASP and the LRS for the stars in common, which shows an overall good agreement. We notice that, even though the peak of the distribution stays in the 1:1 line, the overall shape is slightly asymmetric, in the sense that LRS ages seem slightly underestimated with respect to MRS. We have investigated the possible causes of this and we find that there is a small systematic overestimate in the metallicity of the LRS sample with respect to the MRS one (~0.06 dex in mean), which could explain the small difference in age.

The comparison between MRS-LASP and MRS-CNN also shows a slight tendency of ages derived from CNN to be underestimated with respect to LASP (mean value of −271 Myr), particularly for the metal-poor stars. This difference can be correlated with the difference between the CNN and LASP metallicites (see Fig. A.3), and with a difference in the estimation of the atmospheric parameters from the two methods. As explained in Sect. 3, for our age determination, we choose not to use spectroscopically derived atmospheric parameters, but only de-reddened photometry and spectroscopic metallicity. However, SPInS provides posteriors in all parameters of the input grid, including T_eff, which we compare to spectroscopically derived ones in Fig. A.4. The figure shows that LASP estimates seem to be more coherent with SPInS determinations than CNN ones, which have an overall offset of 124 K and a trend towards hot stars.

Fig. 6 Comparison between the age quoted in the literature from studies indicated in Table 1, and the determination from SPInS using the median values of the filtered stars (orange), and the mode (blue).

Fig. 7 Comparison of the ages from the MRS-LASP and LRS samples obtained by SPInS. The number of stars is indicated in the plot, and median quoted errors per age bin are plotted. The top plot is coloured by density, and the bottom plot by the difference in metallicity.

5.1 Comparison with large age catalogues

We compare the resulting median age per star obtained by SPInS of the MRS and LRS samples with ages coming from recent literature catalogues for stars in common in Fig. 8.

The catalogue of Xiang & Rix (2022) contains 247k stars, and they used a similar method as in this study to derive ages and masses from LAMOST DR7 using Yonsei-Yale stellar isochrones (Demarque et al. 2004). The main difference between their approach and ours is the fact they use as inputs spectroscopic estimates of the absolute magnitude in the K band, instead of Gaia photometry, and estimates of T_eff inferred by the data-driven method run in LAMOST DR7 (DD-Payne). The catalogue of Queiroz et al. (2023) using StarHorse for LAMOST DR7 MRS/LRS sub-giant MSTO stars, and contains age estimates as well as other stellar astrophysical parameters for 120 k and 1.3 M stars, respectively. StarHorse is a Bayesian isochrone-fitting code that uses as many input observables as available (including parallaxes, spectroscopic stellar parameters, and multi-band photometry) to compute the likelihood of the observed quantities for a grid of PARSEC 1.2S + COLIBRI stellar models (Marigo et al. 2017). Aside from the canonical priors, like the IMF, StarHorse uses a prior on the interstellar extinction based on a 3D extinction map as well as generous space density, age, and metallicity priors for the Galactic discs, bulge, halo, GC system, and nearby dwarf galaxies. The catalogue obtained by Kordopatis et al. (2023) contains 5 million stars with stellar ages computed in four different ways, depending on different combinations of parameters to project on the PARSEC 1.2S + COLIBRI isochrones. Here we used their final ‘optimal’ ages, which they provide as a combination of the four different projections (see their Sect. 2.6.2). Following the recommendations in Kordopatis et al. (2023) we filter stars whose estimated relative age uncertainty is larger than 50%. Finally, the recent catalogue by Nataf et al. (2024) provides ages for 289 k stars obtained from multi-band photometry with 3D extinction maps and without any information of the metallicity.

We find a general good agreement in the ages from SPInS with the aforementioned catalogues from the literature. The comparison with Xiang & Rix (2022) shows a marked difference for stars between ~2–4 Gyrs where our results are in general 500 Myr younger than those of Xiang & Rix (2022). These are in general hot stars (T_eff > 6000) for which we noticed that the effective temperatures derived with LAMOST (DR8) LASP pipeline and the data-driven approaches DD-Payne/CNN differ systematically, in the sense that DD-Payne always finds a cooler T_eff for hot stars. This bias could tentatively explain the observed difference, since a cooler T_eff for a given absolute magnitude and metallicity will generally provide older ages. However, this is an indirect explanation because SPInS ages do not use LASP T_eff but only information from derredenned photometry. It is also possible that this difference comes from the two different sets of stellar evolutionary models used. Overall though, we obtain a mean difference and scatter of −0.06/1.12 Gyr with Xiang & Rix (2022) for the LRS. The comparison with Queiroz et al. (2023) is the one that provides the most stars in common in the two samples and shows a good agreement, with a mean difference and scatter of −0.28/1.87 Gyr. In particular, it does not show the previously mentioned artefact for young ages. The comparison with Kordopatis et al. (2023) is more scattered and shows horizontal stripes, probably due to the mild pixelisation of the GSP-Spec results (Recio-Blanco et al. 2023), or alternatively to a different type of interpolation along stellar tracks. The mean difference and scatter are −0.55/2.63 Gyr, though the scatter decreases to 2.2 Gyr if we only keep stars with GSPSpec flags equal to 0 (except for noise, allowed to be 0 or 1). The larger scatter with this catalogue is well explained by the difference in the metallicity between the two catalogues (see Fig. A.5), highlighting the importance of the quality of the metallicity in the age estimates. The comparison with Nataf et al. (2024) shows a very good agreement with their ages, with only faint structures outside the 1:1 line for the case of the youngest stars, which are visible in the LRS sample. We highlight, however, that their estimations of the metallicity have a significant offset of ~0.15 dex with respect to LAMOST determinations (see Fig. A.5), as they already cite in their original paper. Surprisingly, this does not seem to have a big effect in the comparison with their ages, unlike the case of Kordopatis et al. (2023). Xiang & Rix (2022) and Queiroz et al. (2023) metallicities are very similar to ours, as seen in Fig. A.5, because they come from LAMOST even though previous data releases.

In general, we find very little systematics between the literature datasets and the results of SPInS, except for the youngest stars present in the Xiang & Rix (2022) dataset. We expect then similar offsets in these literature catalogues as those found in the present study using clusters (Sect. 4.2, Table 1).

Fig. 8 Comparison of the ages from the MRS (top) and LRS (bottom) samples obtained by SPInS with literature estimates from Xiang & Rix (2022), StarHorse obtained in Queiroz et al. (2023), Kordopatis et al. (2023) and Nataf et al. (2024) (from left to right). The number of stars is indicated in each plot, and median quoted errors per age bin are plotted.

5.2 Kinematics

We have used the 6D phase space information provided by Gaia DR3 with the radial velocities provided by the LAMOST DR8 catalogues (MRS and LRS) to compute action-angle variables and orbital parameters of the two samples with the galpy package (Bovy 2015). Galactocentric cartesian positions (X, Y, Z) and cylindrical positions and velocities (R, ϕ, Z, V_R, V_ϕ, V_Z) were computed assuming the solar values: (X, Y, Z)_⊙ = (−8.34, 0, 0.027) kpc from the center of the Galaxy, an azimuthal velocity at the solar radius of 240 kms⁻¹, and the velocity of the Sun with respect to the local standard of rest (LSR) as (−11.1, 12.24, 7.25) kms⁻¹ (Schönrich et al. 2010). Uncertainties were computed using the quoted individual uncertainties and covariance matrix with the galpy utils package.

Action-angle variables and orbital parameters were obtained using the axisymmetric potential MWPotential2014 implemented in galpy. The integration of the orbits was done up to 1 Gyr with 1000 steps, to obtain eccentricity, apocentre, pericentre, maximum distance from the Galactic plane Z_max, and guiding radius R_g. The obtained kinematical parameters will be used in the next subsections.

6 Discussion: Age–abundance–kinematics relations

We have used our age estimates for the MRS and LRS samples to investigate the age–[Fe/H] and age–[α/Fe] relation in the solar neighbourhood. We use the orbital parameters computed in Sect. 5.2 to dissect the age-abundance relations.

Following the validation tests performed on open and globular clusters in Sect. 4, we have applied a filter to the age catalogue only keeping stars with δAge/Age < 0.2 and δAge < 500 Myr. These filters reduce the size of the two samples to 20822 stars in the MRS (19 911 for the MRS-CNN) and 132 211 in the LRS. They particularly remove a large fraction of stars at ages < 1 Gyr, and a group of stars with solar metallicities with very old inferred ages (most probably low main-sequence stars). We also filter the stars that give multiple solutions in SPInS, detected as explained in Sect. A.1, which remove 3589 and 24 490 additional stars in the MRS and LRS samples, respectively.

Fig. 9 Age-metallicity relation obtained from the MRS (left using LASP abundances, and middle using CNN abundances), and LRS (right) samples. The colour corresponds to the normalised density of stars at each age bin.

Fig. 10 Age-[Fe/H] relation obtained with the LRS sample, isolating the [α/Fe]-enhanced thin (left with 16 800 stars) and thick (right with 7623 stars) discs. The colour corresponds to the normalised density of stars at each age bin (only bins with more than 1 star are shown).

6.1 Age-metallicity relation

We plot the age–[Fe/H] relation (AMR) of our samples in Fig. 9 coloured by the stellar density. The histograms have been normalised for each age bin by its maximum value. The plots for the MRS and the LRS show a hint of two tight sequences that are consistent overall with the already known bimodality of the AMR (Nissen et al. 2020; Jofré 2021; Xiang & Rix 2022; Anders et al. 2023) (see also discussions in Sahlholdt et al. 2022 and Cerqui et al., in preparation). For the LRS, we see a horizontal stripe at [Fe/H] between 0 and −0.1, which causes an underdensity of stars in this region. We have not found the cause of this, but it is most probably an artefact or a bias in the LAMOST LRS abundance pipeline. For the MRS the overall picture is similar when we use the ages obtained using the CNN (middle panel) and LASP abundances (left panel). The thin disc sequence seems more prominent when using LASP abundances.

In Fig. 10, we dissect the AMR of the LRS using simultaneous cuts in angular momentum and [α/Fe] abundances, as done in Xiang & Rix (2022). A low-[α/Fe] branch corresponding to the chemical thin disc is clearly visible up to ~6 Gyr for stars with L_z ≥ 1.1⁸ and [α/Fe] ≤ 0.1. The sequence seems to extend up to 9 Gyr but shows a tail of stars at higher metallicities for ages older than 6 Gyr, which possibly corresponds to a population of radial migrators. On the other hand, selecting stars with L_z ≤ 0.85 and [α/Fe] ≥ 0.2 yields a tight sequence starting at ~9 Gyr corresponding to the [α/Fe]-enhanced disc, or thick disc, which has a steeper slope in the AMR. Similar selections also isolate the two sequences in the MRS sample but with a much smaller number of stars. Between 6 and 9 Gyr, particularly when selecting the intermediate cuts with L_z and [α/Fe], we observe a mixture of the two populations, without a clear split of the two sequences.

The age distribution of the low-α disc shows two distinct peaks at 3 and 6 Gyr with a lower density ‘gap’ around 4 Gyr. We have checked that the gap is still present when we use the sample without filtering the multiple solution stars, to ensure that it is not an artifact of removing stars in a particular region of the CMD. Similar gaps in the young disc are observed in other recent samples of ages, but they are not always located at the same age (e.g. Sahlholdt et al. 2022; Anders et al. 2023; Queiroz et al. 2023). There is a clear relation between age and metallicity in the young peak, and then a more blurred dependence is seen in the second peak after the gap, with a hint of a possible change in the slope. The multimodal age distribution coupled with the tight metallicity dependence also imprints two very close peaks in the metallicity distribution, as seen in the projected histogram. We do not find in our AMR the V shape found by Xiang & Rix (2022), who showed an additional branch with a positive slope. For the high-α disc we obtained a more restricted age distribution with a unique peak at around 11 Gyr, and a long asymmetric tail towards younger ages, which could be due to young α-rich stars (e.g. Cerqui et al. 2023). To better understand these dependences an additional and careful analysis of the selection function is needed.

We plot in Fig. 11 a 2D histogram of the AMR coloured by the mean guiding radius (R_q) per bin for the full LRS and MRSLASP samples (no filter in L_z or [α/Fe]). The non-normalised AMR better highlights the few stars in our sample with metallicities lower than −1, which form a vertical structure particularly visible in the LRS, without any clear dependence in age, and which present small R_g. We also see a dependence of the mean guiding radius along the high-α disc sequence. Overall, the dependence of the AMR with guiding radius is consistent with the general understanding of the Galactic radial metallicity profile coupled with radial mixing (e.g. Haywood 2008; Schönrich & Binney 2009; Minchev et al. 2014; Hayden et al. 2015, 2020; Frankel et al. 2018; Haywood et al. 2024). Even though our sample is made of stars in the very local solar neighbourhood, we can see the effect of radial migration, which, as time goes by, brings stars born at different galactocentric radii to the solar radius. Thus, at a given age we see a clear gradient in guiding radius for the young (<10 Gyr) AMR region, where low metallicity stars have a larger guiding radius than high metallicity stars.

Fig. 11 Non-normalised AMR for the LRS (top) and MRS-LASP (bottom) coloured by the mean guiding radius per bin.

6.2 [α/Fe]–age relation

We plot the relation between [α/Fe] and age in Fig. 12 for the MRS. In the left panels, the colour represents the stellar density, and in the middle and right panels, we plot a 2D histogram weighted by the mean [Fe/H] value per bin and the mean R_g per bin, respectively. We show the relation when using the [Fe/H] and [α/Fe] from LASP abundances in the top plots, and in the bottom plots using the CNN values (notice the change of vertical scale). The usage of the two sets of abundances shows significant differences in the [α/Fe]–age relations. In the case of LASP, we see a slow increase in the [α/Fe]-age relation for stars younger than ~9–10 Gyr. A high-α sequence with a much higher slope starts to be visible at 10 Gyr, which is usually linked to the chemically defined thick disc. In the case of CNN abundances, the thin disc shows a much flatter relation up to ~9–10 Gyr, with a much steeper increase in [α/Fe] for the thick disc sequence.

We notice a blob of stars with old ages (>12 Gyr) that seem to have constant and even decreasing [α/Fe], and that give the impression of a decreasing [α/Fe] versus age sequence for the oldest stars in the LASP case. In the case of CNN, this blob stays at roughly constant values of [α/Fe], as the end of the thick disc. This feature corresponds essentially to a group of stars at relatively low metallicities ([Fe/H] < −0.6) that, in the LASP case, seem to have decreasing [α/Fe] abundances in the [Fe/H] versus [α/Fe] plane. This feature is also seen in other age samples using LAMOST abundances such as Queiroz et al. (2023). In fact, in this region, there are few stars, but normalising the histogram at each age bin enhances this feature. We have not found a clear reason for this, it is possible that the values of the overall α abundances for low metallicity stars in LAMOST suffer from some biases due to the intrinsic difficulty of deriving abundances at low metallicities.

In the middle panels of Fig. 12, we see the imprint of the AMR (Fig. 9) as a colour gradient along the young chemically defined thin disc sequence. In the case of the CNN abundances, we see a clearer stratification of the metallicities across the thin disc. The plot when using the LASP abundances shows a significant group of intermediate-age stars (7–10 Gyr), with high [Fe/H] (>0.1) and low-[α/Fe] abundances. In the right-hand plot, we see how these stars tend to have slightly smaller guiding radii with respect to the rest of the stars at the same age. With a mean R_g ~ 6.7 kpc, most of these stars are most likely to be radial migrators coming from the inner disc. This group of stars is not clearly identified when using the CNN abundances. However, in the case of the CNN abundances, we see a clearer stratification of the thin disc with guiding radius (bottom right plot), showing the effects of radial migration to the [α/Fe] versus age sequence.

Overall, there are remarkable differences between the LASP and the CNN values in the [α/Fe]-age relation, which can lead to different conclusions on the enrichment history of the thick and thin discs. The difference comes mainly from the [α/Fe] values, which have a very poor correlation among the two different determinations. This is a clear example of the importance of the precision and accuracy of the atmospheric parameters and chemical abundances obtained from spectroscopic surveys across the full metallicity range.

6.3 Age–velocity dispersion relation

The relation between the dispersion in the vertical velocity v_z and age in the solar neighbourhood is often called the age-velocity relation (AVR) and has been studied for decades using different tracers (Wielen 1985; Nordström et al. 2004; Casagrande et al. 2011). It is generally described as a power law, though its slope is debated.

We have computed the dispersion in v_z in ~25 bins of age for the filtered MRS and LRS. We have additionally applied a cut on the uncertainty of v_z of <2 kms⁻¹ to avoid outliers. The width of the bins is around 500 Myr, but were slightly tuned for the MRS for old ages to avoid having bins with very few stars. We plot the age-velocity relation (AVR) obtained for the two samples in Fig. 13 coloured by the number of stars, relative to the bin with the maximum number. We have computed the uncertainties in σ_z sampling the uncertainties of the cylindrical velocities 1000 times. Obtained uncertainties are smaller than <0.1 kms⁻¹ in all bins, and are thus smaller than the point size.

We overplot the AVR obtained for a sample of open clusters in the solar neighbourhood obtained by Tarricq et al. (2021), which cover the age range of [0,2.5] Gyr. We overplot with a dotted blue line the power-law fit that they performed, with a slope value of β_OC = 0.19 ± 0.03.

The two youngest age bins (<2 Gyr) of both the MRS and the LRS seem to indicate a flattening of the σ_z at around 12 kms⁻¹ or even an increase in the dispersion, behaviour that is not reproduced by the open cluster sample. The statistics of these bins are smaller than the rest, particularly in the MRS, which count fewer than 100 stars. Our age determination in this range of ages gives slightly overestimated values, as seen for the youngest clusters in Fig. 5. Additionally, taking into account the validation tests done in Sect. 4, it is also possible that the youngest bins have contamination due to BSs, which would then have typical kinematics of older stars, causing an increase in σ_z. Overall, we consider that the two youngest bins are most probably not representative of the true AVR of the Galaxy.

As done by Anders et al. (2023) we perform a power-law fit to our two samples in the age range 2 < Age < 9 Gyr, in orange and green, for the MRS and LRS respectively. We limited the age range to avoid the two youngest age bins and the oldest regions where we see a departure from a simple power law. The uncertainties of the obtained coefficients are extremely small, of the order of 10⁻⁴. We see how the open cluster sample nicely matches the power law fitted to the field stars, and complements it in the younger age range. The oldest point of the open cluster sample, though still compatible with our trend, has a remarkably lower value. This decrease is probably related to the low statistics of the open cluster populations at ages >2 Gyr (reflected by the large uncertainty on the dispersion), which has been attributed to a higher destruction rate of old open clusters. As a direct consequence, the obtained slope with the open clusters is smaller compared to the field stars.

On the old end, a departure from the power-law fit is highlighted by a steepening of the σ_z–age relation at 9–10 Gyr. This has been found in the literature and has been attributed to disc stars that have been kinematically heated due to the Gaia-Sausage-Enceladus merger event (Di Matteo et al. 2019; Belokurov et al. 2020).

Fig. 12 [α/Fe] vs age relation obtained for the MRS sample using LASP abundances (top), and CNN abundances (bottom). In the left panel, the colour corresponds to the normalised density of stars per age bin, in the middle and right panels we colour each bin by the mean value of [Fe/H] and Rg per bin values, respectively.

Fig. 13 Age-velocity relation for the MRS and LRS samples coloured by the relative number of stars in each point (divided by the maximum). Power law fitting the 2 < Age < 9 Gyr range are plotted in orange and green dashed lines. A comparison with the AVR obtained for clusters in the solar neighbourhood (Tarricq et al. 2021) is shown in blue dots and their power-law fits in a blue dashed line.

7 Conclusions

Reliable stellar ages, coupled with 6D phase space information and chemical abundances, are essential to understand the evolution of the Milky Way. In this paper, we have exploited the possibility of using the absolute CMD from Gaia DR3 photometry only, coupled with spectroscopic estimates of the [Fe/H] and [α/Fe] from the LAMOST DR8 medium- and low-resolution samples.

We have used the public SPInS code (Lebreton & Reese 2020) to obtain age estimates and reliable uncertainties of individual MSTO stars and subgiant stars. We implemented new tools in the code, which are now public, to obtain an automatic evaluation of the convergence of the MCMC sampler (see Sect. 2.2), and an option to automatically detect multiple peaks in the case of a multi-modal posterior distribution (see Sect. A.1).

We tested our strategy on a sample of 4374 stars in 14 stellar clusters, including the old metal-poor globular cluster NGC 6397. With this sample, we were able to investigate the possible systematic and random uncertainties that can affect the age distribution of a sample of field stars selected with the same criteria and analysed with the same code. We evaluated the effect of the presence of main-sequence stars in the sample, which biases the cluster ages by overestimating them. These stars can mostly be removed by filtering our uncertainties to 20% and with a hard cut at 500 Myr for the oldest clusters. We also evaluated the effect of the presence of non-single stars, particularly BSs and unresolved binaries, which, unsurprisingly, underestimate and overestimate the cluster ages, respectively. Statistically, in a sample of field stars selected in the same way, the non-single stars would represent around 11% of the sample.

We see how the combination of the posterior distribution of the selected cluster stars yields age distributions with prominent peaks at the cluster’s quoted literature ages, which include ages obtained from isochrone fitting, asteroseismology, eclipsing binaries, and the white dwarf cooling sequences. We applied the same strategy to an exquisite sample of field stars with very small photometric uncertainties (in general smaller than 0.04 in M_G, and of the order of 0.01 in G_BP − G_RP). This sample has radial velocities and [Fe/H] and [α/Fe] abundances that come from the LAMOST DR8 MRS sample (35 096 stars that have LASP abundances and 34 779 stars that have CNN abundances) and/or the LRS sample (243 768 stars). By combining the obtained ages with the 6D phase space and the LAMOST chemical abundances, we were able to investigate the following points:

• The age-metallicity relation (AMR) with the LRS and MRS samples, obtaining a consistent picture with its already-known bimodality;

• We dissected the AMR from the LRS with simple cuts in angular momentum and [α/Fe] abundances, and we were able to isolate a low-[α/Fe] branch visible at young ages up to 6–9 Gyr, and a high-[α/Fe] branch with a steeper slope in the AMR that covers ages older than 9 Gyr;

• While the high-[α/Fe] branch shows a unique peak in the age distribution around 11 Gyr, the low-[α/Fe] branch appears to be much more complex, with at least two peaks at 3 and 6 Gyr. A better understanding of the selection effects is needed to correctly interpret this distribution;

• Using the MRS, we studied the [α/Fe] versus age relation, which shows very different pictures when using LASP or CNN abundances. This highlights the importance of the precision and accuracy of chemical abundances from spectroscopic surveys. We investigated the features in the [α/Fe]=-age relation, coupling it with the [Fe/H] and the guiding radius;

• The LASP [α/Fe] abundances show a slowly increasing trend for stars younger than 10 Gyr and a high-α sequence that starts to be visible at 10 Gyr. We detect a population of radial migrators from the inner disc that is present in this sample for ages between 7–10 Gyr;

• The CNN abundances show a much flatter slope for the thin disc stars, with a steeper thick disc sequence starting at 10 Gyr. A clear stratification of [Fe/H] and R_g is visible across the [α/Fe] range for the thin disc;

• We investigated the AVR in 25 age bins, showing a consistent power law relation in the range of 2–9 Gyr with slopes of 0.44 and 0.55, respectively, for the MRS and LRS samples;

• The simple power law breaks down in our sample for ages lower than 2 Gyr, though the two younger age bins contain fewer stars. Additionally, from the results of our tests using clusters, our method tends to overestimate the ages of younger stars, making this age range unreliable;

• Without these two youngest bins, the obtained AVR is compatible with that obtained with a sample of open clusters previously analysed in the literature;

• The AVR also breaks down at ages older than 9 Gyr, significantly increasing the vertical velocity dispersion; this is possibly related to the Gaia-Sausage-Enceladus merger event.

Data availability

Full catalog is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/692/A243

Acknowledgements

We warmly thank Rosine Lallement for her help and advice when using the extinction maps, and Santi Casissi for his advice when using the BaSTI stellar evolutionary tracks. This work was partially supported by the Spanish MICIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe” by the “European Union” through grant PID2021-122842OB-C21, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. FA acknowledges the grant RYC2021-031683-I funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR. This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/web/gaia/dpac/consortium). We acknowledge the Gaia Project Scientist Support Team and the Gaia DPAC. Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences. This research made extensive use of the SIMBAD database, and the VizieR catalogue access tool operated at the CDS, Strasbourg, France, and of NASA Astrophysics Data System Bibliographic Services.

Appendix A Additional tests and figures

A.1 Detecting multiple solutions in SPInS

In the following section, we explain the details of the algorithm that detects multiple solutions in SPInS output, which has now been implemented in the public SPInS version. This algorithm can be activated by the user, if needed, using the configuration file.

Implemented algorithm

The algorithm’s primary objective is to identify and isolate distinct solutions within multimodal distributions in stellar parameter inference, especially for datasets produced using the SPInS and AIMS (Reese 2016) programs. Hierarchical Density-Based Spatial Clustering of Applications with Noise (HDBSCAN Campello et al. 2013) was chosen for its ability to find groups of varied densities without requiring a fixed number of groups, making it very versatile and appropriate for our requirements. Input data for HDBSCAN in the case of SPInS include the model parameters of the MCMC samples, namely the initial mass log₁₀(M/M_⊙), metallicity [M/H], and Age (in Myrs). These parameters are essential for the grouping process, while the log-likelihood values (lnP) are used after the grouping process to identify the optimal parameters for each group. The data is standardised using StandardScaler, which removes the mean from the MCMC samples and scales them to unit variance.

HDBSCAN requires two primary parameters⁹ to be set: ‘min_cluster_size’ and ‘min_samples’, as well as several hyperparameters such as ‘cluster_selection_method’, a distance threshold ‘cluster_selection_epsilon’, and the ‘allow_single_cluster’ parameter that allows the detection of a single cluster. After testing the algorithm in the specific case of SPInS results, we have chosen to use the ‘cluster_selection_method’=’leaf’, ‘cluster_selection_epsilon’=0.3, ‘allow_single_cluster’=True, and ‘min_cluster_size’ = ‘min_samples’. We have seen that the detection worked efficiently when setting ‘min_samples’ to 1% or 2% of the total number of samples, with a small dependence on the particular case. Therefore the algorithm performs two runs with HDBSCAN using, respectively, the two values of ‘min_samples’, thus yielding two sets of labels that specify to which group each MCMC sample belongs. The algorithm then compares the number of groups created by each run and chooses the one that produces more groups. If both runs generate the same non-zero number of groups, the sizes of the smallest groups are compared, and the run with the most samples in the smallest group is selected.

HDBSCAN labels points as noise if they do not fit well into the identified groups, but at the same time, these points can include useful data points with significant log-likelihood values (ln P). To guarantee that the best log-likelihood values are correctly identified, we reassigned such noise points to the relevant groups. To achieve this, first, noise points (labelled -1) are identified. If there are more than two unique groups (considering all excluding the noise points), the algorithm proceeds as follows: i) the group centres are determined by calculating the mean of the data points within each group; ii) the distances between each noise point and all group centres are determined; iii) if the distance to the nearest group centre is less than the minimum intergroup distance¹⁰, the noise points are reassigned to that group; iv) if there are fewer than two unique groups, all points, including noise points, are assigned to a single group. This method ensures that significant data points, which may have been misclassified as noise due to their distance from dense regions, are accurately reassigned to appropriate groups.

Fig. A.1 Triangle plot of the resulting PDF obtained by SPInS of a turn-off star of NGC 188, projected in the 3D space: log(M0/M⊙), [M/H], Age. The two groups of solutions identified by our algorithm are shown in red and magenta, whereas background noise is plotted in green.

Results obtained with clusters

We detail here the results of applying the previous algorithm to the results of the open and globular clusters. We have visually checked the performance of the algorithm for this smaller sample of stars (compared to the full sample of field stars Sect. 3.3), obtaining overall satisfactory results in the identification of multiple peaks.

From the total sample of 4374 cluster stars (Sect. 4.1) we find 813 stars with a multi-modal solution, such as the one depicted in Fig. A.1 for a turn-off star of NGC 188. Of these, most of them have two solutions, while 42 stars have three solutions, and only 13 stars have more than three solutions.

We plot in Fig. A.2 the CMD of the stars with multiple solutions in the three intermediate-age clusters, for which we find the largest proportion of multi-modal distributions. In the top panels, the colour represents the difference between the age of the most probable PDF peak and the reference age taken as a mean of the literature ages listed in Table 1. In the bottom panels, we plot the ratio of the probabilities of the most probable peak versus the secondary peak. BaSTI stellar tracks of around solar metallicity (similar to the cluster’s metallicities) are overplotted to guide the eye.

Fig. A.2 CMD of the three intermediate-age clusters for which we find the most multiple solutions. Stars for which multiple peaks have been identified are plotted using coloured symbols according to the difference between the age of the most probable age peak and the mean literature age (top panels), and according to the probability ratio of the youngest to the oldest peak (bottom panels). Grey symbols depict stars with a single peak. BaSTI stellar tracks of metallicity −0.08 are shown using orange lines to guide the eye.

From Fig. A.2 we can see that the stars with multi-modal distributions are in general located along the turn-off loop, which is present around solar metallicities for initial masses larger than ~1M_⊙. For smaller masses, the loop is not present any more and thus we do not find stars with multiple peaks in the low main sequence. We also see clear gradients in the age difference of the main peak with respect to the cluster’s real age (top panel, particularly for NGC 6819). Finally, we see how the ratios of the probability of the youngest versus older peak gets more important in certain regions around the loop, particularly in the leftmost side of the MSTO for NGC 6819.

A.2 Additional figures

We show in Fig. A.3 the comparison between SPInS ages for the MRS obtained from LASP and CNN abundances, coloured by the mean difference in the metallicities. A small trend towards metal-poor stars is present (mean value of -271 Myr), which correlates with the differences in the derived metallicities.

Figure A.4 shows the T_eff comparison between the spectroscopic determination provided for the MRS by LASP/CNN, and the median value of the T_eff posterior distribution obtained by SPInS. The obtained mean offsets and scatters are 23/166 K for LASP and 124/149 K for CNN, highlighting a discrepancy between the atmospheric parameters obtained by CNN with respect to the photometrically driven ones obtained by SPInS.

We show in Fig. A.5 the comparison between literature and SPInS ages for the MRS and LRS coloured by the mean difference in the assumed metallicity between the two studies. A running median is shown in red dashed lines. See Fig. 8 for the same plots coloured by stellar density. We see that the largest discrepancies between the derived ages are usually explained by the difference in the assumed metallicity, except for the catalogue of Nataf et al. (2024), for which an overall offset of around 0.15 dex is obtained and no difference in age is present.

Fig. A.3 Age difference between LASP and CNN ages coloured by mean metallicity difference. Running median and contour lines are plotted in red.

Fig. A.4 Spectroscopic determination of Teff provided for the MRS by LASP (top) CNN (bottom), compared with the values obtained from the posterior distributions of SPInS.

Fig. A.5 Comparison of the ages obtained for the MRS (top) and LRS (bottom) sample from SPInS with literature estimations estimates from Xiang & Rix (2022), from StarHorse obtained in Queiroz et al. (2023), Kordopatis et al. (2023) and Nataf et al. (2024) (from left to right), coloured by mean difference in the metallicity from the two catalogues. The number of stars in common is indicated in each panel. A running median is shown in red dashed lines.

References

All Tables

All Figures

	Fig. 1 Distribution of the initial sample of 352 k stars in the LAMOST DR8 LRS sample obtained from the procedure indicated in Sect. 3. Top: galactic X–Y and X–Z distributions coloured by stellar density. The Sun is at (0, 0). Bottom: intrinsic CMD coloured by density (left), and binned with colours representing mean [Fe/H] (right).
In the text

Fig. 2 Intrinsic CMD of the initial sample of LRS (top) and MRS (bottom) stars, divided in four [Fe/H] bins, coloured by the density of stars (logarithmic scale). The number of stars is plotted in each panel. The grey lines represent the cuts done to select MSTO and SGB, as is explained in Sect. 3.3. To guide the eye, we additionally plot in blue two isochrones representative of each [Fe/H] bin at two different ages.

In the text

	Fig. 3 Distribution of uncertainties in absolute magnitude and colour for the selection of MSTO and SG stars in the MRS (bottom) and LRS (top) samples.
In the text

Fig. 4 SPInS results for three of the analysed clusters: NGC 2516 (0.2 Gyr), NGC 3532 (0.4 Gyr), NGC 6819 (2.2 Gyr) and NGC 6397 (12.6 Gyr). Top: intrinsic CMDs of the member stars (grey) and the selection of stars analysed by SPInS coloured by the obtained age (from the median value of the PDFs). Bottom: age histograms of all analysed stars in grey; and a filtered subsample is shown in blue, selected as: relative uncertainties in age better than 20% (and better than 500 Myr), excluding BSs, and excluding stars in the equal-mass binary sequence. Literature ages are marked with vertical dashed orange lines.

In the text

	Fig. 5 Comparison between the age quoted in the literature from studies indicated in Table 1 (diamond symbols), and the age distribution per cluster obtained with SPInS in the form of violin plots. The horizontal black line depicts the mode of the age distribution. We also plot with filled circles and error bars the age estimate coming from the median age per star instead of the full PDF.
In the text

	Fig. 6 Comparison between the age quoted in the literature from studies indicated in Table 1, and the determination from SPInS using the median values of the filtered stars (orange), and the mode (blue).
In the text

	Fig. 7 Comparison of the ages from the MRS-LASP and LRS samples obtained by SPInS. The number of stars is indicated in the plot, and median quoted errors per age bin are plotted. The top plot is coloured by density, and the bottom plot by the difference in metallicity.
In the text

	Fig. 8 Comparison of the ages from the MRS (top) and LRS (bottom) samples obtained by SPInS with literature estimates from Xiang & Rix (2022), StarHorse obtained in Queiroz et al. (2023), Kordopatis et al. (2023) and Nataf et al. (2024) (from left to right). The number of stars is indicated in each plot, and median quoted errors per age bin are plotted.
In the text

	Fig. 9 Age-metallicity relation obtained from the MRS (left using LASP abundances, and middle using CNN abundances), and LRS (right) samples. The colour corresponds to the normalised density of stars at each age bin.
In the text

	Fig. 10 Age-[Fe/H] relation obtained with the LRS sample, isolating the [α/Fe]-enhanced thin (left with 16 800 stars) and thick (right with 7623 stars) discs. The colour corresponds to the normalised density of stars at each age bin (only bins with more than 1 star are shown).
In the text

	Fig. 11 Non-normalised AMR for the LRS (top) and MRS-LASP (bottom) coloured by the mean guiding radius per bin.
In the text

	Fig. 12 [α/Fe] vs age relation obtained for the MRS sample using LASP abundances (top), and CNN abundances (bottom). In the left panel, the colour corresponds to the normalised density of stars per age bin, in the middle and right panels we colour each bin by the mean value of [Fe/H] and Rg per bin values, respectively.
In the text

	Fig. 13 Age-velocity relation for the MRS and LRS samples coloured by the relative number of stars in each point (divided by the maximum). Power law fitting the 2 < Age < 9 Gyr range are plotted in orange and green dashed lines. A comparison with the AVR obtained for clusters in the solar neighbourhood (Tarricq et al. 2021) is shown in blue dots and their power-law fits in a blue dashed line.
In the text

	Fig. A.1 Triangle plot of the resulting PDF obtained by SPInS of a turn-off star of NGC 188, projected in the 3D space: log(M0/M⊙), [M/H], Age. The two groups of solutions identified by our algorithm are shown in red and magenta, whereas background noise is plotted in green.
In the text

Fig. A.2 CMD of the three intermediate-age clusters for which we find the most multiple solutions. Stars for which multiple peaks have been identified are plotted using coloured symbols according to the difference between the age of the most probable age peak and the mean literature age (top panels), and according to the probability ratio of the youngest to the oldest peak (bottom panels). Grey symbols depict stars with a single peak. BaSTI stellar tracks of metallicity −0.08 are shown using orange lines to guide the eye.

In the text

	Fig. A.3 Age difference between LASP and CNN ages coloured by mean metallicity difference. Running median and contour lines are plotted in red.
In the text

	Fig. A.4 Spectroscopic determination of Teff provided for the MRS by LASP (top) CNN (bottom), compared with the values obtained from the posterior distributions of SPInS.
In the text

Fig. A.5 Comparison of the ages obtained for the MRS (top) and LRS (bottom) sample from SPInS with literature estimations estimates from Xiang & Rix (2022), from StarHorse obtained in Queiroz et al. (2023), Kordopatis et al. (2023) and Nataf et al. (2024) (from left to right), coloured by mean difference in the metallicity from the two catalogues. The number of stars in common is indicated in each panel. A running median is shown in red dashed lines.

In the text