Issue 
A&A
Volume 604, August 2017



Article Number  A108  
Number of page(s)  20  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201630090  
Published online  22 August 2017 
A Unified tool to estimate Distances, Ages, and Masses (UniDAM) from spectrophotometric data^{⋆}
^{1} Max Planck Institute for Solar System Research, JustusvonLiebigWeg 3, 37077 Göttingen, Germany
email: mints@mps.mpg.de
^{2} Stellar Astrophysics Centre, Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, 8000 Aarhus C, Denmark
Received: 18 November 2016
Accepted: 19 April 2017
Context. Galactic archaeology, the study of the formation and evolution of the Milky Way by reconstructing its past from its current constituents, requires precise and accurate knowledge of stellar parameters for as many stars as possible. To achieve this, a number of large spectroscopic surveys have been undertaken and are still ongoing.
Aims. So far consortia carrying out the different spectroscopic surveys have used different tools to determine stellar parameters of stars from their derived effective temperatures (T_{eff}), surface gravities (log g), and metallicities ([Fe/H]); the parameters can be combined with photometric, astrometric, interferometric, or asteroseismic information. Here we aim to homogenise the stellar characterisation by applying a unified tool to a large set of publicly available spectrophotometric data.
Methods. We used spectroscopic data from a variety of large surveys combined with infrared photometry from 2MASS and AllWISE and compared these in a Bayesian manner with PARSEC isochrones to derive probability density functions (PDFs) for stellar masses, ages, and distances. We treated PDFs of preheliumcore burning, heliumcore burning, and post heliumcore burning solutions as well as different peaks in multimodal PDFs (i.e. each unimodal subPDF) of the different evolutionary phases separately.
Results. For over 2.5 million stars we report mass, age, and distance estimates for each evolutionary phase and unimodal subPDF. We report Gaussian, skewed, Gaussian, truncated Gaussian, modified truncated exponential distribution or truncated Student’s tdistribution functions to represent each subPDF, allowing us to reconstruct detailed PDFs. Comparisons with stellar parameter estimates from the literature show good agreement within uncertainties.
Conclusions. We present UniDAM, the unified tool applicable to spectrophotometric data of different surveys, to obtain a homogenised set of stellar parameters.
Key words: stars: distances / Galaxy: stellar content / stars: fundamental parameters
The unified tool and the tables with results are available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/604/A108
© ESO, 2017
1. Introduction
The Milky Way Galaxy is a unique object to test our understanding of stellar evolution, galaxy formation, and cosmology. For this test a detailed map of our Galaxy, including bulge, disk, halo, spiral structure, and streams formed by recent mergers, is required. Through an analysis of the Galaxy, we can learn how our Galaxy has formed, evolved, and how it interacts with its surroundings. To build such a map we need to find the distribution of stars in their positions, velocities, chemical compositions, and ages throughout the Galaxy. These parameters can be measured with different kinds of observations, such as astrometry, photometry, spectroscopy, and asteroseismology.
Astrometric observations provide stellar positions, proper motions and, through parallaxes, distances. These kind of data have been available for decades (see e.g. Woolley et al. 1970; Gliese & Jahreiß 1995; Perryman et al. 1997). In the nearest future Gaia (Gaia Collaboration 2016) will vastly increase the precision and amount of such information. The first Gaia data release (Lindegren et al. 2016) already provides proper motions and parallaxes for about two million stars, although the precision of parallaxes in this sample limits their application (see e.g. Stassun & Torres 2017). In the next data releases Gaia will provide highprecision parallaxes and proper motions for hundreds of millions of stars, vastly increasing our knowledge of the Galaxy.
Another rich source of data is spectroscopy, which can provide radial velocities, chemical compositions as well as the effective temperature T_{eff} and the surface gravity log g. These data can be used to derive stellar ages and distances (see below). A growing number of large spectroscopic surveys, such as RAdial Velocity Experiment (RAVE; Kunder et al. 2017), Large Sky Area MultiObject Fibre Spectroscopic Telescope (LAMOST) surveys (Luo et al. 2015), Apache Point Observatory Galactic Evolution Experiment (APOGEE; Ahn et al. 2014), Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), and GaiaESO (Gilmore et al. 2012) provide rich spectroscopic information for millions of stars.
Photometric surveys can be used in two ways. Stromgren or WashingtonDDO51 photometry can be used to estimate stellar parameters such as the effective temperature T_{eff}, surface gravity log g, and metallicity [ Fe / H ] (see Casagrande et al. 2011), or to separate giants from dwarfs for spectroscopic followup (Majewski et al. 2000). Otherwise, broadband photometry is commonly used as a supplement to spectroscopic data to infer stellar distances.
Asteroseismology is a relatively young and very promising method of exploring stars. For lowmass dwarfs, subgiants, and red giant stars asteroseismology can provide a direct measure of mean density and surface gravity. The surface gravities measured by asteroseismic methods have much higher precision than spectroscopic methods. In case the effective temperature T_{eff} or luminosity L are also measured it is possible to obtain stellar mass and radius from the asteroseismic observables. When compared with models stellar ages can also be determined from asteroseismology. COnvection ROtation and planetary Transits (CoRoT; Baglin et al. 2006), Kepler (Borucki et al. 2010), and K2 (Howell et al. 2014) observations provide such asteroseismic data. These datasets provide highprecision data on small patches of the sky and also have proved to be a perfect sample for the calibration of large spectroscopic surveys (see e.g. Hawkins et al. 2016; Wang et al. 2016b; Tayar et al. 2017). Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) and PLAnetary Transits and Oscillations of stars (PLATO; Rauer et al. 2014) space missions are scheduled to be launched in 2018 and 2024, respectively, and will vastly increase the number of stars with asteroseismic data in the coming years.
Fig. 1 Probability distribution function for log(age) (top row), mass (middle row), and distance modulus (bottom row) for two stars from our catalogue. Different evolutionary stages are indicated with different colours (see legend). The stages are mainsequence stars and giants ascending the red giant branch (precorehelium burning, stage I), corehelium burning stars (stage II) and asymptotic giant branch stars (postcorehelium burning, stage III). The histograms represent model distributions and solid lines indicate our fits for individual USPDFs. 
Stellar ages and distances remain among the most challenging parameters to measure. A comprehensive list of age determination methods is given in Soderblom (2010). For a number of stars ages can be derived from asteroseismic observations (Silva Aguirre & Serenelli 2016) or from carbon and nitrogen abundances (Martig et al. 2016). When these data are not available, a typical approach is to compare the parameters directly derived from spectroscopic measurements, which we designate as observed parameters, such as the effective temperature T_{eff}, surface gravity log g, and metallicity [ Fe / H ], to a grid of stellar models. A model or a set of models that have their parameters close to observed parameters give estimates of ages, masses, and absolute magnitudes M_{λ} of a star. Then by comparing the absolute magnitudes to visible magnitudes m_{λ} from photometric surveys we can estimate distances to stars. An overview of this approach is given by Jørgensen & Lindegren (2005) and Burnett et al. (2011).
Proper application of this approach requires some care. First, the transformation from observed (T_{eff},log g, and [ Fe / H ]) to stellar parameters (age and mass) and distance is often degenerate, with the same observables giving two or more possible combinations of stellar parameters. This degeneracy can in some rare cases be resolved when additional observables are available, for example from asteroseismology. Second, the interstellar extinction needs to be accounted for in distance estimations. Extinction values can be taken from external sources or can be derived from observables. Both ways have their advantages and disadvantages. We discuss this in Sect. 4. Third, observed parameters have their uncertainties and correlations that have to be propagated to uncertainties in stellar parameters.
In the literature a number of methods based on the comparison of observed parameters from spectroscopic and photometric surveys with models to estimate distances and other stellar parameters were proposed and used recently. Here we briefly discuss some of them, and how they deal with the issues stated above.
GCS.The GenevaCopenhagen Survey (GCS; Casagrande et al. 2011) team exploited the advantage of having Hipparcos parallaxes for the majority of their objects; this facilitated the calculation of absolute magnitudes for each star. Casagrande et al. (2011) used a Bag of Stellar Tracks and Isochrones (BASTI; Pietrinferni et al. 2009, and references therein) and PAdova and TRieste Stellar Evolution Code (PARSEC; Bressan et al. 2012) isochrones to select models that have T_{eff}, absolute Johnson V magnitude and metallicity close to the observed ones for each star. Applying a Bayesian scheme described in Jørgensen & Lindegren (2005) to selected models, Casagrande et al. derived masses and ages of stars. They used a flat prior on ages and a Salpeter initial mass function (IMF) as a prior for masses.
RAVE.There is a series of papers on distance estimations for stars in the RAVE survey (DR5 is described in Kunder et al. 2017). Breddels et al. (2010) proposed a method for distance estimation for RAVE stars based on a comparison of observed T_{eff},log g, metal abundance [ M / H ] , and colour (J − K_{s}) with Y^{2} models (Demarque et al. 2004). For each star 5000 realisations of observed parameters were sampled from a Gaussian distribution with dispersions equal to the measured uncertainties and for each realisation a closest model was selected. These authors took an average of the model parameters measured in all realisations to derive an absolute J magnitude of the star M_{J}. The difference between the derived absolute magnitude and visible J magnitude from Two Micron All Sky Survey (2MASS) gives a distance. Extinction was ignored in this work. This approach is limited by the fact that it does not take into account the inhomogeneity of models in the T_{eff} − log g plane, effectively increasing the weight for short evolutionary stages and decreasing it for longer ones. This issue was solved in Zwitter et al. (2010) by weighting models with a weight proportional to age and mass range represented by each model. A likelihood depending on the difference between observed and model T_{eff} and log g was also added. Other important changes were applied, including a change from Y^{2} to PARSEC (Bressan et al. 2012) isochrones, the addition of a prior on mass (assuming Chabrier 2003, IMF), and the application of a volume correction. Zwitter et al. calculated an absolute J magnitude as a weighted mean of the absolute magnitudes derived from luminosities of the models. The difference between the visible J magnitude from 2MASS and the absolute magnitude gives the distance modulus for each star. As in Breddels et al. (2010), extinction was ignored.
Binney et al. (2014) further developed the above method by adding priors from the Galactic structure; they provide priors on age, metallicity, and positions from halo, thin, and thick disk models. A kinematic correction (see Schönrich et al. 2012) was also applied. Extinction was included into distance calculations. An exponential prior on the value of ln(A_{V}) was imposed with the extinction value at infinity A_{V∞}(b,l) taken from Schlegel et al. (1998). The extinction at a given distance was calculated as , where ρ(s) is the density of extincting material along the line of sight, taken from the model of the Galaxy (see the Eq. (10) in Binney et al. 2014). This is so far the most advanced method and it was applied with minor modifications to LAMOST data as well (see below). Distance moduli (but not ages and masses) were recalculated with the same method for RAVE DR5. Binney et al. (2014) solved the problem of multimodal probability distribution functions (PDFs) for the distance modulus by fitting a Gaussian mixture model to it with up to three Gaussians. This approach works fine in most cases. However, as we illustrate below in Fig. 1 it cannot be applied to mass and log(age) PDFs because they can be skewed or truncated, a shape which is hard to fit with a small set of Gaussians. Truncated shapes of the PDF arise from a limited range of allowed masses and log(age)s. For Binney et al. (2014) limits are imposed by an age prior; see their Eqs. (3)–(5).
APOKASC.Rodrigues et al. (2014) applied Bayesian methods to estimate distances and extinctions for approximately 2000 red giant stars from the joint APOGEE and Kepler Asteroseismic Science Consortium (APOKASC) sample (Pinsonneault et al. 2014), which is a part of APOGEE (Ahn et al. 2014), covering the Kepler field of view. They supplemented spectroscopic parameters T_{eff} and [ M / H ] with asteroseismic data from Kepler. As alluded to before, from asteroseismic values Δν and ν_{max} and knowing T_{eff}, it is possible to derive an estimate of stellar radius R and mass M, using scaling relations from Kjeldsen & Bedding (1995), i.e.
This puts more constraints on stellar models, thus increasing the precision of stellar parameters and distance determinations. Using PARSEC isochrones, Rodrigues et al. (2014) built PDFs for stellar parameters (mass, radius, and surface gravity) and stellar absolute magnitudes. The latter were then combined with photometric data from Sloan Digital Sky Survey (SDSS), 2MASS, and Widefield Infrared Survey Explorer (WISE) to be converted to the PDFs of distance modulus μ_{d} and extinction A_{K}. The mode and 68% confidence intervals of the PDFs were calculated for both distance and extinction. Rodrigues et al. (2014) noted that over onethird of stars in their sample have bimodal PDFs. Bimodal PDFs were treated in the same way as singlepeaked PDFs. Using the mode allows one to select the highest peak of the PDF and other peaks only show themselves by broadening of confidence intervals. Only distance estimates were published by Rodrigues et al. (2014).
LAMOST.The LAMOST team is also working on estimating distances to stars from spectroscopic data. First and second public data releases of the project include spectral properties for about one and two million stars, respectively (Luo et al. 2015).
Carlin et al. (2015) used a Bayesian approach to derive distances from LAMOST DR1 data combined with 2MASS photometry. They used the Dartmouth Stellar Evolution Database (Dotter et al. 2008) and a Bayesian technique similar to that by Burnett et al. (2011) to derive the PDF of the absolute magnitude for each star. This was then converted to distances using 2MASS photometry. Interstellar extinction was ignored in this work. Carlin et al. (2015) performed a comparison with RAVE distances to test their method. The derived distances are systematically smaller by 12% than those derived by Zwitter et al. (2010) with 16% spread. Given the precision of the LAMOST data, they derived distance uncertainties to be on the order of 40%.
Wang et al. (2016a) applied the Bayesian approach from Binney et al. (2014) to derive parallaxes and extinctions for LAMOST data. Again, 2MASS photometry was used. The reported uncertainty in parallax is about 20% for dwarf stars and 40% for giants. Kinematic correction (see Schönrich et al. 2012) was applied using PPMXL (Roeser et al. 2010) and UCAC4 (Zacharias et al. 2013) data. Data from Carlin et al. (2015) and Wang et al. (2016a) are not yet publicly available.
Distances are provided in the LAMOST Galactic AntiCentre project data release (Yuan et al. 2015). These data include spectroscopic measurements of T_{eff},log g, and [ Fe / H ] and photometry from 2MASS and Xuyi Schmidt Telescope Photometric Survey (Zhang et al. 2014). In their work, Yuan et al. (2015) applied two different methods to get distances. In the first method, which they call “empirical”, stars are divided into four groups (OB stars, giants, and two groups of dwarfs). For each group absolute magnitudes were calculated using a thirdorder polynomial of T_{eff},log g, and [ Fe / H ]. Polynomials were derived by fitting data from the parts of the medium resolution INT Library of Empirical Spectra (MILES) library (SánchezBlázquez et al. 2006) corresponding to each group. The precision of the obtained distance modulus is about for GKM giants and for other groups. A second, “isochrone” distance estimate was derived using the isochrones of Dartmouth Stellar Evolution Database (Dotter et al. 2008). For each star a model with closest values of T_{eff},log g, and [ Fe / H ] was selected from the database. A difference between the visible magnitudes of the star and absolute magnitudes for the closest model provides distance. For both methods extinction values were derived from the LAMOST data using the starpairs method, which is described in Yuan et al. (2014). The “isochrone” method provides distances that are about 5 percent lower than those derived by the “empirical” method.
Studies listed above use similar methods, but the implementation can vary, leading to different results even for the same input data. Moreover, while distances are typically calculated, mass and, most importantly, age estimates are less common. The amount of complementary spectroscopic data available in different surveys calls for a more unified approach. In this paper we present a Unified tool to estimate Distances, Ages, and Masses from spectrophotometric data (hereafter UniDAM). There are two major points in which we differ from studies listed above.
First, whereas most of the previously published studies were dedicated to data from a single survey, we processed data from several large surveys with one tool. For some surveys no data on distances, masses, and ages are publicly available to date. For others our results are consistent with previously published studies with the advantage that our catalogue was produced with the same method, isochrones, and priors on parameters for all surveys. Thus all differences in results for different surveys can be attributed to systematic differences in parameters determined in the spectroscopic surveys. We provide more details on spectroscopic surveys used in Sect. 2. Another advantage of using many surveys simultaneously comes from the fact that different surveys probe different parts of the Galaxy because of different observing strategies and locations of telescopes. Therefore we do not simply increase the statistics, but have a more complete coverage of the Galaxy.
Second, we try to lift the degeneracy of the transformation from observed to stellar parameters by representing PDFs as sums of unimodal functions (unimodal subPDFs or USPDF) for each evolutionary stage. Thus we separate out physically different solutions. This allows us to increase the precision of stellar parameters for each solution.
2. Data samples used
We used observable parameters from a set of publicly available spectroscopic surveys in our work. All surveys were crossmatched with 2MASS (Skrutskie et al. 2006) and AllWISE (Cutri et al. 2014) to get the infrared photometry. We used only “clean” photometry that is only bands that are not affected by low photometric quality, contamination, or confusion. This was achieved by taking only bands with 2MASS quality flag (Qfl) set to ’A’ and AllWISE bands with the contamination and confusion flag (ccf) set to zero and photometric quality flag (qph) set to ’A’. We also requested that the reported uncertainty in magnitude has a positive value. Table 1 summarises properties of the spectroscopic surveys from which we extracted our input data. We discuss some of them below, focusing on parameters for each survey, which we added or modified for our purposes.
Spectroscopic surveys from which data were used.
2.1. APOGEE and APOKASC
We used APOGEE data from SDSS DR12 (Alam et al. 2015) and DR13 (SDSS Collaboration et al. 2017). We kept only those stars that belong to the Main Survey Targets^{1} and have their temperatures, gravities, and metallicities measured. Both DR12 and DR13 were used, as they differ mainly in spectroscopic calibration, and it is interesting to test how that influences the estimates of age, mass, and distance. We include our results for DR13 data in our final catalogue, whereas results for DR12 data are provided as a separate table.
We use as a separate input survey the APOKASC sample (Pinsonneault et al. 2014), although it is in this context just a subset of APOGEE. Therefore the result for this sample is not included in our final catalogue. These data were used to compare the results of Rodrigues et al. (2014) with the prospect of the inclusion of asteroseismic data (see Sect. 5.2.2).
2.2. LAMOST
The second public data release of the LAMOST project (Luo et al. 2015) contains spectral parameters for over 2 million stars. However, the uncertainties in the stellar parameters reported, i.e. 170 K in T_{eff}, 0.5 dex in log g and 0.2 dex in [ Fe / H ], are too high for these data to be used reliably for the model fitting. Therefore we decided not to use the main LAMOST dataset. We focused instead on the LAMOST Galactic AntiCenter (LAMOSTGAC) project second data release (Xiang et al. 2017). This data release contains spectral parameters for about onethird of a million stars in the direction of the Galactic anticenter in its main sample. The bright sample contains over a million stars from a larger area. A different processing pipeline was used by LAMOSTGAC team, which resulted in substantially lower parameter uncertainties of 115 K in T_{eff}, 0.2 dex in log g, and 0.13 dex in [ Fe / H ].
An additional dataset derived from LAMOST DR2 data was prepared with The Cannon tool (Ness et al. 2015; Ho et al. 2017). This tool allows the transfer of parameters from highresolution APOGEE spectra to LAMOST data using stars observed by both surveys for the calibration. This method transfers APOGEE uncertainties in the measured parameters to the LAMOST data, improving the precision of obtained parameters. To account for calibration uncertainties we added in quadrature the median APOGEE absolute uncertainties to the formal uncertainties reported by The Cannon. Another benefit of The Cannon tool is that it measures the value of [ α/ Fe ], which is not provided by LAMOST. The Cannon tool was only calibrated for giant stars, which are available from APOGEELAMOST overlap and therefore the LAMOSTCANNON sample contains only giant stars.
2.3. RAVE surveys
The fifth data release of the RAVE project (Kunder et al. 2017) contains spectral parameters for almost half a million stars. This release contains a flag indicating whether the fitting algorithm has converged, but it turns out that even for stars with this flag set to zero (indicating that the fit converged) there are clear concentrations of values of effective temperatures and gravities towards grid points. This feature is known to the community (see Binney et al. 2014). We added in the output catalogue a flag that indicates if log T_{eff} is within 0.01 dex of a grid point or if log g or [ Fe / H ] are on a grid point. About onethird of the stars are affected by clustering around grid points.
The RAVEon (Casey et al. 2017) is a product of processing of original RAVE spectra with The Cannon tool. Calibration set was constructed from the overlap of RAVE with APOGEE giants and K2/EPIC survey Huber et al. (2016). In addition to T_{eff},log g, and [ Fe / H ], the output of The Cannon tool also contains the value of [ α/ Fe ] and abundances for several chemical elements. The RAVEon data refer to exactly the same stars as the main RAVE survey, but the reported stellar parameters might be slightly different and the quoted uncertainties are smaller, therefore we chose to use RAVEon in our catalogue, providing results for the main RAVE survey in a separate table.
2.4. GenevaCopenhagen survey
GenevaCopenhagen survey (GCS) is the only nonspectroscopic survey used in this work. GCS is a photometric survey, which contains T_{eff},log g, and [ Fe / H ] derived using Stromgren photometry. We used GCS reanalysed data published by Casagrande et al. (2011). We exclude 15% of the stars for which no estimates of T_{eff},log g, and [ Fe / H ] are provided or for which no photometry was present. The latter was mainly because a number of GCS sources are too bright for 2MASS.
2.5. SEGUE
We used SEGUE data from SDSS DR12 (Yanny et al. 2009) with internal uncertainties from the SDSS database. We add in quadrature the internal and systematic uncertainties derived by Allende Prieto et al. (2008) of 130 K in T_{eff}, 0.21 dex in log g and 0.11 dex in [ Fe / H ]. The SEGUE survey is based on SDSS photometry, which is deeper than 2MASS, therefore for about onehalf of SEGUE targets no 2MASS or AllWISE photometry is available or the photometry is very uncertain. We do not use such stars in our work.
2.6. GaiaESO
For the GaiaESO survey, the data release 2 (Gilmore et al. 2012) was used. For nearly half of its nearly 15 000 spectra T_{eff},log g and [ Fe / H ] are available. So we used approximately 7 000 sources from this survey.
2.7. AMBRE
Atmospheric Parameters and Chemical Abundances from Stellar Spectra (AMBRE; Worley et al. 2016) project released parameters extracted from the automatic analysis of the ESO spectral data archives for over 4500 observations (over 2000 sources). No photometry or positional information are provided in the project data, so we attempted to get this information using target names. With the SIMBAD service (Wenger et al. 2000) we obtained positions for nearly 1500 sources, having a total of 3400 observations in the AMBRE survey.
2.8. GALAH
Martell et al. (2017) describe the GALactic Archaeology with HERMES (GALAH) survey first data release. Stellar parameters were derived for 2576 GALAH stars with the Spectroscopy Made Easy (SME) tool (Valenti & Piskunov 2012). These data were used as a training sample for The Cannon tool, which was then used to derive stellar parameters for the rest of the survey. Martell et al. (2017) provide typical uncertainties of the Cannon tool used to derive spectral parameters and internal precision of the SME. We added them in quadrature to get uncertainties of 108 K in T_{eff}, 0.3 dex in log g and 0.11 dex in [ Fe / H ].
2.9. Mock survey
In addition to real survey data we also created a mock survey to test our UniDAM tool. In this case we have full control on both the input parameters for our tool and the desired output parameters of the star. We produced mock surveys by sampling a number of models from PARSEC isochrones Bressan et al. (2012, see Sect. 3).
We stress that the choice of models was aimed at covering model parameter space. So our mock survey does not resemble observed stellar surveys, which are typically magnitudelimited, nor a physical distribution of stars in masses and ages. We motivate our choice by the need to study the behaviour of our tool over a large parameter range. We chose isochrones with 8 different metallicities and 20 ages, which we selected at random. From each isochrone we randomly selected 20 models (with mass below 4 M_{⊙}). We used T_{eff},log g, [ Fe / H ] as well as 2MASS and AllWISE magnitudes for each selected model. Highmass stars were excluded because of their rarity. We took absolute magnitudes from PARSEC models as our “observed” magnitudes, thus setting the distance to 10 pc and extinction to zero. Parameter uncertainties were taken to be 100 K for T_{eff}, 0.1 dex for log g, and [ Fe / H ], and for each magnitude m_{λ}, which is similar to uncertainties in real spectroscopic and photometric surveys.
We prepared four mock surveys. In the first survey we took spectral and photometric values as provided by the PARSEC models. In the second survey we perturbed photometric parameters with random Gaussian noise, while keeping original spectroscopic parameters. In the third we perturbed spectral parameters with random Gaussian noise, while keeping original photometry. In the last survey all parameters were perturbed. Perturbation spread was always taken to be equal to the chosen parameter uncertainties. This allows us to control how uncertainties in observations influence our results.
3. Isochrones
We used PARSEC 1.2S isochrones (Bressan et al. 2012), which provide a large sample of models covering a wide range of stellar parameters. These data include effective temperatures, surface gravities, radii, and absolute photometric magnitudes for a models covering large ranges in metallicities, ages, and masses. We selected nearly three million models that cover the following ranges:

