© C. Babusiaux et al. 2020

1. Introduction

Uncovering the Galactic structure within the Galactic plane is a challenging issue due to the mix between stars and dust at different distances, the dust affecting the light of the stars, which becomes fainter and redder.

Several methods have now been developed to draw 3D extinction maps: fully model-based methods (e.g. Drimmel & Spergel 2001), methods that use a stellar distribution model but derive a non-parametric 3D extinction map (Marshall et al. 2006; Chen et al. 2013; Schultheis et al. 2014 using the Besançon model), methods based on the distribution of stars near the main-sequence turn-off (Gontcharov 2017) or, the most common, using individual stellar distance and extinction estimates that are then inverted (Arenou et al. 1992; Vergely et al. 2010; Lallement et al. 2019; Berry et al. 2012; Chen et al. 2014, 2019; Hanson & Bailer-Jones 2014; Rezaei Kh et al. 2017; Anders et al. 2019). Green et al. (2014) sample the full probability density function of distance and reddening for individual stars, derived on main-sequence star’s broad band photometry, to build their 3D extinction map, taking into account the survey selection function. Sale (2012) used a full hierarchical model to handle simultaneously the mean-distance-extinction relationship for a sightline and the individual stellar properties.

To derive stellar density distributions, most methods are parametric (e.g. Drimmel & Spergel 2001; Reylé et al. 2009). Non-parametric stellar density models have been derived up to now when the extinction could be handled independently, for example assuming that most of the extinction occurs in the foreground of the structure under study: this is suitable at high galactic latitudes (e.g. de Jong et al. 2008) and for the bulge structure outside of the Galactic plane (e.g. López-Corredoira et al. 2000; Wegg & Gerhard 2013, both using deconvolution methods). Other methods specifically studying the bar structure have searched for the red clump position using a magnitude that is independent of extinction (e.g. Stanek et al. 1994; Babusiaux & Gilmore 2005; Nishiyama et al. 2005; Cabrera-Lavers et al. 2008; Wegg et al. 2015).

Here we wish to work within the Galactic plane and derive the non-parametric distribution of both the extinction and the stellar density at the same time. This is the first time this has been attempted in the Galactic disc. For this, we use a Bayesian deconvolution method (Richardson 1972; Lucy 1974) using all the stellar information available within a given line of sight. We present here the algorithm we developed, FEDReD (Field Extinction – Distribution Relation Deconvolver). It is designed to work using one reference catalogue on which the completeness model will be based and any other survey that can provide complementary information on the stars observed. In the description and applications presented here, we chose to use near-infrared surveys such as 2MASS (Skrutskie et al. 2006) or the UKIDSS Galactic Plane Survey (Lucas et al. 2008) as reference catalogues as they can probe large distances in high extinction fields, and Gaia data (Gaia Collaboration 2016, 2018b) as complementary information. In Sect. 2 we present the method, in Sect. 3 the H-R diagram (HRD) priors we constructed, and in Sect. 4 the results in both a simulated field and selected test fields.

2. Method

We present here in detail the Bayesian deconvolution method. In practice it is composed in two main steps: in the first one we treat each star individually, extracting its join probability of having a given extinction and distance, and in the second one we combine those probabilities into the full line of sight extinction and distance probability distribution.

2.1. Bayesian deconvolution

We wish to derive the probability distribution P(A₀, D) that gives the probability of a star having both an interstellar extinction A₀ (extinction at 550 nm, which is roughly the centre of the V band, e.g. Bailer-Jones 2011) and a distance D along a given line of sight. The variation of the extinction with distance is given by P(A₀|D) and the stellar density distribution along the line of sight is given by P(D).

We first assume that we are observing all the N stars along a given line of sight, each star observed O_j having several observables (here as a minimum near-IR (NIR) magnitudes and potentially optical magnitudes and parallax). We wish to derive

$$\begin{array}{c}\hfill P(A_0,D)={\displaystyle \sum _{j=1}^{N}P(A_0,D|O_j)P(O_j).}\end{array}$$(1)

The sum is discrete instead of the usual integral as we are observing a finite number of stars.

Using Bayes’ theorem we get

$$\begin{array}{cc}& P(A_0,D|O_j)=\frac{P(O_j|A_0,D)P(A_0,D)}{P(O_j)}\hfill \\ \hfill & =\frac{P(O_j|A_0,D)P(A_0,D)}{{\int }_{(A_0,D)}P(O_j|A_0,D)P(A_0,D)d{A}_{0}dD}.\hfill \end{array}$$(2)

Following the well-known Richardson-Lucy deconvolution algorithm, we can estimate P(A₀, D) by iteratively computing h(A₀, D) (Lucy 1974, with ξ = A, D and x_n = O_j):

$$\begin{array}{c}\hfill h_{k+1}(A_0,D)=\frac{1}{N}{\displaystyle \sum _j\frac{P(O_j|A_0,D)h_k(A_0,D)}{{\int }_{(A_0,D)}P(O_j|A_0,D)h_k(A_0,D)d{A}_{0}dD}.}\end{array}$$(3)

The initial values h₀(A₀, D) are discussed in Sect. 2.4.

However we do not observe all the stars, but we can model (Sect. 2.3) the selection function (S) through a model of the completeness of our near-infrared data: P(S|m_J, m_H, m_K). We can then iteratively compute P(A₀, D|S):

$$\begin{array}{c}\hfill P(A_0,D|S)={\displaystyle \sum _{j}P(A_0,D|O_j,S)P(O_j|S),}\end{array}$$(4)

with O_j still being an observed star with at least NIR magnitude observables and P(O_j|S) being the probability of an observed star being in the selected sample. Similarly to Eq. (2), we have

$$\begin{array}{c}\hfill P(A_0,D|O_j,S)=\frac{P(O_j|A_0,D,S)P(A_0,D|S)}{{\int }_{(A_0,D)}P(O_j|A_0,D,S)P(A_0,D|S)d{A}_{0}dD}\end{array}$$(5)

with

$$\begin{array}{c}\hfill P(O_j|A_0,D,S)=\frac{P(S|O_j,A_0,D)P(O_j|A_0,D)}{P(S|A_0,D)}.\end{array}$$(6)

As the observed star is actually observed, we have P(S|O_j) = 1. We therefore estimate P(A₀, D|S) by iteratively computing h(A₀, D|S):

