Ensemble age inversions for large spectroscopic surveys

TL;DR

This paper introduces an age inversion method that reconstructs the true stellar age distribution from probability density functions, enabling better understanding of galactic evolution from large spectroscopic surveys.

Contribution

The paper presents a novel age inversion technique that models the log(age) distribution as a weighted sum of mono-age populations, validated on simulated and real survey data.

Findings

Successfully recovers true age distributions from simulated data.

Effectively reproduces age-metallicity trends in large surveys.

Applicable to samples with over 1,000 stars.

Abstract

Galactic astrophysics is now in the process of building a multi-dimensional map of the Galaxy. For such a map, stellar ages are the essential ingredient. Ages are however measured only indirectly by comparing observational data with models. It is often difficult to provide a single age value for a given star, as several non-overlapping solutions are possible. We aim at recovering the underlying log(age) distribution from the measured log(age) probability density function for an arbitrary set of stars. We build an age inversion method, namely, we represent the measured log(age) probability density function as a weighted sum of probability density functions of mono-age populations. Weights in that sum give the underlying log(age) distribution. Mono-age populations are simulated so that the distribution of stars on the log g-[Fe/H] plane is close to that of the observed sample. We tested the age inversion method on simulated data, demonstrating that it is capable of properly recovering the true log(age) distribution for a large (N > 103) sample of stars. The method was further applied to large public spectroscopic surveys. For RAVE-on, LAMOST and APOGEE we also applied age inversion to mono-metallicity samples, successfully recovering age-metallicity trends present in higher-precision APOGEE data and chemical evolution models. We conclude that applying an age inversion method as presented in this work is necessary to recover the underlying age distribution of a large (N > 103 ) set of stars. These age distributions can be used to explore for instance age-metallicity relations.

Equations

$$N({\tau }_b)=\sum _{i:{\tau }_i\in ({\tau }_b-\Delta ,{\tau }_b+\Delta )}{w}_i,$$

$$N_{unc}(\tau )=\sum _i{w}_i\mathcal{N}(\tau |{\tau }_i,{\sigma }_{\tau ,i}),$$

$$N_{fit}(\tau )=\sum _i{w}_i{F}_i(\tau ,p_i),$$

$$C(\tau )={\int }_0^{\infty }P(\tau ,\tilde{\tau })N(\tilde{\tau })d\tilde{\tau }.$$

$$C_j=\sum _i{P}_{j,i}{N}_i.$$

$$L_0=-\sum _{j=1}^n(C_j-\sum _{i=1}^n{P}_{ji}{N}_i{)}^2.$$

$$L=-\sum _{j=1}^n(C_j-\sum _{i=1}^n{P}_{ji}{N}_i{)}^2-\lambda \sum _{i=1}^n{\left(\frac{d{N}_i}{d\tau }\right)}^2,$$

$$\frac{d{N}_i}{d\tau }=\mathcal{T}\text{@vec}N,$$

$$\mathcal{T}=\left(\begin{array}{cccccccc}-1.5& 2.& -0.5& 0& \cdots & \cdots & \cdots & 0\\ 0& -1.5& 2.& -0.5& 0& \cdots & \cdots & 0\\ 0& 0& -1.5& 2.& -0.5& 0& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \cdots & 0& 0& -1.5& 2.& -0.5& 0\\ 0& \cdots & \cdots & 0& 0& -1.5& 2.& -0.5\\ 0& \cdots & \cdots & 0& 0.5& -2& 1.5& 0\\ 0& \cdots & \cdots & \cdots & 0& 0.5& -2& 1.5\end{array}\right)$$

$$f^{\prime }(x_i)=-\frac{3}{2}f(x_i)+2f(x_{i+1})-\frac{1}{2}f(x_{i+2}),$$

$$f^{\prime }(x_i)=\frac{1}{2}f(x_{i-2})-2f(x_{i-1})+\frac{3}{2}f(x_i)$$

$$\text{@vec}{C}^{\prime }=P^{\prime }\text{@vec}N,$$

$$P^{\prime }$$

$$\text{@vec}{C}^{\prime }$$

$$\text{@vec}{C}_{q_0}^{\prime }=\text{@vec}{P}_{p_0}^{\prime }\text{@vec}N,$$

$$R_{m,n,p_0,q_0}={\Vert \begin{array}{c}\text{@vec}{C}_n-\text{@vec}{P}_m\text{@vec}{N}_{p_0,q_0}\end{array}\Vert }_2.$$

$$Q_{p_0,q_0}=\frac{1}{24}\left(\sum _p\sum _q{R}_{m,n,p_0,q_0}-R_{p_0,q_0,p_0,q_0}\right).$$

