FEDReD I: 3D extinction and stellar maps by Bayesian deconvolution

Context. While Gaia enables to probe in great detail the extended local neighbourhood, the thin disk structure at larger distances remains sparsely explored. Aims. We aim here to build a non-parametric 3D model of the thin disc structures handling both the extinction and the stellar density simultaneously. Methods. We developed a Bayesian deconvolution method in two dimensions: extinction and distance. It uses a reference catalogue which completeness information defines the selection function. It is designed so that any complementary information from other catalogues can be added. It has also been designed to be robust to outliers, frequent in crowded fields, and differential extinction. The prior information is designed to be minimal: only a reference H-R diagram. We derived for this an empirical H-R diagram of the thin disk using Gaia DR2 data and synthetic isochrone-based H-R diagrams can also be used. Results. We validated the method on simulations and real fields using 2MASS and UKIDSS data complemented by Gaia DR2 photometry and parallaxes. We detail the results of two test fields: a 2MASS field centred around the NGC 4815 open cluster which shows an over-density of both extinction and stellar density at the cluster distance, and a UKIDSS field at $l=10^{\circ}$ where we recover the position of the Galactic bar.


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 through the extinction.
Several methods have now been developed to draw 3D extinction maps. Fully model based (e.g. Drimmel & Spergel 2001), using a stellar distribution model but deriving a nonparametric 3D extinction map (Marshall et al. 2006;Chen et al. 2013;Schultheis et al. 2014, using the Besançon model), the distribution of stars near the main-sequence turn-off (Gontcharov 2017) or, the most common, using individual stellar distance and extinction estimates which are then inverted (Arenou et al. 1992;Vergely et al. 2010;Lallement et al. 2019;Berry et al. 2012;Chen et al. 2014;Hanson & Bailer-Jones 2014;Rezaei Kh. et al. 2017;Chen et al. 2019;Anders et al. 2019). Green et al. (2014) samples the full probability density function of distance, 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) uses a full hierarchical model to handle simultaneously the mean-distanceextinction 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, e.g. assuming that most of the extinction occurs in the foreground of the structure under-study: 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 been searching for the red clump position using a magnitude 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 is attempted in the Galactic disk. 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 which can provide complementary information on the stars observed. We choose in the description and applications presented here to use nearinfrared 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 et al. 2016;Gaia Collaboration et al. 2018b) as complementary information. We present in Sect. 2 the method, in Sect. 3 the H-R diagram (HRD) priors we constructed and in Sect. 4 results in both a simulated field and selected test fields.

Bayesian deconvolution
We wish to derive the probability distribution P(A 0 , D) which gives the probability of a star to have both an interstellar extinction A 0 (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. P(A 0 |D) gives the variation of the extinction with distance and P(D) gives the stellar density distribution along the line of sight.
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 NIR magnitudes and potentially optical magnitudes and parallax). What we which to derive is The sum is discrete instead of the usual integral as we are observing a finite number of stars.
We have through Bayes' theorem: Following the well-known Richardson-Lucy deconvolution algorithm, we can estimate P(A 0 , D) by iteratively computing h(A 0 , D) (Lucy 1974, with The initial values h 0 (A 0 , 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 ). What we can compute iteratively is then in fact P(A 0 , D|S ): still with O j being an observed star with at minimum NIR magnitudes observables. P(O j |S ) is the probability of an observed star to be in the selected sample. Similarly to Eq. 2 we have The observed star being actually observed, P(S |O j ) = 1. We therefore estimate P(A 0 , D|S ) by iteratively computing h(A 0 , D|S ): with ψ k (A 0 , D) = P(O j |A 0 , D) h k (A 0 , D|S )/P(S |A 0 , 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 0 , D|S ) which we will note in the followingP(A 0 , D|S ) = h K (A 0 , D). We can then retrieve an estimate of P(A 0 , D) witĥ 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-match between observations in different filters and catalogues can be numerous. Those false observations imply to use a robust method to deriveP(A 0 |D) while they should not impact significantlyP(D).

Individual probabilities
To derive each star's probability P(O j |A 0 , 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 0 , has an apparent magnitude m i : 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 0 through the formalism of Danielski et al. (2018), using the same coefficients as Lallement et al. (2019): 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.
To compute P(m|m i ) and P(˜ | ) 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 0 , D) for a thin 2-D 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 0 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.
Article number, page 2 of 12 C. Babusiaux et al.: FEDReD I: 3D extinction and stellar maps by Bayesian deconvolution 2.3. Selection function model P(S |A 0 , D) As previously, we compute the probability to be selected using the H-R diagram prior: We adopted the following model for the completeness of the surveys which is generic enough to reproduce simply typical completeness curves: P(S |m) is the probability for a star to be observed in a given photometric band given its magnitude true m. α, β and m * are three parameters defining 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 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 2 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 not valid anymore when a large stellar density is present near the end of the completeness survey, e.g. 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.
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 to be selected through the cumulative distribution function of the magnitude errors σ m derived with Eq. 14: For UKIDSS the global selection probability is then the product of Eq. 13 and 15.