$$\begin{array}{c}\hfill h_{k+1}(A_0,D|S)=\frac{1}{N}{\displaystyle \sum _j\frac{{\psi }_k(A_0,D)}{{\int }_{A_0,D}{\psi }_k(A_0,D)d{A}_{0}dD}}\end{array}$$(7)

with ψ_k(A₀, D) = P(O_j|A₀, D) h_k(A₀, D|S)/P(S|A₀, D), all the observed stars contributing with the same weight to the selected sample. At the last iteration K (see Sect. 2.5 for the convergence criteria), we have an estimate of P(A₀, D|S), which we will denote in the following as $\widehat{P}(A_0,D|S)=h_K(A_0,D)$. We can then retrieve an estimate of P(A₀, D) with

$$\begin{array}{c}\hfill \widehat{P}(A_0,D)\propto \frac{\widehat{P}(A_0,D|S)}{P(S|A_0,D)}.\end{array}$$(8)

We assume here that all sources in the catalogue are real stars, which is unfortunately not the case, in particular in crowded fields where false detections and false cross-matches between observations in different filters and catalogues can be numerous. Those false observations require the use of a robust method to derive $\widehat{P}(A_0|D),$ while they should not significantly impact $\widehat{P}(D)$.

2.2. Individual probabilities P(O_j|A₀, D)

To derive each star’s probability P(O_j|A₀, D), we compare its observables to the properties of all points of an intrinsic HRD, either isochrone-based or empirical (see Sect. 3.2 for details on the HRD). A given HRD point i with an absolute magnitude M_i at a distance D and with an extinction A₀, has an apparent magnitude m_i of

$$\begin{array}{c}\hfill m_i=M_i+5logD-5+k_m(i,A_0)A_0,\end{array}$$(9)

where k_m is the extinction coefficient in the given m photometric band. We take into account the fact that k_m is actually a function of the intrinsic colour of the star (known through i) and of the extinction itself A₀ through the formalism of Danielski et al. (2018), using the same coefficients as Lallement et al. (2019):

$$\begin{array}{c}\hfill k_m(i,A_0)=a_1+a_2{X}_i+a_3{X}_i^2+a_4{X}_i^3+a_5{A}_0+a_6{A}_0^2+a_7{X}_i{A}_0,\end{array}$$(10)

where X_i is by default the G − K colour of the HRD point i. If the HRD is based on isochrones, X_i can be chosen to be the stellar temperature. We can then derive

$$\begin{array}{cc}& P(O_j|A_0,D)={\displaystyle \sum _{i}P(O_j|A_0,D,i)P(i)}\hfill \\ \hfill & ={\displaystyle \sum _i\prod _{m\in \{J,H,K\}}P(\stackrel{\sim }{m}|m_i)P(\stackrel{\sim }{\varpi }|\varpi )P(i)}\hfill \end{array}$$(11)

To compute $P(\stackrel{\sim }{m}|m_i)$ and $P(\stackrel{\sim }{\varpi }|\varpi )$ we assume Gaussian observational errors on the magnitudes for the NIR surveys, on the flux for Gaia, and on the parallax.

We derive P(O_j|A₀, D) for a thin 2D grid of distances and extinction. The distance being computed using the magnitudes, we do not use a constant step in distance but in distance modulus μ with a step of 0.05 mag, corresponding to the typical photometric error of the input catalogues. We therefore work in dμ = 5/log(10) dD/D. Similarly, we choose a step in extinction A₀ of 0.05 mag. The grid typically extends from 0.1 to 30 kpc in distance and from 0 to 30 mag in extinction, although this can be adapted to the field of view and the survey to optimise the computation time. Illustrations of the results for different stellar types are provided in Fig. A.1. It shows that the information is mostly carried by red clump stars and that the Gaia parallax and/or photometry is needed to differentiate a red clump star from a red dwarf.

2.3. Selection function model P(S|A₀, D)

As previously, we compute the probability of being selected using the H-R diagram prior:

$$\begin{array}{c}\hfill P(S|A_0,D)={\displaystyle \sum _{i}P(S|A_0,D,i)P(i)=\sum _i\prod _{m\in \{J,H,K\}}P(S|m_i)P(i).}\end{array}$$(12)

We adopted the following model for the completeness of the surveys, which is generic enough to reproduce typical completeness curves:

$$\begin{array}{c}\hfill P(S|m)\propto \{\begin{array}{cc}1-exp[-{\left(\frac{m^{\ast }-m}{\beta }\right)}^{\alpha }]\hfill & \mathrm{if}m\le m^{\ast }\hfill \\ 0\hfill & \mathrm{if}m>m^{\ast }\hfill \end{array}\end{array}$$(13)

where P(S|m) is the probability of a star being observed in a given photometric band given its magnitude true m. The three parameters α, β, and m^* define the completeness function, which depends on the survey and on the crowding of the field of view. We used simulations to derive the parameters of this completeness function for the different surveys. We simulated a few typical fields with the Marshall et al. (2006) extinction model, the Fux (1999) stellar density distribution, and by modelling the errors from the observed catalogues. We fitted the parameters α, β, and m^* of Eq. (13) using the observed J, H, and K magnitude distributions P(m|S) and the simulated ones P(m), and solving P(S|m) ∝ P(m|S)/P(m) on their cumulative distribution functions. We found that α = 2 and β = 2 were globally appropriate for the UKIDSS survey and a sharper curve with α = 10, β = 1 for the 2MASS survey, in agreement with the studies of Lucas et al. (2008) for UKIDSS and Skrutskie et al. (2006) for 2MASS. We found that for both surveys, a reasonable approximation of m^* can be obtained by adding two magnitudes to the maximum of the observed magnitude distribution P(m|S). We note that this approximation is only valid when the distribution of stars is relatively smooth. For example, it is no longer valid when a large stellar density is present near the end of the completeness survey, typically in fields dominated by the bar feature. In those fields, parameters must be adjusted either through simulations, or better, through direct image completeness tests (e.g. Surot et al. 2019).

For the UKIDSS data, photometric errors can go very high so we restricted the data used to stars with photometric errors lower than 0.1 mag. This means in practice restricting UKIDSS photometry to be roughly within J <  19, H <  18, K <  17. We take this truncation into account in our selection function model. For this, we first model the photometric errors in the band m by a fit on the observables:

$$\begin{array}{c}\hfill \sigma (m)=a+b{e}^{\mathit{cm}}.\end{array}$$(14)

As the errors in the UKIDSS survey are a direct function of the observed magnitude, we derive from this the magnitude m_σ corresponding to our truncation on σ. We then have the probability of a theoretical star of magnitude m being selected through the cumulative distribution function of the magnitude errors σ_m derived with Eq. (14):

$$\begin{array}{c}\hfill P(S|m)=P(\stackrel{\sim }{m}<m_{\sigma }|m,\sigma (m)).\end{array}$$(15)

For UKIDSS the global selection probability is then the product of Eqs. (13) and (15).

2.4. Initial values

The construction of the initial value of the iteration, h₀(A₀, D|S), needs two initial conditions: P₀(D) and P₀(A₀|D). Then it is simply computed as

$$\begin{array}{c}\hfill h_0(A_0,D|S)\propto P(S|A,D)P_0(A_0|D)P_0(D).\end{array}$$(16)

2.4.1. Initial distance distribution P₀(D)

We only take into account the cone effect on a constant underlying stellar density profile: P₀(D) ∝ D². We also tested using a disc exponential profile start but that does not influence our results at all and we therefore stay with the simple flat start.

2.4.2. Initial extinction as a function of distance P₀(A₀|D)

We also use a flat start here: P₀(A₀|D) = 1. However, a number of degenerate solutions occur due to the confusion between giants and red dwarfs (see Fig. A.1), especially in crowded fields with high extinctions where there are not enough Gaia DR2 stars to constrain the solution at small distances (typically below 1 kpc but depending on the cone aperture chosen). To avoid this, we tested whether adding an extra, simple, local prior could help; for example, this stops the extinction from being too high locally: P₀(A₀|D) = 0 for A₀ >  10D, with D the distance in kiloparsecs. The local map of Lallement et al. (2019) confirms that this is a very safe prior. However, the algorithm is robust enough so that it is not needed in practice, at least in well-behaved fields. For UKIDSS fields, where the red clump information only starts at relatively large distances and with very little Gaia information, adding this simple local prior is sometimes not enough and using a prior based on 2MASS data over a larger field of view is needed. However such a prior should be robust enough to differential extinction in order to be used safely.

2.5. Convergence

Assessing the convergence of such a deconvolution is always tricky. We decided to stop the iterations when the convergence rate slows down. With

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_k={\displaystyle \sum _{A_0,D}{[h_k(A_0,D|S)-h_{k-1}(A_0,D|S)]}^2}\end{array}$$(17)

we derive the convergence criteria

$$\begin{array}{c}\hfill cr=\left|\frac{{\mathrm{\Delta }}_{k+1}-{\mathrm{\Delta }}_k}{{\mathrm{\Delta }}_k}\right|<0.05\end{array}$$(18)

The 0.05 threshold for the convergence criteria is somewhat arbitrary but was tested on both simulations and real data to enable us to reach the final shape of the $\widehat{P}(A_0|D)$ relation without introducing too much noise into the overall $\widehat{P}(A_0,D)$. We limited the number of iterations to 50, which was reached in a few crowded areas only. We checked on those areas that the algorithm is robust enough to recover $\widehat{P}(A_0|D)$, although the resulting matrix $\widehat{P}(A_0,D)$ is quite noisy.

2.6. Deriving A₀(D)

We derive $\widehat{P}(A_0|D)$ from $\widehat{P}(A_0,D)$ obtained at the end of the deconvolution process with

$$\begin{array}{c}\hfill \widehat{P}(A_0|D)\propto \frac{\widehat{P}(A_0,D)}{\widehat{P}(D)}=\frac{\widehat{P}(A_0,D)}{{\int }_{A_0}\widehat{P}(A_0,D)d{A}_0}.\end{array}$$(19)

We now search for the relation A₀(D) corresponding to the maximum of the probability $\widehat{P}(A_0|D)$ with the physical constraint that A₀(D) should increase with D. For this, we randomly generate 10 000 monte-carlo solutions of A₀(D) (called MCS hereafter) according to the following algorithm. We randomly select a distance D_l within our working area following the probability weights $\widehat{P}(D|S)={\int }_{A_0}\widehat{P}(A_0,D|S)d{A}_0$. For each distance D_l we randomly select a corresponding extinction within the possible values allowed by the increasing constraint and the extinctions already assigned to the previous distances. This random selection of A₀(D_l) is done using the probability weights defined by $\widehat{P}(A_0|D_l)$. We initiate the generation by setting A₀(0) = 0 and A₀(D_max) = A_max. We compute the total log-likelihood of a MCS by $log(\mathcal{L})=\sum _{l}log(\widehat{P}(A_0(D_l),D_l))$. We then select the best 1000 MCSs. We re-generate 10 000 random MCSs but this time further constraining the solutions to be within the extinction envelop of the 1000 best MCSs for each distance. If the log likelihood of the newly generated solution is better than the worst of the MCSs, the new solution replaces it. Finally a median cubic spline fit with an increasing constraint (Ng & Maechler 2007) is applied on the final 1000 MCSs. Those solutions are illustrated in Fig. 1 for our default simulation (Sect. 4.1). The 68% confidence interval (equivalent to 1σ for a normal distribution) is derived from the quantiles of the MCS distributions.

Fig. 1. Resulting $\widehat{P}(A_0|D)$ from the first deconvolution. The green line is the real relation of the simulation. The black line is the resulting relation A0(D) (obtainted through a constrained median cubic spline fit on the MCSs). The red dotted lines are the minimum and maximum envelop of the MCSs. For the second deconvolution, P0(A0, D) is initialised with a zero probability outside of this red envelop.

Fig. 2. Resulting $\widehat{P}(A_0|D)$ from the second deconvolution. The black line is the final A(D) estimated relation. Dotted lines corresponds to the 1σ CI. The green line is the real relation of the simulation. The vertical lines correspond to Dmin and Dmax, and show the valid distance range derived from the red clump star saturation magnitude and completeness limit respectively.

As illustrated in Appendix A, red clump stars provide the strongest constraints on the joint distance-extinction distribution thanks to their small intrinsic dispersion in absolute magnitude and colour. To avoid noise induced by degenerate solutions, we apply the previously described procedure to select the best MCSs on the distance interval where red clump stars are expected to provide information. To do so, we used the red clump absolute magnitudes of Ruiz-Dern et al. (2018) and considered when the red clump is expected to saturate to set the minimum distance and to reach the completeness limit to set the maximum distance. For the first generation of MCSs, the log-likelihoods are only computed within the red clump distance range derived assuming no extinction. For the second generation the distance range is updated using the maximum and minimum extinctions of the best MCS. The resulting MCS results are considered valid within the final distance interval [D_min,D_max] still defined by the expected red clump saturation and completeness limit, using this time the extinction provided by the 1σ confidence interval of the derived $\widehat{P}(A_0|D)$.

A second deconvolution process is done using the first one as the initial probability distribution, that is setting P₀(A₀, D) with a Gaussian distribution according to the MCSs quantiles, zero outside the MCSs envelop, and a flat probability outside D_min and D_max just taking into account the maximum and minimum (respectively) extinction of the last valid distance. This new starting distribution removes in particular the degenerate solutions seen at small distances with large extinction due to the confusion between giants and red dwarfs. This second step is needed mainly to derive accurately the distance distribution, since the first pass already effectively derives the A₀(D) relation. A last restriction on the definition of our estimated distance validity range for our results is the addition of a constraint on D_min and D_max to exclude distances too far away or too close within our cone angle to have enough stars: $\widehat{P}(D<D_{\min }|S)>0.01$ and $\widehat{P}(D<D_{\max }|S)<0.99$. This last restriction is needed in particular for anti-centre areas with a low intrinsic stellar density at large distances.

2.7. Deriving $\widehat{\mathit{P}}(\mathit{D})$

We simply compute $\widehat{P}(D)={\int }_{A_0}\widehat{P}(A_0,D)d{A}_0$. However the result is quite noisy, as usual for Richardson-Lucy deconvolution, and recovering its error bar is not obvious. We chose here to estimate the confidence interval using a simple bootstrap method on the second deconvolution. For this, we bootstrap the input stars and the first deconvolution prior. For the latter, we bootstrap the MCSs, we select a random one, and we put as prior a flat distribution within the 2σ interval defined by the MCS centred around this random MCS. We then remove the cone effect to recover $\rho (D)\propto \widehat{P}(D)/D^2$. The result can be seen in Fig. 3 for two mock catalogues of the same field with 2MASS and UKIDSS photometric properties (see Sect. 4.1). We only fit ρ(D) within the valid distance range as defined previously, that is where red clump stars are within the rough completeness regime in all three NIR bands. The original P(D) is recovered in all cases within two sigma but is quite noisy, in particular at small distances.

Fig. 3. 1σ confidence interval of the derived stellar density ρ(D) as obtained by bootstrapping. We show here the results for two mock catalogues of the same area at l = 10°, dark grey: 2MASS-like, light grey: UKIDSS-like. The green line indicates the input simulation stellar density.

3. H-R diagram priors

We implemented two different H-R diagram priors, an empirical one based on Gaia observations and a theoretical one based on isochrones. One advantage of the Gaia empirical HRD is that it does not need initial mass function (IMF), metallicity, or age priors, and it takes into account naturally the presence of binaries. However the theoretical HRD based on isochrones is still useful if we want to add other constraints than parallax and photometry, for example spectroscopic information, or if we want to test the impact of variations of the HRD within the Galaxy. Another motivation for implementing an empirical HRD is the known mismatch between the atmosphere models used in the isochrones and the observed intrinsic colour-colour relations, in particular for cool stars (e.g. Aringer et al. 2016; Ruiz-Dern et al. 2018).

3.1. Gaia empirical HRD

We use by default an empirical HRD based on Gaia DR2. We restrict ourselves to a distance above the plane |Z| < 50 pc as we are looking here only in the Galactic plane. This value is a trade-off between having enough stars to sample the giant branch and staying as close as possible to the Galactic plane, that is within a relatively homogeneous stellar population mixture. We select only low extinction stars to build the empirical HRD to avoid adding uncertainties associated with a local extinction map, but as a consequence this HRD can only be used in regions with relatively high extinction so that our extinction residuals become negligible. This is the case for the Galactic disc fields for which this HRD has been built. We applied the same astrometric and photometry filters as in Appendix B of Gaia Collaboration (2018a) with the exception of the photometric flux error one¹. Instead we used a sharp limit in magnitude, G <  20, G_RP <  19, and G_BP <  18, to have a simple selection function model. To select low extinction stars we also used the 3D extinction map of Capitanio et al. (2017), but to get enough red giants close to the plane we used the rather large limit of E(B − V) < 0.05 mag. We use this sample to build a Hess Diagram on a grid of G_BP − G_RP colour (step 0.01 mag) and M_G (step 0.02 mag), that is P(M_G, G_BP − G_RP|S′), with S′ being the HRD stars selection function. To correct for the selection function we derive

$$\begin{array}{c}\hfill P(M_G,G_{\mathit{BP}}-G_{\mathit{RP}})\propto \frac{P(M_G,G_{\mathit{BP}}-G_{\mathit{RP}}|S^{\prime })}{P(S^{\prime }|M_G,G_{\mathit{BP}}-G_{\mathit{RP}})}.\end{array}$$(20)

As previously, we call i one of those HRD point (M_G, G_BP − G_RP). We compute the probability of an HRD point being in our selection:

$$\begin{array}{c}\hfill P(S^{\prime }|i)={\displaystyle \int \int {\int }_{d,l,b}P(S^{\prime }|i,d,l,b)P(d,l,b)dddldb}\end{array}$$(21)

We assume an homogeneous sky distribution, which is a good enough approximation for our work within the |Z| < 50 pc constraint of our local sample. We therefore have P(d, l, b) = d²cos(b). We move our integral in parallax space instead of distance as this is our observable:

$$\begin{array}{c}\hfill P(S^{\prime }|i)={\displaystyle \int \int {\int }_{\varpi ,l,b}P(S^{\prime }|i,\varpi ,l,b)cos(b)/{\varpi }^{4}d\varpi dldb.}\end{array}$$(22)

We take into account the selection on the parallax relative uncertainty of 10%, the |Z| < 50 pc constraint and the distance borders of the E(B − V) < 0.05 contours. The parallax uncertainties are assumed to depend on the magnitude and we do not take into account the second order dependency on the colour nor on sky position. Consequently,

$$\begin{array}{cc}& P(S^{\prime }|i,\varpi ,l,b)=P(\stackrel{\sim }{\varpi }/{\sigma }_{\varpi }>10|M_G,\varpi )\hfill \\ \hfill & P(E(B-V)<0.05|\varpi ,l,b)\hfill \\ \hfill & P(|Z|<50|\varpi ,l,b)\hfill \\ \hfill & P(G<20,G_{\mathit{RP}}<19,G_{\mathit{BP}}<18|\varpi ,i).\hfill \end{array}$$(23)

We model σ_ϖ as a function of G = M_G − 5 − 5 log(ϖ) by fitting a random sample representative of the Gaia data: σ_ϖ = 0.023 + exp(0.828G − 16.9) for G >  13 and σ_ϖ = 0.04 for G <  13.

$$\begin{array}{c}\hfill P(\stackrel{\sim }{\varpi }/{\sigma }_{\varpi }>10|M_G,\varpi )=1-P(\stackrel{\sim }{\varpi }<10{\sigma }_{\varpi }|\varpi ,{\sigma }_{\varpi }(G)),\end{array}$$(24)

the last term being determined by the cumulative distribution function of the Gaussian centred on ϖ with dispersion σ_ϖ.

The E(B − V) < 0.05 constraint corresponds to d <  d_max(l, b), so

$$\begin{array}{c}\hfill P(E(B-V)<0.05|\varpi ,l,b)=1-P(\stackrel{\sim }{\varpi }<1/d_{\max }|\varpi ,d_{\max }(l,b)).\end{array}$$(25)

The |Z| < 50 pc constraint corresponds to ϖ >  |sinb|/0.05. We therefore have

$$\begin{array}{c}\hfill P(|Z|<50|\varpi ,l,b)=1-P(\stackrel{\sim }{\varpi }<\frac{|sinb|}{0.05}|\varpi ,b)\end{array}$$(26)

The maximum distance probed during the integration is set by the extinction constraint, which corresponds to 1 kpc.

The photometric bands used in this study, M_X = M_BP, M_RP, M_J, M_H, M_K are added to the (M_G,G_BP − G_RP) HRD using colour-colour relations M_G − M_X as a function of G_BP − G_RP. We chose a seventh order polynomial model for those relations and ensured that we did not extrapolate them. The calibrations are done using stars with extinction lower than E(B − V) < 0.01 mag, photometric errors less than 2% in the G band, and 5% in G_BP and G_RP, applying the photometric excess flux filter (Evans et al. 2018), G >  6 to avoid saturation, G_BP <  18 to avoid background subtraction issues (Evans et al. 2018), 2MASS photometric quality “AAA”, and photometric errors in J, H, and K smaller than 0.05 mag. As there are not enough very red stars with our strict criteria, we increased the extinction criterion to E(B − V) < 0.015 mag for G_BP − G_RP >  4. We apply three different calibrations for (i) the red giants with M_G <  4.5 mag and M_G <  −5.5 + 9(G_BP − G_RP) (ii) the white dwarfs with M_G <  10 + 2.6(G_BP − G_RP), and (iii) the dwarfs in between. To ensure the continuity between those calibrations, the red giant calibration has been derived with all stars with M_G <  2.5, and the white dwarf relations have been derived using both the dwarf and the white dwarfs sample. As the different colour-colour relations are correlated, we fitted them simultaneously² within each population.

For a better modelling of the bottom of the HRD, quite important for the pollution of nearby red dwarfs in our colour-magnitude diagrams (CMD), we applied the same procedure on a Gaia-2MASS HRD, P(M_G, J − K), without constraint on G_BP nor G_RP, using J − K as the reference colour and selecting stars with J <  15.8 and K <  14.3 (which corresponds to the > 99% completeness of the 2MASS catalogue (Skrutskie et al. 2006)). For the colour-colour relations of the faintest red dwarfs, we had to use the Phoenix relations (Baraffe et al. 2015) to derive the G_BP photometry for J − K >  1.5. We merged the results of both HRDs at M_G = 6.5. Figure 4 shows the difference between both HRD densities for the bottom of the main sequence: at M_G <  11 the higher resolution of Gaia leads to more intrinsically faint stars being observed than 2MASS, while for fainter absolute magnitudes the G_BP quality criteria leads to too significant incompleteness in the Gaia data to be properly modelled. To confirm our interpretation of the differences we show on the same plot the HRD densities obtained with the theoretical HRD described in Sect. 3.2. However the exact density of low mass dwarfs has no implication for this work, which concentrates on more distant stars; they only need to be present with the correct colour-colour relations to avoid them ending as strong outliers.

Fig. 4. Relative stellar densities along the HRD P(MG) as derived from Gaia data only (P(MG, GBP − GRP), black dots), Gaia and 2MASS (P(MG, J − K), red stars), and Padova isochrones (blue triangles). The different densities are normalized at MG = 6.5.

Fig. 5. Gaia DR2 empirical H-R diagram. Here we overlay some stars at different evolutionary stages for which we study their individual P(O|A0, D) in Appendix A.

3.2. Theoretical H-R diagram prior

To build a theoretical HRD, we use the PARSEC isochrones (Bressan et al. 2012) with a step of 0.1 Gyr in age between [0.1, 13.4] and a step of 0.05 dex in [M/H] between [−2.15, 0.5]. Each isochrone point i, corresponding to a metallicity [M/H]_i, age τ_i, and mass ℳ_i, has a weight associated to it P(i) according to the IMF P(ℳ), an age distribution P(τ), and an age-metallicity relation (AMR) P([M/H]|τ). In other words

$$\begin{array}{c}\hfill P(i)=P([\mathrm{M}/\mathrm{H}],\tau ,\mathcal{M})=P(\mathcal{M})P(\tau )P([\mathrm{M}/\mathrm{H}]|\tau )\end{array}$$(27)

We use the Chabrier (2001) log-normal IMF (integrated over the mass interval between isochrone points), the Rocha-Pinto et al. (2000) AMR, and a constant star formation rate (SFR). Here we are observing in the Galactic plane, so we can easily compute a rough correction to the relative number of stars of each age due to the different scale height H_z of the populations as a function of age. Indeed if we assume that the density distribution at each age τ can be modelled under the assumption of an isothermal disc, the solution of the Jeans equation is then a sech² profile:

$$\begin{array}{c}\hfill \rho (z)={\rho }_0{\mathrm{sech}}^2\left(\frac{z}{2H_z}\right).\end{array}$$(28)

If we integrate over z and assume that all the populations have the same radial density profile, we have

$$\begin{array}{c}\hfill \mathrm{\Sigma }={\displaystyle {\int }_z}\rho (z)dz=4{\rho }_0{H}_z\end{array}$$(29)

and therefore a constant SFR corresponds to a local density

$$\begin{array}{c}\hfill P(\tau )\propto {\rho }_0\propto 1/H_z.\end{array}$$(30)

We are using here the default Hz of Trilegal 1.7 (Girardi et al. 2005): H_z = 0.095(1 + τ/5.55)^1.6666. The resulting prior HRD is shown in Fig. 6.

Fig. 6. Alternative H-R diagram to Fig. 5 based on Padova isochrones.

4. Tests and results

We made extensive tests on our algorithm, first using simulations, then on real 2MASS and UKIDSS data combined with Gaia DR2. To check our results on real fields, we looked at a few fields where we knew what to expect, like the ones described below, and at several ones presenting different issues (e.g. high crowding, low stellar density towards the anti-centre, convergence issues). For those, we checked how well our derived A₀(D) function enables us to recover the red clump track. We also checked that both 2MASS and UKIDSS provided consistent results within the uncertainties.

For 2MASS we selected stars with good photometric quality flags (A, B, C, or D). Following the prescription of Lucas et al. (2008), we corrected the errors provided with the UKIDSS catalogue by

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{cor}}=\sqrt{{(1.2\sigma )}^2+0.{02}^2},\end{array}$$(31)