$$Q_{fit}(\log\lambda)=\left\{\begin{array}[]{l}a+c(\log\lambda-\log b_{1}),\textrm{if}\,\lambda\leq b_{1}\\ a+d(\log\lambda-\log b_{1}),\textrm{if}\,b_{1}<\lambda\leq b_{2}\\ a+d(\log b_{2}-\log b_{1}),\textrm{if}\,\lambda>b_{2}\end{array}\right..$$

$$\tau ([Fe/H])=\frac{{\int }_{{\tau }_{\min }}^{{\tau }_{\max }}\mathrm{C}(\tau ,[Fe/H])\mathrm{d}\tau }{{\tau }_{\max }-{\tau }_{\min }},$$

$$[Fe/H](\tau )=\frac{{\int }_{[Fe/H{]}_{\min }}^{[Fe/H{]}_{\max }}\mathrm{C}(\tau ,[Fe/H])\mathrm{d}[Fe/H]}{[Fe/H{]}_{\max }-[Fe/H{]}_{\min }}.$$

Full text

*Aims. *We aim at recovering the underlying log(age) distribution from the measured log(age) probability density function for an arbitrary set of stars.

*Methods. *We build an age inversion method, namely, we represent the measured log(age) probability density function as a weighted sum of probability density functions of mono-age populations. Weights in that sum give the underlying log(age) distribution. Mono-age populations are simulated so that the distribution of stars on the $\log g$-[Fe/H] plane is close to that of the observed sample.

*Results. *We tested the age inversion method on simulated data, demonstrating that it is capable of properly recovering the true log(age) distribution for a large ($N>10^{3}$) sample of stars. The method was further applied to large public spectroscopic surveys. For RAVE-on, LAMOST and APOGEE we also applied age inversion to mono-metallicity samples, successfully recovering age-metallicity trends present in higher-precision APOGEE data and chemical evolution models.

*Conclusions. *We conclude that applying an age inversion method as presented in this work is necessary to recover the underlying age distribution of a large ($N>10^{3}$) set of stars. These age distributions can be used to explore for instance age-metallicity relations.

1 Introduction

Stellar ages and distances are important ingredients to compose a reliable picture of the Galaxy. In combination with chemical abundances and kinematics they help us recover the history of star formation and satellite accretion for the Galaxy from its formation to present day. Large astrometric surveys, of which Gaia (Perryman et al., 2001) with its over $10^{9}$ objects is now the dominating one, provide us with kinematic information, such as positions, proper motions and parallaxes. To obtain radial velocities and physical properties of the stars, such as temperatures, surface gravities and chemical compositions, spectroscopic data can be used. Modern spectroscopic surveys provide such data for millions of stars.

While distances, chemical abundances and kinematics are typically contained in survey results, ages are not. This is because they are not related directly to any observable parameter and therefore estimating ages is more complicated and in most cases involves models of stellar evolution. Soderblom (2010) gives a large overview of different methods used to derive stellar ages. One of the methods listed there, the so-called “isochrone matching”, is useful to derive both ages and distances, and is widely applied to spectroscopic data. This method is based on comparison of quantities measured spectroscopically (like the effective temperature $T_{\rm{eff}}$, surface gravity $\log g$ and metallicity [Fe/H]) to a set of models. A subset of models with parameters close to the observed ones gives an estimate of the age and luminosity of the star. The luminosity combined with visible magnitudes from photometric surveys and extinction values gives an estimate of the distance, independent of the astrometric parallax.

The recent work of Minchev et al. (2019) argued that the lack of age information in Galactic Archaeology can lead to severe misinterpretation of the Milky Way formation and evolution. The authors showed that a number of chemo-kinematical relations used to study the Milky Way are plagued by a phenomenon known as Yule-Simpson’s paradox, which has the effect of erasing or completely reversing the trends seen in mono-age populations when age (or birth radius) is marginalized over.

In Mints & Hekker (2017, hereafter Paper 1) we introduced the implementation of this approach named Unified tool for Distance, Age and Mass estimations (UniDAM111http://www2.mps.mpg.de/homes/mints/unidam.html). This tool uses PARSEC models (Bressan et al., 2012) and infrared photometry from 2MASS (Skrutskie et al., 2006) and AllWISE (Cutri, R. M. et al., 2014) surveys to produce probability density functions (PDFs) in distance, log(age) and mass for a given star. A further extension of UniDAM that includes the use of Gaia parallaxes was presented in Mints & Hekker (2018) and Mints (2018). The output of UniDAM contains for a given star one or several solutions, with each solution having a unimodal probability density functions (PDF) in log(age), mass and distance (labelled as unimodal sub-PDF, or USPDF). We report for each solution, along with mean, median and mode values, the standard deviation and confidence intervals, also a label indicating the type of best-fitting unimodal function used to fit the USPDF and parameters of that function. This is done in order to overcome problems that are inherent to the isochrone matching method: non-Gaussianity and multi-modality of the produced PDFs. There are cases when age PDFs are close to being Gaussian, and thus ages can be inferred with high precision. This occurs when high quality data are used and a specific subset of stars is analysed, like it was done, for example, in Tucci Maia et al. (2016) with high-resolution spectra of solar twins and in Wu et al. (2017) with spectroscopic and asteroseismic data for main sequence turn-off stars. In a general case like a large spectroscopic survey, derivation of an age value from an age PDF is often more problematic. This is illustrated in Fig. 1, where we show examples of log(age) PDFs for several stars. Some log(age) PDFs are in fact close to exponential functions, which corresponds to flat PDFs in linear ages. For such stars we can put little or no constraints on their age. An important consequence of these problems is the fact that it is difficult or even impossible to provide for a given star a single age estimate that will be both accurate and precise. In the literature mean values (see, for example Feuillet et al., 2016; Queiroz et al., 2018) or modes (Xiang et al., 2017) of PDFs were used as age proxies. The implicit assumption is that these proxies are unbiased – at least in a statistical sense. In this work, we argue that this assumption is not always valid, and we provide a way of reconstructing the age distribution of a stellar population.

2 PDF quality

For our task we will be using UniDAM results presented in Mints & Hekker (2017). Fits to log(age) PDFs produced by UniDAM allow us to quickly reconstruct these PDFs for any subset of the survey. Here, we first want to show that these fits give a reliable representation of log(age) PDFs. We do that by comparing for each survey the stacked PDF for all stars with the sum of PDF fits produced by UniDAM. In this section, we analyse log(age) distributions to show that using USPDF fits from our catalogue have an advantage over using mean values alone or mean values and uncertainties. This is illustrated in Fig. 2, which shows several representations of log(age) distribution, for every survey processed by UniDAM:

• •

Histogram of mean (blue line in Fig. 2), median (purple line) or mode (brown line) values of USPDFs for all stars:

[TABLE]

where $\tau_{b}$ are centres of bins in log(age), $\Delta$ is the bin half-width and $w_{i}$ is the weight of the USPDF $i$222In UniDAM, USPDF weight is the fraction of the total PDF contained in the USPDF.. These representations thus retain only one parameter per USPDF, and therefore become very noisy on small stellar samples.

• •

Histograms of mean (orange line), median or mode values of USPDFs that are smoothed with log(age) uncertainties as:

[TABLE]

where $\sigma_{\tau,i}$ is the uncertainty in log(age) for USPDF $i$, $\mathcal{N}(\tau|m,v)$ is the PDF of the normal distribution with mean $m$ and standard deviation $v$ and the summation is done over all USPDFs for all stars in the survey. This representation is much smoother than histograms of mean, median or mode values, although it retains only two parameters ($\tau_{i}$ and $\sigma_{\tau,i}$) for each USPDF. Smoothed distribution of median and mode values are not shown in Fig. 2 to avoid overloading the plot.

• •

Stacked fits to USPDFs (green dashed line in Fig. 2):

[TABLE]

where $F_{i}$ is the best-fitting function for a solution’s PDF in log(age) and $p_{i}$ are the parameters of the fit. Again, the summation is done over all solutions (USPDFs) for all stars in the survey.

• •

“Full PDF” – the sum of all PDFs for all stars in the survey (red line in Fig. 2) – this sum is produced by UniDAM by directly adding up all PDFs for each solution in the fitting process. For each solution the detailed PDFs are however not stored, as it will require much more space than already used and will be difficult to work with.

The last two representations save more detailed information about log(age) PDFs than histograms and smoothed histograms of single parameters, and hence are more helpful in recovering real underlying age distributions.

Fig. 2 shows that stacked fits to USPDFs $N_{fit}(\tau)$ (see Equation 3) are close to the full PDF, which confirms that fits are a reliable proxy for the full PDF. This means, that for a subset of survey stars one can use stacked fits to USPDFs for stars in this subset in place of the full PDF for the same subset. The advantage is that UniDAM output allows to reconstruct the log(age) PDF for an arbitrary subset of stars without re-running the fit and without storing full PDF data for each star, which would take about 30 times more space. For that reason, in our subsequent analysis we use $N_{fit}(\tau)$.

From Fig. 2, it is also clear that different representation of the age distribution have a substantially different shape for most surveys. Thus, the choice of the representation might affect further conclusions based on the age distribution of stars. It is however unclear which distribution is closest to the real age distribution. As we show below with the use of simulated data for which the age distribution is known, none of them is in fact a good representation of the real age distribution, so none can be reliably used as an age tracer, at least not directly.

To further illustrate the problem, we build an artificial survey using stellar isochrone models with a constant age $\tau_{0}$ and feed it into UniDAM. The real distribution in log(age) will be a delta function $F(\tau)=\delta(\tau-\tau_{0})$. However, the log(age) PDF produced by UniDAM will not be a delta function. A few examples of these log(age) PDFs as produced by UniDAM are shown in Fig. 3. There we show log(age) PDFs for mono-age populations simulated in a way to mimic APOGEE (Majewski et al., 2017) and RAVE-on (Casey et al., 2016) surveys (see Section 4 for details on how the simulation were done). Strong difference in PDF shapes between APOGEE- and RAVE-on-based populations are due to different fractions of main-sequence stars in APOGEE and RAVE-on surveys: main sequence stars contribute more to the high-age part of the log(age) PDF. In most cases, log(age) PDFs have a peak near the population age, however the distributions are clearly non-Gaussian: they are very broad and in some cases show more than one peak. This is not a mistake, but an inherent property of the Bayesian isochrone fitting method used in UniDAM. In fact, almost any other isochrone fitting method will have the same property: isochrones overlap on the Hertzsprung–Russell diagram, resulting in the degenerate relation between observed spectroscopic parameters and physical properties of stars.

We can try to reconstruct the real log(age) distribution from the observed log(age) PDF for a given stellar population. To do that, we need to simulate a set of mono-age mock catalogues (hereafter MAMC) and produce log(age) PDFs for each catalogue in this set. If we have an observed mono-age population (for example, an open cluster), we can compare its log(age) PDFs to the log(age) PDFs for MAMCs. The age for MAMC with the closest PDF will be the estimate of the age of the observed population. If the population consists of stars with a set of ages $\tau_{1},\tau_{2},...\tau_{n}$, its log(age) PDF will be the sum of log(age) PDFs for MAMCs with same ages $\tau_{1},\tau_{2},...\tau_{n}$. We can formulate an inverse problem and try for an arbitrary observed population to represent its log(age) PDF as a linear combination of the MAMC log(age) PDFs, with coefficients as a function of log(age) giving an estimate of the true age distribution in the observed population. The details of this approach are explained below.

3 Inversion method

3.1 Main equation

The idea of the method proposed here is to try to represent the stacked log(age) PDF or analogous function for a given set of stars with the linear combination of log(age) PDFs for mono-age mock catalogues (MAMC). The method of MAMC construction is described below in Section 4. The important point here is that MAMC should be constructed in such a way that the distribution of its stars in $\log g$ and metallicity $[\rm{Fe/H}]$ will be close to the distribution in these parameters for the considered set of stars. In our approach MAMC with log(age) of $\tilde{\tau}$ contains stars with their spectrophotometric parameters taken from PARSEC model isochrones with the same log(age) $\tilde{\tau}$. Uncertainties were assigned to be close to observational ones (see more on how the uncertainties were chosen in Section 4). Because of the uncertainties in parameters, UniDAM will give for each star a log(age) PDF which will be different from a delta function $\delta(\tau-\tilde{\tau})$. In some cases, especially for stars on the main sequence, the PDF in log(age) will not have its peak around $\tilde{\tau}$, but will rather be close to an exponential function with its peak at the largest value of log(age) covered by models (see Fig. 1). The exponential distribution in log(age) is equivalent to a flat distribution in linear age, which means that age is poorly constrained or unconstrained in these cases. A stacked PDF of a MAMC, similarly, does not necessarily have its maximum at $\tilde{\tau}$. In Fig. 3 we show PDFs for several MAMC constructed to have their $\log g$ and $[\rm{Fe/H}]$ distributions similar to those for APOGEE (Majewski et al., 2017) and RAVE-on (Casey et al., 2016) surveys.

We designate the log(age) PDF, as produced by UniDAM for a given survey, as $C(\tau)$. We aim at estimating the underlying log(age) distribution of stars $N(\tilde{\tau})$ for the same survey. For the MAMC with log(age) $\tilde{\tau}$ we designate the log(age) PDF as $P(\tau,\tilde{\tau})$. If the sizes of the survey and all MAMCs are infinite, the following equation should hold:

[TABLE]

Thus, to get the unknown age distribution $N$ from the observed PDF $C$ we need to solve the integral equation 4. We note however that in reality Equation 4 holds only approximately, because the number of stars in the survey and in each MAMC is finite, thus both $C(\tau)$ and $P(\tau,\tilde{\tau})$ are subject to stochastic deviations. These deviations limit our ability to find the solution for Equation 4.

Let us consider as an example a case where the survey contains several mono-age populations (for example, a set of open clusters). This implies that $N(\tilde{\tau})=\sum_{i}N_{i}\delta(\tilde{\tau}-\tau_{i})$, with $N_{i}$ being the fractional weight of the $i$-th population, which has log(age) $\tau_{i}$. Following Equation 4, we get that $C(\tau)=\sum_{i}N_{i}P(\tau,\tau_{i})$, in other words, log(age) PDF for such a survey is a weighted sum of PDFs of mono-age populations. In fact, the PDF from UniDAM $C(\tau)$ is defined over the grid of log(ages) $(\tau_{i},i=1,..n)$. Similarly, MAMC PDFs $P(\tau,\tilde{\tau})$ are defined only for $\tau\in\tau_{i}$. Thus we can substitute the function $P(\tau,\tilde{\tau})$ with a quadratic matrix of the size $n\times n$ and functions $N(\tilde{\tau})$ and $C(\tau)$ with vectors $\@vec{N}(\tilde{\tau})$ and $\@vec{C}(\tau)$ of the size $n$, and Equation 4 can be rewritten as a system of linear equations:

[TABLE]

There is an obvious non-negativity constraint for this system: $N_{i}\geq 0$ for all values of $i$. Therefore, similarly to the example provided above, we try to represent the survey population as a superposition of $i$ mono-age populations with ages $\tau_{i}$.

3.2 Solving the main equation

In the system of linear equations (Equation 5), $C_{j}$ are taken from the survey log(age) PDF and $P_{j,i}$ comes from MAMC log(age) PDFs. We want to solve this system of linear equations for $N_{i}$, bearing in mind the non-negativity constraint. It is possible to find a solution by means of non-negative least squares method (NNLS, see Lawson & Hanson, 1995), which maximizes the following function:

[TABLE]

This will typically lead to a result in which only several components of $\@vec{N}$ will be non-zero, which is not physical – we expect a rather smooth log(age) distribution.

In order to get a smooth result we use regularized likelihood maximization. We chose a Tikhonov regularization, that favours smoother solutions by adding a sum of squares of the solution’s first derivatives333We performed tests with a version of Tikhonov regularization where the sum of squares of the solution values, rather than derivatives, is used, which resulted in very similar however slightly noisier results.. We thus maximize the following function:

[TABLE]

where $\lambda$ is a regularization parameter. Solutions with narrow spikes will produce large absolute values of $\frac{dN_{i}}{d\tau}$ and will have smaller values of likelihood function $L$, as compared to smoother solutions, even if the latter produce larger differences between $C_{j}$ and $\sum_{i=1}^{n}P_{ji}N_{i}$.

We can use finite difference formulas to calculate $\frac{dN_{i}}{d\tau}$ from $N_{i}$, and rewrite:

[TABLE]

where, $\mathcal{T}$ is a Toeplitz matrix representation of the first derivative:

[TABLE]

In this matrix, the $i$-th row represents coefficients for the 2nd order forward formula for the first derivative:

[TABLE]

with an exception of the last two rows, that represent similar coefficients for the backward formula:

[TABLE]

Thus, Equation 5 can be rewritten as:

[TABLE]

where:

[TABLE]

In this notation, maximization of $L$ from Equation 7 is equivalent to minimization of $\begin{Vmatrix}\@vec{C}^{\prime}-P^{\prime}\@vec{N}\end{Vmatrix}_{2}$, with a constraint that all components of $\@vec{N}$ are non-negative. Again, we can use NNLS to obtain the optimal vector $\@vec{N}$. The solution can be characterised by the sum of residuals $R=\begin{Vmatrix}\@vec{C}-P\@vec{N}\end{Vmatrix}_{2}$, which indicates, how close the log(age) PDF for the derived age distribution is to $C(\tau)$.

3.3 Searching for optimal $\lambda$ value

The problem of the regularization approach is that the parameter $\lambda$ has to be properly chosen in order to obtain the correct solution of the problem. Two extreme cases are $\lambda=0$ and $\lambda=\infty$. In the first case no regularization is active, and Equation 12 is reduced to Equation 5. In that case the solution is most precise however not smooth. If $\lambda=\infty$ the regularization dominates the solution, such that $N(\tau)$ becomes a constant, equal to the mean value of $C(\tau)$ – this is the smoothest possible solution, which is very imprecise. In this section we describe an empirical method of obtaining the value of $\lambda$ that gives the result we consider to be optimal.

First of all, we note that both $P$ and $\@vec{C}$ are subject to statistical variations due to the limited number of sources in both the observed survey and in the simulated MAMC. Because each star can contribute to a wide range of log(age) values, variations of the PDF at different values of log(age) $\tau$ are not independent. In order to account for this effect properly, we chose to produce five realizations of MAMC (differing only by random number algorithm seed) and to split the observed survey into five parts. This gives us five $P$ matrices $P_{p},p=1,2,..5$ and five $\@vec{C}$ vectors $\@vec{C}_{q},q=1,2,..5$. Using this information, we can make use of the cross-validation technique, requiring that the solution obtained with one combination of $p$ and $q$ values should be good for all other combinations.

For a given value of regularization parameter $\lambda$ we can thus build a set of 25 equations similar to Equation 12, with all possible combinations of $\@vec{P}_{p}$ and $\@vec{C}_{q}$. Let us now focus on the combination with $p=p_{0}$ and $q=q_{0}$. For this combination, we write Equation 12 as:

[TABLE]

and designate a solution of this equation as $\@vec{N}_{p_{0},q_{0}}$. We can further build a matrix of residuals with respect to all combinations of $\@vec{C}$ and $\@vec{P}$:

[TABLE]

As our quality parameter we take the following expression:

[TABLE]

All $R_{m,n,p_{0},q_{0}}$ and ${Q_{p_{0},q_{0}}}$ are functions of the regularization parameter $\lambda$. We show an example of $R_{m,n,p_{0},q_{0}}(\lambda)$ and $Q_{p_{0},q_{0}}(\lambda)$ functions in Fig. 4. $R_{p_{0},q_{0},p_{0},q_{0}}$ decreases as $\lambda$ decreases – with less regularization it is possible to have more accurate solution of Equation 12. For values of $R_{m,n,p_{0},q_{0}}$ for $m\neq p_{0}$ or $n\neq q_{0}$ the decreasing trend with decreasing $\lambda$ stops or even reverses at some point. This happens because solutions $\@vec{N}_{p_{0},q_{0}}$ of Equation 12 that are accurate for the combination $\@vec{P}_{p_{0}}$ and $\@vec{C}_{q_{0}}$ become too “specialized”, and are not as accurate for other combinations of $\@vec{P}_{p}$ and $\@vec{C}_{q}$. As an optimal value of $\lambda$ we will take the point where the rapid decrease of $Q_{p_{0},q_{0}}$ (hereafter labelled simply as $Q$) with decreasing $\lambda$ stops or slows down. To find this point, we first fit a following piecewise linear function to the $Q$-function in log-log space:

[TABLE]

Function $Q_{fit}$ has three linear segments, with transitions at values of $\lambda=b_{1}$ and $\lambda=b_{2}$. Slopes are $c$ in the first segment (close to zero in Fig. 4), $d$ in the second segment and zero in the third segment. Parameter $a$ represents the value of $Q_{fit}(b_{1})$. This choice of the representation for the piecewise linear function is motivated by the ease of the initial guess for the fitted parameters. All five parameters ($a,b_{1},b_{2},c,d$) are fitted simultaneously.

As the first estimate for the optimal regularization parameter value we take the first turning point of the fitted function $\lambda_{est}=b_{1}$. In many cases the turn of the $Q$ function is not as sharp as the one of the fitted function, and $\lambda_{est}$ is an overestimation of the optimal value of $\lambda$. To correct for that, we take as our final estimate $\lambda_{f}$ the first point to the left of $\lambda_{est}$ where $Q$ and $Q_{fit}$ intersect, as shown in Fig. 4. This estimate in some cases tends to produce a noisier “over-fitted” result, while the $\lambda_{est}$ tends to produce smoother “under-fitted” one. On the average, however, $\lambda_{f}$ produces better results than $\lambda_{est}$. Our choice of the $\lambda_{f}$ is further validated with tests described below in Section 5.

We can repeat the above procedure for every possible combination of $p_{0}$ and $q_{0}$, thus obtaining 25 different values of $\lambda_{f}$ and corresponding solutions $\@vec{N}_{p_{0},q_{0}}$ of Equation 15. As the final solution we can take the mean of $\@vec{N}_{p_{0},q_{0}}$, and as a measure of the uncertainty – their standard deviation. For a more detailed uncertainty analysis see Section 5.3.

4 Mono-age mocks construction

In order to properly recover the age distribution for a given survey, which we call a base survey, we need to construct MAMCs in such a way, that the distribution of stars in each MAMC in physical parameters and their uncertainties will be close to that of the base survey, at the same time retaining age information. To achieve that, we follow the procedure below.

We build MAMC using PARSEC models (Bressan et al., 2012) to simulate stars. In order to mimic the random observational scatter we take instead of each model a sample of 25 models with random normally distributed perturbations added in $T_{\rm{eff}},\log g$ and $[\rm{Fe/H}]$. To assign a proper amplitude of perturbations, we randomly select uncertainties $\sigma_{P}=\{\sigma_{T},\sigma_{\log g},\sigma_{[\rm{Fe/H}]}\}$ from the base survey data for stars that have $P=\{T_{\rm{eff}},\log g\,\textrm{and}\,[\rm{Fe/H}]\}$ close to the model, such that the difference between physical parameters for the model and for the base survey star is smaller than the mean uncertainty in the respective parameter for the base survey $P_{model}-P_{star}<\textrm{mean}(\sigma_{P})$. This provides a way to reconstruct the scatter and the systematic variations of $\sigma_{P}$ across the parameter space. Selected uncertainties are not only used as perturbation amplitudes, they are also assigned as “observational” uncertainties for models. The photometric magnitudes of all models are perturbed with a Gaussian noise with a scale of $0.^{m}025$, to reproduce typical photometric uncertainties in 2MASS and AllWISE.

The set of models is then truncated in the $T_{\rm{eff}},\log g,[\rm{Fe/H}]$ space to the footprint of the base survey in that space, models outside of that footprint are excluded. In order to build MAMC from that set of models, we sample a pre-defined number of models $n_{\textrm{mock}}$ as follows. We bin stars from the base survey in $\log g-[\rm{Fe/H}]$ space. The fraction of models sampled for a given MAMC from each bin is equal to the fraction of base survey stars in that bin. Within each bin, models are selected randomly with probabilities proportional to the fraction of the initial mass function represented by that model.

We cannot use all three physical parameters ($T_{\rm{eff}},\log g,[\rm{Fe/H}]$) for binning, as in this case the resulting MAMC will be too close to the base survey itself, retaining almost no information about the MAMC age value. On the other hand, using just one parameter will produce MAMC that will not resemble the base survey, making inversion impossible. We decided to use as a first parameter for the binning $[\rm{Fe/H}]$, as it is the only parameter out of the three that is independent of stellar mass, and as such is expected to be the same for all stars of the same origin. The second parameter is $\log g$, as it helps to distinguish main sequence stars and giants, that have very different contribution the log(age) PDF. Such distinction is impossible to make with $T_{\rm{eff}}$.

The above procedure is repeated for each value of log(age) $\tau$ over the considered grid: $6.61\leq\tau\leq 10.13$ dex with a step of $0.02$ dex or age between approximately $4\cdot 10^{6}$ and $13.5\cdot 10^{9}$ years. As a result we get a set of MAMCs that is fed into UniDAM to obtain a set of log(age) PDFs $P_{j,i}$.

5 Tests with mock data

In this section, we use simulated data to validate the choice of the optimum regularization parameter $\lambda$ and to determine the influence of the survey size $n_{\textrm{sim}}$ and MAMC size $n_{\textrm{mock}}$ on the precision and the accuracy of the inversion result. We simulate a survey with pre-defined age distribution $\@vec{N}_{input}$ as a concatenation of MAMC simulated in a way described in Section 4, with the size of $i$-th MAMC defined by $N_{input,i}$. We than obtain with UniDAM the log(age) PDF for the simulated survey and apply the inversion as described in the Section 3.2. We can than compare the result of the inversion with the input age distribution, to get an estimate of the accuracy and precision of the method. Below we describe the process in detail.

5.1 Input distributions

For this work we chose six input age distributions, aiming to emulate critical as well as more common cases. Every distribution was generated using one of the three base surveys: APOGEE DR14 (Majewski et al., 2017), LAMOST DR3 (Luo et al., 2015) and RAVE-on (Casey et al., 2016). In Fig. 5 we show the input log(age) distributions and PDFs produced by UniDAM for the three base surveys. “Age” and “box” inputs implement a spike and step function. Such functions are hard to reproduce with our method, as the uncertainties in age determination tend to smooth out the distribution. “Block” and “bump” inputs show the opposite case of smooth slowly changing PDFs. Hence, we expect that the inversion will work best for them. “Combined” and “wave” inputs represent intermediate cases. Interestingly, log(age) PDFs for “bump” and “wave” simulations have very similar shapes for all three base surveys, despite the very different input log(age) distributions. Note that the choice of base survey affects the log(age) PDFs generated by UniDAM. Most importantly, the log(age) PDF is typically not close to the true age distribution, with the only exception of the APOGEE result for the “bump” distribution.

5.2 Test results

In Figs. 6 and 7 we show examples of test results, with different input distributions and varying $n_{sim}$ and $n_{mock}$. There we compare input distributions (blue lines) with age inversion results (red lines), and “observed” log(age) PDFs (dashed black lines) with log(age) PDFs as predicted by age inversion (grey lines, one for each combination of simulated survey and MAMC realizations).

Figs. 6 and 7 illustrate that the result of the inversion is almost insensitive to $n_{\textrm{mock}}$ for small values of $n_{\textrm{sim}}$, as expected (see discussion in Section 5.3). On the other hand, for $n_{\textrm{sim}}\geq 25000$, MAMC size $n_{\textrm{mock}}$ plays a larger role. It is also clear that an increase in $n_{\textrm{sim}}$ and $n_{\textrm{mock}}$ increases the sensitivity of the inversion to rapid changes in log(age) distribution. For example, in the case of “wave” test, $n_{\textrm{mock}}=250$ and $n_{\textrm{sim}}=25000$ is enough to properly recover the second peak of log(age) distribution at $\tau\approx 9.6$ dex, though we need to increase $n_{\textrm{mock}}$ to 5000 in order to properly recover the first, much narrower peak at $\tau\approx 8.9$ dex. For all simulations, “observed” (black lines) and predicted (grey lines) log(age) PDFs become almost indistinguishable for $n_{\textrm{mock}}=5000$ and $n_{\textrm{sim}}=25000$.

5.3 Uncertainty analysis

An important part of every scientific result is a proper uncertainty. We can use our tests to provide a way to estimate uncertainties for age inversion results.

In Section 3.3 we gave a method of obtaining the optimal value of regularization parameter $\lambda$. This parameter and a corresponding solution of Equation 12 can be obtained for every combination of the five survey log(age) PDFs $\@vec{C}_{q}$ and five MAMC log(age) PDF sets $\@vec{P}_{p}$. Hence we can have 25 different solutions of Equation 15, with there differences being due to statistical variations in $\@vec{P}$ and $\@vec{C}$. Variations between these solutions can be used to estimate the uncertainty of the solutions. Because of the smoothness of the solution, which is imposed by the regularization, and because of the similarity between log(age) PDFs for MAMCs with similar log(age) values, differences between solutions at different values of log(age) $\tau$ are highly correlated. For these reasons, the distribution of differences between the inversion result and the true underlying age distribution will not be a Gaussian but a heavy-tailed distribution.

In the case of the tests that we are performing, we also use the knowledge of the true age distribution to find the relation between the variation between solutions $\sigma(\tau)$ and the true error $\delta(\tau)$ (the true error is defined as the absolute difference between input and output log(age) distributions). The fraction $F(\omega)$ of values of $\tau$ with $\delta(\tau)<\omega\sigma(\tau)$ can be calculated for each result. In Fig. 8 we show distributions of $F(\omega)$ for $\omega=1$ and $\omega=3$. If $\sigma(\tau)$ is a correctly determined Gaussian uncertainty, than $F(1)=0.68$ and $F(3)=0.997$, as expected for 1- and 3-sigma confidence intervals. Measured $F(1)$ values for all our tests show a broad distribution with a mean of 0.67 and median of 0.685 – very close to expectations. Measured $F(3)$ values also have a broad distribution, with a mean of 0.93 and median of 0.96 – lower than the expected value, making $3\sigma(\tau)$ effectively a 2-sigma rather than 3-sigma confidence interval estimate. This is because uncertainties have a heavy-tail non-Gaussian distribution. There are no clear trends visible for $F(\omega)$ values with $n_{\textrm{sim}}$ and $n_{\textrm{mock}}$, at least within the considered ranges of these parameters.

We note here, that the scatter of $F(\omega)$ around the median values is high, which means that values of $\sigma(\tau)$ can be seen only as an approximation of the uncertainty.

One would expect the fractional uncertainty values to decrease as $(n_{\textrm{sim}}n_{\textrm{mock}})^{-1/2}$. In reality, there relation is not as strong, because of complex correlations between the real age distribution and the corresponding log(age) PDF. Our estimate shows that the fractional uncertainty scales as $n_{\textrm{sim}}^{-0.4}n_{\textrm{mock}}^{-0.15}$. This means that one has to increase $n_{\textrm{mock}}$ by an order of magnitude to improve the result by about $30\%$. Importantly, for a given $n_{\textrm{mock}}$ there seem to exist a maximum value of $n_{\textrm{sim}}$, above which the fractional uncertainty does not decrease with $n_{\textrm{sim}}$ at all. This is caused by the fact that in that regime fractional uncertainty becomes dominated by variations between realizations of MAMC, and not by the survey. The opposite is also true – for a given $n_{\textrm{sim}}$ there is a maximum value of $n_{\textrm{mock}}$, above which the fractional uncertainty does not decrease, as it is dominated by survey PDF variations. The exact value of the maximum of $n_{\textrm{mock}}$ for a given $n_{\textrm{sim}}$ depends on the base survey and on the shape of the underlying age distribution and can vary by over an order of magnitude. The maximum considered value $n_{\textrm{mock}}=5\,000$ is sufficient for surveys as large as $n_{\textrm{sim}}=50\,000$. This is still below typical sizes of spectroscopic surveys. However, if we want to slice the survey into parts to trace, for example, age-metallicity relations, as is done below in Section 6.2, the size of each part will be on the order of $n_{\textrm{sim}}=50\,000$ or even smaller. Even more importantly, all these measurements are done for artificial data, where systematic uncertainties are zero by definition, as simulated and mock catalogues are created from the same set of models. For real data, as we will show below, systematics will likely be the dominating source of uncertainty.

6 Real data applications

6.1 Inversion of full surveys

We apply the inversion method described above to several large spectroscopic surveys, namely APOGEE (Majewski et al., 2017), GALAH (Martell et al., 2016), RAVE-on (Casey et al., 2016) and LAMOST (Luo et al., 2015). The results are presented in Fig. 9. These plots indicate that all surveys have a bimodal age distribution, although the location, width and relative amplitudes for the two modes are different. There are several possible explanations for this fact. The first possibility is that the trend is real and is caused by the target selection in the survey, with stars of different metallicity being observed in a different parts of the Galaxy. The second possibility is that one of the implicit assumptions of the method does not hold, namely, that in a given sample age does not depend on $\log g$ or $[\rm{Fe/H}]$. Within the metallicity bin for a magnitude limited survey, observed stars with lower values of $\log g$, and thus with higher luminosities, are on the average located at larger distances than those with higher $\log g$, and thus with lower luminosities. Therefore age distribution might vary over the observed $\log g$ range. This may cause the inversion method to give incorrect results. It remains unclear however, which of the effects dominates.

Last but not least, PARSEC models and spectroscopic measurements can have a systematic offset between them, which can cause a complex systematic age bias, which cannot be accounted for in modelling.

Possible systematic effects manifest themself also in the difference between the inversion result and log(age) PDFs of the survey, indicated in Fig. 9 with black dashed and grey solid lines. This difference is considerably larger than the one we obtain for the simulated data, where the systematic offset is zero, even though the survey size $n_{\textrm{sim}}$ is larger than those considered in the simulation.

6.2 Inversion of mono-metallicity populations

The limitation of the inversion of the full survey is that the age distribution is assumed to be the same for all distances and all metallicities covered by the survey, which in general does not hold, as more metal poor stars tend to be older. In order to mitigate this problem we made a separate study of RAVE-on, LAMOST (Luo et al., 2015) DR4 and APOGEE stars, splitting them in several metallicity bins and applying age inversion to each bin separately.

In Fig. 10 we present the inputs (log(age) PDFs) and outputs (underlying age distributions) of the age inversion procedure for RAVE-on stars. Inputs (log(age) PDFs) are typically broad and for metal-poor populations show an extra narrow peak at log(age) $\tau\approx 9$, which is associated with red-clump stars. The underlying age distributions computed through the inversion are a lot smoother and narrower, with no secondary peaks. The peak of the underlying age distribution goes from log(age) $\tau=9.37$ for highest metallicity bin $[\rm{Fe/H}]=0.4$ dex to the maximum possible log(age) of $\tau=10.13$ for $[\rm{Fe/H}]<-0.4$ dex. At very low metallicities $[\rm{Fe/H}]\leq-1.2$ a gradual broadening of the underlying age distribution is observed, which is caused by the presence of a lower number of stars in metal-poor bins and thus a reduced age resolution (see test results with low $n_{\textrm{sim}}$ value in Fig. 6 and 7).

For LAMOST, we use a constraint of $T_{\rm{eff}}<7000$K in this work, as there are visible pipeline artefacts in the Hertzsprung–Russell diagram at higher temperatures, which we are not able to simulate. Even with that cut, LAMOST contains about an order of magnitude more stars than RAVE-on or APOGEE. The results of the age inversions for the LAMOST sample are shown in Fig. 11 and for the APOGEE sample in Fig. 12. Results for LAMOST show narrower underlying age distributions than RAVE-on or APOGEE results, likely due to larger statistics and hence higher age resolution.

For LAMOST, two peaks are visible in the underlying age distributions for metallicities $[\rm{Fe/H}]\geq-0.2$ dex: one at around $\tau\approx 9.35$ and one at $\tau\approx 9.83$. The latter peak becomes dominant and shifts to higher $\tau$ values, as $[\rm{Fe/H}]$ increases. Similarly, for APOGEE a secondary peak is visible in the underlying age distributions for metallicities $[\rm{Fe/H}]\geq-0.2$ dex at around $\tau\approx 8.9$. In both cases the result is that the mean age for bins with $[\rm{Fe/H}]\geq+0.2$ dex is higher than that for solar-like metallicities. This does not seem to be physical – we expect more metal rich stars to be systematically younger. This might be attributed to the limitations of the method listed above in the Section 6.1. An alternative explanation is that the effect is physical, and is caused by the fact that metal rich ($[\rm{Fe/H}]\geq+0.2$ dex) originate in the inner part of the Galaxy, and it takes a certain amount of time for the migration process to bring them to the solar vicinity, where they can be observed (see Minchev et al., 2018). Thus a small deficit of young metal rich stars can be expected.

It is very likely that all three effects are acting simultaneously. For better results, selection effects (like those presented in Mints & Hekker (2019)) have to be taken into account. Binning the data in distance or galactic latitude is also desirable. Both tasks are beyond the scope of the current work.

Overall, the underlying age distributions obtained through the inversion allows to trace age-metallicity trends that are much less obvious if log(age) PDFs are considered. It also allows to remove structures like low-age spikes, that do not seem to be real and are likely due to observational scatter in physical parameters of older stars.

6.3 Metallicity-age and age-metallicity relations

Results from Section 6.2 allow us to study metallicity-age and age-metallicity relations for RAVE-on, LAMOST and APOGEE. Both relations can be defined through a two-dimensional distribution function $C(\tau,[\rm{Fe/H}])=C_{[\rm{Fe/H}]}(\tau)N([\rm{Fe/H}])$, where $C_{[\rm{Fe/H}]}(\tau)$ is the log(age) PDF in the bin with metallicity $[\rm{Fe/H}]$, and $N([\rm{Fe/H}])$ is the number of stars in the survey as a function of metallicity. From that we can define log(age) as a function of metallicity as:

[TABLE]

and similarly metallicity as a function of log(age):

[TABLE]

Here, $(\tau_{min},\tau_{max})$ is the considered range in log(age) and $([\rm{Fe/H}]_{min},[\rm{Fe/H}]_{max})$ is the range in metallicity.

In Figs. 13, 14 and 15 we show $\tau([\rm{Fe/H}])$ with blue solid and $[\rm{Fe/H}](\tau)$ with blue dotted lines. In these plots we will continue to use solid lines for $\tau([\rm{Fe/H}])$ relations and dotted lines for $[\rm{Fe/H}](\tau)$ relations. These two lines are nearly orthogonal, which is the consequence of the $C(\tau,[\rm{Fe/H}])$ distribution being very broad.

We can now replace log(age) PDFs $C(\tau,[\rm{Fe/H}])$ with the underlying age distribution from the inversion $N(\tau,[\rm{Fe/H}])$ in Equation 19 and 20. The results are shown in Figs. 13, 14 and 15 with red solid and dotted lines. They are a lot closer to each other, which is a sign of a much tighter $C(\tau,[\rm{Fe/H}])$ function.

These data were compared to the results of the recent work by Feuillet et al. (2018), who presented the analysis of a sample of 721 nearby red giant stars selected from APOGEE (Majewski et al., 2017). All these stars are closer than 400 pc and have reliable TGAS parallaxes, which allows to determine log(age) to 0.07 dex precision. This sample was used, among other applications, to derive mean ages as a function of metallicity, shown in Figs. 13, 14 and 15 with orange solid line. We note that the data presented in Feuillet et al. (2018), as well as our results for $\tau([\rm{Fe/H}])$ give age distribution and mean age in each metallicity bin, thus measuring age as a function of metallicity $\tau([\rm{Fe/H}])$. At the same time, models of chemical evolution typically focus on the metallicity as a function of age $[\rm{Fe/H}](\tau)$ (see for example Fig. 4 in Minchev et al., 2013), which is strictly speaking a different function. Functions $\tau([\rm{Fe/H}])$ and $[\rm{Fe/H}](\tau)$ can be close to each other only if the two-dimensional distribution $C(\tau,[\rm{Fe/H}])$ is tight. Hence in a general case care must be taken in comparing $\tau([\rm{Fe/H}])$ and $[\rm{Fe/H}](\tau)$, as it might lead to wrong conclusions.

If we compare the age-metallicity distributions from log(age) PDFs $C(\tau,[\rm{Fe/H}])$ with those obtained from the underlying age distributions computed through the inversion, $N(\tau,[\rm{Fe/H}])$, we see that the latter are much closer to Feuillet et al. (2018), than the former. The small systematic offset can be attributed to the fact that, as opposed to the solar neighbourhood APOGEE subsample used by Feuillet et al. (2018), surveys used here contain more thick-disk and halo stars that are systematically older.

We can compare our age-metallicity and metallicity-age trends with those predicted by chemical evolution models. We use the model described in Minchev et al. (2013). We added an additional smoothing to the two-dimensional age-metallicity distribution predicted by this model. The smoothing scale was chosen to be close to the typical uncertainty in age (2 Gyrs) and metallicity (0.1 dex). We then calculated $\tau([\rm{Fe/H}])$ and $[\rm{Fe/H}](\tau)$ for this distribution and show it in Figs. 13, 14 with black solid and 15 with black dotted lines. Note that mean age $\tau$ starts to increase as a function of $[\rm{Fe/H}]$ for $[\rm{Fe/H}]>0$ for the model data as well as LAMOST and APOGEE inversion results. This is very likely caused by the absence of young metal-rich stars in the solar vicinity, which are formed in the inner Galaxy and need time to migrate outwards. Mean ages for metal poor stars are also systematically larger for inversion results than those predicted by models. This might be caused by the fact that UniDAM allows for stellar ages up to the age of the Universe ($\approx 13.5\,$Gyrs), while the maximum stellar age in the Minchev et al. (2013) model is $11.175\,$Gyrs (before the smoothing was applied). In general, on all three plots $\tau([\rm{Fe/H}])$ and $[\rm{Fe/H}](\tau)$ from the model and from age inversion results are close to each other.

7 Summary and outlook

In this work we present age inversions – a method of revealing the underlying distribution of a large ($N>10^{3}$) ensembles of stars in log(age) from their cumulative log(age) PDF produced by UniDAM. This allows us to remove biases inherent to isochrone fitting. The method was tested on simulated data showing that it allows us to reconstruct the underlying age distribution of stars in the simulated sample. The inversion results are a lot closer to the real underlying age distribution than the stacked log(age) PDF or the distribution of mean PDF values.

The method was further applied to data produced by UniDAM for real surveys, deriving different age distributions for different surveys. This is expected, as surveys use different observational strategies and focus at different samples of stars. Systematic offsets between survey data and PARSEC models used in this work and survey pipeline artefacts are the main limitations of the method. Due to this we are limited in our ability to generate mono-age populations that will be distributed closely to the survey stars, which is needed for a reliable age inversion.

We also apply age inversion method to mono-metallicity slices of RAVE-on, LAMOST and APOGEE surveys. This allows us to trace how the age distribution changes as a function of metallicity, and to reconstruct both metallicity-age and age-metallicity relations for RAVE-on, LAMOST and APOGEE samples, successfully removing artefacts of the isochrone fitting. We obtain results similar to those published in Feuillet et al. (2018) for a much smaller high-precision APOGEE sample and to those predicted by chemical evolution models. The number of stars in both RAVE-on and LAMOST surveys, at least for metallicities close to solar, is high enough to make further division in, for example, galactic latitude, distance or $\alpha$-abundance bins. This is considered beyond the scope of the current work.

In the future, systematic offsets between data and models are expected to decrease as the uncertainties (both random and systematic) of observations decrease and as models improve. Furthermore, the future spectroscopic surveys such as 4MOST (de Jong et al., 2016) and WEAVE (Dalton et al., 2014) will provide data for millions of stars. The age inversion method presented here will be suited to provide the underlying age distribution of such large survey data sets as a whole as well as as functions of, for example, metallicity, location in the Galaxy and $\alpha$-elements abundances.

One can be tempted to use the vast amount of data contained in Gaia to derive age distributions from photometry alone, like it was done, for example, in Dolphin (2013). However, it might be difficult given the lack of metallicity information in photometric data, which is essential for unambiguous age determination. In any case, such kind of study will differ substantially from that presented in this work, and the new method to solve Equation 4 might be needed.

It is possible to make use of Gaia parallax data to improve the precision of ages derived from spectrophotometric data, as it was done in Mints & Hekker (2018) and Mints (2018). However, this will require a more complex simulation strategy that will take into account proper parallax uncertainty distribution. Apart from that, Gaia parallax zero-point offset (see, for example, Leung & Bovy, 2019) have to be treated properly, as it can inflict a systematic age bias.

Acknowledgements

Authors thank the SAGE group at Max Planck Institute for Solar System Research for many fruitful discussions. Authors thank the anonymous referee for a detailed report with many useful suggestions, which helped us to improve the manuscript substantially.

The research leading to the presented results has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement (No 338251, StellarAges).

This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013). This research made use of matplotlib, a Python library for publication quality graphics (Hunter, 2007). This research made use of SciPy (Scipy team, 2001). This research made use of TOPCAT, an interactive graphical viewer and editor for tabular data (Taylor, 2005). Funding for RAVE has been provided by: the Australian Astronomical Observatory; the Leibniz-Institut fuer Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory; and the Universities of Groningen, Heidelberg and Sydney. The RAVE web site is at https://www.rave-survey.org. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University. 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. The GALAH survey is based on observations made at the Australian Astronomical Observatory, under programmes A/2013B/13, A/2014A/25, A/2015A/19, A/2017A/18. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present.