Issue 
A&A
Volume 629, September 2019



Article Number  A127  
Number of page(s)  14  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201935864  
Published online  16 September 2019 
Ensemble age inversions for large spectroscopic surveys
^{1}
LeibnizInstitut für Astrophysik Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany
email: amints@aip.de, minzastro@gmail.com
^{2}
Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
^{3}
Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, 8000 Aarhus C, Denmark
Received:
10
May
2019
Accepted:
10
August
2019
Context. Galactic astrophysics is now in the process of building a multidimensional map of the Galaxy. For such a map, stellar ages are an essential ingredient. Ages are measured only indirectly however, by comparing observational data with models. It is often difficult to provide a single age value for a given star, as several nonoverlapping solutions are possible.
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 monoage populations. Weights in that sum give the underlying log(age) distribution. Monoage 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 RAVEon, LAMOST and APOGEE we also applied age inversion to monometallicity samples, successfully recovering age–metallicity trends present in higherprecision 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 age–metallicity relations, for instance.
Key words: methods: statistical / stars: fundamental parameters / Galaxy: stellar content
© ESO 2019
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 the present day. Large astrometric surveys, of which Gaia (Perryman et al. 2001) with its morethan 10^{9} objects is now the dominating one, provide us with kinematic information such as positions, proper motions, and parallaxes. Spectroscopic data can be use to obtain radial velocities and the physical properties of stars such as temperatures, surface gravities, and chemical compositions. 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 by these latter authors, the socalled isochrone matching, can be used to derive both ages and distances, and is widely applied to spectroscopic data. The isochrone matching method is based on comparison of quantities measured spectroscopically (like the effective temperature T_{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 a star. The luminosity combined with visible magnitudes from photometric surveys and extinction values gives an estimate of the distance, independent of the astrometric parallax.
Recent work by Minchev et al. (2019) argues 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 chemokinematical relations used to study the Milky Way are plagued by a phenomenon known as YuleSimpson’s paradox, which has the effect of erasing or completely reversing the trends seen in monoage populations when marginalizing over age (or birth radius).
In Mints & Hekker (2017; hereafter Paper 1) we introduced the implementation of an approach named Unified tool for Distance, Age, and Mass estimations (UniDAM^{1}). This tool uses PARSEC models (Bressan et al. 2012) and infrared photometry from 2MASS (Skrutskie et al. 2006) and AllWISE (Cutri 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 PDF in log(age), mass, and distance (labelled as unimodal subPDF, or USPDF). We report for each solution, along with mean, median, and mode values, the standard deviation and confidence intervals, as well as a label indicating the type of bestfitting 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 isochronematching method: nonGaussianity and multimodality of the produced PDFs. There are cases where age PDFs are close to being Gaussian, and thus ages can be inferred with high precision. This occurs when highquality data are used and a specific subset of stars is analysed, as was done, for example, by Tucci Maia et al. (2016) with highresolution spectra of solar twins and by Wu et al. (2017) with spectroscopic and asteroseismic data for main sequence turnoff stars. In a general case, as in 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 close to exponential functions, which correspond 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, e.g. 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.
Fig. 1. Examples of log(age) PDF for various stars. We note that PDFs often show nonGaussian features: heavy tails, truncation, multiple modes. 
Fig. 2. Different representations of log(age) distribution for public spectroscopic surveys. See Sect. 2 for a discussion. 
2. Quality of probability density functions
For our task we use 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 confer 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:
$$\begin{array}{c}\hfill N({\tau}_{\mathrm{b}})={\displaystyle \sum _{i:{\tau}_{i}\in ({\tau}_{\mathrm{b}}\mathrm{\Delta},{\tau}_{\mathrm{b}}+\mathrm{\Delta})}}{w}_{i},\end{array}$$(1)
where τ_{b} are the centres of bins in log(age), Δ is the bin halfwidth, and w_{i} is the weight of the USPDF i^{2}. 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
$$\begin{array}{c}\hfill {N}_{\mathrm{unc}}(\tau )={\displaystyle \sum _{i}}{w}_{i}\mathcal{N}(\tau {\tau}_{i},{\sigma}_{\tau ,i}),\end{array}$$(2)
where σ_{τ, i} is the uncertainty in log(age) for USPDF i, 𝒩(τ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 (τ_{i} and σ_{τ, i}) for each USPDF. The smoothed distribution of median and mode values is not shown in Fig. 2 to avoid overloading the plot.
– Stacked fits to USPDFs (green dashed line in Fig. 2):
$$\begin{array}{c}\hfill {N}_{\mathrm{fit}}(\tau )={\displaystyle \sum _{i}}{w}_{i}{F}_{i}(\tau ,{p}_{i}),\end{array}$$(3)
where F_{i} is the bestfitting function for the PDF of a solution 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 this would require much more space than already used and would be difficult to work with.
These 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.
Figure 2 shows that stacked fits to USPDFs N_{fit}(τ) (see Eq. (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 of doing so is that UniDAM output allows us to reconstruct the log(age) PDF for an arbitrary subset of stars without rerunning 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}(τ).
From Fig. 2, it is also clear that a different representation of the age distribution has a substantially different shape for most surveys. Therefore, 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 are in fact a good representation of the real age distribution, so none of them 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 τ_{0} and feed it into UniDAM. The real distribution in log(age) will be a delta function F(τ) = δ(τ − τ_{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 monoage populations simulated in such a way as to mimic APOGEE (Majewski et al. 2017) and RAVEon (Casey et al. 2017) surveys (see Sect. 4 for details on how the simulations were done). The strong difference in PDF shapes between APOGEE and RAVEonbased populations are due to the different fractions of main sequence stars in APOGEE and RAVEon surveys: main sequence stars contribute more to the highage part of the log(age) PDF. In most cases, log(age) PDFs have a peak near the population age, however the distributions are clearly nonGaussian; 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.
Fig. 3. Examples of log(age) PDFs derived with UniDAM (same as red lines in Fig. 2) for simulated monoage populations with log g and [Fe/H] distributions taken from APOGEE (left) and RAVEon (right). Vertical lines indicate the true age for each population. 
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 monoage mock catalogues (MAMCs) and produce log(age) PDFs for each catalogue in this set. If we have an observed monoage population (e.g. an open cluster), we can compare its log(age) PDFs to the log(age) PDFs for MAMCs. The age for the 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 τ_{1}, τ_{2}, …τ_{n}, its log(age) PDF will be the sum of log(age) PDFs for MAMCs with the same ages τ_{1}, τ_{2}, …τ_{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 aim 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 monoage mock catalogues (MAMC). The method of MAMC construction is described below in Sect. 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 [Fe/H] is close to the distribution in these parameters for the considered set of stars. In our approach, MAMC with log(age) of $\stackrel{\sim}{\tau}$ contains stars with their spectrophotometric parameters taken from PARSEC model isochrones with the same log(age) $\stackrel{\sim}{\tau}$. Uncertainties were assigned to be close to observational ones (see more on how the uncertainties were chosen in Sect. 4). Because of the uncertainties in parameters, UniDAM will give for each star a log(age) PDF which is different from a delta function $\delta (\tau \stackrel{\sim}{\tau})$. In some cases, especially for stars on the main sequence, the PDF in log(age) will not have its peak around $\stackrel{\sim}{\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 $\stackrel{\sim}{\tau}$. In Fig. 3 we show PDFs for several MAMC constructed to have their log g and [Fe/H] distributions similar to those for APOGEE (Majewski et al. 2017) and RAVEon (Casey et al. 2017) surveys.
We designate the log(age) PDF, as produced by UniDAM for a given survey, as C(τ). We aim at estimating the underlying log(age) distribution of stars $N(\stackrel{\sim}{\tau})$ for the same survey. For the MAMC with log(age) $\stackrel{\sim}{\tau}$, we designate the log(age) PDF as $P(\tau ,\stackrel{\sim}{\tau})$. If the sizes of the survey and all MAMCs are infinite, the following equation should hold:
$$\begin{array}{c}\hfill C(\tau )={{\displaystyle \int}}_{0}^{\infty}P(\tau ,\stackrel{\sim}{\tau})N(\stackrel{\sim}{\tau})\mathrm{d}\stackrel{\sim}{\tau}.\end{array}$$(4)
Thus, to obtain the unknown age distribution N from the observed PDF C we need to solve the integral Eq. (4). We note however that in reality Eq. (4) holds only approximately, because the number of stars in the survey and in each MAMC is finite, and therefore both C(τ) and $P(\tau ,\stackrel{\sim}{\tau})$ are subject to stochastic deviations. These deviations limit our ability to find the solution for Eq. (4).
Let us consider as an example a case where the survey contains several monoage populations (e.g. a set of open clusters). This implies that $N(\stackrel{\sim}{\tau})=\sum _{i}{N}_{i}\delta (\stackrel{\sim}{\tau}{\tau}_{i})$, with N_{i} being the fractional weight of the ith population, which has log(age) τ_{i}. Following Eq. (4), we obtain that C(τ) = ∑_{i}N_{i}P(τ, τ_{i}), in other words, log(age) PDF for such a survey is a weighted sum of PDFs of monoage populations. In fact, the PDF from UniDAM C(τ) is defined over the grid of log(ages) (τ_{i}, i = 1, …n). Similarly, MAMC PDFs $P(\tau ,\stackrel{\sim}{\tau})$ are defined only for τ ∈ τ_{i}. Thus we can substitute the function $P(\tau ,\stackrel{\sim}{\tau})$ with a quadratic matrix of the size n × n and functions $N(\stackrel{\sim}{\tau})$ and C(τ) with vectors $\mathit{N}(\stackrel{\sim}{\tau})$ and C(τ) of the size n, and Eq. (4) can be rewritten as a system of linear equations:
$$\begin{array}{c}\hfill {C}_{j}={\displaystyle \sum _{i}}{P}_{j,i}{N}_{i}.\end{array}$$(5)
There is an obvious nonnegativity constraint for this system: N_{i} ≥ 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 monoage populations with ages τ_{i}.
3.2. Solving the main equation
In the system of linear equations (Eq. (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 nonnegativity constraint. It is possible to find a solution by means of the nonnegative least squares method (NNLS; see Lawson & Hanson 1995), which maximises the following function:
$$\begin{array}{c}\hfill {L}_{0}={\displaystyle \sum _{j=1}^{n}}{({C}_{j}{\displaystyle \sum _{i=1}^{n}}{P}_{\mathit{ji}}{N}_{i})}^{2}.\end{array}$$(6)
This will typically lead to a result in which only a few of the components of N will be nonzero, 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 first derivatives of the solution^{3}. We thus maximise the following function:
$$\begin{array}{c}\hfill L={\displaystyle \sum _{j=1}^{n}}{{\displaystyle ({C}_{j}\sum _{i=1}^{n}{P}_{\mathit{ji}}{N}_{i})}}^{2}\lambda \sum _{i=1}^{n}{\left(\frac{\mathrm{d}{N}_{i}}{\mathrm{d}\tau}\right)}^{2},\end{array}$$(7)
where λ is a regularization parameter. Solutions with narrow spikes will produce large absolute values of $\frac{\mathrm{d}{N}_{i}}{\mathrm{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}_{\mathit{ji}}{N}_{i}$.
We can use finite difference formulas to calculate $\frac{\mathrm{d}{N}_{i}}{\mathrm{d}\tau}$ from N_{i}, and rewrite:
$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}_{i}}{\mathrm{d}\tau}=\mathcal{T}\mathit{N},\end{array}$$(8)
where, 𝒯 is a Toeplitz matrix representation of the first derivative:
$$\begin{array}{c}\hfill \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\\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \phantom{\rule{0.166667em}{0ex}}& \vdots \\ \vdots & \phantom{\rule{0.166667em}{0ex}}& \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 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)\end{array}$$(9)
In this matrix, the ith row represents coefficients for the secondorder forward formula for the first derivative:
$$\begin{array}{c}\hfill f\prime ({x}_{i})=\frac{3}{2}f({x}_{i})+2f({x}_{i+1})\frac{1}{2}f({x}_{i+2}),\end{array}$$(10)
with the exception of the last two rows, which represent similar coefficients for the backward formula:
$$\begin{array}{c}\hfill f\prime ({x}_{i})=\frac{1}{2}f({x}_{i2})2f({x}_{i1})+\frac{3}{2}f({x}_{i}).\end{array}$$(11)
Thus, Eq. (5) can be rewritten as
$$\begin{array}{c}\hfill \mathit{C}\prime =P\prime \mathit{N},\end{array}$$(12)
where
$$\begin{array}{cc}\hfill P\prime & =\left(\begin{array}{c}P\\ \lambda \mathcal{T}\end{array}\right)\hfill \end{array}$$(13)
$$\begin{array}{cc}\hfill \mathit{C}\prime & ={(\mathit{C}0\cdots 0)}^{T}.\hfill \end{array}$$(14)
In this notation, maximization of L from Eq. (7) is equivalent to minimization of ${\parallel \begin{array}{c}\mathit{C}\prime P\prime \mathit{N}\end{array}\parallel}_{2}$, with a constraint that all components of N are nonnegative. Again, we can use NNLS to obtain the optimal vector N. The solution can be characterised by the sum of residuals $R={\parallel \begin{array}{c}\mathit{C}P\mathit{N}\end{array}\parallel}_{2}$, which indicates how close the log(age) PDF for the derived age distribution is to C(τ).
3.3. Searching for optimal λ value
The problem of the regularization approach is that the parameter λ has to be properly chosen in order to obtain the correct solution of the problem. Two extreme cases are λ = 0 and λ = ∞. In the first case no regularization is active, and Eq. (12) is reduced to Eq. (5). In that case the solution is the most precise, however it is not smooth. If λ = ∞, the regularization dominates the solution, such that N(τ) becomes a constant, equal to the mean value of C(τ) – this is the smoothest possible solution, which is very imprecise. In this section we describe an empirical method for obtaining the value of λ that gives the result we consider to be optimal.
First of all, we note that both P and 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 the PDF in a wide range of log(age) values, variations of the PDF at different values of log(age) τ 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 C vectors C_{q}, q = 1, 2, …5. Using this information, we can make use of the crossvalidation 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 the regularization parameter λ we can thus build a set of 25 equations similar to Eq. (12), with all possible combinations of P_{p} and C_{q}. Let us now focus on the combination with p = p_{0} and q = q_{0}. For this combination, we write Eq. (12) as:
$$\begin{array}{c}\hfill \mathit{C}{\prime}_{{q}_{0}}=\mathit{P}{\prime}_{{p}_{0}}\mathit{N},\end{array}$$(15)
and designate a solution of this equation as N_{p0, q0}. We can further build a matrix of residuals with respect to all combinations of C and P:
$$\begin{array}{c}\hfill {R}_{m,n,{p}_{0},{q}_{0}}={\parallel \begin{array}{c}{\mathit{C}}_{n}{\mathit{P}}_{m}{\mathit{N}}_{{p}_{0},{q}_{0}}\end{array}\parallel}_{2}.\end{array}$$(16)
We take the following expression as our quality parameter:
$$\begin{array}{c}\hfill {Q}_{{p}_{0},{q}_{0}}=\frac{1}{24}({\displaystyle \sum _{p}\sum _{q}}{R}_{m,n,{p}_{0},{q}_{0}}{R}_{{p}_{0},{q}_{0},{p}_{0},{q}_{0}}).\end{array}$$(17)
All R_{m, n, p0, q0} and Q_{p0, q0} are functions of the regularization parameter λ. We show an example of R_{m, n, p0, q0}(λ) and Q_{p0, q0}(λ) functions in Fig. 4. R_{p0, q0, p0, q0} decreases as λ decreases – with less regularization it is possible to obtain a more accurate solution of Eq. (12). For values of R_{m, n, p0, q0} for m ≠ p_{0} or n ≠ q_{0} the decreasing trend with decreasing λ stops or even reverses at some point. This happens because solutions N_{p0, q0} of Eq. (12) that are accurate for the combination P_{p0} and C_{q0} become too “specialised”, and are not as accurate for other combinations of P_{p} and C_{q}. As an optimal value of λ we take the point where the rapid decrease of Q_{p0, q0} (hereafter labelled simply as Q) with decreasing λ stops or slows down. To find this point, we first fit the following piecewise linear function to the Qfunction in loglog space:
$$\begin{array}{c}\hfill {Q}_{\mathrm{fit}}(log\lambda )=\{\begin{array}{c}a+c(log\lambda log{b}_{1}),\phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.166667em}{0ex}}\lambda \le {b}_{1}\hfill \\ a+d(log\lambda log{b}_{1}),\phantom{\rule{1em}{0ex}}\mathrm{if}\phantom{\rule{0.166667em}{0ex}}{b}_{1}<\lambda \le {b}_{2}\hfill \\ a+d(log{b}_{2}log{b}_{1}),\phantom{\rule{1em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\mathrm{if}\phantom{\rule{0.166667em}{0ex}}\lambda >{b}_{2}\hfill \end{array}.\end{array}$$(18)
Fig. 4. Illustration for the choice of optimal λ. The plot shows all residual functions R_{m, n} (grey lines, see Eq. (16)) and highlights the R_{p0, q0} in red and Q (Eq. (17)) in blue. The fitted piecewise linear function Q_{fit} is shown in orange (see Eq. (18)). Vertical lines indicate the turning points of the fitted piecewise linear function Q_{fit} (orange dashed lines) and the finally adopted λ value at the intersection point of Q and Q_{fit} (blue dashed line). 
Function Q_{fit} has three linear segments, with transitions at values of λ = b_{1} and λ = 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 λ_{est} = b_{1}. In many cases the turn of the Q function is not as sharp as that of the fitted function, and λ_{est} is an overestimation of the optimal value of λ. To correct for this, we take the first point to the left of λ_{est}, where Q and Q_{fit} intersect, as our final estimate, λ_{f}, as shown in Fig. 4. This estimate in some cases tends to produce a noisier “overfitted” result, while the λ_{est} tends to produce smoother “underfitted” one. On average, however, λ_{f} produces better results than λ_{est}. Our choice of λ_{f} is further validated with tests described below in Sect. 5.
We can repeat the above procedure for every possible combination of p_{0} and q_{0}, thus obtaining 25 different values of λ_{f} and corresponding solutions N_{p0, q0} of Eq. (15). As the final solution we can take the mean of N_{p0, q0}, and as a measure of the uncertainty – their standard deviation. For a more detailed uncertainty analysis see Sect. 5.3.
Fig. 5. Input log(age) distributions and log(age) PDFs produced with UniDAM for different base surveys. Plots are for six cases considered in this work. For the “age” input the true PDF was rescaled to make other lines visible. 
Fig. 6. Results of the age inversion for box test data. Upper panels are for the smaller MAMC (n_{mock} = 250), and lower panels are for the larger MAMC (n_{mock} = 5000). Left panels are for the smaller simulated survey (n_{sim} = 1000), and right panels are for the larger simulated survey (n_{sim} = 25 000). The solid blue line shows the input log(age) distribution, and the solid red line is the result of the inversion with 68 (dark shading) and 95 (light shading) percent confidence intervals. Input log(age) PDFs for five realizations of the simulated catalogue are plotted with grey lines and the black dashed line shows the log(age) PDF inferred from the inversion result. The inset in each panel shows the crossvalidation curves (red) and piecewiselinear fit (blue) for one realization of the simulated catalogue (see Sect. 3.3 for the description of the crossvalidation curve). 
4. Monoage mocks construction
In order to properly recover the age distribution for a given survey, which we refer to as the 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 is close to that of the base survey, and at the same time retain age information. To achieve this, we follow the procedure below.
We build MAMCs using PARSEC models (Bressan et al. 2012) to simulate stars. In order to mimic the random observational scatter, instead of taking each model, we take a sample of 25 models with random normally distributed perturbations added in T_{eff}, log g, and [Fe/H]. To assign a proper amplitude of perturbations, we randomly select uncertainties σ_{P} = {σ_{T}, σ_{log g}, σ_{[Fe/H]}} from the base survey data for stars that have P = {T_{eff}, log g and [Fe/H]} close to the model, such that the difference between the physical parameters for the model and those for the base survey star is smaller than the mean uncertainty in the respective parameter for the base survey P_{model} − P_{star} < mean(σ_{P}). This provides a way to reconstruct the scatter and the systematic variations of σ_{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\stackrel{\phantom{\rule{0.333333em}{0ex}}\text{m}}{.}025$, to reproduce typical photometric uncertainties in 2MASS and AllWISE.
The set of models is then truncated in the T_{eff}, log g, [Fe/H] space to the footprint of the base survey in that space, and models outside of that footprint are excluded. In order to build MAMC from that set of models, we sample a predefined number of models n_{mock} as follows. We bin stars from the base survey in log g − [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_{eff}, log g, [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 an MAMC that will not resemble the base survey, making inversion impossible. We decided to use [Fe/H] as a first parameter for the binning, 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 between main sequence stars and giants, which have very different contributions to the log(age) PDF. Such a distinction is impossible to make with T_{eff}.
The above procedure is repeated for each value of log(age) τ over the considered grid: 6.61 ≤ τ ≤ 10.13 dex with a step of 0.02 dex or age between approximately 4 × 10^{6} and 13.5 × 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 λ and to determine the influence of the survey size n_{sim} and MAMC size n_{mock} on the precision and the accuracy of the inversion result. We simulate a survey with predefined age distribution N_{input} as a concatenation of MAMC simulated in a way described in Sect. 4, with the size of ith MAMC defined by N_{input, i}. We then obtain with UniDAM the log(age) PDF for the simulated survey and apply the inversion as described in Sect. 3.2. We can then 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 RAVEon (Casey et al. 2017). 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, and therefore 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. We 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 being 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).
Figures 6 and 7 illustrate that the result of the inversion is almost insensitive to n_{mock} for small values of n_{sim}, as expected (see discussion in Sect. 5.3). On the other hand, for n_{sim} ≥ 25 000, MAMC size n_{mock} plays a larger role. It is also clear that an increase in n_{sim} and n_{mock} increases the sensitivity of the inversion to rapid changes in log(age) distribution. For example, in the case of the wave test, n_{mock} = 250 and n_{sim} = 25 000 is enough to properly recover the second peak of log(age) distribution at τ ≈ 9.6 dex, though we need to increase n_{mock} to 5000 in order to properly recover the first, much narrower peak at τ ≈ 8.9 dex. For all simulations, observed (black lines) and predicted (grey lines) log(age) PDFs become almost indistinguishable for n_{mock} = 5000 and n_{sim} = 25 000.
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 Sect. 3.3 we give a method for obtaining the optimal value of the regularization parameter λ. This parameter and a corresponding solution of Eq. (12) can be obtained for every combination of the five survey log(age) PDFs C_{q} and five MAMC log(age) PDF sets P_{p}. Hence we can have 25 different solutions for Eq. (15), with the differences between them being due to statistical variations in P and 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) τ 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 heavytailed distribution.
Fig. 8. Distributions of the fraction F(ω) of values of τ with true error δ(τ) < ωσ(τ), where σ(τ) is the estimated uncertainty. Shown are distributions for ω = 1 and ω = 3 for simulated data. Vertical lines and numbers represent median values of each distribution. 
Fig. 9. Results of age inversion for several surveys. The solid red line is the computed underlying age distribution with 68 (dark shading) and 95 (light shading) percent confidence intervals. Input log(age) PDFs for five parts of each survey are plotted with grey lines and the black dashed line shows the log(age) PDF inferred from the inversion result. The inset in each panel shows the crossvalidation curves (red) and piecewiselinear fit (blue) for one part of the survey (see Sect. 3.3 for the description of the crossvalidation curve). 
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 σ(τ) and the true error δ(τ) (the true error is defined as the absolute difference between input and output log(age) distributions). The fraction F(ω) of values of τ with δ(τ) < ωσ(τ) can be calculated for each result. In Fig. 8 we show distributions of F(ω) for ω = 1 and ω = 3. If σ(τ) is a correctly determined Gaussian uncertainty, then F(1) = 0.68 and F(3) = 0.997, as expected for 1 and 3σ 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σ(τ) effectively a 2sigma rather than 3sigma confidence interval estimate. This is because uncertainties have a heavytail nonGaussian distribution. There are no clear trends visible for F(ω) values with n_{sim} and n_{mock}, at least within the considered ranges of these parameters.
We note here that the scatter of F(ω) around the median values is high, which means that values of σ(τ) can be seen only as an approximation of the uncertainty.
One would expect the fractional uncertainty values to decrease as (n_{sim}n_{mock})^{−1/2}. In reality, their relation is not as strong due to complex correlations between the real age distribution and the corresponding log(age) PDF. Our estimate shows that the fractional uncertainty scales as ${n}_{\mathrm{sim}}^{0.4}{n}_{\mathrm{mock}}^{0.15}$. This means that one has to increase n_{mock} by an order of magnitude to improve the result by about 30%. Importantly, for a given n_{mock} there seems to exist a maximum value of n_{sim}, above which the fractional uncertainty does not decrease with n_{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_{sim} there is a maximum value of n_{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_{mock} for a given n_{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_{mock} = 5000 is sufficient for surveys as large as n_{sim} = 50 000. This is still below typical sizes of spectroscopic surveys. However, if we want to slice the survey into parts, to trace age–metallicity relations, for example, as is done below in Sect. 6.2, the size of each part will be on the order of n_{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 show below, systematic errors 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. 2017), RAVEon (Casey et al. 2017), 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 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 [Fe/H]. Within the metallicity bin for a magnitudelimited 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 also manifest themselves 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_{sim} is larger than those considered in the simulation.
6.2. Inversion of monometallicity 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 RAVEon, LAMOST (Luo et al. 2015), DR4, and APOGEE stars, splitting them into 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 RAVEon stars. Inputs (log(age) PDFs) are typically broad and for metalpoor populations show an extra narrow peak at log(age) τ ≈ 9, which is associated with redclump 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) τ = 9.37 for the highest metallicity bin [Fe/H] = 0.4 dex to the maximum possible log(age) of τ = 10.13 for [Fe/H] < −0.4 dex. At very low metallicities, [Fe/H] ≤ −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 metalpoor bins and thus a reduced age resolution (see test results with low n_{sim} values in Figs. 6 and 7).
Fig. 10. Results of age inversion for RAVEon metallicity slices: log(age) PDFs as produced by UniDAM (left) and underlying age distribution computed through the inversion (right). Plots for various metallicities are offset in vertical direction for visual purposes, with numbers between plots indicating [Fe/H] values. 
For LAMOST, we use a constraint of T_{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 RAVEon 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 RAVEon 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 [Fe/H] ≥ −0.2 dex: one at around τ ≈ 9.35 and one at τ ≈ 9.83. The latter peak becomes dominant and shifts to higher τ values, as [Fe/H] increases. Similarly, for APOGEE a secondary peak is visible in the underlying age distribution for metallicities [Fe/H] ≥ −0.2 dex at around τ ≈ 8.9. In both cases, the result is that the mean age for bins with [Fe/H] ≥ +0.2 dex is higher than that for solarlike metallicities. This does not seem to be physical – we expect more metalrich stars to be systematically younger. This might be attributed to the limitations of the method listed above in the Sect. 6.1. An alternative explanation is that the effect is physical, and is caused by the fact that metal rich ([Fe/H] ≥ +0.2 dex) stars 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 metalrich 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 us to trace age–metallicity trends that are much less obvious if log(age) PDFs are considered. It also allows structures like lowage spikes to be removed, which 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 Sect. 6.2 allow us to study metallicity–age and age–metallicity relations for RAVEon LAMOST and APOGEE. Both relations can be defined through a twodimensional distribution function C(τ, [Fe/H]) = C_{[Fe/H]}(τ)N([Fe/H]), where C_{[Fe/H]}(τ) is the log(age) PDF in the bin with metallicity [Fe/H], and N([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:
$$\begin{array}{c}\hfill \tau ([\mathrm{Fe}/\mathrm{H}])=\frac{{\int}_{{\tau}_{\mathrm{min}}}^{{\tau}_{\mathrm{max}}}C(\tau ,[\mathrm{Fe}/\mathrm{H}])\mathrm{d}\phantom{\rule{0.166667em}{0ex}}\tau}{{\tau}_{\mathrm{max}}{\tau}_{\mathrm{min}}},\end{array}$$(19)
and similarly metallicity as a function of log(age):
$$\begin{array}{c}\hfill [\mathrm{Fe}/\mathrm{H}](\tau )=\frac{{\int}_{{[\mathrm{Fe}/\mathrm{H}]}_{\mathrm{min}}}^{{[\mathrm{Fe}/\mathrm{H}]}_{\mathrm{max}}}C(\tau ,[\mathrm{Fe}/\mathrm{H}])\mathrm{d}\phantom{\rule{0.166667em}{0ex}}[\mathrm{Fe}/\mathrm{H}]}{{[\mathrm{Fe}/\mathrm{H}]}_{\mathrm{max}}{[\mathrm{Fe}/\mathrm{H}]}_{\mathrm{min}}}.\end{array}$$(20)
Here, (τ_{min}, τ_{max}) is the considered range in log(age) and ([Fe/H]_{min}, [Fe/H]_{max}) is the range in metallicity.
In Figs. 13–15 we show τ([Fe/H]) with blue solid and [Fe/H](τ) with blue dotted lines. In these plots, we continue to use solid lines for τ([Fe/H]) relations and dotted lines for [Fe/H](τ) relations. These two lines are nearly orthogonal, which is the consequence of the C(τ, [Fe/H]) distribution being very broad.
Fig. 13. Age–metallicity (τ([Fe/H])) relations for nearby APOGEE stars (orange line, Feuillet et al. 2018) and for RAVEon stars (solid lines, this work). Blue lines are values derived from log(age) PDFs, red lines are values derived from the underlying age distributions obtained from inversion, and black lines are values derived from Minchev et al. (2013) chemical evolution models. Dotted lines show mean metallicity as a function of log(age) for the same data. 
We can now replace log(age) PDFs C(τ, [Fe/H]) with the underlying age distribution from the inversion N(τ, [Fe/H]) in Eqs. (19) and (20). The results are shown in Figs. 13–15 with red solid and dotted lines. These two red lines are a lot closer to each other, which is a sign of a much tighter C(τ, [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 the determination of 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–15 with orange solid lines. We note that the data presented in Feuillet et al. (2018), as well as our results for τ([Fe/H]), give age distribution and mean age in each metallicity bin, thus measuring age as a function of metallicity τ([Fe/H]). At the same time, models of chemical evolution typically focus on the metallicity as a function of age [Fe/H](τ) (see for example Fig. 4 in Minchev et al. 2013), which is strictly speaking a different function. Functions τ([Fe/H]) and [Fe/H](τ) can be close to each other only if the twodimensional distribution C(τ, [Fe/H]) is tight. Hence in a general case care must be taken in comparing τ([Fe/H]) and [Fe/H](τ), as this comparison might lead to inaccurate conclusions.
If we compare the age–metallicity distributions from log(age) PDFs C(τ, [Fe/H]) with those obtained from the underlying age distributions computed through the inversion, N(τ, [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 thickdisk and halo stars that are systematically older.
We can compare our age–metallicity and metallicity–age trends with those predicted by chemical evolution models. For this comparison, we use the model described in Minchev et al. (2013). We added an additional smoothing to the twodimensional age–metallicity distribution predicted by this model. The smoothing scale was chosen to be close to the typical uncertainty in age (2 Gyr) and metallicity (0.1 dex). We then calculated τ([Fe/H]) and [Fe/H](τ) for this distribution and show it in Figs. 13 and 14 with black solid lines and Fig. 15 with black dotted lines. We note that mean age τ starts to increase as a function of [Fe/H] for [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 metalrich stars in the solar vicinity, which are formed in the inner Galaxy and need time to migrate outwards. Mean ages for metalpoor stars are also systematically larger for inversion results than those predicted by models. This might be caused by the fact that UniDAM accommodates stellar ages up to the age of the Universe (≈13.5 Gyr), while the maximum stellar age in the Minchev et al. (2013) model is 11.175 Gyr (before the smoothing was applied). In general, on all three plots τ([Fe/H]) and [Fe/H](τ) 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 for revealing the underlying distribution of a large (N > 10^{3}) ensemble of stars in log(age) from their cumulative log(age) PDFs 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 on 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 monoage populations that will be distributed closely to the survey stars, which is needed for a reliable age inversion.
We also apply the age inversion method to monometallicity slices of RAVEon, 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 RAVEon, 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, highprecision APOGEE sample and to those predicted by chemical evolution models. The number of stars in both RAVEon and LAMOST surveys, at least for metallicities close to solar, is high enough to make further divisions in, for example, galactic latitude, distance, and α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 providing the underlying age distribution of such large survey data sets as a whole and as functions of, for example, metallicity, location in the Galaxy, and αelements abundances.
One could be tempted to use the vast amount of data contained in Gaia to derive age distributions from photometry alone, as was done, for example, in Dolphin (2013). However, this might be difficult given the lack of metallicity information in photometric data, which is essential for unambiguous age determination. In any case, such a study would differ substantially from the present one, and the new method to solve Eq. (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 was done in Mints & Hekker (2018) and Mints (2018). However, this would require a more complex simulation strategy that takes into account proper parallax uncertainty distribution. In addition, Gaia parallax zeropoint offset (see, e.g. Leung & Bovy 2019) would have to be treated properly, as it can inflict a systematic age bias.
Acknowledgments
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 communitydeveloped core Python package for Astronomy (Astropy Collaboration 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 LeibnizInstitut 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 (ERCStG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST0908326); 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 and 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.ravesurvey.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. SDSSIV acknowledges support and resources from the Center for HighPerformance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSSIV 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, HarvardSmithsonian 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), MaxPlanckInstitut für Astronomie (MPIA Heidelberg), MaxPlanckInstitut für Astrophysik (MPA Garching), MaxPlanckInstitut 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 MultiObject 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.
References
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Casey, A. R., Hawkins, K., Hogg, D. W., et al. 2017, ApJ, 840, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Cutri, R. M., Wright, E. L., Conrow, T., et al. 2014, VizieR Online Data Catalog: II/328 [Google Scholar]
 Dalton, G., Trager, S., Abrams, D. C., et al. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, SPIE Conf. Ser., 9147, 91470L [Google Scholar]
 de Jong, R. S., Barden, S. C., BellidoTirado, O., et al. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, SPIE Conf. Ser., 9908, 99081O [Google Scholar]
 Dolphin, A. E. 2013, ApJ, 775, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Feuillet, D. K., Bovy, J., Holtzman, J., et al. 2016, ApJ, 817, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Feuillet, D. K., Bovy, J., Holtzman, J., et al. 2018, MNRAS, 477, 2326 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Lawson, C. L., & Hanson, R. J. 1995, Solving Least Squares Problems (Philadelphia: SIAM) [Google Scholar]
 Leung, H. W., & Bovy, J. 2019, MNRAS, 489, 2079 [NASA ADS] [CrossRef] [Google Scholar]
 Luo, A.L., Zhao, Y.H., Zhao, G., et al. 2015, Res. Astron. Astrophys., 15, 1095 [NASA ADS] [CrossRef] [Google Scholar]
 Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Martell, S., Sharma, S., Buder, S., et al. 2017, MNRAS, 465, 3203 [NASA ADS] [CrossRef] [Google Scholar]
 Minchev, I., Chiappini, C., & Martig, M. 2013, A&A, 558, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minchev, I., Anders, F., RecioBlanco, A., et al. 2018, MNRAS, 481, 1645 [NASA ADS] [CrossRef] [Google Scholar]
 Minchev, I., Matijevic, G., Hogg, D. W., et al. 2019, MNRAS, 487, 3946 [NASA ADS] [Google Scholar]
 Mints, A. 2018, ArXiv eprints [arXiv:1805.01640] [Google Scholar]
 Mints, A., & Hekker, S. 2017, A&A, 604, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mints, A., & Hekker, S. 2018, A&A, 618, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mints, A., & Hekker, S. 2019, A&A, 621, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Queiroz, A. B. A., Anders, F., Santiago, B. X., et al. 2018, MNRAS, 76, 2556 [NASA ADS] [CrossRef] [Google Scholar]
 Scipy Team 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/ [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Soderblom, D. R. 2010, ARA&A, 48, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, M. B. 2005, in Astronomical Data Analysis Software and Systems XIV, eds. P. Shopbell, M. Britton, & R. Ebert, ASP Conf. Ser., 347 [Google Scholar]
 Tucci Maia, M., Ramírez, I., Meléndez, J., et al. 2016, A&A, 590, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wu, Y.Q., Xiang, M.S., Zhang, X.F., et al. 2017, Res. Astron. Astrophys., 17, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Xiang, M., Liu, X., Shi, J., et al. 2017, ApJS, 232, 2 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Examples of log(age) PDF for various stars. We note that PDFs often show nonGaussian features: heavy tails, truncation, multiple modes. 

In the text 
Fig. 2. Different representations of log(age) distribution for public spectroscopic surveys. See Sect. 2 for a discussion. 

In the text 
Fig. 3. Examples of log(age) PDFs derived with UniDAM (same as red lines in Fig. 2) for simulated monoage populations with log g and [Fe/H] distributions taken from APOGEE (left) and RAVEon (right). Vertical lines indicate the true age for each population. 

In the text 
Fig. 4. Illustration for the choice of optimal λ. The plot shows all residual functions R_{m, n} (grey lines, see Eq. (16)) and highlights the R_{p0, q0} in red and Q (Eq. (17)) in blue. The fitted piecewise linear function Q_{fit} is shown in orange (see Eq. (18)). Vertical lines indicate the turning points of the fitted piecewise linear function Q_{fit} (orange dashed lines) and the finally adopted λ value at the intersection point of Q and Q_{fit} (blue dashed line). 

In the text 
Fig. 5. Input log(age) distributions and log(age) PDFs produced with UniDAM for different base surveys. Plots are for six cases considered in this work. For the “age” input the true PDF was rescaled to make other lines visible. 

In the text 
Fig. 6. Results of the age inversion for box test data. Upper panels are for the smaller MAMC (n_{mock} = 250), and lower panels are for the larger MAMC (n_{mock} = 5000). Left panels are for the smaller simulated survey (n_{sim} = 1000), and right panels are for the larger simulated survey (n_{sim} = 25 000). The solid blue line shows the input log(age) distribution, and the solid red line is the result of the inversion with 68 (dark shading) and 95 (light shading) percent confidence intervals. Input log(age) PDFs for five realizations of the simulated catalogue are plotted with grey lines and the black dashed line shows the log(age) PDF inferred from the inversion result. The inset in each panel shows the crossvalidation curves (red) and piecewiselinear fit (blue) for one realization of the simulated catalogue (see Sect. 3.3 for the description of the crossvalidation curve). 

In the text 
Fig. 7. Same as Fig. 6, but for wave input. Line styles and colours are as in Fig. 6. 

In the text 
Fig. 8. Distributions of the fraction F(ω) of values of τ with true error δ(τ) < ωσ(τ), where σ(τ) is the estimated uncertainty. Shown are distributions for ω = 1 and ω = 3 for simulated data. Vertical lines and numbers represent median values of each distribution. 

In the text 
Fig. 9. Results of age inversion for several surveys. The solid red line is the computed underlying age distribution with 68 (dark shading) and 95 (light shading) percent confidence intervals. Input log(age) PDFs for five parts of each survey are plotted with grey lines and the black dashed line shows the log(age) PDF inferred from the inversion result. The inset in each panel shows the crossvalidation curves (red) and piecewiselinear fit (blue) for one part of the survey (see Sect. 3.3 for the description of the crossvalidation curve). 

In the text 
Fig. 10. Results of age inversion for RAVEon metallicity slices: log(age) PDFs as produced by UniDAM (left) and underlying age distribution computed through the inversion (right). Plots for various metallicities are offset in vertical direction for visual purposes, with numbers between plots indicating [Fe/H] values. 

In the text 
Fig. 11. Same as Fig. 10, but for the LAMOST survey. 

In the text 
Fig. 12. Same as Fig. 10, but for the APOGEE survey. 

In the text 
Fig. 13. Age–metallicity (τ([Fe/H])) relations for nearby APOGEE stars (orange line, Feuillet et al. 2018) and for RAVEon stars (solid lines, this work). Blue lines are values derived from log(age) PDFs, red lines are values derived from the underlying age distributions obtained from inversion, and black lines are values derived from Minchev et al. (2013) chemical evolution models. Dotted lines show mean metallicity as a function of log(age) for the same data. 

In the text 
Fig. 14. Same as Fig. 13, now for LAMOST stars. 

In the text 
Fig. 15. Same as Fig. 13, now for APOGEE stars. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.