and we only selected stars with a photometric error lower than 0.1 mag. For the cross-match with Gaia DR2, we used the cross-match with 2MASS provided within Gaia DR2 (Marrese et al. 2019) and a simple cross-match within a radius of 0.15″ for UKIDSS. We applied the same Gaia photometric and astrometric filters as detailed in Hottier et al. (2020): phot_bp_rp_excess_factor > 1.3 + 0.06 (bp_rp)², G_BP >  18, astrometric_chi2_al/(astrometric_n_good_obs_al-5) < 1.44 max(1, exp(−0.4*(G-19.5))), ϖ + 3σ_ϖ <  0. We took into account the 3 mmag mag⁻¹ drift of the G band, we added quadratically 10 mmag to the photometric uncertainties to take into account the systematics, and we corrected the parallax from the −0.03 mas zero point.

4.1. Simulation

We tested our procedure on a simulation, as illustrated in Figs. 1 and 2, corresponding to either 2MASS or UKIDSS observations towards l = 10°. We simulated stars with intrinsic stellar properties randomly taken from the Hess diagram described in Sect. 3.1 and placed them along the line of sight following the Fux (1999) model stellar distance distribution and the cone effect. The extinction density is assumed to be proportional to the Fux (1999) model gas density in this direction and the proportion factor is simply derived assuming an integrated extinction along the line of sight of $A_0^{\infty }$ = 32 mag. An intrinsic dispersion in the extinction is added using a log-normal distribution with σ_A = 0.05. The photometric and parallax errors of 2MASS, UKIDSS, and Gaia are added assuming a simple increase in the errors with magnitude following a fit of Eq. (14) on catalogue data. The completeness is then simulated following Eq. (13) with α = 10, β = 1, $m_K^{*}=14.1$, $m_H^{*}=14.8$, $m_J^{*}=16.6$ for 2MASS and α = β = 2 and $m_K^{*}=19$, $m_H^{*}=19.5$, $m_J^{*}=22$ for UKIDSS. GaiaG photometry and parallaxes are kept only if they satisfy the same completeness model like Eq. (13) with $m_G^{*}=20.7$ and G_BP and G_RP photometry with $m_{G_{\mathit{BP}}}^{\ast }=20.9$ and $m_{G_{\mathit{RP}}}^{\ast }=19.5$ (the exact values are not important as they do not enter the catalogue completeness model, but just allow us to take into account that Gaia information is not present for the faintest nor reddest stars). In this way, we have built mock catalogues of about 4000 stars satisfying our photometric criteria. Ten percent of the UKIDSS mock catalogue has Gaia parallax information compared to 50% for the 2MASS one. The UKIDSS mock catalogue is represented in Fig. 7 using the magnitude independent of extinction K_JK (e.g. Babusiaux & Gilmore 2005):