Initial values
The construction of the initial value of the iteration, h 0 (A 0 , D|S ), needs two initial conditions: P 0 (D) and P 0 (A 0 |D). Then it is simply computed as: We only take into account the cone effect on a constant underlying stellar density profile: P 0 (D) ∝ D 2 . We tested also using a disc exponential profile start but that does not influence at all our results and we therefore stay with the simple flat start.

P 0 (A 0 |D)
Here also we use a flat start : P 0 (A 0 |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 distance (typically below 1 kpc but depending on the cone aperture chosen). To avoid this, we tested that adding an extra simple local prior could help, e.g. that the extinction cannot be too high locally: P(A 0 |D) = 0 for A 0 > 10D, with D the distance in kpc. The local map of Lallement et al. (2019) confirms that this is a very safe prior. But 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 starts only at relatively large distances and with very few 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.

Convergence
Assessing the convergence of such a deconvolution is always tricky. We decided to stop the iterations when the convergence rate slows down.
The 0.05 threshold for the convergence criteria is somewhat arbitrary but was tested on both simulations and real data to enable to reach the final shape of theP(A 0 |D) relation without introducing too much noise in the overallP(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, although the resulting ma-trixP(A 0 , D) is quite noisy, the algorithm to recoverP(A 0 |D) (Section 2.6) is robust enough to cope with these areas.

Deriving A 0 (D)
We deriveP(A 0 |D) fromP(A 0 , D) obtained at the end of the deconvolution process witĥ Article number, page 3 of 12 A&A proofs: manuscript no. chloe We now search for the relation A 0 (D) corresponding to the maximum of the probabilityP(A 0 |D) with the physical constraint that A 0 (D) should increase with D. For this, we randomly generate 10 000 monte-carlo solutions of A 0 (D) (called MCS hereafter) according to the following algorithm. We randomly select a distance D l within our working area following the probability weightsP(D|S ) = A 0P A 0 , D|S ) dA 0 . For each distance D l we randomly select a corresponding extinction within the possible values allowed by the increasing constrain and the extinctions already assigned to the previous distances. This random selection of A 0 (D l ) is done using the probability weights defined bŷ P(A 0 |D l ). We initiate the generation by setting A 0 (0) = 0 and A 0 (D max ) = A max . We compute the total log-likelihood of a MCS by log(L) = l log(P(A 0 (D l ), D l )). We then select the best 1 000 MCSs. We re-generate 10 000 random MCSs but this time further constraining the solutions to be within the extinction envelop of the 1 000 best MCSs for each distance. If the log likelihood of the new generated solution is better than the worse of the MCSs, the new solution replaces it. Finally a median cubic spline fit with an increasing constrain (Ng & Maechler 2007) is applied on the final 1 000 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 MCSs distributions.
As illustrated in appendix A, red clump stars are providing the strongest constraints on the distance/extinction distribution thanks to their small intrinsic dispersion in absolute magnitude and colour. To avoid noise induced by degenerate solutions, we therefore 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 computed only 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 derivedP(A 0 |D).
A second deconvolution process is done using the first one as the initial probability distribution, e.g. setting P 0 (A 0 , D) with a Gaussian distribution according to the MCSs quantiles, zero outside the MCSs envelop, and a flat probability outside D min /D max just taking into account the maximum / 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. Such a second step is needed mainly to derive accurately the distance distribution, the first pass deriving already well the extinction/distance 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 /D max to exclude distances too far away or too close within our cone angle to have enough stars:P(D < D min | S ) > 0.01 and P(D < D max | S ) < 0.99. This last restriction is needed in particular for anticentre areas with a low intrinsic stellar density at large distance.  The green line is the real relation of the simulation. The vertical lines correspond to D min and D max , e.g. show the valid distance range derived from the red clump star saturation magnitude and completeness limit respectively.

DerivingP(D)
We simply computeP(D) = A 0P (A 0 , D) dA 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 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 ρ(D) ∝P(D)/D 2 . The result can be seen in Fig. 3   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, e.g. 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 2 sigma but quite noisy, in particular at small distances.

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 IMF, metallicity nor age priors and 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, e.g. spectroscopic information, or if we want to test the impact of variations of the HRD within the Galaxy. An other 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).

