A Unified tool to estimate Distances, Ages and Masses (UniDAM) from spectrophotometric data

TL;DR

UniDAM is a unified Bayesian tool that estimates stellar distances, ages, and masses from spectrophotometric data, harmonizing results across multiple large surveys for over 2.5 million stars.

Contribution

This paper introduces UniDAM, a homogenized Bayesian method for deriving stellar parameters from diverse spectrophotometric data sets, improving consistency across surveys.

Findings

Estimated parameters for over 2.5 million stars.

Good agreement with literature stellar parameters.

Detailed PDFs for different evolutionary phases.

Abstract

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. So far consortia carrying out the different spectroscopic surveys have used different tools to determine stellar parameters of stars from their derived effective temperatures (Teff), surface gravities (log g) and metallicities ([Fe/H]) possibly 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. We use spectroscopic data from a variety of large surveys combined with infra-red photometry from 2MASS and AllWISE and compare these in a Bayesian manner with PARSEC isochrones to derive probability density functions (PDFs) for stellar masses, ages and distances. We treat PDFs of pre-helium-core burning, helium-core burning and post helium-core burning solutions as well as different peaks in multi-modal PDFs (i.e. each unimodal sub-PDF) of the different evolutionary phases separately. For over 2.5 million stars we report mass, age and distance estimate for each evolutionary phase and unimodal sub-PDF. We report Gaussian, skewed Gaussian, truncated Gaussian, modified truncated exponential distribution or truncated Student's t-distribution functions to represent each sub-PDF, allowing to reconstruct detailed PDFs. Comparisons with stellar parameter estimates from the literature show good agreement within uncertainties. We present UniDAM - the unified tool applicable to spectrophotometric data of different surveys to obtain a homogenised set of stellar parameters.

Figures

Tables

Table 1. Table 1: Spectroscopic surveys from which data were used (see Section 2 for more information on surveys). First column contains an abbreviated designation of the survey with a star symbol (*) indicating surveys that were included in the final catalogue (other surveys are available as separate tables). The second column lists the approximate number of observations that were used in our work. The third column is the resolution of the spectra from which T eff , log ⁡ g subscript 𝑇 eff 𝑔 T_{\rm{eff}},\log g and [ Fe / H ] delimited-[] Fe H [\rm{Fe/H}] were derived. The fourth, fifth, and sixth columns list median uncertainties in T eff subscript 𝑇 eff T_{\rm{eff}} , log ⁡ g 𝑔 \log g and [ Fe / H ] delimited-[] Fe H [\rm{Fe/H}] . The last column shows the reference.

Survey	N sources	Resolution	$ℳ\left(\Delta T\right)$ (K)	$ℳ\left(\Delta \log g\right)$ (dex)	$ℳ\left(\Delta \left[\mathrm{Fe}/H\right]\right)$ (dex)	Reference
APOGEE (DR12)	$\mathrm{88\hspace{0.17em}000}$	$\mathrm{22\hspace{0.17em}500}$	$91.5$	0.11	0.03	1
APOGEE (DR13)*	$\mathrm{89\hspace{0.17em}000}$	$\mathrm{22\hspace{0.17em}500}$	$91.5$	0.11	0.03	2
APOKASC	$\mathrm{2\hspace{0.17em}000}$	$\mathrm{22\hspace{0.17em}500}$	$91.5$	0.11	0.03	3
LAMOST-GAC (Main sample)*	$\mathrm{368\hspace{0.17em}000}$	$\mathrm{1\hspace{0.17em}800}$	$115$	0.19	0.15	4
LAMOST-GAC (Bright sample)*	$\mathrm{1\hspace{0.17em}075\hspace{0.17em}000}$	$\mathrm{1\hspace{0.17em}800}$	$100$	0.15	0.13	5
LAMOST-CANNON*	$\mathrm{450\hspace{0.17em}000}$	$\mathrm{1\hspace{0.17em}800}$	$96.3$	0.13	0.05	6
RAVE (DR5)	$\mathrm{450\hspace{0.17em}000}$	$\mathrm{7\hspace{0.17em}500}$	$92$	0.20	0.10	7
RAVE-on*	$\mathrm{450\hspace{0.17em}000}$	$\mathrm{7\hspace{0.17em}500}$	$85$	0.14	0.07	8
GCS*	$\mathrm{13\hspace{0.17em}800}$	$\mathrm{20\hspace{0.17em}000}$	$80$	0.10	0.10	9
SEGUE*	$\mathrm{277\hspace{0.17em}500}$	$\mathrm{2\hspace{0.17em}000}$	$145$	0.26	0.13	10
Gaia-ESO (DR2)*	$\mathrm{7\hspace{0.17em}000}$	$\mathrm{16\hspace{0.17em}000}$	$50$	0.10	0.07	11
AMBRE*	$\mathrm{3\hspace{0.17em}400}$	$\mathrm{16\hspace{0.17em}000}$	$120$	0.20	0.10	12
GALAH (DR1)*	$\mathrm{28\hspace{0.17em}000}$	$\mathrm{10\hspace{0.17em}700}$	$108$	0.30	0.11	13
Mock	$4\times \mathrm{8\hspace{0.17em}000}$	-	$100$	0.10	0.10	-