$$\begin{array}{c}\hfill K_{\mathit{JK}}=K-\frac{k_K}{k_J-k_K}(J-K),\end{array}$$(32)

Fig. 7. Colour-magnitude diagram of the UKIDSS mock catalogue built from the Fux model stellar and gas particle distributions towards l = 10°. The magnitude independent of extinction KJK is used. The corresponding distance and extinction for a red clump star on this diagram are indicated on the right and top axis.

with k_J and k_K the extinction coefficients in the J and K bands respectively.

Figure 8 shows the final results of the deconvolution on the UKIDSS mock catalogue. We see that the bootstrap confidence interval is smaller than the one derived directly from the full P(A₀|D) in Fig. 2, which also takes into account the intrinsic dispersion of the extinction as well as the deconvolution artefacts and is therefore the confidence interval to be used. The residuals are within the 1σ confidence interval (with a dispersion of 0.55 mag) but are correlated by the fact that we impose a continuously increasing fit and the deconvolution artefacts. For the density the result is within the 2σ bootstrap confidence interval. Here again the residuals are correlated and correspond to an error of about 20% on the density estimation.

Fig. 8. Results of the deconvolution of the simulation for A0(D) (top) and ρ(D) (bottom) within the [Dmin,Dmax] distance range. The black line is the deconvolution result. The light grey area in the top panel corresponds to the 1σ confidence interval of A0(D) derived from the full P(A0|D) (Fig. 2) while the darker grey area in both panels shows the 1σ confidence interval derived from the bootstrap. The green line is the input relation used in the simulation. Red dashed line: isochrones HRD. Blue dotted line: Cardelli et al. (1989) extinction law. Orange dot-dashed line: assuming completeness parameters α = 10 and β = 1.