The 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 tradeoff between having enough stars to sample the giant branch and staying as close as possible to the Galactic plane, e.g. within a relatively homogeneous stellar population mixture. We select only low extinction stars to build the empirical HRD to avoid adding uncertainties associated with any 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, which is the case for the Galactic disk fields for which this HRD has been built. We applied the same astrometric and photometry filters as in Appendix B of Gaia Collaboration et al. (2018a) with the exception of the photometric flux error one 1 . 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 an Hess Diagram on a grid of G BP − G RP colour (step 0.01 mag) and M G (step 0.02 mag), e.g. P(M G , G BP − G RP |S ), with S the HRD stars selection function. To correct for the selection function we derive As previously, we call i one of those HRD point (M G , G BP −G RP ).
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 2 cos(b). We move our integral in parallax space instead of distance as this is our observable: 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, 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.
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 The |Z| < 50 pc constraint corresponds to > | sin b|/0.05.
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 7 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 2 within each population.
For a better modelling of the bottom of the HRD, quite important for the pollution of nearby red dwarfs in our CMDs, 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 strong 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 on this work which concentrate on more distant stars, they only need to be present with the correct colour-colour relations to avoid them ending as strong outliers.

Theoretical H-R diagram prior
To build a theoretical HRD, we use the PARSEC isochrones (Bressan et al. 2012   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 SFR. We are observing here 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 2 profile: If we integrate over z and assume that all the populations have the same radial density profile, we have and therefore a constant SFR corresponds to a local density 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.

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, as the ones described below, and at several ones presenting different issues (high crowding, low stellar density towards the anti-centre, convergence issues...). For those, we checked how well our derived A 0 (D) function permit 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 correct the errors provided with the UKIDSS catalogue by the following: and we selected only 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.

Simulation
We tested our procedure on a simulation, as illustrated in Fig. 1 and Fig. 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 =32 mag. An intrinsic dispersion in the extinction is added using a log-normal distribution with σ A =0.05. 2MASS, UKIDSS and Gaia photometric and parallax errors are added assuming a simple increase of the parallax 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. Gaia G photometry and parallaxes are kept only if it satisfy the same completeness model as Eq. 13 with m * G = 20.7 and G BP and G RP photometry with m * G BP = 20.9 and m * G RP = 19.5 (the exact values are not important as they do not enter the catalogue completeness model, but just allows us to take into account that Gaia information is not present for the faintest / reddest stars). We build this way mock catalogues of about 4 000 stars satisfying our photometric criteria. 10% 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): with k J and k K the extinction coefficients in the J and K bands respectively. Figure 8 shows in dark and grey 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 0 |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 continuous 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.
We tested the influence of our choices of a number of parameters on the simulation.
To test the HRD prior's influence, we processed our default simulation, done with the Gaia empirical HRD, using the isochrone-based HRD described in Sec. 3.2. We see in Fig. 8 that the A 0 (D) relation is reasonably well recovered although slightly shifted. The bar overdensity is still visible in the ρ(D) distribution but is noisier.
Article number, page 7 of 12 A&A proofs: manuscript no. chloe 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 change sensibly the results. We also tested changing the α and β parameters to 10 and 1 respectively (e.g. the 2MASS sharper values) for the UKIDSS simulation. Fig. 8 shows that this affect, as expected, only ρ(D) and not A 0 (D).

Field NGC 4815
We looked at the FEDReD capabilities in the field around NGC 4815 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 having Gaia parallaxes. This field is complex for FEDReD as it suffers from differential extinction requesting a very small field of view and therefore a small number of stars, and has the presence of a cluster which will differ 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 compare 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 2 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 from 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 ISM 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 Article number, page 8 of 12 C. Babusiaux et al.: FEDReD I: 3D extinction and stellar maps by Bayesian deconvolution be at 3.5±0.1 kpc, which 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.

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, e.g. 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 by the 2MASS one for stars brighter than J = 13.25, H = 12.75, 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 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 disk 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 large distance spread. This large dispersion could be due to us seeing both the disk end and the bar, 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 disk. The extension of the method to other longitudes to constrain the bar/disk interface will be presented in a forthcoming work.

Conclusion
We presented here a Bayesian deconvolution method, FEDReD, allowing us to derive at the same time the extinction distribution and stellar density maps 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 disk using 2MASS and Gaia DR2 is presented in Hottier et al. (2020). Applications to UKIDSS and VVV 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 disk (metallicity gradient, changing ratio thin/thick disk) 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 for some stars of the field of view can be implemented such as spectroscopic and asteroseismology data.