Table 2. Table 2: PARSEC model columns (as named in the output of http://stev.oapd.inaf.it/cgi-bin/cmd_2.7 ).

Column	Description	Unit
Z	Metallicity	-
log(age/yr)	Age	log(years)
M_ini	Initial mass	$M_{\odot }$
M_act	Actual mass	$M_{\odot }$
logL/Lo	Luminosity	-
logTe	Effective temperature	-
logG	Gravity	-
…	Set of absolute magnitudes (see text)	mag
int_IMF	Value of the cumulative IMF function $F\left(M_{ini}\right)$	-
stage	Evolutionary stage	-

Table 3. Table 3: Values of the extinction coefficients R λ subscript 𝑅 𝜆 R_{\lambda} and C λ subscript 𝐶 𝜆 C_{\lambda} used for 2MASS and AllWISE photometry. Values were taken from Yuan et al. ( 2013 )

Band	Value of $R_{\lambda }$	Value of $C_{\lambda }$
J	0.720	2.35
H	0.460	1.5
K	0.306	1
W1	0.180	0.59
W2	0.160	0.52

Table 4. Table 4: Measures in the output table. For items denoted by a cross ( † † \dagger ) multiple columns are provided; see section 4.3.2 for details. More than one USPDF (i.e. more than one line) can be reported for each stage of each star.

Column name	Units	Description
id		Unique ID of the star from the input data
stage		Stage number (I, II or III)
uspdf_priority		Priority order of a given USPDF (starting from 0)
uspdf_weight		Weight $V_m$ of a given USPDF
total_uspdfs		Number of USPDF with $V_m>0.03$
p_best		Probability for a best-fitting model (see Section 4.2)
p_sed		$p$-value from ${\chi }^2$ SED fit (see Section 4.2)
quality		Quality flag (see Section 6.1)
distance_modulus$†$	mag	Distance modulus ${\mu }_d$
distance$†$	kpc	Distance $d$
parallax$†$	mas	Parallax $\pi$
extinction$†$	mag	Extinction $A_K$ in 2MASS K-band
mass$†$	$M_{\odot }$	Mass
age$†$	log(yr)	Logarithm of age $\tau$
Distance modulus - logarithm of age relation:
dm_age_slope		Slope of the relation
dm_age_intercept		Intercept of the relation
dm_age_scatter		Scatter of the relation
Logarithm of age - Distance modulus relation:
age_dm_slope		Slope of the relation
age_dm_intercept		Intercept of the relation
age_dm_scatter		Scatter of the relation
Distance modulus - logarithm of mass relation:
dm_mass_slope		Slope of the relation
dm_mass_intercept		Intercept of the relation
dm_mass_scatter		Scatter of the relation

Table 5. Table 5: Results of mock-survey comparison. ℳ ​ ( t ) ℳ 𝑡 \mathcal{M}\!\left(t\right) is the median value of t 𝑡 t ; Δ ​ X = X 0 − X Δ 𝑋 subscript 𝑋 0 𝑋 \Delta X=X_{0}-X is the difference between input mock survey value X 0 subscript 𝑋 0 X_{0} and our estimate X 𝑋 X ; and σ X subscript 𝜎 𝑋 \sigma_{X} is the uncertainty of our estimate. O b ​ e ​ s ​ t subscript 𝑂 𝑏 𝑒 𝑠 𝑡 O_{best} and O a ​ l ​ l subscript 𝑂 𝑎 𝑙 𝑙 O_{all} are outlier fractions calculated using the highest weight USPDF and all USPDFs, respectively. See text for details.

Parameter	$ℳ\left(\frac{{\sigma }_X}{X}\right)$	$ℳ\left(\frac{\Delta X}{X}\right)$	$M.A.D.$	$O_{best}$	$O_{all}$
No perturbation
mass	0.17	-0.02	0.10	0.06	0.01
age	0.03	0.01	0.02	0.06	0.00
distance	0.13	-0.02	0.08	0.06	0.00
Photometry perturbation
mass	0.17	-0.00	0.10	0.04	0.01
age	0.03	0.00	0.02	0.04	0.00
distance	0.13	0.01	0.08	0.03	0.00
Spectral parameters perturbation
mass	0.17	-0.00	0.12	0.06	0.02
age	0.03	0.00	0.02	0.06	0.01
distance	0.13	0.00	0.11	0.05	0.01
Photometry and spectral parameters perturbation
mass	0.17	-0.01	0.12	0.06	0.02
age	0.03	0.01	0.02	0.06	0.01
distance	0.13	0.00	0.11	0.06	0.01

Table 6. Table 6: Differences between values provided in the literature and our results. ℳ ​ ( t ) ℳ 𝑡 \mathcal{M}\!\left(t\right) is the median value of t 𝑡 t . Δ ​ X = X 0 − X Δ 𝑋 subscript 𝑋 0 𝑋 \Delta X=X_{0}-X is the difference between input value X 0 subscript 𝑋 0 X_{0} and our estimate X 𝑋 X , σ X subscript 𝜎 𝑋 \sigma_{X} is the uncertainty of our estimate. M . A . D . formulae-sequence M A D \rm{M.A.D.} is the median absolute deviation of Δ ​ X X Δ 𝑋 𝑋 \frac{\Delta X}{X} , O b ​ e ​ s ​ t subscript 𝑂 𝑏 𝑒 𝑠 𝑡 O_{best} is the outlier fraction calculated using highest-weight USPDF. For RAVE comparison is done for two publications: (1) – Zwitter et al. ( 2010 ) , (2) – Binney et al. ( 2014 ) . See text for details.

Survey	Value	$ℳ\left(\frac{{\sigma }_X}{X}\right)$	$ℳ\left(\frac{\Delta X}{X}\right)$	$M.A.D.$	$O_{best}$
GCS	$\pi$	$0.13$	$-0.065$	$0.03$	0.03
	$M$	$0.064$	$0.024$	$0.024$	0.02
	$\tau$	$0.02$	$-0.000$	$0.005$	0.05
APOKASC	$d$	$0.11$	$-0.017$	$0.05$	0.15
RAVE (1)	${\mu }_{d,\text{Z}}$	$0.044$	$-0.011$	$0.050$	0.20
RAVE (2)	${\mu }_{d,\text{B}}$	$0.044$	$-0.03$	$0.055$	0.25
	$M$	$0.15$	$-0.016$	$0.12$	0.17
	$\tau$	$0.027$	$0.009$	$0.018$	0.13
LAMOST-GAC	$d_{\text{emp}}$	$0.20$	$0.1$	$0.08$	0.03
(main sample)	$d_{\text{iso}}$	$0.20$	$0.02$	$0.063$	0.05

Table 7. Table 7: Value of ℳ ​ ( Δ ​ X X ) ℳ Δ 𝑋 𝑋 \mathcal{M}\!\left(\frac{\Delta X}{X}\right) for measures without and with the volume correction applied.

Survey	Value	$ℳ\left(\frac{\Delta X}{X}\right)$ without volume correction	$ℳ\left(\frac{\Delta X}{X}\right)$ with volume correction
GCS	$\pi$	$-0.091$	$-0.059$
	$M$	$0.027$	$0.015$
	$\tau$	$-0.000$	$0.000$
APOKASC	$d$	$-0.001$	$-0.021$

Equations

$$\Delta \nu$$

$${\nu }_{\max }$$

$$w_j=w_{\mathrm{age},\mathrm{j}}{w}_{Z,j}{w}_{\mathrm{mass},\mathrm{j}}.$$

$$P(\text{X}|\text{O})=\frac{P(\text{X})P(\text{O}|\text{X})}{P(\text{O})}\propto P(\text{X})P(\text{O}|\text{X}).$$

$$P(\text{X})=F_{\text{IMF}}(M).$$

$$P(\text{O}|\text{X})=P(\text{O}|{\text{X}}_m,{\text{X}}_p)=e^{-L_{iso}-L_{sed}}.$$

$$L_{iso}$$

$$=\sum _{i\in (T_{\mathrm{eff}},\log g,[Fe/H])}\frac{({\mathcal{I}}_i(\text{X})-O_i{)}^2}{2{\sigma }_{O,i}^2}.$$

$$m_{\lambda }={\mu }_d+M_{\lambda }+A_{\lambda }={\mu }_d+M_{\lambda }+C_{\lambda }{A}_K,$$

$$L_{sed}=\sum _{\lambda }\frac{(m_{\lambda }-M_{\lambda }-C_{\lambda }{A}_K-{\mu }_d{)}^2}{2{\sigma }_{m_{\lambda }}^2}-V_{corr}({\mu }_d),$$

$$V_{corr}({\mu }_d)=\log d^2=\log 10^{2(0.2{\mu }_d+1)}=(0.4{\mu }_d+2)\log 10.$$

$$P(\tau )=\text{iiiint}P(\text{X})P(\text{O}|\text{X})dMd[Fe/H]\mathrm{d}{\text{X}}_{\mathrm{p}},$$

$$P(\text{O}|\text{X})=\sum _{j}P(\text{O}|{\text{X}}_{m,j},{\text{X}}_p){\delta }_{\tau ,j}{\delta }_{M,j}{\delta }_{[Fe/H],\mathrm{j}}{w}_j,$$

$$P(\tau )=\sum _j\text{iiiint}$$

$$\times w_{j}dMd[Fe/H]\mathrm{d}{\text{X}}_{\mathrm{p}},$$

$$\{\begin{array}{ll}\frac{\partial L_{sed}}{\partial {\mu }_d}& =0\\ \frac{\partial L_{sed}}{\partial A_K}& =0\end{array}.$$

$$\{\begin{array}{ll}{\sum }_{\lambda }\frac{1}{{\sigma }_{m_{\lambda }}^2}{\mu }_d+{\sum }_{\lambda }\frac{C_{\lambda }}{{\sigma }_{m_{\lambda }}^2}{A}_K& ={\sum }_{\lambda }\frac{m_{\lambda }-M_{\lambda }}{{\sigma }_{m_{\lambda }}^2}-0.4\log 10\\ {\sum }_{\lambda }\frac{C_{\lambda }}{{\sigma }_{m_{\lambda }}^2}{\mu }_d+{\sum }_{\lambda }\frac{C_{\lambda }^2}{{\sigma }_{m_{\lambda }}^2}{A}_K& ={\sum }_{\lambda }\frac{C_{\lambda }(m_{\lambda }-M_{\lambda })}{{\sigma }_{m_{\lambda }}^2},\end{array}$$

$$H=\left|\begin{array}{cc}\frac{{\partial }^2{L}_{sed}}{\partial {\mu }_d^2}& \frac{{\partial }^2{L}_{sed}}{\partial {\mu }_d\partial A_K}\\ \frac{{\partial }^2{L}_{sed}}{\partial {\mu }_d\partial A_K}& \frac{{\partial }^2{L}_{sed}}{\partial A_K^2}\end{array}\right|=\left|\begin{array}{cc}{\sum }_{\lambda }\frac{1}{{\sigma }_{m_{\lambda }}^2}& {\sum }_{\lambda }\frac{C_{\lambda }}{{\sigma }_{m_{\lambda }}^2}\\ {\sum }_{\lambda }\frac{C_{\lambda }}{{\sigma }_{m_{\lambda }}^2}& {\sum }_{\lambda }\frac{C_{\lambda }^2}{{\sigma }_{m_{\lambda }}^2}\end{array}\right|.$$

$$P(\tau )$$

$$\times \delta ({\text{X}}_{p,j}-{\text{X}}_p)w_{j}dMd[Fe/H]\mathrm{d}{\text{X}}_{\mathrm{p}}$$

$$=\sum _{j}P(\text{O}|{\text{X}}_{m,j},{\text{X}}_{p,j})F_{IMF}(M_j)\delta ({\tau }_j-\tau )w_j,$$

$$P({\mu }_d)$$

$$P(A_K)$$

$$W_j=P(\text{O}|{\text{X}}_j)F_{IMF}(M_j)w_j=e^{-L_{iso}-L_{sed}}{w}_j{F}_{IMF}(M_j).$$

$$V_m=\frac{{\sum }_{j\in m}{W}_j}{\sum W_j}.$$

$$D_{KL}=\sum _i{H}_i\log \frac{H_i}{F_i},$$

$${\tau }^{\prime }$$

$$M^{\prime }$$

$$Y_{i,m}=\frac{{\sum }_{j\in m}{W}_j{Y}_{i,j}}{{\sum }_{j\in m}{W}_j},$$

$${\sigma }_{Y_{i,m}}=\sqrt{\frac{{\sum }_{j\in m}{W}_j(Y_{i,j}-Y_{i,m}{)}^2}{{\sum }_{j\in m}{W}_j}}.$$

$$f(x,\alpha ,\sigma ,a,b)=\frac{\varphi (\xi )}{\sigma Z},$$

$$f(x,\alpha ,\sigma ,s)=\frac{2}{\sigma }\varphi (\xi )\Phi (s\xi ).$$

$$f_t(x,\nu )=\frac{\Gamma \left(\frac{\nu +1}{2}\right)}{\sqrt{\nu \pi }\Gamma \left(\frac{\nu }{2}\right)}{\left(1+\frac{x^2}{\nu }\right)}^{-\frac{\nu +1}{2}},$$

$$f_{\text{exp}}(x,\alpha ,\sigma ,a,b)=\{\begin{array}{ll}C{e}^{-\frac{|x-\alpha |}{\sigma }},& \text{if}a<x<b\\ 0,& \text{otherwise.}\end{array}$$

Full text

*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_{\rm 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 pre-helium-core burning, helium-core burning, and post helium-core burning solutions as well as different peaks in multimodal PDFs (i.e. each unimodal sub-PDF) 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 sub-PDF. We report Gaussian, skewed, Gaussian, truncated Gaussian, modified truncated exponential distribution or truncated Student’s t-distribution functions to represent each sub-PDF, 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.

.

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ß, 1991; Perryman et al., 1997). In the nearest future Gaia (Gaia Collaboration et al., 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, 2016). In the next data releases Gaia will provide high-precision 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_{\rm{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., 2016), Large Sky Area Multi-Object 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 Gaia-ESO (Gilmore, 2015) provide rich spectroscopic information for millions of stars.

Photometric surveys can be used in two ways. Stromgren or Washington-DDO51 photometry can be used to estimate stellar parameters such as the effective temperature $T_{\rm{eff}}$, surface gravity $\log g$, and metallicity $[\rm{Fe/H}]$ (see Casagrande et al., 2011), or to separate giants from dwarfs for spectroscopic follow-up (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 low-mass 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_{\rm{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 high-precision 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, and Tayar et al. 2017, accepted). 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.

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_{\rm{eff}}$, surface gravity $\log g$, and metallicity $[\rm{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_{\lambda}$ of a star. Then by comparing the absolute magnitudes to visible magnitudes $m_{\lambda}$ 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_{\rm{eff}},\log g$, and $[\rm{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 Section 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 Geneva-Copenhagen 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_{\rm{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., 2016). Breddels et al. (2010) proposed a method for distance estimation for RAVE stars based on a comparison of observed $T_{\rm{eff}},\log g$, metal abundance $[\rm{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_{\rm{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_{\rm{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\infty}(b,l)$ taken from Schlegel et al. (1998). The extinction at a given distance was calculated as $A_{Vprior}(b,l,s)=A_{V\infty}(b,l)\int_{0}^{s}\rho(s)ds/\int_{0}^{\infty}\rho(s)ds$, where $\rho(s)$ is the density of extincting material along the line of sight, taken from the model of the Galaxy (see the Equation 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 Figure 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 Equations 3, 4, and 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_{\rm{eff}}$ and $[\rm{M/H}]$ with asteroseismic data from Kepler. As alluded to before, from asteroseismic values $\Delta\nu$ and $\nu_{\rm{max}}$ and knowing $T_{\rm{eff}}$, it is possible to derive an estimate of stellar radius $R$ and mass $M$, using scaling relations from Kjeldsen & Bedding (1995), i.e.

[TABLE]

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 Wide-field Infrared Survey Explorer (WISE) to be converted to the PDFs of distance modulus $\mu_{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 one-third of stars in their sample have bimodal PDFs. Bimodal PDFs were treated in the same way as single-peaked 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 Anti-Centre project data release (Yuan et al., 2015). These data include spectroscopic measurements of $T_{\rm{eff}},\log g$, and $[\rm{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 third-order polynomial of $T_{\rm{eff}},\log g$, and $[\rm{Fe/H}]$. Polynomials were derived by fitting data from the parts of the medium resolution INT Library of Empirical Spectra (MILES) library (Sánchez-Blázquez et al., 2006) corresponding to each group. The precision of the obtained distance modulus is about $0.^{m}65$ for GKM giants and $0.^{m}3$ 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_{\rm{eff}},\log g$, and $[\rm{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 star-pairs 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 (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 Section 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 sub-PDFs 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 cross-matched 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.

2.1 APOGEE and APOKASC

We used APOGEE data from SDSS DR12 (Alam et al., 2015) and DR13 (SDSS Collaboration et al., 2016). We kept only those stars that belong to the Main Survey Targets111see APOGEE target selection description at http://www.sdss.org/dr12/irspec/targets/ 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 Section 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_{\rm{eff}}$, $0.5\,$dex in $\log g$ and $0.2\,$dex in $[\rm{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 Anti-Center (LAMOST-GAC) project second data release (Xiang et al., 2017). This data release contains spectral parameters for about one-third of a million stars in the direction of the Galactic anti-center in its main sample. The bright sample contains over a million stars from a larger area. A different processing pipeline was used by LAMOST-GAC team, which resulted in substantially lower parameter uncertainties of $115\,$K in $T_{\rm{eff}}$, $0.2\,$dex in $\log g,$ and $0.13\,$dex in $[\rm{Fe/H}]$.

An additional dataset derived from LAMOST DR2 data was prepared with The Cannon tool (Ness et al., 2015; Ho et al., 2016). This tool allows the transfer of parameters from high-resolution 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 $[\alpha/\rm{Fe}]$, which is not provided by LAMOST. The The Cannon tool was only calibrated for giant stars, which are available from APOGEE-LAMOST overlap and therefore the LAMOST-CANNON sample contains only giant stars.

2.3 RAVE surveys

The fifth data release of the RAVE project (Kunder et al., 2016) 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_{\rm{eff}}$ is within $0.01\,$dex of a grid point or if $\log g$ or $[\rm{Fe/H}]$ are on a grid point. About one-third of the stars are affected by clustering around grid points.

The RAVE-on (Casey et al., 2016) 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_{\rm{eff}},\log g$, and $[\rm{Fe/H}]$, the output of The Cannon tool also contains the value of $[\alpha/\rm{Fe}]$ and abundances for several chemical elements. The RAVE-on 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 RAVE-on in our catalogue, providing results for the main RAVE survey in a separate table.

2.4 Geneva-Copenhagen survey

Geneva-Copenhagen survey (GCS) is the only non-spectroscopic survey used in this work. GCS is a photometric survey, which contains $T_{\rm{eff}},\log g$, and $[\rm{Fe/H}]$ derived using Stromgren photometry. We used GCS re-analysed data published by Casagrande et al. (2011). We exclude $15\%$ of the stars for which no estimates of $T_{\rm{eff}},\log g$, and $[\rm{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_{\rm{eff}}$, $0.21\,$dex in $\log g$ and $0.11\,$dex in $[\rm{Fe/H}]$. The SEGUE survey is based on SDSS photometry, which is deeper than 2MASS, therefore for about one-half 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 Gaia-ESO

For the Gaia-ESO survey, the data release 2 (Gilmore, 2015) was used. For nearly half of its nearly $15\,000$ spectra $T_{\rm{eff}},\log g$ and $[\rm{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 $4\,500$ observations (over $2\,000$ 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 $1\,500$ sources, having a total of $3\,400$ observations in the AMBRE survey.

2.8 GALAH

(Martell et al., 2016) 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. (2016) 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_{\rm{eff}}$, $0.3\,$dex in $\log g$ and $0.11\,$dex in $[\rm{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 Section 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 magnitude-limited, 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 $\rm{M}_{\odot}$). We used $T_{\rm{eff}},\log g,[\rm{Fe/H}]$ as well as 2MASS and AllWISE magnitudes for each selected model. High-mass 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_{\rm{eff}}$, $0.1\,$dex for $\log g$, and $[\rm{Fe/H}]$, and $0^{m}.03$ for each magnitude $m_{\lambda}$, 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 $[\rm{Fe/H}]$

• •

$6.6:10.13$ in log(age) $\tau$, corresponding to $4\cdot 10^{6}:13.5\cdot 10^{9}$ years

• •

$0.09:67\,\rm{M}_{\odot}$ 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,

[TABLE]

The PARSEC isochrone models are calculated for the bin mid-points 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 $\tau$. 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_{\rm{age}}$ the time span represented by the isochrone. It is calculated as $w_{\rm{age},j}=(10^{\tau_{j+1}}-10^{\tau_{j-1}})/2$, so time span range is defined by mid-points between isochrones in age.

Observations provide $[\rm{M/H}]$ or $[\rm{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 $[\rm{Fe/H}]$ is smaller than the mean uncertainty $\sigma_{[\rm{Fe/H}]}$ of iron abundance measure at a given $[\rm{Fe/H}]$ in the most precise input data, i.e. APOGEE data. Typically, uncertainties in $[\rm{Fe/H}]$ are smaller for metal-rich stars than for metal-poor stars with $\sigma_{[\rm{Fe/H}]}\propto Z^{-0.15}$. Therefore $\Delta_{Z}$ – the bin width in $Z$ – is roughly $\Delta_{Z}\propto Z\sigma_{[\rm{Fe/H}]}\propto 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 $w^{\prime}_{Z}$ proportional to the width of the bin in $[\rm{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_{\rm{eff}}$ and $\log g$, stellar parameters like age, mass, and luminosity are changing slowly with $[\rm{Fe/H}]$ and $Z$, so the variations in weights are second-order effects. Thus $w_{Z}$ is the second-order 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_{\rm{mass}}$ the width of the bin in mass.

One of nine evolutionary stages is assigned to every PARSEC model. We grouped these stages into main-sequence stars and giants ascending the red giant branch (pre-core-helium burning; stage I), core-helium burning stars (stage II), and asymptotic giant branch stars (post-core-helium burning; stage III). These stage labels were used to separate models with different internal structures.

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 $\textbf{O}=(T_{\rm{eff}},\log g,[\rm{Fe/H}],m_{\lambda},\sigma_{T_{\rm{eff}}},\sigma_{\log g},\sigma_{[\rm{Fe/H}]},\sigma_{m_{\lambda}})$. Here, $m_{\lambda}$ indicates visible magnitudes in several photometric bands and $\sigma_{x}$ is the uncertainty of the parameter $x$. These values were taken from surveys listed in Table 1. When $\alpha$ element abundances were available, the metallicity $[\rm{Fe/H}]$ was corrected with the relation $[\rm{Fe/H}]=[\rm{Fe/H}]_{0}+\log(1.+0.638[\alpha/\rm{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, $\textbf{X}_{m}=(\tau,M,[\rm{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, $\textbf{X}_{p}=(\mu_{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}\geq 0$ (see discussion in section 4.4). The full output parameter vector is then $\textbf{X}=\textbf{X}_{m}\cup\textbf{X}_{p}$.

The probability of having parameters X with given observables O can be expressed via Bayesian formula as

[TABLE]

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 Section 5.2.3. We used a mass prior based on the IMF $F_{\textrm{IMF}}$ from Kroupa & Weidner (2003). Therefore,

[TABLE]

Isochrones give us for each $\textbf{X}_{m}$ a new vector $\textbf{O'}=(T_{\rm{eff}},\log g,[\rm{Fe/H}],M_{\lambda})$, where $M_{\lambda}$ indicates absolute magnitudes in several photometric bands. So we can define a function $\mathcal{I}$ as $\textbf{O'}=\mathcal{I}(\textbf{X}_{m})$. Noticeably, $[\rm{Fe/H}]$ is contained in both O and $\textbf{X}_{m}$, so $\mathcal{I}_{[\rm{Fe/H}]}(\textbf{X}_{m})\equiv[\rm{Fe/H}]$. We express $P(\textbf{O}|\textbf{X})$ using two log-likelihoods

[TABLE]

Here, $L_{iso}$ is a measure of the separation between observed spectral parameters $T_{\rm{eff}},\log g$, and $[\rm{Fe/H}]$ and those predicted by model parameters $\textbf{X}_{m}$. Assuming Gaussian uncertainties in O, we can write

[TABLE]

We use Equation 7 for the log-likelihood 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_{\lambda}$ obtained from photometric surveys) and that predicted by isochrone model for X. Visible magnitudes $m_{\lambda}$ come from 2MASS and AllWISE. These magnitudes are related to the absolute magnitudes $M_{\lambda}$ in $\textbf{O}^{\prime}$ by the following relation:

[TABLE]

where the extinction in band $\lambda$ is defined as $A_{\lambda}=\frac{R_{\lambda}}{R_{K}}A_{K}=C_{\lambda}A_{K}$, with the extinction coefficients $R_{\lambda}$ 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 $\mu_{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 three-dimensional 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 three-dimensional map by Green et al. (2015) covers a large fraction of the sky, but a full-sky three-dimensional 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.

For $L_{sed}$ we use the following expression:

[TABLE]

where the summation is carried out over all bands, for which photometry is available for a given star and $V_{corr}(\mu_{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 section 5.3. Using the relation between distance modulus and distance ($d=10^{0.2\mu_{d}+1}$), we can write

[TABLE]

We use both $L_{iso}$ and $L_{sed}$ in $P(\textbf{O}|\textbf{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 multi-dimensional PDF of output parameters, $P(\textbf{X}|\textbf{O})$, defined in Equation 4, over all other parameters. For example, for log(age) $\tau$ one has to calculate

[TABLE]

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(\textbf{O}|\textbf{X})$ with the sum of delta functions

[TABLE]

where we write for brevity $\delta_{\tau,j}=\delta(\tau_{j}-\tau)$, $\delta_{[\rm{Fe/H}],j}=\delta([\rm{Fe/H}]_{j}-[\rm{Fe/H}])$, $\delta_{M,j}=\delta(M_{j}-M)$. Here we have to use volumes represented by each model $w_{j}$ from section 3, which reflect the volume of the parameter space represented by the model. The summation is carried out over all models and $\textbf{X}_{m,j}=(\tau_{j},M_{j},[\rm{Fe/H}]_{j})$ is a vector of parameters of the model $j$. Therefore we can write, using equations 5 and 12,

[TABLE]

with again the summation carried out over all models. We need to keep integration over $\textbf{X}_{p}$, because both $\mu_{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(\textbf{O}|\textbf{X})$ is very small. We chose a threshold of $L_{iso}<8$. In this case $T_{\rm{eff}},\log g$, and $[\rm{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 three-dimensional 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 $\mu_{d}$ and $A_{K}$. Therefore, for a given model $j$, $P(\textbf{O}|\textbf{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:

[TABLE]

This is equivalent to the following set of equations:

[TABLE]

which is solved for $\mu_{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 $\mu_{d}$ from the first part of LABEL:eq:system. 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 Section 4.2).

The covariance matrix of $P(\textbf{O}|\textbf{X}_{p,j})$ is exactly the inverse Hessian matrix $H$ of $L_{sed}$ , i.e.

[TABLE]

It is important to note that $H$ depends only on $C_{\lambda}$, which are constants and photometric uncertainties $\sigma_{m_{\lambda}}$, thus $H$ has a constant value for a given star.

The width of the $L_{sed}$ distribution in $\mu_{d}$ for a given model is thus $\Delta_{\mu_{d}}=\sqrt{H^{-1}_{0,0}}$, which is of the order of $\sigma_{m}$ and is about an order of magnitude smaller than a typical uncertainties in $\mu_{d}$ that we derive. This is not true for the extinction, but we are not focused on derivation of high-quality 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 $\mu_{d}$ and $A_{K}$, but we ignore it here. This is justified because $\Delta_{\mu_{d}}$ is approximately one order of magnitude smaller than a typical uncertainties in $\mu_{d}$ for a star, so a correlation between derived $\mu_{d}$ and $A_{K}$ in a two-dimensional PDF $P(\mu_{d},A_{K})$ for a given star is dominated by scatter in $\mu_{d}$ and $A_{K}$ for models used to build the PDF, rather than by correlation of $\mu_{d}$ and $A_{K}$ for each model.

As far as $P(\textbf{O}|\textbf{X}_{p})$ is a bivariate Gaussian function, the integral over $d\textbf{X}_{p}$ in Equation 13 is exactly the value of $P(\textbf{O}|\textbf{X}_{p})$ at the location of its maximum, derived in LABEL:eq:system. So we can replace the integral in Equation 13 with a delta function

[TABLE]

where $\textbf{X}_{p,j}$ is the solution of LABEL:eq:system for model $j$. A similar equation can be written for $P(M)$.

For $\mu_{d}$ and $A_{K}$ each model contributes to the PDF a Gaussian summand with a width of $\Delta_{\mu_{d}}$ and $\Delta_{A_{K}}=\sqrt{H^{-1}_{1,1}}$, respectively. We can correct for using delta function in place of a bivariate Gaussian for $P(\textbf{O}|\textbf{X}_{p})$ by adding a Gaussian smoothing multiplier with the corresponding width

[TABLE]

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 $\chi^{2}$ distribution to get the $p$ values from our log-likelihoods. We characterise our model quality by the $p$ value corresponding to the value of $\chi^{2}=L_{iso}+L^{\prime}_{sed}$ for the model with the highest $P(\textbf{O}|\textbf{X}_{m,j},\textbf{X}_{p,j})$. Here we use $L^{\prime}_{sed}=L_{sed}+V_{corr}(\mu_{d})$. We added $V_{corr}(\mu_{d})$, thus removing the volume correction. This is necessary because $V_{corr}(\mu_{d})$ adds a non-$\chi^{2}$ summand, that depends on $\mu_{d}$. Thus, the $\chi^{2}$ value used to compute the $p$ value is not the lowest possible value, but that corresponding to the model with the highest $P(\textbf{O}|\textbf{X}_{m,j},\textbf{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 Section 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 chi-square value $\chi^{2}=-L^{\prime}_{sed}$. 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_{\rm{eff}}$ inconsistent with the SED from photometry. We flagged data with $p_{sed}<0.1$ (see Section 6.1).

4.3 Calculating final values

From Equation 17 we can define a weight for each model $j$ as

[TABLE]

The PDF in each parameter can thus be calculated as a distribution of parameters for models, with weights $W_{j}$. For the PDF in $\mu_{d}$ and $A_{K}$ we smooth histograms with Gaussian kernel of width $\Delta_{\mu_{d}}$ or $\Delta_{A}$, respectively (see Equation 18).

4.3.1 Determination of unimodal sub-PDFs

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 non-unique. For a given combination of $T_{\rm{eff}},\log g$, and $[\rm{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_{\rm{eff}},\log g$, and $[\rm{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 Figure 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 Section 3, i.e. pre-core-helium burning, core-helium burning, and post-core-helium burning stars (plotted in Figure 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_{\rm{eff}},\log g,[\rm{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 sub-groups of models within one stage that are physically different. A good example is stage I, which contains both main-sequence 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 Figures 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 Appendix 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 Section 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

[TABLE]

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 $\tau$, distance modulus $\mu_{d}$, distance $d$, parallax $\pi$, 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 non-symmetric. 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 t-distribution (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

[TABLE]

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\rm{M}_{\odot}$ if we consider the full range of metallicities. This limit produces sharp cut-offs in histograms, making truncated functions a natural choice.

A modified truncated exponential distribution in log-age is equivalent to a flat distribution in ages, thus such a fit indicates that age is poorly constrained. Such PDFs are typical for main-sequence 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: $\tau=a_{1}\mu_{d}+b_{1}$ (distance modulus versus logarithm of age; red lines on left panels of Figure 2), $\mu_{d}=a_{2}\tau+b_{2}$ (logarithm of age versus distance modulus; blue lines on left panels of Figure 2), and $M=10^{a_{3}\mu_{d}+b_{3}}$ (distance modulus versus logarithm of mass; red lines on right panels of Figure 2). An illustration of correlations is given in Figure 2, where we show the two-dimensional PDFs for several stars and our fits to them. From the right-hand 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 main-sequence stars and lower giant branch stars. The shape of two-dimensional PDF can be quite complex, like in panels a and c of Figure 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 Figure 2. Correlations between distance modulus, log(age), and log(mass) can be used, for example, if the new distance estimate $\mu^{\prime}_{d}$ 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:

[TABLE]

In the future work we will show applications of these relations, which are beyond the scope of this work.

Summing up, we chose to provide for each USPDF $m$ and each stellar parameter designated as $Y_{i,m}$, where $Y_{i}\in(M,\tau,\mu_{d},d,\pi,A_{K})$ the following quantities:

• •

A weighted mean (catalogue column suffix _mean) of model values $Y_{i,j}$ for all models,

[TABLE]

where the summation is carried out over all models within USPDF $m.$

• •

A weighted standard deviation (suffix _err),

[TABLE]

• •

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 t-distribution, ¡¡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 t-distribution 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 one-sigma uncertainties) or 99.73% (for three-sigma 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 Figure 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 t-distribution to fit the histogram for the higher age solution. The small USPDF visible for mass around $1\rm{M}_{\odot}$ and log(age) of $9.3$ for stage I was excluded because its weight $V_{m}$ is below the accepted $0.03$ threshold.

4.4 Role of age and extinction cuts

We chose to impose hard cuts on log(age) $\tau\leq 10.13$ and extinction $A_{K}\geq 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 $\tau=8$ and one with $\tau=10.7$, where all other parameters are equal (see Figure 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 Figure 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) $\tau\leq 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 Figure 3). As a result, we get a “realistic” first solution with $\tau=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 t-distribution. 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.

We now consider the case when only one solution is available for a given star, which is $\tau=10.7$. We can use such values of $\tau$ as an indication of a problem in spectroscopic parameters, photometry, or in isochrones. If the cut in log(age) $\tau\leq 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 $\tau\leq 10.13$ represent observables and thus have $p_{best}$ close to unity; see section 4.2. This means that we still have a reliable log(age) PDF for $\tau\leq 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 $\tau=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 best-fitting 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 Section 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 Section 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 $\mathcal{M}\!\left(\frac{\sigma_{X}}{X}\right)$ ; this is an internal precision measure.

• •

Median relative deviation (median bias) of derived value $\mathcal{M}\!\left(\frac{\Delta X}{X}\right)=\mathcal{M}\!\left(\frac{X_{0}-X}{X}\right)$; this shows whether the values that we calculate are systematically offset with respect to the input.

• •

Median absolute relative deviation of derived value $\rm{M.A.D.}=\mathcal{M}\!\left(\left|\frac{\Delta X}{X}-\mathcal{M}\!\left(\frac{\Delta X}{X}\right)\right|\right)$; 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 three-sigma 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 three-sigma 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 $\sigma_{MAD}$ relates to standard deviation $\sigma$ as $\sigma\approx 1.4826\,\sigma_{MAD}$, therefore we expected median absolute relative deviation to be approximately equal to two-thirds 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 three-sigma 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 three-sigma confidence intervals (i.e. $O_{all}\gtrsim 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.

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 Figure 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 LAMOST-GAC 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.

5.2.1 GCS parallaxes, masses, and ages

The GCS (Casagrande et al., 2011, and references therein) mainly covers nearby main-sequence stars. The big advantage is that for most of them parallaxes were measured by HIPPARCOS. The three top rows of Table 6 and Figure 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 $\mu_{d,\textrm{Z}}$ by Zwitter et al. (2010) is around $-0.01$. As compared to the Binney et al. (2014) results ($\mu_{d,\textrm{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 $\mu_{d,\textrm{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 Figure 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 Figure 4).

We show in Figure 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 Hertzsprung-Russell diagram. We chose RAVE because it contains estimates for both distance modulus and log(age) and because it covers both main-sequence and giant stars. There is clearly a good agreement in both distance modulus and log(age) for the main-sequence stars and large fraction of giant branch, including the red clump. A disagreement for pre-main-sequence 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.

5.2.4 LAMOST-GAC 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 LAMOST-GAC 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 Equation 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 main-sequence 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 main-sequence 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 ($\approx 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 ($\approx 0.15$). The logarithm of age estimates, which are anti-correlated with distance, increased, but only by a fraction of $0.005$ (log(age) fractional uncertainties $\approx 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 section 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 Section 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 Figure 6. This figure shows medians of different quantities in each bin on the Hertzsprung-Russell (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 Figure 6 show uncertainties in measured log(age)s and masses. Log-ages 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 best-fitting 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

• highest-weight 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 main-sequence 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 Figure 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_{\rm{eff}},\log g$, and $[\rm{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, accepted).

An advantage of our approach is that we represent multi-peaked 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.

Acknowledgements

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 Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

This publication makes use of data products from the Wide-field 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 Leibniz-Institut fuer Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; The Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science 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.rave-survey.org.

Funding for SDSS-III 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 SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III 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

Appendix A Splitting the PDF into unimodal sub-PDFs

Here we describe a method used to split complex PDFs into a set of unimodal sub-PDFs.

First, we produced a histogram for logarithm of stellar masses of models of the same evolutionary stage (weighted by $W_{j}$, see Equation 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\leq j\leq i+n\}$ for a local minimum and $h_{i}=max\{h_{j},i-n\leq j\leq 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\rm{M}_{\odot}$, but gives a smoother histogram for high-mass 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 $\phi(x)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}x^{2}\right)$ as the PDF of the standard normal distribution and $\Phi(x)$ as it’s cumulative distribution. This PDF is designated as a Gaussian in the text.

Introducing $\xi=\frac{x-\alpha}{\sigma}$ and $Z=\Phi(\frac{b-\alpha}{\sigma})-\Phi(\frac{a-\alpha}{\sigma})$, we have for a truncated Gaussian

[TABLE]

if $a<x<b$ and $f(x;\alpha,\sigma,a,b)=0$ otherwise. Here $\alpha$ is a location, $\sigma$ - scale, $a,b$ - lower and upper limits.

For skewed Gaussian with a shape parameter $s$

[TABLE]

We use the definition of the truncated Student’s t-distribution

[TABLE]

where $\nu$ is the “number of degrees of freedom” (which can be an arbitrary real number here). Again, $\Phi_{t}(x,\nu)$ is the cumulative distribution function.

A modified truncated exponential distribution with lower and upper limits $a$ and $b$, respectively is defined as

[TABLE]

Here, $C$ is the normalisation constant, so that $\int_{a}^{b}f_{\textrm{exp}}(x,\alpha,\sigma,a,b)=1$.

For a truncated Student’s t-distribution with lower and upper limits $a$ and $b$, respectively, we define $\xi=\frac{x-\alpha}{\sigma}$ and $Z_{t}=\Phi_{t}(\frac{b-\alpha}{\sigma},\nu)-\Phi_{t}(\frac{a-\alpha}{\sigma},\nu)$. Then for the PDF we have

[TABLE]