We tested the influence of our choices of a number of parameters on the simulation. To test the influence of the HRD prior, we processed our default simulation, done with the Gaia empirical HRD, using the isochrone-based HRD described in Sect. 3.2. We see in Fig. 8 that the A₀(D) relation is reasonably well recovered although slightly shifted. The bar overdensity is still visible in the ρ(D) distribution but is noisier.

We tested the influence of the extinction law adopted by processing our default simulation, constructed with the Fitzpatrick & Massa (2007) extinction law, assuming in FEDReD the Cardelli et al. (1989) extinction law. We see in Fig. 8 that the results are quite similar to the HRD change.

Concerning our completeness model, FEDReD estimates magnitude limits slightly different from the input ones through the estimation using the maximum of the observed magnitude distribution, but they are still within 0.4 mag of the input ones. We checked that providing the exact input completeness values did not noticeably change the results. We also tested changing the α and β parameters to 10 and 1 respectively (e.g. the 2MASS sharper values) for the UKIDSS simulation. Figure 8 shows that this only affects, as expected, ρ(D) and not A₀(D).

4.2. Field NGC 4815

We looked at the FEDReD capabilities in the field around NGC 4815, which has been studied in extinction with the Gaia-ESO Survey UVES observations by Puspitarini et al. (2015). We used 711 2MASS stars located in an area of 0.1°x0.1° around l = 306.6°, b = −2.1°; 87% of those stars have Gaia parallaxes. This field is complex for FEDReD as it suffers from differential extinction requiring a very small field of view and therefore a small number of stars, and has the presence of a cluster that differs from the mix of age and metallicities of our empirical HRD. To compare our results in Fig. 9 with the ones of Puspitarini et al. (2015), we updated the distances of the latter with the Gaia DR2 distances using the inverse of the parallax. We also compared our results with the maps of Marshall et al. (2006) and Lallement et al. (2019), our results being in between both maps, with a better agreement with Puspitarini et al. (2015). This field is indicated as having a convergence issue in Green et al. (2019). The updated Gaia DR2 distances confirm that the two stars with lower extinction are foreground stars, as suspected by Puspitarini et al. (2015). We do not recover the same shape at small distances as Lallement et al. (2019), which can be due either to a too relaxed definition of D_min on our side, considering the very few stars present in our small field of view to drive the solution, or to the too big resolution of the Lallement et al. (2019) map for this specific area. We confirm that a dust cloud is present at the cluster distance. We also confirm that the extinction continues to increase beyond the cluster, in phase with the higher velocity interstellar medium structures seen in the HI data and not detected in the stars studied by Puspitarini et al. (2015). The extinction is likely to continue to increase beyond our distance limit as we do not reach the total extinction of 4.4 mag indicated by the map of Schlegel et al. (1998). Concerning the stellar density, we recover the overdensity linked to the presence of the cluster, which we estimate to be at 3.5 ± 0.1 kpc. This is consistent with the results of Cantat-Gaudin et al. (2018). Using the isochrone HRD instead of the empirical Gaia one leads to consistent results within the uncertainties.