10^{4}:0.06 dex in metallicity (Z), corresponding to − 2.2:0.6 dex in [ Fe / H ];

6.6:10.13 in log(age) τ, corresponding to 4 × 10^{6}:13.5 × 10^{9} yr;

0.09:67 M_{⊙} in mass.
The density of models varies within the ranges indicated above. The reason for this is that isochrones are designed to reproduce details of stellar evolution. Therefore, there are a relatively large number of models covering some rapid stages of evolution and there are a lower number of models for stages of slow evolution. To account for this, we introduced a value w_{j} that is a measure of the volume of the parameter space (metallicity, age, and mass) represented by each model. Otherwise we would be biased towards rare evolutionary stages.
We calculated w_{j} for each model as a product of width of the bin in each dimension represented by the model, (3)The PARSEC isochrone models are calculated for the bin midpoints and they are not equal to the average model in each bin, which makes the binning somewhat arbitrary.
The PARSEC isochrones are equally spaced in log(age)s τ. Therefore the density of models with lower ages is higher than that of models with higher ages. This has to be compensated for to avoid a bias towards lower ages. We took for the age bin width w_{age} the time span represented by the isochrone. It is calculated as w_{age,j} = (10^{τj + 1} − 10^{τj − 1}) / 2, so time span range is defined by midpoints between isochrones in age.
Observations provide [ M / H ] or [ Fe / H ], which are proportional to the logarithm of Z. We created a grid in Z, ranging from 10^{4} to 0.05, such that the spacing between the values of [ Fe / H ] is smaller than the mean uncertainty σ_{[ Fe / H ]} of iron abundance measure at a given [ Fe / H ] in the most precise input data, i.e. APOGEE data. Typically, uncertainties in [ Fe / H ] are smaller for metalrich stars than for metalpoor stars with σ_{[ Fe / H ]} ∝ Z^{0.15}. Therefore Δ_{Z} – the bin width in Z – is roughly Δ_{Z} ∝ Zσ_{[ Fe / H ]} ∝ Z^{1 − 0.15} = Z^{0.85}. The width of the bin in Z was used for w_{Z}, thus ensuring a flat prior in Z. To check the impact of this, we performed tests with proportional to the width of the bin in [ Fe / H ]. These tests showed that this difference has little impact on our results. This is caused by the fact that for given values of T_{eff} and log g, stellar parameters like age, mass, and luminosity are changing slowly with [ Fe / H ] and Z, so the variations in weights are secondorder effects. Thus w_{Z} is the secondorder effect, but we kept it to keep the flat prior in physical quantity Z.
Masses for models were selected by the PARSEC algorithm to track the shape of the isochrone as well as possible. This results in more models in more curved parts of the isochrone. Such an approach produces heavily inhomogeneous coverage of the mass range, which has to be corrected for. We used for w_{mass} the width of the bin in mass.
One of nine evolutionary stages is assigned to every PARSEC model. We grouped these stages into mainsequence stars and giants ascending the red giant branch (precorehelium burning; stage I), corehelium burning stars (stage II), and asymptotic giant branch stars (postcorehelium burning; stage III). These stage labels were used to separate models with different internal structures.
PARSEC model columns.
For each model we used basic physical information (see Table 2) and 2MASS and AllWISE absolute magnitudes that were derived from the luminosities. Other magnitudes are often available, but we did not use them for the reasons discussed below.
4. Methodology
The method used in our tool is similar to the Bayesian method described in Rodrigues et al. (2014). We introduced the vector O for input (“observed”) parameters and their uncertainties O = (T_{eff},log g, [ Fe / H ] ,m_{λ},σ_{Teff},σ_{log g},σ_{[ Fe / H ]},σ_{mλ}). Here, m_{λ} indicates visible magnitudes in several photometric bands and σ_{x} is the uncertainty of the parameter x. These values were taken from surveys listed in Table 1. When α element abundances were available, the metallicity [ Fe / H ] was corrected with the relation [ Fe / H ] = [ Fe / H ] _{0} + log (1. + 0.638 [ α/ M ]) (see Salaris et al. 1993). Additional input parameters for each star can be used in O, for example masses and radii derived from asteroseismic data or parallaxes from Gaia.
We used two vectors for output parameters. The first vector, X_{m} = (τ,M, [ Fe / H ]) represents stellar model parameters log(age), mass, and metallicity. These parameters are taken from isochrone models and therefore have discrete values. We always refer to the actual rather than initial stellar mass because this quantity can be measured from other data, for example from asteroseismic quantities or from binary orbital solutions. The second vector, X_{p} = (μ_{d},A_{K}), where subscript p stands for photometry, represents distance modulus and extinction; these parameters can formally have any value, but we set a physically motivated limit A_{K} ≥ 0 (see discussion in Sect. 4.4). The full output parameter vector is then X = X_{m} ∪ X_{p}.
The probability of having parameters X with given observables O can be expressed via Bayesian formula as (4)We used flat priors on age (in linear scale) and metallicity Z, which means a star formation rate that is constant in time and is independent of Z. The quantitative effect of different priors is described below in Sect. 5.2.3. We used a mass prior based on the IMF F_{IMF} from Kroupa & Weidner (2003). Therefore, (5)Isochrones give us for each X_{m} a new vector O′ = (T_{eff},log g, [ Fe / H ] ,M_{λ}), where M_{λ} indicates absolute magnitudes in several photometric bands. So we can define a function ℐ as O′ = ℐ(X_{m}). Noticeably, [ Fe / H ] is contained in both O and X_{m}, so ℐ_{[ Fe / H ]}(X_{m}) ≡ [ Fe / H ]. We express P(O  X) using two loglikelihoods (6)Here, L_{iso} is a measure of the separation between observed spectral parameters T_{eff},log g, and [ Fe / H ] and those predicted by model parameters X_{m}. Assuming Gaussian uncertainties in O, we can write (7)We use Eq. (7) for the loglikelihood as in most cases spectroscopic surveys do not provide information about the correlations of the uncertainties on the different parameters. When this information or the PDFs of the spectroscopic parameters are known, this information can be included in L_{iso}.
L_{sed} is a measure of similarity of the observed spectral energy distribution (SED) (set of m_{λ} obtained from photometric surveys) and that predicted by isochrone model for X. Visible magnitudes m_{λ} come from 2MASS and AllWISE. These magnitudes are related to the absolute magnitudes M_{λ} in O′ by the following relation: (8)where the extinction in band λ is defined as , with the extinction coefficients R_{λ} taken from Yuan et al. (2013) and summarised in Table 3. Therefore to compare observed visible and model absolute magnitudes we need to know the distance modulus μ_{d} and the extinction A_{K} in the direction of the star. The latter can be obtained from extinction maps or by comparing magnitudes in different bands. We chose the second method and calculated the extinction value for each star using photometric infrared data. This allowed us to take into account variations of extinction that might occur on scales smaller than the typical map resolution. We were also able to calculate extinction for any position on the sky and any distance, whereas a detailed threedimensional extinction maps are created only for the nearest kiloparsec (Gontcharov 2012) or for the Galactic plane (Sale et al. 2014), which is not sufficient for our purpose. A more recent threedimensional map by Green et al. (2015) covers a large fraction of the sky, but a fullsky threedimensional extinction map is still not available.
We use 2MASS and AllWISE data as infrared bands are much less affected by interstellar extinction. Extinction in optical bands is generally higher than in the infrared and can have higher spectral variations between different points on the sky (see a discussion in Majewski et al. 2011). By using infrared data alone we increased the precision of our distance estimates at a cost of decreasing the precision for the extinction estimate. As far as we focus on distances, this seems a fair trade.
Values of the extinction coefficients R_{λ} and C_{λ} used for 2MASS and AllWISE photometry.
For L_{sed} we use the following expression: (9)where the summation is carried out over all bands, for which photometry is available for a given star and V_{corr}(μ_{d}) is a volume correction. We introduce volume correction to compensate for the fact that with a given field of view we probe larger space volume at larger distances than at smaller distances. See a discussion of the effect of volume correction in Sect. 5.3. Using the relation between distance modulus and distance (d = 10^{0.2μd + 1}), we can write (10)We use both L_{iso} and L_{sed} in P(O  X). Therefore on top of the spectroscopic parameters we utilise additional information, namely, the SED of the star. The drawback is that this also brings in the systematic errors of both stellar spectra modelling in PARSEC and possible large errors in photometry in the case of a mismatch between spectroscopic and photometric surveys.
4.1. Probability distribution functions
In order to get the PDF in each parameter, we need to marginalise a multidimensional PDF of output parameters, P(X  O), defined in Eq. (4), over all other parameters. For example, for log(age) τ one has to calculate (11)with the integral taken over the whole parameter space. In practice, we have a discrete sample of models from isochrones. So we can replace P(O  X) with the sum of delta functions (12)where we write for brevity δ_{τ,j} = δ(τ_{j} − τ), δ_{[ Fe / H ] ,j} = δ( [ Fe / H ] _{j} − [ Fe / H ]), δ_{M,j} = δ(M_{j} − M). Here we have to use volumes represented by each model w_{j} from Sect. 3, which reflect the volume of the parameter space represented by the model. The summation is carried out over all models and X_{m,j} = (τ_{j},M_{j}, [ Fe / H ] _{j}) is a vector of parameters of the model j. Therefore we can write, using Eqs. (5) and (12), (13)with again the summation carried out over all models. We need to keep integration over X_{p}, because both μ_{d} and A_{K} are continuous values.
We can make two important simplifications here. First, there is no need to sum over all models because for most of them P(O  X) is very small. We chose a threshold of L_{iso}< 8. In this case T_{eff},log g, and [ Fe / H ] are within 4 sigma uncertainties from observed values. We verified that increasing the threshold from 8 to 12.5 (or going from a combined 4 sigma to 5 sigma uncertainty threshold) leads to marginal changes in the output; parameter estimates change by more than 3% for only less than 2 percent of the stars. Because models are selected from threedimensional space, this comes at a cost of doubling the number of models to be considered. A decreasing the likelihood threshold however leads to more significant changes in the resulting parameters.
Second, L_{sed} is a quadratic form in μ_{d} and A_{K}. Therefore, for a given model j, P(O  X_{p,j}) = exp( − L_{sed}) is a bivariate Gaussian distribution. The location of the maximum of this distribution can be found by solving the system of equations as follows: (14)This is equivalent to the following set of equations: (15)which is solved for μ_{d} and A_{K}. If A_{K}< 0, which is not physical, or if only one magnitude is available we set A_{K} = 0 and obtained μ_{d} from the first part of Eq. (15). By doing so we increase L_{sed} for a given model, which decreases the contribution of this model to the PDFs. In some cases A_{K} is zero for all models for a star. This indicates either that the extinction for this star is statistically indistinguishable from zero or that there is a mismatch between spectral and photometric data, which results in using visible magnitudes from a different star. In the first case we still produce a reliable distance estimate, while in the second case the obtained L_{sed} is high and the quality of the result is low (see Sect. 4.2).
The covariance matrix of P(O  X_{p,j}) is exactly the inverse Hessian matrix H of L_{sed} , i.e. (16)It is important to note that H depends only on C_{λ}, which are constants and photometric uncertainties σ_{mλ}, thus H has a constant value for a given star.
The width of the L_{sed} distribution in μ_{d} for a given model is thus , which is of the order of σ_{m} and is about an order of magnitude smaller than a typical uncertainties in μ_{d} that we derive. This is not true for the extinction, but we are not focused on derivation of highquality extinctions. Moreover, tests show that the error we bring into mean extinction values by this simplification is small. Furthermore there is an obvious correlation between μ_{d} and A_{K}, but we ignore it here. This is justified because Δ_{μd} is approximately one order of magnitude smaller than a typical uncertainties in μ_{d} for a star, so a correlation between derived μ_{d} and A_{K} in a twodimensional PDF P(μ_{d},A_{K}) for a given star is dominated by scatter in μ_{d} and A_{K} for models used to build the PDF, rather than by correlation of μ_{d} and A_{K} for each model.
As far as P(O  X_{p}) is a bivariate Gaussian function, the integral over dX_{p} in Eq. (13) is exactly the value of P(O  X_{p}) at the location of its maximum, derived in Eq. (15). So we can replace the integral in Eq. (13) with a delta function (17)where X_{p,j} is the solution of Eq. (15) for model j. A similar equation can be written for P(M).
For μ_{d} and A_{K} each model contributes to the PDF a Gaussian summand with a width of Δ_{μd} and , respectively. We can correct for using delta function in place of a bivariate Gaussian for P(O  X_{p}) by adding a Gaussian smoothing multiplier with the corresponding width
4.2. Quality of model fit to the data
It is important to quantify how well a set of models represents the observed parameters of a star. To accomplish this, we used the χ^{2} distribution to get the p values from our loglikelihoods. We characterise our model quality by the p value corresponding to the value of for the model with the highest P(O  X_{m,j},X_{p,j}). Here we use . We added V_{corr}(μ_{d}), thus removing the volume correction. This is necessary because V_{corr}(μ_{d}) adds a nonχ^{2} summand, that depends on μ_{d}. Thus, the χ^{2} value used to compute the p value is not the lowest possible value, but that corresponding to the model with the highest P(O  X_{m,j},X_{p,j}), thus the value that is closest to observables. The number of degrees of freedom in this case equals the number of observables, which is the number of available magnitudes plus three (for the temperature, surface gravity, and metallicity dimensions). This p value is designated as p_{best}. Low values of p_{best} can be caused either by observables falling out of the range covered by the models or by inconsistencies between observed stellar parameters and the observed SED. We flagged data with p_{best}< 0.1 (see Sect. 6.1).
In addition to p_{best}, we quantify how good our models are at representing the observed SED. We use the same model used for p_{best} and report as p_{sed} the p value corresponding to the chisquare value . The number of degrees of freedom in this case is equal to the number of available magnitudes. Low values of p_{sed} might be caused, for example, by a mismatch between spectral and photometric data, which results in using visible magnitudes from a different star. Another possible reason is a problematic spectral parameter estimation, which makes T_{eff} inconsistent with the SED from photometry. We flagged data with p_{sed}< 0.1 (see Sect. 6.1).
4.3. Calculating final values
From Eq. (17) we can define a weight for each model j as (19)The PDF in each parameter can thus be calculated as a distribution of parameters for models, with weights W_{j}. For the PDF in μ_{d} and A_{K} we smooth histograms with Gaussian kernel of width Δ_{μd} or Δ_{A}, respectively (see Eq. (18).
4.3.1. Determination of unimodal subPDFs
For each combination of stellar mass, age, and metallicity, PARSEC models provide a single combination of effective temperature and surface gravity. The transition from the effective temperature, surface gravity, and metallicity to stellar mass and age is however nonunique. For a given combination of T_{eff},log g, and [ Fe / H ] with their uncertainties it is possible to find more than one corresponding model with different combinations of age and mass. For example, red clump stars and red giant stars can have similar spectral parameters T_{eff},log g, and [ Fe / H ], but different ages, masses, and internal structures. Therefore, distributions in ages, masses, and absolute magnitudes (and thus distances) are in some cases different from Gaussian. This is illustrated in Fig. 1, where we show typical PDFs in log(age), mass, and distance modulus for two stars. Some of the distributions shown are multimodal with two or more peaks. Reporting mean values and standard deviations do not capture that properly. Mode values, as used by Rodrigues et al. (2014), give the value of the highest peak only. Full distributions can be provided, but they are often considered too complex for further analysis. We suggest an intermediate solution: split PDFs into several USPDFs with each of these described by a unimodal function, assuming that this represents a group of models with similar stellar structure.
We split all models in three evolutionary stages, described in Sect. 3, i.e. precorehelium burning, corehelium burning, and postcorehelium burning stars (plotted in Fig. 1 with red, blue, and yellow, respectively). Splitting our results this way has a benefit in case of overlapping isochrones from different evolutionary stages for a given T_{eff},log g, [ Fe / H ] combination. Without this split, we would combine values for substantially different evolutionary stages that are not physical.
Splitting in evolutionary stages is not enough, as due to curvature of isochrones we can have subgroups of models within one stage that are physically different. A good example is stage I, which contains both mainsequence and giant stars. This results in multimodal distributions of models in space of stellar parameters. An example of such a situation is given in Jørgensen & Lindegren (2005; see their Figs. 1 and 2). To split a multimodal distribution into several unimodal distributions we applied an additional empirically derived routine to our PDFs, which is described in Sect. A. This routine works in the vast majority of cases. Those cases in which our splitting of the PDFs breaks down typically have too few models to produce a histogram (such cases are given the quality flag “N”, see Sect. 6.1).
The overall weight V_{m} for a USPDF m is defined as a ratio of the sum of weights of models within the USPDF and the sum of weights of all models (20)
4.3.2. Output values
We provide output for each USPDF such that it is possible to reproduce the PDF in each parameter we are interested in mass M, logarithm of age τ, distance modulus μ_{d}, distance d, parallax π, and extinction A_{K}. Even for a unimodal distribution the mean, median, or mode might be a poor estimate for the value of interest in case the distribution is nonsymmetric. Values of mode and median might produce less bias but should be used with care as they are not proper moments of the PDF and many statistical methods rely on moments. To provide a simple representation of each USPDF (for each parameter), we fit them with a Gaussian, a skewed Gaussian, a truncated Gaussian, a modified truncated exponential distribution (MTED) and a truncated Student’s tdistribution (see definitions of these functions in the Appendix B). For the truncated functions the upper and lower truncation limits were not fitted, but were set to upper and lower limits of the considered USPDF. We selected the function that gives the lowest symmetric Kullback–Leibler divergence value, which is the measure of the information gain (21)where H_{i} are histogram counts and F_{i} are fitted function values.
Truncated functions were taken because we have a natural upper limit on the age of the star, which is the age of the Universe and therefore there is a lower limit on the mass of a star that left the main sequence, which is approximately 0.7 M_{⊙} if we consider the full range of metallicities. This limit produces sharp cutoffs in histograms, making truncated functions a natural choice.
A modified truncated exponential distribution in logage is equivalent to a flat distribution in ages, thus such a fit indicates that age is poorly constrained. Such PDFs are typical for mainsequence stars, where, as expected, it is hard to constrain age from spectrophotometric data.
In rare (less than 1%) cases in which the fit did not converge for all five fitting functions. This is primarily caused by long tails in the distributions or insufficient data for a proper fit. In this case we reported only mean and standard deviation of the data as fit parameters.
An important property of our result is that values of distance, mass, and log(age) for models are strongly correlated. We used the fact that distance modulus, log(age), and logarithm of mass have a nearly linear correlation within every USPDF in most cases. We report coefficients (slope a and intercept b) of a weighted linear fit and a scatter around it for three relations for each USPDF: τ = a_{1}μ_{d} + b_{1} (distance modulus versus logarithm of age; red lines on left panels of Fig. 2), μ_{d} = a_{2}τ + b_{2} (logarithm of age versus distance modulus; blue lines on left panels of Fig. 2), and M = 10^{a3μd + b3} (distance modulus versus logarithm of mass; red lines on right panels of Fig. 2). An illustration of correlations is given in Fig. 2, where we show the twodimensional PDFs for several stars and our fits to them. From the righthand side panels it is clear that the relation of distance modulus and logarithm of mass is close to linear; the mass is plotted in linear scale, thus our fits are not straight lines. The relation of distance modulus and log(age) is weak for the mainsequence stars and lower giant branch stars. The shape of twodimensional PDF can be quite complex, like in panels a and c of Fig. 2. For these cases the scatter is large and our relations does not work. For giant stars these correlations are much more pronounced, as can be seen in panels e and g of Fig. 2. Correlations between distance modulus, log(age), and log(mass) can be used, for example, if the new distance estimate is obtained from some external source (like Gaia) for a star in our catalogue. We verified that our estimates for mass and log(age) can then be corrected for by the value of the slope times the difference between the externally determined distance modulus and our estimate as follows: In the future work we will show applications of these relations, which are beyond the scope of this work.
Fig. 2 Twodimensional PDFs for typical lower mainsequence (a) and b)), lower giant branch (c) and d)), red clump (e) and f)), and upper giant branch (g) and h)) stars. The left column shows distance modulus and log(age) PDFs and the right column shows distance modulus and mass PDFs. The red line shows τ(μ_{d}) and log M(μ_{d}) fits and the blue line shows – μ_{d}(τ) fit. Shading indicates PDF values (in log scale). Panels e) and f) have two lines of each colour because we detected two USPDFs for this star. 
Summing up, we chose to provide for each USPDF m and each stellar parameter designated as Y_{i,m}, where Y_{i} ∈ (M,τ,μ_{d},d,π,A_{K}) the following quantities:

A weighted mean (catalogue column suffix _mean) of modelvalues Y_{i,j} for all models, (24)where the summation is carried out over all models within USPDF m.

A mode of the USPDF (suffix _mode).

A weighted median value (suffix _median).

A character indicating which fitting function was chosen: «G» for Gaussian, «S» for skewed Gaussian, «T» for truncated Gaussian, «L» for MTED, «P» for truncated Student’s tdistribution, «E» if the fit failed for all five functions, and «N» if there was not enough data for a fit (suffix _fit).

Parameters for a chosen fit (suffix _par). The first two values are location and shape for the chosen best fitting function. For a Gaussian function, by definition, location parameter is equal to the mean value and shape parameter to the variance. If the chosen function is a skewed Gaussian then the third value is the skew value. If the chosen function is a truncated Gaussian or MTED, then third and fourth values are lower and upper limits. If the chosen function is Student’s tdistribution then the third value is the number of degrees of freedom and the fourth and fifth values are lower and upper limits.

One and three sigma confidence intervals. These are defined as a region containing 68.27% (for onesigma uncertainties) or 99.73% (for threesigma uncertainties) of the USPDF, positioned to minimise its span. By construction, such a confidence interval always includes the mode value. For a Gaussian distribution this is equivalent to a range centred on the mean value with width of one or three standard deviations. (suffixes _low_1sigma, _up_1sigma, _low_3sigma, _up_3sigma).
We report all USPDFs with weights V_{m} higher than 0.03. Integer priority values starting from 0 were assigned to each USPDF in order of decreasing weights V_{m}. We list all measures provided for our catalogue in Table 4.
In Fig. 1 we show examples of different USPDFs and fits to them. For the first star (left column) three different evolutionary stages are possible. For the evolutionary stage I a truncated Gaussian is required to fit log(age) distribution and for the evolutionary stage II a skewed Gaussian is needed to fit distribution in mass and distance modulus. For the second star (right column) mainly stage II is possible, but the distribution for this stage can be split in two parts. We need a truncated Student’s tdistribution to fit the histogram for the higher age solution. The small USPDF visible for mass around 1M_{⊙} and log(age) of 9.3 for stage I was excluded because its weight V_{m} is below the accepted 0.03 threshold.
Measures in the output table.
4.4. Role of age and extinction cuts
We chose to impose hard cuts on log(age) τ ≤ 10.13 and extinction A_{K} ≥ 0. This is not an obvious choice, so we justify it here.
We first consider a star for which two equally good solutions are possible: one with τ = 8 and one with τ = 10.7, where all other parameters are equal (see Fig. 3). If we do not use the hard cut on ages, both solutions are reported with good quality flags and equal weights W_{1,2} = 0.5 (blue lines in Fig. 3). The unphysical age value for the second solution might be used as a sign of a problem either with data or with models. It is therefore likely that this solution or even both solutions for this star will be dropped from further study. If we, on the other hand, use the cut in log(age) τ ≤ 10.13, we will keep both solutions, but their weights will change. Only the tail of the USPDF for the second solution will be retained. Because the sum of USPDF weight still has to be unity, the weight of the first solution will increase. We get for this example W_{1} = 0.89 and W_{2} = 0.11 (red lines in Fig. 3). As a result, we get a “realistic” first solution with τ = 8 and retain a part of the second solution. For the second solution, and, in general, for all cases when part of the USPDF is cut away, the mean of the USPDF is a poor measure of log(age) and it is biased towards lower values. But in such cases the USPDF is fitted with either a MTED, a truncated Gaussian or a truncated Student’s tdistribution. So instead of a solution with high weight and correct, but unphysical, mean log(age), we get a solution with lower weight and biased value of the mean with a proper fitting function.
Fig. 3 Sketch of PDF in log(age) for a star with two possible ages without (blue lines) and with (red lines) hard upper limit on log(age). See text for details. 
We now consider the case when only one solution is available for a given star, which is τ = 10.7. We can use such values of τ as an indication of a problem in spectroscopic parameters, photometry, or in isochrones. If the cut in log(age) τ ≤ 10.13 is applied, there are several possibilities. In some cases the PDF in log(age) follows an exponential distribution, which means that the age is poorly constrained for a given star and extending the range of possible ages does not improve this.
An other possibility is that models that have τ ≤ 10.13 represent observables and thus have p_{best} close to unity; see Sect. 4.2. This means that we still have a reliable log(age) PDF for τ ≤ 10.13, but the mean, mode, and median values might be biased. Without the age cut, the mean, mode, and median log(age) values will be above τ = 10.13, and such solution will likely be excluded from further analysis, despite the fact that a fraction of it is reliable.
In yet another case the value for a bestfitting model probability p_{best} will be small, which will indicate potential problems in either the data or with the models. Such cases will be flagged as unreliable with and without age cut.
The same arguments as discussed above are applicable to the cut in extinction. Moreover, a cut in extinction has minor influence on the result, as extinction values are typically very small, i.e. about 10 times smaller than the derived uncertainty in the distance modulus. We verified that negative extinctions typically arise for faint stars, for which photometric uncertainties are large. In the vast majority of cases for which the derived value of extinction is negative, the value is still consistent with zero within uncertainties.
5. Tests: comparison with other measurements
We tested our UniDAM tool in two ways. We first applied it to mock surveys (see Sect. 2.9). By doing that we checked the accuracy and performance of the tool and explored the effect of random perturbations added to input values. Then we proceeded with comparing our parameter measurements for real stars with those obtained by other groups and presented in the literature. The aim of these exercises was to check the quality of our estimates compared to results obtained by the consortia of the different surveys and the sensitivity of our results to priors.
5.1. Mock survey
We ran our UniDAM tool on all four mock surveys, described in Sect. 2.9. Knowing the input values allowed us to evaluate which of the reported measures, that is mean, median, or mode, is the best proxy. In agreement with Jørgensen & Lindegren (2005) we find that the mode is less biased than the mean or median, but produces slightly more outliers.
Mean, mode, and median values show similar qualitative patterns, so we used only mean output values of the highest weight USPDF for comparison. We compared derived mean values X with input values X_{0}. We considered several measures of interest as follows:

Median fractional uncertainty of derived value; this is an internal precision measure.

Median relative deviation (median bias) of derived value ; this shows whether the values that we calculate are systematically offset with respect to the input.

Median absolute relative deviation of derived value ; this shows how scattered our derived values are with respect to input. Median absolute deviation is a much better estimate of scatter than standard deviation in the presence of outliers (Leys et al. 2013).

Outlier fraction rate O; this is the fraction of stars for which the input value X_{0} lies outside the threesigma confidence interval. We use two values: O_{best} is calculated using only highest weight USPDF, whereas O_{all} is calculated using all USPDFs (i.e. X_{0} lies outside the threesigma confidence intervals of all reported USPDFs).
Measures for all four mock surveys are listed in Table 5. For a normally distributed random variable, median absolute deviation σ_{MAD} relates to standard deviation σ as σ ≈ 1.4826 σ_{MAD}, therefore we expected median absolute relative deviation to be approximately equal to twothirds of the median fractional uncertainty. As can be seen from the second and fourth columns of Table 5, this is the case for our data, which means that the distribution of the offsets is close to normal.
Outliers are expected for two reasons. First, a model with the “correct” (i.e. input) values might not belong to the highest weight USPDF. This is revealed by O_{best}. We detected that in about one to seven percent of the cases input parameters are better recovered with USPDF that has second (or even the third) priority. This happens primarily in the upper part of the giant branch, where isochrone overlap is highest. This is inherent to the method; we seek a model that most likely represents the data. If the mock star was taken to be in some short phase of its evolution (thus having low model weight w_{j}), chances are high that we assign a highest weight USPDF not to this phase but to a much longer phase, which has similar observables. Second, because the threesigma range includes by definition 99.7% of the data, we expect at least a fraction of 0.003 of stars for which no USPDF recovers the “correct” values within threesigma confidence intervals (i.e. O_{all} ≳ 0.003) due to our random perturbations added to input values. In fact, this fraction is slightly higher, of the order of 0.01–0.02. We checked that this is caused by a combination of both the perturbations of input values and models selected for the mock sample that is in a very short phase of evolution. In the latter case USPDFs might be pulled away from the “correct” solution by nearby models with higher weights. Another possible case in which there might be no “correct” solution found is when the mock star is located on the edge of the parameter space covered by models.
In cases where perturbations were added, the values of the bias, median absolute relative deviation, and outlier fractions increase. This increase is most prominent in the outlier fractions; this is caused by the fact that sometimes even a small perturbation of input parameters might change the priorities of USPDFs, thus changing the parameters of the highest weight USPDF by a large value.
Results of mocksurvey comparison.
We also tested if distance modulus or parallax provide a better estimate of distance than distance value itself. We find that this is not the case, and all three estimates give very similar precision, with distance itself showing a slightly smaller fraction of outliers. This contradicts the statement of Binney et al. (2014) that “the most reliable distance indicator is the expectation of parallax”. This might be because Binney et al. (2014) compared their derived values with parallaxes from Hipparcos, and comparing parallaxes with parallaxes is likely less biased. We nevertheless provide all three estimates for each star in the output catalogue.
5.2. Literature
We compared our results with values available in the literature. Results of this comparison are shown in Fig. 4 and in Table 6, and are discussed here. The aim of this comparison is to show that our results are consistent with previous studies. In most cases these previous studies were based on similar data and methods. So the differences that appear are primarily due to different models used and differences in details of the method implementation. The exceptions are GCS parallaxes that are coming from Hipparcos and APOKASC distances, derived with asteroseismic values of log g, which are more precise than spectroscopic values. In both cases our results are consistent with published data.
We verified that our extinction estimates are consistent with those provided in surveys. In fact, the differences between our extinction estimates and those in surveys are comparable to differences between extinctions derived with different methods for the same survey, for example, in the LAMOSTGAC data (Xiang et al. 2017). We do not provide a detailed analysis here, as the derivation of precise extinctions is beyond the scope of this work.
Fig. 4 Distributions of relative difference (left) and difference in units of uncertainty (right) for values provided in the literature X_{0} (see legend) and our results X. Circles indicate mean values, star symbols indicate median values. For RAVE survey values from Binney et al. (2014) were used for comparison. 
Differences between values provided in the literature and our results.
5.2.1. GCS parallaxes, masses, and ages
The GCS (Casagrande et al. 2011, and references therein) mainly covers nearby mainsequence stars. The big advantage is that for most of them parallaxes were measured by Hipparcos. The three top rows of Table 6 and Fig. 4 show differences between parallaxes, masses, and log(age)s from GCS and our estimate. We detected a small bias in parallaxes and masses but a negligible bias in log(age)s. Median absolute deviations are three times lower than fractional uncertainties, which means that our method is consistent with GCS results.
5.2.2. Distances of APOKASC red giants
Rodrigues et al. (2014) have determined distances for about 2000 red giant stars from the APOKASC sample. Our distances are less precise than those of Rodrigues et al. (2014) because we do not include asteroseismic data. This test therefore helps us estimate the quality of distance estimations we make for the whole APOGEE sample, as stellar parameters in APOGEE DR13 were calibrated with the use of asteroseismic data from APOKASC. We predicted slightly larger distances (− 0.017 relative offset), but both bias and scatter are well below the mean fractional uncertainty (0.11) of our derived distances. The origin of the bias is likely the difference in how the distance value is calculated when distance PDF is multimodal. This is supported by the fact that for stars with unimodal PDFs we got a relative distance bias of less than 0.0025.
5.2.3. RAVE stars
We ran the UniDAM tool on RAVE DR4 data (Kordopatis et al. 2013) and compared these findings with the results of Zwitter et al. (2010) for distance moduli and Binney et al. (2014) for distance moduli, log(age)s, and masses. We used DR4 data here, as these were used by Zwitter et al. (2010) and Binney et al. (2014). The relative difference between our distance estimates and μ_{d,Z} by Zwitter et al. (2010) is around − 0.01. As compared to the Binney et al. (2014) results (μ_{d,B}), our distance moduli have a relative difference of − 0.03. The median absolute deviations are large in both cases and are comparable to or larger than the mean relative uncertainties of our values.
The reason for a larger difference for μ_{d,B} is that Binney et al. (2014) use strong priors on distances, metallicities, and ages coming from a model of the Galaxy. These priors are decreasing functions of distance from the Galactic centre and from the Galactic plane. Therefore they decrease with distance from the Sun for the majority of directions probed by RAVE. The prior that decreases with distance results in smaller estimates for stellar distances and thus slightly smaller masses and larger ages as compared to our results. Difference in log(age) are further enhanced due to age priors used by Binney et al. (2014).
As can be seen in panels c and e of Fig. 4, distributions of differences between our results for log(age)s and masses and Binney et al. (2014) results are bimodal, with a secondary peak at approximately − 0.75 in relative mass difference and 0.2 in relative log(age) difference. The same can be seen in panel d. In panel f the second peak is out of the plotted range. This second peak contains about 12% of stars and is caused by a difference in the evolutionary stages accepted in Binney et al. (2014) and by our UniDAM tool. A similar pattern but with a much smaller secondary peak can be seen with data from the mock survey (black histogram in panels c and d of Fig. 4).
We show in Fig. 5 the distributions of the median difference between our results and the Binney et al. (2014) results for distance modulus and log(age) on the HertzsprungRussell diagram. We chose RAVE because it contains estimates for both distance modulus and log(age) and because it covers both mainsequence and giant stars. There is clearly a good agreement in both distance modulus and log(age) for the mainsequence stars and large fraction of giant branch, including the red clump. A disagreement for premainsequence stars and large and hot (thus most massive) giants is primarily due to difference in the models and priors used. Similar plots can be produced for other datasets, revealing similar patterns.
Fig. 5 HertzsprungRussell diagrams of RAVE data showing colour differences between our results and RAVE results (Binney et al. 2014) for distance moduli (top panel) and log(age)s (bottom panel). 
Value of for measures without and with the volume correction applied.
Fig. 6 Median values across the catalogue: a) median uncertainty in log(age); b) median uncertainty in mass; c) p_{best}; d) p_{sed}; e) median weight of the highest weight USPDF; and f) median number of USPDFs with weight V> 0.03 per star. 
5.2.4. LAMOSTGAC distances
We compared our results with two distance estimates provided in Luo et al. (2015). Our values are systematically smaller by a fraction of 0.1 as compared to their “empirical” estimates based on the MILES library. We have much better agreement with estimates based on isochrones from Dartmouth Stellar Evolution Database (0.02 fractional difference). Luo et al. (2015) do not provide uncertainties for their distance estimates. Relative uncertainties of our distance estimates for LAMOSTGAC are higher than estimates build on data from other surveys due to the higher uncertainties in spectral parameters, which lead to a fractional uncertainties on our distances of 0.2.
5.3. Effect of the volume correction
We ran tests to see how much the use of the volume correction (see Eq. (19)) affects our results. Volume correction can be seen as a distant prior that ensures constant number density. In general, if the distance prior is a decreasing (increasing) function, the resulting distance is smaller (larger) than in the case of a flat prior. The size of this effect depends on the relative variation of the prior function within the uncertainty range of the parameter. We chose two datasets for the test: GCS and APOKASC giants (Rodrigues et al. 2014). The GCS dataset contains primarily mainsequence stars with distances derived from Hipparcos parallaxes. These parallaxes are in most cases more precise than our distance measurements. Distances of APOKASC giants were derived using asteroseismic data, and therefore should also be more precise than our measurements. We ran our UniDAM tool with and without volume correction. We selected only USPDFs with highest weight for analysis in each case. We then explored how parallaxes, distances, masses, and log(age)s were affected by volume correction. For multimodal cases it is important that the use of the volume correction might change the relative weights of USPDFs, so that priorities might also change. The result of our experiment is that in 7% cases for APOGEE the assigned evolutionary stage changed when we applied the volume correction. This did not happen for GCS as PDFs are unimodal in most cases for mainsequence stars in that survey. By removing volume correction we decreased distance estimations in both datasets by a fraction of 0.032 if the assigned evolutionary stage did not change. This is well below the median relative distance uncertainties that we have (≈ 0.13). The mass estimates are correlated with distance, and decreased by a fraction of 0.03, again, this is well below relative uncertainties in mass that we find (≈ 0.15). The logarithm of age estimates, which are anticorrelated with distance, increased, but only by a fraction of 0.005 (log(age) fractional uncertainties ≈ 0.03). So the conclusion here is that the volume correction has a measurable and well understood effect on measured parameters, but this effect is smaller than our typical parameter uncertainties. This effect is systematic and has to be taken into account when comparing with results obtained with distance priors; see for example Sect. 5.2.3. However, we expect the contribution from the (unknown) systematic uncertainties of spectroscopic measurements to be at least as high as the influence of the volume correction.
We also compare how the volume correction affects the agreement between our measurements and data from the literature, which is described above in Sect. 5.2. The effect of volume correction is summarised in the Table 7. For the GCS sample there is a clear advantage of using the volume correction. Without volume correction our parallax estimates are lower than those in GCS by a fraction 0.091. If we use the volume correction, our parallax estimates increase on average, which improves the agreement (fractional difference of 0.059; see the first row of Table 7). The same applies for log(age) and mass estimates. As for the APOKASC sample, there seems to be an opposite result, as we seem to overestimate the distance compared to Rodrigues et al. (2014). This is likely caused by the fact that distances provided in Rodrigues et al. (2014) are in fact modes of probability density function. If we use modes instead of means for our USPDFs for distances, than we get a relative difference of less than 10^{3} if the volume correction is included and a relative difference around 0.014 if there is no volume correction used.
The effect of the volume correction increases gradually with increasing distances, where our distance uncertainty is larger. This is caused by an increase in the relative variation of the value of the volume correction within the distance uncertainty. The effect of volume correction is approximately proportional to a square of the uncertainty in distance modulus. If the assigned evolutionary stage changes, the estimates of distance, mass, and log(age) can change by a large amount, sometimes by more than 50%.
6. Stellar parameters catalogue
We provide a catalogue of stellar distances, masses, and log(age)s determined with the UniDAM tool described in this manuscript. Our catalogue contains over 3.8 million rows (one row for each USPDF) for over 2.5 million stars. We summarise some properties of this catalogue in Fig. 6. This figure shows medians of different quantities in each bin on the HertzsprungRussell (HR) diagram. Data from all input spectroscopic surveys have been used to produce this figure. Quantification of differences between spectroscopic data from different surveys and effects of incompleteness and selection are beyond the scope of this paper and will be addressed in future work. Here, we are interested in a qualitative description of how the quality of our estimates vary in different parts of the HR diagram. Panels a and b of Fig. 6 show uncertainties in measured log(age)s and masses. Logages are best constrained in the upper part of the main sequence, where parameters change fast as a star leaves the main sequence. On the contrary, masses are much better constrained on the main sequence.
Panels c and d show median values of p_{best} (the probability for a bestfitting model) and p_{sed} (a measure of how good we reproduce SED with our model). Patterns on both panels are similar with worse results close to the edges of the region covered by the PARSEC isochrones and additionally for p_{sed} between the main sequence and giant branch.
Panels e and f show median values for the weights V_{0} of the highest weight USPDF and total number of USPDFs with V_{i}> 0.03. The patterns are nearly inverse: on the main sequence we typically have only one USPDF with weight equals to unity, whereas on the giant branch the number of USPDFs increases and the weights of the highest weight USPDF decreases. It is important that for the giant branch we typically have two or more USPDFs, therefore using just the one with the highest weight is insufficient; the best solution is to use all USPDFs with their relative weight taken into account.
6.1. Quality flags
The output catalogue contains a quality column, which indicates how reliable data contained in each row are. Values have been assigned as follows:

1
– single PDF;

A
– highestweight USPDF has power of 0.9 or more;

B
– 1st and 2nd priority USPDFs together have power of 0.9 or more;

C
– 1st, 2nd, and 3rd priority USPDFs together have power of 0.9 or more;

D
– 1st, 2nd, and 3rd priority USPDFs together have power of less than 0.9;

L
– low power USPDF (between 0.03 and 0.1);

E
– USPDF has p_{sed}< 0.1 (possibly bad photometry);

X
– highest weight USPDF has p_{best}< 0.1 (likely off the model grid);

N
– USPDF has less than 10 models (unreliable result).
Although the quality value provides some information on the quality of the parameter estimation, it is not recommended to select stars based on that value alone (apart from removing unreliable results with values E, N, or X), because the quality value varies heavily over the HR diagram: for mainsequence stars the quality is in most cases 1 or A, whereas for giants quality B, C, or even D are much more common. This is illustrated by the distribution of the number of USPDFs in panel f of Fig. 6. There are 2% cases where a highest weight USPDF has quality E, 4.2% cases with quality X, and less than 0.01% cases with quality N.
7. Discussion and conclusions
We provide a catalogue of distances, log(age)s, and masses for over 2.5 million stars. This number will increase as new data is made available, for example new data releases for surveys already included, or data from new surveys. Gaia data will be of high value and can be used as an independent test of our distances or as a parallax prior. In the latter case it should improve our extinction, mass, and log(age) estimates considerably.
In the current version of our UniDAM tool we use infrared magnitudes, T_{eff},log g, and [ Fe / H ] as inputs to derive distances, log(age)s, and masses of stars. The tool was also successfully used to derive temperatures for a APOKASC sample, with inputs being surface gravities, and masses derived from seismic information and spectroscopic metallicities (Tayar et al. 2017).
An advantage of our approach is that we represent multipeaked PDFs for parameters with a sum of unimodal distributions. Additionally we provide parameters of fits representing each distribution and the correlations between distance modulus, log(age), and mass. Therefore our catalogue contains not only mean values and uncertainties, but detailed information on PDFs. This allows us to apply more sophisticated analysis to the dataset to reveal both global and local structures in the Galaxy.
The next step will be to add proper motion data, thus obtaining all six dimensions of stellar positions and velocities. Combination of positions and velocities with ages, metallicities, and (where available) chemical abundances will open up new possibilities to study Galactic structure. Furthermore, it is important to get a correct estimate of the selection function, as this might affect results not only quantitatively, but also qualitatively, as was shown by Bovy et al. (2012). We intend to produce a selection function for our catalogue and then proceed to study Galactic structure on large and small scales.
See APOGEE target selection description at http://www.sdss.org/dr12/irspec/targets/
Acknowledgments
Authors thank the anonymous referee for a detailed report with many useful suggestions. It 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 publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. 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. This publication makes use of data products from the Widefield Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, and NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology. WISE and NEOWISE are funded by the National Aeronautics and Space Administration. 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 website is at https://www.ravesurvey.org. Funding for SDSSIII has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSSIII website is http://www.sdss3.org/. SDSSIII is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSSIII Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
References
 Ahn, C. P., Alexandroff, R., Allen de Prieto, C., et al. 2014, ApJS, 211, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Alam, S., Albareti, F. D., Allen de Prieto, C., et al. 2015, ApJS, 219, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Allen de Prieto, C., Sivarani, T., Beers, T. C., et al. 2008, AJ, 136, 2070 [NASA ADS] [CrossRef] [Google Scholar]
 Baglin, A., Auvergne, M., Barge, P., et al. 2006, in The CoRoT Mission PreLaunch Status – Stellar Seismology and Planet Finding, eds. M. Fridlund, A. Baglin, J. Lochard, & L. Conroy, ESA SP, 1306, 33 [Google Scholar]
 Binney, J., Burnett, B., Kordopatis, G., et al. 2014, MNRAS, 437, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bovy, J., Rix, H.W., & Hogg, D. W. 2012, ApJ, 751, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Breddels, M. A., Smith, M. C., Helmi, A., et al. 2010, A&A, 511, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Burnett, B., Binney, J., Sharma, S., et al. 2011, A&A, 532, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Carlin, J. L., Liu, C., Newberg, H. J., et al. 2015, AJ, 150, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Casagrande, L., Schönrich, R., Asplund, M., et al. 2011, A&A, 530, A138 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casey, A. R., Hawkins, K., Hogg, D. W., et al. 2017, ApJ, 840, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G. 2003, PASP, 115, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Cutri, R. M., et al. 2014, VizieR Online Data Catalog: II/328 [Google Scholar]
 Demarque, P., Woo, J.H., Kim, Y.C., & Yi, S. K. 2004, ApJS, 155, 667 [NASA ADS] [CrossRef] [Google Scholar]
 Dotter, A., Chaboyer, B., Jevremović, D., et al. 2008, ApJS, 178, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2016, A&A, 595, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gilmore, G., Randich, S., Asplund, M., et al. 2012, The Messenger, 147, 25 [NASA ADS] [Google Scholar]
 Gliese, W., & Jahreiß, H. 1995, VizieR Online Data Catalog: V/70A [Google Scholar]
 Gontcharov, G. A. 2012, Astron. Lett., 38, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Hawkins, K., Masseron, T., Jofré, P., et al. 2016, A&A, 594, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ho, A. Y. Q., Ness, M. K., Hogg, D. W., et al. 2017, ApJ, 836, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Huber, D., Bryson, S. T., Haas, M. R., et al. 2016, ApJS, 224, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Jørgensen, B. R., & Lindegren, L. 2005, A&A, 436, 127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kjeldsen, H., & Bedding, T. R. 1995, A&A, 293, 87 [NASA ADS] [Google Scholar]
 Kordopatis, G., Gilmore, G., Steinmetz, M., et al. 2013, AJ, 146, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., & Weidner, C. 2003, ApJ, 598, 1076 [NASA ADS] [CrossRef] [Google Scholar]
 Kunder, A., Kordopatis, G., Steinmetz, M., et al. 2017, AJ, 153, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Leys, C., Ley, C., Klein, O., Bernard, P., & Licata, L. 2013, J. Exp. Social. Psychol., 49, 764 [Google Scholar]
 Lindegren, L., Lammers, U., Bastian, U., et al. 2016, A&A, 595, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luo, A.L., Zhao, Y.H., Zhao, G., et al. 2015, RA&A, 15, 1095 [Google Scholar]
 Majewski, S. R., Ostheimer, J. C., Kunkel, W. E., & Patterson, R. J. 2000, AJ, 120, 2550 [NASA ADS] [CrossRef] [Google Scholar]
 Majewski, S. R., Zasowski, G., & Nidever, D. L. 2011, ApJ, 739, 25 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Martell, S., Sharma, S., Buder, S., et al. 2017, MNRAS, 465, 3203 [NASA ADS] [CrossRef] [Google Scholar]
 Martig, M., Fouesneau, M., Rix, H.W., et al. 2016, MNRAS, 456, 3655 [NASA ADS] [CrossRef] [Google Scholar]
 Ness, M., Hogg, D. W., Rix, H.W., Ho, A. Y. Q., & Zasowski, G. 2015, ApJ, 808, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49 [NASA ADS] [Google Scholar]
 Pietrinferni, A., Cassisi, S., Salaris, M., Percival, S., & Ferguson, J. W. 2009, ApJ, 697, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Pinsonneault, M. H., Elsworth, Y., Epstein, C., et al. 2014, ApJS, 215, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 9143, 914320 [Google Scholar]
 Rodrigues, T. S., Girardi, L., Miglio, A., et al. 2014, MNRAS, 445, 2758 [NASA ADS] [CrossRef] [Google Scholar]
 Roeser, S., Demleitner, M., & Schilbach, E. 2010, AJ, 139, 2440 [NASA ADS] [CrossRef] [Google Scholar]
 Salaris, M., Chieffi, A., & Straniero, O. 1993, ApJ, 414, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Sale, S. E., Drew, J. E., Barentsen, G., et al. 2014, MNRAS, 443, 2907 [NASA ADS] [CrossRef] [Google Scholar]
 SánchezBlázquez, P., Peletier, R. F., JiménezVicente, J., et al. 2006, MNRAS, 371, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., Binney, J., & Asplund, M. 2012, MNRAS, 420, 1281 [NASA ADS] [CrossRef] [Google Scholar]
 SDSS Collaboration, Albareti, F. D., Allende Prieto, C., et al. 2017, ApJS, submitted [arXiv:1608.02013] [Google Scholar]
 SilvaAguirre, V., & Serenelli, A. M. 2016, Astron. Nachr., 337, 823 [NASA ADS] [CrossRef] [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]
 Stassun, K. G., & Torres, G. 2017, ApJ, submitted [arXiv:1609.05390] [Google Scholar]
 Tayar, J., Somers, G., Pinsonneault, M. H., et al. 2017, ApJ, 840, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Valenti, J. A., & Piskunov, N. 2012, Astrophysics Source Code Library [record ascl: 1202.013] [Google Scholar]
 Wang, J., Shi, J., Zhao, Y., et al. 2016a, MNRAS, 456, 672 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., Wang, W., Wu, Y., et al. 2016b, AJ, 152, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Woolley, R., Epps, E. A., Penston, M. J., & Pocock, S. B. 1970, Royal Observatory Annals, 5 [Google Scholar]
 Worley, C. C., de Laverny, P., RecioBlanco, A., Hill, V., & Bijaoui, A. 2016, A&A, 591, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Xiang, M., Liu, X., Yuan, H., et al. 2017, MNRAS, 467, 1890 [NASA ADS] [Google Scholar]
 Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, H. B., Liu, X. W., & Xiang, M. S. 2013, MNRAS, 430, 2188 [NASA ADS] [CrossRef] [Google Scholar]
 Yuan, H.B., Liu, X.W., Xiang, M.S., et al. 2014, in Setting the scene for Gaia and LAMOST, eds. S. Feltzing, G. Zhao, N. A. Walton, & P. Whitelock, IAU Symp., 298, 240 [Google Scholar]
 Yuan, H. B., Liu, X. W., Huo, Z. Y., et al. 2015, MNRAS, 448, 855 [NASA ADS] [CrossRef] [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhang, H.H., Liu, X.W., Yuan, H.B., et al. 2014, RA&A, 14, 456 [Google Scholar]
 Zwitter, T., Matijevič, G., Breddels, M. A., et al. 2010, A&A, 522, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Splitting the PDF into unimodal subPDFs
Here we describe a method used to split complex PDFs into a set of unimodal subPDFs.
First, we produced a histogram for logarithm of stellar masses of models of the same evolutionary stage (weighted by W_{j}, see Eq. (19)). Then, we detected local minima and maxima of this histogram. Local minima (or maxima) are defined as locations of bins that have lower (higher) value h_{i} than all other bins within the window: i.e. h_{i} = min { h_{j},i − n ≤ j ≤ i + n } for a local minimum and h_{i} = max { h_{j},i − n ≤ j ≤ i + n } for a local maximum. Window size n was taken to be 3 for maxima and 2 for minima. Differences in window sizes are caused by the need to locate minima with high precision and to avoid too many maxima in noisy data. Formally, it is possible to have more than one local minimum between two local maxima; we split only by the lowest of them in this case. We split the sample at positions of local minima that are lower than 0.75 times the value of the smallest of the two enclosing maxima. We thus could have one or several USPDFs, for each evolutionary stage.
We chose the histogram in logarithm of mass to split the multimodal PDFs as the logarithmic scale is close to linear one around 1 M_{⊙}, but gives a smoother histogram for highmass stars. Mass is a better choice to split the PDF because values of log(age)s are quantised by construction of the isochrones, and distances are much less sensitive to evolutionary stage.
Appendix B: Distributions used in the paper
Here we give definitions for some functions used in the paper.
We define as the PDF of the standard normal distribution and Φ(x) as it’s cumulative distribution. This PDF is designated as a Gaussian in the text.
Introducing and , we have for a truncated Gaussian (B.1)if a<x<b and f(x;α,σ,a,b) = 0 otherwise. Here α is a location, σ – scale, a,b  lower and upper limits.
For skewed Gaussian with a shape parameter s(B.2)We use the definition of the truncated Student’s tdistribution(B.3)where ν is the “number of degrees of freedom” (which can be an arbitrary real number here). Again, Φ_{t}(x,ν) is the cumulative distribution function.
A modified truncated exponential distribution with lower and upper limits a and b, respectively is defined as (B.4)Here, C is the normalisation constant, so that .
For a truncated Student’s tdistribution with lower and upper limits a and b, respectively, we define and . Then for the PDF we have (B.5)
All Tables
Values of the extinction coefficients R_{λ} and C_{λ} used for 2MASS and AllWISE photometry.
All Figures
Fig. 1 Probability distribution function for log(age) (top row), mass (middle row), and distance modulus (bottom row) for two stars from our catalogue. Different evolutionary stages are indicated with different colours (see legend). The stages are mainsequence stars and giants ascending the red giant branch (precorehelium burning, stage I), corehelium burning stars (stage II) and asymptotic giant branch stars (postcorehelium burning, stage III). The histograms represent model distributions and solid lines indicate our fits for individual USPDFs. 

In the text 
Fig. 2 Twodimensional PDFs for typical lower mainsequence (a) and b)), lower giant branch (c) and d)), red clump (e) and f)), and upper giant branch (g) and h)) stars. The left column shows distance modulus and log(age) PDFs and the right column shows distance modulus and mass PDFs. The red line shows τ(μ_{d}) and log M(μ_{d}) fits and the blue line shows – μ_{d}(τ) fit. Shading indicates PDF values (in log scale). Panels e) and f) have two lines of each colour because we detected two USPDFs for this star. 

In the text 
Fig. 3 Sketch of PDF in log(age) for a star with two possible ages without (blue lines) and with (red lines) hard upper limit on log(age). See text for details. 

In the text 
Fig. 4 Distributions of relative difference (left) and difference in units of uncertainty (right) for values provided in the literature X_{0} (see legend) and our results X. Circles indicate mean values, star symbols indicate median values. For RAVE survey values from Binney et al. (2014) were used for comparison. 

In the text 
Fig. 5 HertzsprungRussell diagrams of RAVE data showing colour differences between our results and RAVE results (Binney et al. 2014) for distance moduli (top panel) and log(age)s (bottom panel). 

In the text 
Fig. 6 Median values across the catalogue: a) median uncertainty in log(age); b) median uncertainty in mass; c) p_{best}; d) p_{sed}; e) median weight of the highest weight USPDF; and f) median number of USPDFs with weight V> 0.03 per star. 

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.