Fig. 9. Field of NGC 4815. Top: 2MASS CMD. In green the red clump track corresponding to our results with its 1σ confidence interval, in blue the Marshall et al. (2006) results and in red the Lallement et al. (2019) ones. Middle: extinction, 1σ confidence interval in grey (see Fig. 8). Dotted line: FEDReD result using the isochrone HRD. Black points: Puspitarini et al. (2015) updated with the Gaia DR2 distances, in red for members according to Friel et al. (2014). Bottom: stellar density. The cluster distance (Cantat-Gaudin et al. 2018) is indicated in red.

4.3. Field 9P

We looked at the capability of FEDReD to detect the bar signature using the field l = 9.55° ,b = −0.09° studied in detail in Babusiaux & Gilmore (2005) with CIRSI near-infrared photometry and in Babusiaux et al. (2014) with GIRAFFE spectroscopy. We took an area of 0.16°   ×  0.16° leading to about 10 000 UKIDSS stars with a photometric uncertainty lower than 0.1 mag in J, H, and K, that is up to J = 19, H = 17.9, and K = 16.8 mag. To improve the convergence at small distances, we completed UKIDSS with 2MASS photometry and we replaced the UKIDSS photometry with the 2MASS one for stars brighter than J = 13.25, H = 12.75, and K = 12.0, following Lucas et al. (2008). For this we derived and applied colour-colour calibrations on well-behaved stars of both surveys following Hodgkin et al. (2009). We used the isochrones generated with the UKIDSS filters and we transformed the empirical HRD from 2MASS to the UKIDSS photometric system using the transformations of Hodgkin et al. (2009). We used the same extinction coefficients as previously (i.e. Lallement et al. 2019). In this field, the over-density in stellar counts due to the Galactic bar occurs brighter than the completeness limit. We checked through simulations that our default way to estimate the completeness parameters presented in Sect. 2.3 is indeed adapted to this field.

The results are presented in Fig. 10. Both the empirical HRD and the synthetic one give consistent results within the uncertainties. Our results are barely overlapping in distance with the ones of Lallement et al. (2019) but are consistent within the uncertainties. We find a higher extinction than Marshall et al. (2006), more in agreement with the results of Babusiaux et al. (2014). The red clump track, clearly visible in the CMD, is well recovered. We confirm the increase in extinction in the disc up to the bar location seen in Babusiaux et al. (2014) and see the decrease of the extinction material afterwards. We see the location of the bar-driven overdensity at 4.9 ± 0.2 kpc, which is consistent with the value of Babusiaux & Gilmore (2005). We also confirm the bar’s large distance spread. This large dispersion could be due to us seeing both the disc end and the bar, which are too close to be separated. The main increase in extinction seems to be slightly in front of the density peak, in agreement with what would be expected if we are seeing here the bar close to reaching the disc. The extension of the method to other longitudes to constrain the bar/disc interface will be presented in a forthcoming work.

Fig. 10. Field l = 9.6° ,b = 0°. Top: UKIDSS CMD. In green the red clump track corresponding to our results with its 1σ confidence interval, in blue the Marshall et al. (2006) results, in green the Lallement et al. (2019) ones. Middle: extinction, 1σ confidence interval in grey (see Fig. 8). Dotted line: FEDReD result using the isochrone HRD. Thin black line: isocontours of the spectroscopic sample results of Babusiaux et al. (2014). Bottom: stellar density. The bar distance determined by Babusiaux & Gilmore (2005) is indicated with a red line and the distance spread in light red.

5. Conclusion

We presented here a Bayesian deconvolution method, FEDReD, allowing us to derive the extinction distribution and stellar density maps at the same time, taking into account the incompleteness of the surveys. We showed the performances of the algorithm on simulated data and on two test fields, one using 2MASS data centred around NGC 4815 and another using UKIDSS data towards the galactic bar at l = 10°. The first full application of the method to construct an extinction map of the Galactic disc using 2MASS and Gaia DR2 is presented in Hottier et al. (2020). Applications to UKIDSS and VVV (Surot et al. 2019) data are underway.

FEDReD is quite robust to differential extinction for its extinction derivation part, since it converges towards the median extinction behaviour. It is however important to select an homogeneous extinction behaviour to recover correctly the density distribution. We have seen towards NGC 4815 that it can work with a rather limited number of stars. Using the Gaia DR2 empirical HRD provides an accurate description of the local HRD, which we have shown to work well towards different parts of the disc. Still, variations of the HRD within the disc (metallicity gradient, changing ratio of thin/thick disc) can be preferred and implemented easily within FEDReD using the isochrone module. FEDReD has been designed to be flexible in its observable inputs so that any other knowledge of stars in the field of view can be implemented, such as spectroscopic and asteroseismology data.

Acknowledgments

We thank the referee for useful comments that helped to improve the paper. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work makes use of data products from the 2MASS, 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. This work is based in part on data obtained as part of the UKIRT Infrared Deep Sky Survey.

References

Appendix A: Individual P(O_j|A₀, D)

Fig. A.1. Individual P(Oj|A0, D)P0(D) for various stellar types indicated in Fig. 5 using as observables J, H, K photometry only, displayed with a square-root colour scale. The prior on stellar density is chosen here to be uniform, e.g. containing only the cone effect: P0(D) ∝ D2. All stars are located at 4 kpc with an extinction A0 = 3 mag, with the exception of the red dwarf, which is located at 0.1 kpc without extinction. The real position of the star is indicated by a white point. We see in this plot that the information is mostly carried by the red clump stars and that the Gaia parallax and/or photometry is needed to differentiate a red clump star from a red dwarf.

All Figures

Fig. 1. Resulting $\widehat{P}(A_0|D)$ from the first deconvolution. The green line is the real relation of the simulation. The black line is the resulting relation A0(D) (obtainted through a constrained median cubic spline fit on the MCSs). The red dotted lines are the minimum and maximum envelop of the MCSs. For the second deconvolution, P0(A0, D) is initialised with a zero probability outside of this red envelop.

In the text

	Fig. 2. Resulting $\widehat{P}(A_0|D)$ from the second deconvolution. The black line is the final A(D) estimated relation. Dotted lines corresponds to the 1σ CI. The green line is the real relation of the simulation. The vertical lines correspond to Dmin and Dmax, and show the valid distance range derived from the red clump star saturation magnitude and completeness limit respectively.
In the text

	Fig. 3. 1σ confidence interval of the derived stellar density ρ(D) as obtained by bootstrapping. We show here the results for two mock catalogues of the same area at l = 10°, dark grey: 2MASS-like, light grey: UKIDSS-like. The green line indicates the input simulation stellar density.
In the text

	Fig. 4. Relative stellar densities along the HRD P(MG) as derived from Gaia data only (P(MG, GBP − GRP), black dots), Gaia and 2MASS (P(MG, J − K), red stars), and Padova isochrones (blue triangles). The different densities are normalized at MG = 6.5.
In the text

	Fig. 5. Gaia DR2 empirical H-R diagram. Here we overlay some stars at different evolutionary stages for which we study their individual P(O|A0, D) in Appendix A.
In the text

	Fig. 6. Alternative H-R diagram to Fig. 5 based on Padova isochrones.
In the text

	Fig. 7. Colour-magnitude diagram of the UKIDSS mock catalogue built from the Fux model stellar and gas particle distributions towards l = 10°. The magnitude independent of extinction KJK is used. The corresponding distance and extinction for a red clump star on this diagram are indicated on the right and top axis.
In the text

Fig. 8. Results of the deconvolution of the simulation for A0(D) (top) and ρ(D) (bottom) within the [Dmin,Dmax] distance range. The black line is the deconvolution result. The light grey area in the top panel corresponds to the 1σ confidence interval of A0(D) derived from the full P(A0|D) (Fig. 2) while the darker grey area in both panels shows the 1σ confidence interval derived from the bootstrap. The green line is the input relation used in the simulation. Red dashed line: isochrones HRD. Blue dotted line: Cardelli et al. (1989) extinction law. Orange dot-dashed line: assuming completeness parameters α = 10 and β = 1.

In the text

Fig. 9. Field of NGC 4815. Top: 2MASS CMD. In green the red clump track corresponding to our results with its 1σ confidence interval, in blue the Marshall et al. (2006) results and in red the Lallement et al. (2019) ones. Middle: extinction, 1σ confidence interval in grey (see Fig. 8). Dotted line: FEDReD result using the isochrone HRD. Black points: Puspitarini et al. (2015) updated with the Gaia DR2 distances, in red for members according to Friel et al. (2014). Bottom: stellar density. The cluster distance (Cantat-Gaudin et al. 2018) is indicated in red.

In the text

Fig. 10. Field l = 9.6° ,b = 0°. Top: UKIDSS CMD. In green the red clump track corresponding to our results with its 1σ confidence interval, in blue the Marshall et al. (2006) results, in green the Lallement et al. (2019) ones. Middle: extinction, 1σ confidence interval in grey (see Fig. 8). Dotted line: FEDReD result using the isochrone HRD. Thin black line: isocontours of the spectroscopic sample results of Babusiaux et al. (2014). Bottom: stellar density. The bar distance determined by Babusiaux & Gilmore (2005) is indicated with a red line and the distance spread in light red.

In the text

Fig. A.1. Individual P(Oj|A0, D)P0(D) for various stellar types indicated in Fig. 5 using as observables J, H, K photometry only, displayed with a square-root colour scale. The prior on stellar density is chosen here to be uniform, e.g. containing only the cone effect: P0(D) ∝ D2. All stars are located at 4 kpc with an extinction A0 = 3 mag, with the exception of the red dwarf, which is located at 0.1 kpc without extinction. The real position of the star is indicated by a white point. We see in this plot that the information is mostly carried by the red clump stars and that the Gaia parallax and/or photometry is needed to differentiate a red clump star from a red dwarf.

In the text