© R. Lallement et al. 2019

1 Introduction

Evolutionary models of the Milky Way require measurements of the spatial distribution of massive numbers of stars as well as their dynamical, physical, and chemical properties. These quantities are currently provided by the ESA satellite Gaia (Gaia Collaboration 2016, 2018b) and by ground-based surveys. In this context, 3D maps of the Galactic interstellar (IS) dust are an additional and mandatory tool because they allow dereddening of stellar spectra and suppression of degeneracies between stellar temperature and dust reddening on the one hand, and because they may shed important additional light on the star formation through the interplay between star and IS matter on the other hand. Fortunately, Gaia parallaxes and Gaia photometric data complemented by ground-based photometric data allow the construction of the required IS dust maps in parallel with the stellar studies. 3D mapping is based on the tomographic inversion of distance-limited IS absorption data for large numbers of targets that are distributed in space at known locations: Gaia parallaxes evidently provide the target locations, and Gaia together with ground-based photometric measurements provide estimates of the reddening along each sightline. Independently of these evolutionary aspects, 3D maps of the Galactic dust are a general tool for many various purposes, such as studies of foreground, environment, or background to specific objects, and models of light or particle propagation.

One form of 3D dust mapping is the construction of radial profiles of color excess or extinction, sightline by sightline. The first map produced in this way was based on HIPPARCOS data (Arenou et al. 1992). Later on, maps of radial profiles and their derivatives, that is, reddening or extinction per unit distance, were built by Marshall et al. (2006), Chen et al. (2013), and Schultheis et al. (2014) based on the Besançon model of stellar population synthesis (Robin et al. 2012) and on photometry obtained in the Two Micron All Sky Survey (2MASS), the Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE), and the Vista Variables in the Via Lactea (VVV) survey, respectively. Majewski et al. (2011) derived individual stellar reddening and radial profiles by combining data from Spitzer/GLIMPSE and 2MASS. Berry et al. (2012) used the Sloan Digital Sky Survey (SDSS) and 2MASS to constrain stellar spectral energy distributions (SEDs) and reddening profiles with reddening-free empirical spectra. Following a hierarchical Bayesian method devised by Sale (2012), Sale et al. (2014) used photometric data from the INT Photometric H-Alpha Survey (IPHAS) to simultaneously reconstruct stellar properties and radial extinction profiles at a 10′ resolution for −5 ≤b ≤ 5° and 30 ≤l ≤ 210° up to 5 kpc. Green et al. (2015, 2018) used data from 2MASS and the Panoramic Survey Telescope And Rapid Response System (Pan-STARRS) 1 and a refined Bayesian method to derive radial reddening profiles at very high angular resolution (on the order of 7′) for the sky regions accessible to Pan-STARRS (0 ≤l ≤ 240° at low latitudes). Both Sale et al. (2014) and Green et al. (2015, 2018) used derivatives of extinction profiles to produce 3D maps of the dust density. Kos et al. (2014) used spectrometric and photometric data of the RAdial Velocity Experiment (RAVE) and presented reddening and also diffuse interstellar band (DIB) strength radial profiles for the longitude interval l = 180 to l = 60°. Recently, Chen et al. (2019) used Gaia DR2 parallactic distances and combined photometry from Gaia DR2, 2MASS, and the Wide-field Infrared Survey Explorer (WISE) to derive extinction profiles along the Plane at high angular resolution (6′) by means of a new machine-learning algorithm trained on datasets obtained with the Large Sky Area Multi-Object Fibre Spectroscopic Telescope (LAMOST), the Sloan Extension for Galactic Understanding and Exploration (SEGUE), and the APO Galactic Evolution Experiment (APOGEE).

Full 3D inversions differ from the previous technique by imposing spatial correlations between volume densities of IS matter in all directions, thereby linking adjacent sightlines and deriving the volume density of the differential extinction everywhere in 3D space. The differential extinction at point P in 3D space (also called extinction density) is the amount of extinction per unit distance that stellar light is subject to when it travels in the region enclosing P, here expressed in magnitude per parsec. In the case of the nearby interstellar medium (ISM), 3D inversions following a technique developed by Vergely et al. (2001) have been applied to color excesses from ground-based photometry (Vergely et al. 2010; Lallement et al. 2014) or composite sources combining photometry and DIBs (Capitanio et al. 2017; Lallement et al. 2018). Sale & Magorrian (2014, 2015, 2018) developed a Gaussian field method adapted to realistic multiscale IS matter distributions. Rezaei Kh. et al. (2017, 2018) developed a nonparametric 3D inversion method based on an isotropic Gaussian process to derive the 3D dust distribution. They tested and then applied their new technique on extinction-distance measurements derived from the APOGEE/Kepler Asteroseismic Science Consortium (APOKASC) catalog of Rodrigues et al. (2014) and APOGEE data, respectively.

In Sect. 2 we estimate extinctions based on Gaia and 2MASS photometric measurements and apply various selection criteria and filters. Section 3 details the hierarchical inversion technique, its limitations, and the chosen parameters, and we present the resulting multiresolution map by means of several images of the dust density in various selected planes. In Sect. 4 we compare the dust distribution within the Plane with tracers of star formation and molecular clouds. Finally, in Sect. 5 we discuss perspectives for 3D dust mapping. The appendix shows some comparisons between reddening profiles obtained by integration through our 3D dust distribution and those derived by Green et al. (2018) based on Pan-STARRS and 2MASS photometry.

2 Deriving the extinction

To obtain the extinction, we used the intrinsic color–color relations and modeled the extinction coefficients in the following way. To derive the extinction of individual stars, we combined photometry from Gaia DR2 (Gaia Collaboration 2018b) and 2MASS (Skrutskie et al. 2006). We used the catalog cross-match provided within the Gaia archive (Marrese et al. 2018). We applied the Gaia photometric and astrometric filters as provided in Gaia Collaboration (2018a; photometric uncertainties lower than 2% in G, 5% in G_BP and G_RP, and filters on the photometric excess factor and an equivalent to the astrometric excess noise filter). We did not apply the filter on the number of visibility periods because it removed specific regions of the sky while the most critical outliers (negative parallax and too large parallaxes) are already removed by our other filters. We then selected stars whose relative parallax uncertainties were lower than 20%. For 2MASS we selected stars with a 2MASS photometric quality flag AAA and photometric uncertainties lower than 0.05 mag. Only stars at the top of the Hertzsprung–Russell diagram (HRD) (M_G < 5) were selected to probe large distances and minimize the influence of binaries. The photometry was handled as suggested by Weiler (2018), with two different zero-points, filters, and photometric calibrations for the G _BP photometry for G brighter and fainter than 11 mag and a correction for the G magnitude drift (Arenou et al. 2018; Weiler 2018; Casagrande & VandenBerg 2018) of 3.5 mmag per magnitude. Stars brighter than G < 6 mag were removed to avoid saturation problems and stars fainter than G_BP < 18 were removed to avoid background subtraction problems (Evans et al. 2018; Arenou et al. 2018). After we applied all these criteria, the total number of selected objects was ~27 340 000.

The photometric calibration was performed on low-extinction stars selected to have E(B−V) < 0.1 mag using the previous 3D extinction map of Capitanio et al. (2017). A seven-degree monotone polynomial¹ was adjusted to the color-color relations G−X as a function of G−K, removing strong (ten times the median absolute deviation (MAD) of the residuals) outliers one by one, as illustrated in Fig. 1. We calibrated bright stars using stars with G < 10.5 and applied this calibration to all stars with G < 11 mag. The color range of our calibration is $-1.0<{(G-K)}_0<5.6$.

The extinction coefficients were derived using the Fitzpatrick & Massa (2007) extinction law and the Kurucz SEDs (Castelli & Kurucz 2003), as described in Danielski et al. (2018). The 2MASS transmissions were taken from Cohen et al. (2003), and the Gaia transmissions are those of Weiler (2018). We here selected only the top of the HRD. We therefore adapted the surface gravity of the Kurucz spectra to the temperature with log g = 4 for T _eff > 5250 and log g = − 8.3 + 0.0023 T_eff for cooler stars (simple fit adjusted on APOGEE Majewski et al. 2017 data). We also increased the degree of the polynomial fit of the extinction coefficient model: $$\begin{array}{lll}{k}_{\text{m}}\hfill & =\hfill & a_1+a_{2}C+a_3{C}^2+a_4{C}^3+a_5{A}_0+a_6{A}_0^2+a_7{A}_0^3+a_8{A}_{0}C\hfill \\ \hfill & \hfill & +a_9{A}_0{C}^2+a_{10}C{A}_0^2\hfill \end{array}$$(1)

with $C={(G-K)}_0$. The fit was performed on a grid with a spacing of 250 K in T_eff with 3500 < T_eff < 1000 K and a step that linearly increased by 0.01 mag in A₀ with 0.01 < A₀ < 20 mag. This increase in step size with A₀ was made to ensure the best model of the extinction coefficients at low extinction values, which dominates our sample. The residuals are smaller than 0.6% for A₀ < 10 mag in the G band.

We tested the influence on our results of the choice of Gaia passband used to derive the extinction coefficients. The differences between the revised Gaia DR2 passbands provided by Evans et al. (2018) and those of Weiler (2018) are small. We preselected intrinsically bright stars using the 2MASS K_s magnitude, which is less strongly affected by extinction: $$K_{\text{s}}+5+5{\log }_{10}\left(\frac{\varpi +{\varsigma }_{\varpi }}{1000}\right)<4,$$(2)

where ϖ is the Gaia DR2 parallax. This cut allowed us to retrieve all stars with M_G < 5 up to 5 mag of A₀ extinction and covers the range of most of our extinctions. For each star, we determined extinction A₀ and intrinsic color ${(G-K_{\text{s}})}_0$ and their associated uncertainties through a maximum likelihood estimation², $$\mathcal{L}={\displaystyle \prod _{X}P}(G-X|A_0,{(G-K)}_0),$$(3)

for X = G_BP, G_RP, J, H. To avoid local minima, three different initial values were tested: ${(G-K)}_0=G-K$ (no extinction), ${(G-K)}_0=1.5,$ and ${(G-K)}_0={(G-K)}_{\text{max}}$ (the maximum color of our intrinsic color relation, e.g., 5.6 mag). We used the intrinsic color-color relations to derive ${(G-X)}_0$ from ${(G-K)}_0$. Then $P(G-X|A_0,{(G-K)}_0)=P(G-X|{(G-X)}_0+k[A_0,{(G-K)}_0]A_0)$. We modeled this probability by a Gaussian, quadratically adding the photometric error in the X band and the intrinsic scatter of the color–color relation. Negative values of A₀ were allowed to ensure a Gaussian uncertainty model, which is needed for the inversion method. However, we did not extrapolate the extinction coefficients because the negative values of A₀ were replaced by 0 to derive the extinction coefficients. The typical resulting uncertainties are 0.3 mag in A₀ and 0.2 mag in ${(G-K)}_0$.

A chi-square test was performed to check the validity of the resulting parameters, removing stars with ap-value limit lower than 0.05. We removed stars with uncertainties on the derived A₀ and ${(G-K)}_0$ higher than0.5 and 0.4 mag, respectively, as well as outliers for which A₀ was not compatible with being positive at 3 σ. To ensure that we remained within our calibration interval, we kept only stars with dereddened $M_{{\text{G}}_0}<5$. Inspection of the dereddened HRD showed that the bluest and reddest stars were not correctly recovered. We therefore also removed stars with ${(G-K)}_0<0.5$ and ${(G-K)}_0>2.9$ mag.

We tested our individual extinctions by checking the extinction distribution in low-extinction regions (Fig. 2) and the shape of the dereddened HRD (Fig. 3). We also compared our extinctions with previous work, Fig. 4 showing the comparison with extinctions derived from APOGEE spectroscopic data by Santiago et al. (2016).

Fig. 1 Example of an intrinsic color-color relation fit (top) and its residuals (bottom) for G − GBP vs. G − K. In red we plot the outliers that were removed in the process.

Fig. 2 A0 extinction distribution (mag) in low-extinction selections. In blue we show the local bubble selection (ϖ > 14 mas, corresponding to the first 70 pc), and in green we plot a high-latitude selection (|b| > 80°). The median residuals are smaller than 0.03 mag with a standard deviation of 0.2 mag.

Fig. 3 Dereddened Hess diagram for a subset of stars with a relative parallax uncertainty smaller than 10%. The grayscale corresponds to the square root of the stellar density. The Red Clump is the most prominent feature. Its shape can be compared to Fig. 2 of Ruiz-Dern et al. (2018).

3 Inversion of extinctions

3.1 Description of the inversion technique

The inversion method we used to estimate the local differential extinction (in mag.pc⁻¹ units) is derived from the method presented in Sect. 3 of Lallement et al. (2018). It was initially developed to construct a 3D distribution of local dust that was to be used as a complement to large-scale models in a new version of the Gaia Universe Model Snapshot (GUMS) model (Robin et al. 2015, 2012).

The basic principles of the method are as follows. The treatment is hierarchical from the point of view of the spatial resolution (or structure scale size), with increased resolution at each iteration. For each considered resolution, it is carried out in two steps: a first step consists of estimating the extinction in radial directions by means of a robust Bayesian inversion of individual extinctions, then a second step is a full 3D inversion of the results obtained in the first step to produce differential extinction in 3D space at the resolution corresponding to the current hierarchical level. The covariance functions for the two radial and 3D steps are adapted to the current resolution. At variance with previous works based on a covariance function that was a combination of two different functions and aiming at better representing both dense and more diffuse clouds, here a unique Gaussian covariance function is used in the inversions. The FWHM of the Gaussian kernel is equal to the chosen scale and subsequently decreases in the iterative steps. The hierarchical approach makes it possible to process a very large volume of data with a strongly varying spatial density of targets (i.e., number of stars per unit volume) that enter the inversion, while ensuring coherence of the information at a given scale.

The iterations proceed by progressively refining the solution (or equivalently, the resolution), meaning that this approach makes it possible to manage the larger scales first and then go to the finest scales. Because the data coverage is not homogeneous, the spatial resolution obtained in the final solution varies from one place to another. In order to estimate the extinction on a given spatial scale, we averaged extinction data in boxes organized along radial beams, with a box size d_b corresponding to half of the considered scale size, that is, when we considered the 200 pc scale, we computed the average extinction in boxes of 100 pc size using all the stars contained in each box. The boxes and the decomposition in radial beams are defined in the following way: the beam angle is the angle at which a segment the size of the box size d_b is seen from the Sun at mid-distance between the Sun and the extremity of the computational volume in the given direction. A box at distance d along the beam is a truncated cone whose summit angle is this previously defined angle and whose thickness along the radial direction is d_b. The decomposition of the 4π sphere in beams and the distribution of box centers in the radial radiation are such that boxes overlap by 50% in both radial and azimuthal directions to ensure a consistency between average extinctions of two adjacent boxes. When fewer stars are in a box than a chosen scale-dependent limit, the information from the box is not taken into account in the inversion. Here our series of scales was 500, 400, 300, 200, 100, 50, and 25 pc. The minimum number of stars required in a box depends on the scale and was 10, 10, 10, 10, 5, 3, and 2 for the 500 pc, 400 pc, 300 pc, 200 pc, 100 pc, 50 pc, and 25 pc scales, respectively. Figure 5 shows the decomposition in boxes with volumes similar to those of the mid-distance boxes used for the first and last iterative steps, along the Galactic plane. The number of targets in each bin is color-coded, allowing us to determine which cells were used (lower limit of ten (two) targets for the first (last) iteration). The right panel shows the final resolution achieved in the Plane, which results from the cell sizes and the minimum numbers of targets at the various steps. As expected, the resolution is mainly governed by the radial distance. We note, however, that it is particularly poor in several regions, for instance, for Galactic longitudes 25 or 85 °. We return to this point below. More precisely, for each scale n, the inversion algorithm is as follows:

• 1)

  We start by calculating the limits of the boxes of scale n along beams covering the whole sky. The angular aperture of the beams as the radial sampling along each beam depends on the scale.

• 2)

  Then, we assign boxes to all target stars and calculate the average extinction in each box with a weighting based on the individual errors. Moreover, the standard deviation of the extinction in each box is calculated from the variability of the extinction inside the box, that is, this variability is considered to be an extinction error on the computed average.

• 3)

  For each of the radial beams, the extinction density is inverted using the Bayesian method presented in Vergely et al. (2010). The beams are processed independently of each other. The autocorrelation function is a Gaussian kernel whose correlation length corresponds to the scale we considered. An important point is the choice of the prior: at each stage, the prior distribution is the interpolated 3D extinction density obtained at the previous stage of scale n − 1.

• 4)

  Extinction densities obtained along all different beams are assembled. An error is estimated based on the standard deviation (std) calculated among the nearest-neighbor extinction densities. This is the std divided by the square root of N, where N is the number of neighbors.

• 5)

  The extinction densities from all different beams are then inverted in 3D to update the 3D distribution obtained at the n − 1 scale and produce a 3D distribution at the scale n. The 3D density distribution obtained at the n − 1 scale is used as a prior distribution for the inversion. To do so, the update of the cube is made by an optimal nonlinear interpolation method, which imposes the positivity of the solution and integrates a 3-σ filtering (a data point that after the first iteration deviates from 3σ of the model is discarded). The autocorrelation function used in the 3D inversion is the same as that in step 3.

At the first step n = 1 (or at the first spatial scale of 400 pc), the 3D extinction distribution entered as a prior is a homogeneous model with exponential decay on both sides of the Galactic plane, with a characteristic height of 200 pc. Extinction data were those described in the previous section for which the distance from the solar plane Z is smaller than or equal to 400 pc and the distance along the radial axis X and rotation direction Y is smaller than 3 kpc (~16 130 000 objects). Distances were derived from Gaia DR2 parallaxes, after subtraction of the zero-point −0.03 mas (Arenou et al. 2018; Lindegren et al. 2018).

Fig. 4 Comparison between the extinction A0 derived here and AV determined by Santiago et al. (2016) using APOGEE data.

3.2 Limitations and biases of the inversion

This hierarchical inversion has suppressed or at least strongly attenuated the main limitation of our full 3D mapping technique (used recently by Capitanio et al. 2017; Lallement et al. 2018), namely the use of a unique correlation length. This single length could not be adequate simultaneously for low and high target space density, that is, in general, for small and large distances. When the correlation length was chosen too wide, then some information at small distances was not used. Conversely, when it was chosen too narrow, structures in distant regions were unrealistic. Still, the best 3D spatial resolution associated with the last step is limited by our target space density and does not yet allow us to reach the scale of the filaments revealed by 2D maps or the very high angular resolution of the Green et al. (2018) reddening profiles, as shown by the comparisons with Planck in Fig. 8 (see below). However, in principle, iterative steps can be further continued to achieve a better resolution without introducing artifacts far from the Sun. Searches for the lowest achievable size are in progress and will benefit from future Gaia data releases. This optimal size depends primarily on the spatial density of targets, but also on distance and extinction uncertainties.

A bias of this inversion is the use of distances that are directly obtained from Gaia parallaxes, at all iterative steps. This means that distances are calculated as the inverse of the DR2 parallax, after correction for the zero-point of −0.03 mas. It is well known that using this relationship underestimates actual distances. This has been studied in several ways and was recently reinvestigated by Luri et al. (2018) in the context of Gaia DR2. In our previous maps, most of the inverted cloud complexes were at small distances (e.g., ≤1 kpc), and, more importantly, for the large majority of targets beyond several hundred parsec, the distances were photometric. Therefore, the problem was neglected. With the current dataset, we aim at mapping structures up to 3 kpc, and all distances are parallactic. Biases may therefore be significant, and it is important to check the potential consequences of our simplification.

A first estimate can be made that is only based on the relative error on the parallax: following Arenou & Luri (1999), if errors on the observed parallaxes are Gaussian, then for small relative errors and in the absence of negative or null parallaxes, the expected bias on the distance, that is, the difference between the true distance and the expectation of distance using the inverse of the parallax 1∕ϖ, that is, E(d|ϖ = 1∕ϖ), can be approximated as $$d-E(d|\varpi )\approx \frac{1}{\varpi }\cdot {\left(\frac{\sigma (\varpi )}{\varpi }\right)}^2.$$(4)

We restricted the input source catalog to targets with relative errors on parallaxes lower than 20%, or $\frac{\sigma (\varpi )}{\varpi }\le 0.2$, therefore wecan assume that we are close to the small error regime and that this formula, which applies to the true parallax, can also provide an order of magnitude in the case of a measured parallax. At the outer boundaries of our map (3 kpc along the X, Y axes in the Plane) and for the faintest targets, the bias in distance should therefore be limited to a few hundred parsec.

To estimate the distance bias in a way that takes into account not only the parallax error, but also our target selections, we assumed that our method is equivalent to a weighted mean on the inverse of the parallax and simulated with a Monte Carlo the Gaia errors on the Gaia HRD with M_G > 5 shifted to different distances and extinctions. We applied the ϖ∕σ_ϖ > 5 criterion and simplified the magnitude cut to G < 20. The resulting bias as a function of distance and extinction is presented in Fig. 7. The figure shows that at fixed reddening, the effect of using 1∕ϖ and selecting the sources as we did is to stretch the distance scale. Quantitatively, it is clear from the figure that because most of our extinctions are lower than A₀ = 4, biases should be limited to 200 pc at 3 kpc and to less than 100 pc within 2 kpc.

These estimates imply that for the large majority of our input data and at all distances, biases should be smaller than the resolution currently achieved (see Fig 5), for instance, biases of ~150–200 pc at the outer boundaries of the map are smaller than the 300–500 pc (poor) resolution achieved at those distances. For this reason, we neglected this effect for our mapping here. An a posteriori correction of the maps is under study, but is beyond the scope of this article.

As described in previous work, the main advantage of our methods is the use of spatial correlations in all directions in3D space to achieve a global picture of the matter distribution. The price to pay is the lower limit on the sizes of the inverted structures. This shortcoming is increasingly important at increasing distances, and it is mandatory to understand its impact. A first way to estimate the quality of the 3D distribution is illustrated in Fig. 8. The top map is the latest 2D dust map from Planck, based on the optical thickness of the dust at 353 GHz (Planck Collaboration X 2016). The level of detail is extremely high, as expected. The middle 2D map represents the locations of all our input targets that are distant by more than 1 kpc. The color scale represents the reddening associated with each star. The level of detail reached by the Gaia-2MASS extinction dataset is also quite remarkable. The third 2D map shown at the bottom is obtained through integration of the 3D distribution of differential extinction (taken from the last stage of the hierarchical inversion) in radial directions spaced every 0.25° in longitude and latitude, from the Sun to the limits of our computational volume. It allows us to clearly determine which structures were captured by the whole inversion process, how they were smoothed, and which were smoothed out. All major structures that appear in the first two maps are retrieved, but the level of detail is inferior, as expected, and a few very small filamentary structures are no longer visible in this 2D representation. However, they have some effect on the 3D distribution, which is visible in images of the extinction density in selected planes, for example.

Fig. 5 Left panel: 250 × 250 × 250 pc bins similar in volume to those used in the first step of the inversion. Shown are the cells centered on the Plane. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The colorcode corresponds to the number of target stars per bin, i.e., a quantity proportional to the targetspace density. Only bins with more than ten targets are used at this stage. Middle panel: 12.5 × 12.5 × 12.5 pc bins used in the last step of the inversion. Only bins with more than two targets are used. Right panel: resolution achieved in the Plane based on the hierarchical inversion. In regions of high target space density, essentially close to the Sun, the resolution is 25 pc and is achieved after the last step. In regions of scarce targets, mainly at large distances, the resolution is coarse, may reach 500 pc, and is achieved at the first steps. These regions are not improved during the following steps.

Fig. 6 Inverted 3D distribution of differential extinction in the Galactic plane and as it is obtained after each of the first four iterative steps. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The achieved resolution is300, 200, 100, and 50 pc, respectively.

Fig. 7 Biases on the parallax distance in parsec from a Monte Carlo simulation adapted to our dataset.

Fig. 8 Top: Planck map of the Galactic dust optical thickness at τ353 GHz. Middle: estimated extinction toward all targets located beyond 1 kpc. Bottom: 3D differential extinction integrated from the Sun to the limits of the computational volume.

Fig. 9 Dust density along the Plane, based on the 3D distribution obtained after the last inversion step. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The color scale represents the differential extinction in units of mag per parsec. The dashed black line represents the distance beyond which the final resolution of 25 pc is not achievable due to target scarcity (see text). The dotted white contours correspond to a differential extinction of 0.003 mag.pc−1 and delimit the dense areas.

4 Main features and comparisons with other maps and ISM tracers

Figures 9–13 display images of the extinction density in various horizontal and vertical planes. In each of these figuresa white line separates regions where the target density is above or below the threshold of two per box chosen for the final resolution of25 pc. That is to say, regions beyond this limit have only few target and a resulting poorer resolution and should be considered with more caution. The limit has been smoothed for a better visualization. The horizontal plane containing the Sun (i.e., parallel to the Galactic plane and very close in distance, hereafter called solar plane or Plane) in Fig. 9 is particularly interesting because it covers a 6 × 6 kpc × kpc area and reveals the main nearby features. Figures 10 and 11 display images of the dust density in horizontal planes at various distances from the solar plane, from Z = −300 pc to Z = +300 pc in steps of 50 pc. Figures 12 and 13 display images of the dust densityin vertical planes containing the Sun and oriented along longitudinal axes varying in steps of 10°.

It is interesting to compare the 3D distribution of dust with tracers of star-forming regions and molecular clouds (MCs). In addition to the density maps alone, we therefore also display the locations of several tracers of star formation and MCs, each time superimposed on the density map of the solar plane. We display the dust density in logarithmic scale for clarity, and it is color-coded in a different way, with reduced contrast.

• Figure 14 displays the locations of nearby masers from Reid et al. (2014), Xu et al. (2016) and Honma et al. (2012) that possess parallax distances and are located within 50 pc from the Plane. In addition to the masers, the spiral arms, which were adjusted to the data by Reid et al. (2014) and Hou & Han (2014), are also displayed; here, the Perseus, Local and Carina–Sagittarius arms as well as the Local Spur inferred by Xu et al. (2016).

• Figure 15 shows the projections onto the Plane of the HII regions from the catalogs of Russeil (2003) and Hou & Han (2014). The Russeil (2003) star-formation regions are deduced from Hα, H109α, CO, radio continuum, and absorption lines, and distances are based on parallaxes or photometric or kinematic data, depending on the objects. Additional HII regions compiled by Hou & Han (2014) that correspond to various transitions of H, He, CS, and CO and more recent observations (see references in Hou & Han 2014) are also represented. We restrict the figure to objects with distances estimated from stellar parallaxes and/or stellar photometric determinations that are located within less than 50 pc from the solar plane.

• Figure 16 displays the O stars from the Reed (2003) catalog that have very recently been cross-matched for their distances with Gaia DR2 by Xu et al. (2018) and lie within 50 pc from the Plane.

• Figure 17 displays the Miville-Deschênes et al. (2016) mole- cular clouds derived from CO spectral-cube decomposition and kinematic models, completed by the fraction of the giant molecular clouds listed by Hou & Han (2014) that are located within 3 kpc. Again, only objects closer than 50 pc from the Plane are retained.

We discuss the most prominent dust cloud structures and their degree of agreement with the various tracers, based on the images quoted above. Discussions of specific regions are beyond the scope of this article.

In the solar plane in Figs. 9 and 14, a first prominent region is the well-delineated inner part of the Carina–Sagittarius arm in the fourth quadrant and half of the first quadrant. Its series of compact dust cloud complexes is very well aligned in a l = 45–225° direction, that is, a direction that is significantly tilted from the expected spiral arm direction. This orientation is already partly visible in the maps of Kos et al. (2014) and Chen et al. (2019). In the fourth quadrant the clouds in front of this huge complex are as close as 1 kpc to the Sun. A more distant (~2 kpc and more) and less compact dense dust area is also found beyond the first rank of clouds, from l = 300 to l = 0°. The more diffuse aspect probably partly arises because this region is close to the limits of our maps and is mapped at lower resolution. There is no apparent clear continuity between the two regions in the Plane. However, examining the cloud distribution below and above the Plane is very informative. As shown in Figs. 10 and 11, this second rank of structures is particularly prominent well above the solar plane and is oriented in a homogeneous way in the Z = +100 pc and even Z = +150 pc images, but its orientation is drastically different from that of the large structure close to the Plane, here along a −90, +90° axis. We labeled it upper Sagittarius–Carina in Fig. 14. The closest cloud complex is not visible above the solar plane, but it is clearly visible below the Plane down to −150 pc. We labeled it lower Sagittarius–Carina. Another way to look at this complex structure of dust clouds in Carina–Sagittarius is by inspection of the vertical planes of Fig. 13. The fourth quadrant corresponds to the left parts of the images. Especially along l = 270 and l = 230° (left parts of first and seventh images) lies a complex wavy structure throughout the solar plane that is the counterpart of what is seen in horizontal planes.

The second striking and prominent feature is the more than 2 kpc long series of dense dust clouds towards l ~ 45°. This structure is close to but not exactly radial. We nicknamed it the split. It is located between the Local arm and the Carina–Sagittarius arm. This structure is visible at various degrees in the previous maps of Marshall et al. (2006); Chen et al. (2019); Sale et al. (2014); Kos et al. (2014); Green et al. (2015, 2018), but has never been commented on. Interestingly, it has characteristics that differentiate it from other main regions. Figures 14–16 show almost no match between the split and HII regions, masers, or O-star concentrations. In contrast, Fig. 17 reveals a very large concentration of molecular clouds all along the split. Interestingly, Skowron et al. (2018) also found a strong concentration of Cepheids right along the split that they attributed to a ~30-Myr-old star formation burst.

A third at least 2 kpc long series of compact regions starts close to the Sun and is oriented along l = 75–90°. It corresponds to the well-known Cygnus Rift and can be seen as the main part of the Local arm. Here again, the orientation of the structure varies with distance Z to the mid-plane, as shown in Figs. 10 and 11. Below the mid-plane, the orientation is toward l = ~55°, while at Z ≥ 50 pc, it rotates and reaches l = ~90 °. Interestingly, only the ≥ 1.2 kpc part of this structure is rich in O stars, HII regions, masers, and molecular clouds, while the closer part shows none of these tracers.

A shorter structure of about 0.8 kpc length is also visible in the first quadrant. This structure is anchored in the inner part of the split and lies roughly parallel to the Cygnus Rift, it is labeled Vul in Fig 14. Similarly to the split, it is rich in MCs, but tracers of star formation are absent. Farther away at d ≥1.5 kpc and toward l ~ 60° are the Vul OB4 masers and another elongated structure that could be the foreground counterpart of the local spur of Xu et al. (2016).

In the second and third quadrants lie series of clouds that are approximately aligned along the rotation direction and are located at a distance of about 2 kpc from the l = 90-270 axis. These complexes correspond to the expected location of the Perseus arm. They are distributed in a much less compact way than the structures described previously, suggesting that the outer Perseus is a more diffuse arm. Surprisingly, only very few dust clouds are visible between l = 120° and l = 150°, where a prolongation of Perseus and clouds at about 2 kpc might be expected. The Perseus arm is indeed indicated in this region by several HII regions, O stars, and masers. A potential explanation might be a bias of our mapping caused by the scarcity of targets in this area, which is due to the high opacity of clouds in the foreground. However, the threshold white line in Fig. 9 shows that other regions with the same foreground opacity have well-defined dust clouds. Moreover, no detected MCs lie in this area, suggesting that both dust and molecular gas have been processed as a result of strong star formation.

The inner part of the third quadrant is characterized by the giant cavity centered on l = 240°, which has been clearly defined in our previous 3D maps. This cavity is surrounded by the main star-forming regions of Orion and Vela. It is adjacent and connected to the Orion-Eridanus cavity below the Plane. The latter is better visible in the third vertical plane in Fig 12. At all distances from the Plane, this region of the third quadrant extends between 150 pc to 1.2 kpc. We note that this region is devoid of dust, but is filled with ionizedgas (Lallement 2015), suggesting that most of the dust evaporated during the event that ionized the gas. As we have noted in previous maps, this cavity is aligned with a narrower cavity elongated along l = 60°, that is, a cavity separating the Cygnus Rift and the split. This axis also corresponds to an ionization gradient in the nearby ISM (Wolff et al. 1999).

Fig. 10 Differential extinction in horizontal planes 300, 250, 200, 150, 100, and 50 pc below the Galactic plane. The color coding and the significance of the black line are the same as in Fig. 9.

Fig. 11 Same as Fig. 10 for planes above the Galactic plane. We note that the orientation of the Carina–Sagittarius region in the fourth quadrant, especially at Z = +100, +150 pc, differs by a large angle (~ 45°) from the orientation of the compact closer part shown in Fig. 9.

Fig. 12 Dust distribution in vertical planes containing the Sun and for longitudes spaced every 10°. The quantity represented is the square root of the differential extinction. The color scale is indicated in Fig. 9. The dotted gray line represents the distance beyond which the final resolution of 25 pc cannot be achieved because there are too few targets (see text). The dashed light gray contours correspond to a differential extinction of 0.003 mag.pc−1. Note the wavy structure throughout the Plane for longitudes 200, 220, and 250°.

Fig. 13 Continuation of Fig 12. Note the wavy structure throughout the Plane, especially for longitudes 270 and 330°.

5 Discussion and perspectives

We have used the Gaia and 2MASS photometric data to estimate the extinction toward ~27 million carefully selected target stars with a Gaia DR2 parallax uncertainty below 20%. Approximately 16 million distance-extinction pairs for the closest objects were inverted to produce 3D maps of the dust density within 400 pc from the Plane and in a 6 kpc by 6 kpc region along the Plane. To do so, a new hierarchical inversion algorithm was developed that has theadvantage of adapting the map resolution to the target space density. Following this procedure, the minimum size of the inverted structures varies from 25 pc (a resolution achieved within 1 kpc) to 500 pc (in a few regions thatare distant by more than 3 kpc from the Sun). Tools adapted to the new 3D map will be available online³. They allow drawing an image of the dust distribution in any chosen plane, regardless of whether it contains the Sun, as well as reddening profiles. Images and cumulative extinctions as well as the 3D distribution itself will be available for download. We emphasize that these full 3D inversions have the advantage of showing the relation of cloud complexes to each other in 3D space. As we have cautioned in previous similar works, the absolute values of the integrated extinction between the Sun and each point in the computation volume may underestimate the extinction when sightlines cross cloud cores. This loss of angular resolution is due to the regularization constraints and is the price to pay for using full 3D correlations between neighboring points. Moreover, because of the general bias toward less obscured, hence brighter stars, the differential extinction is likely to be biased low, especially at large distances.

The dust density distribution appears to be dominated by elongated structures with a double system of orientations. The first quadrant is dominated by two elongated structures. One of the two is coincident with the well-known Cygnus Rift, which usually is identified as the Local arm. Its orientation varies from close to the direction of rotation l = 80–90. ° above the mid-plane to about l = 60. ° below the mid-plane. The second structure, here nicknamed the split, is oriented very differently from the direction of rotation and is visible as a long arm segment linking the Perseus and Sagittarius–Carina spiral arms (see, e.g., Hou & Han 2014; Vallée 2017). The latter structure is different from the spur described by Xu et al. (2016), which is mostly farther away from our map and may be visible as the prolongation of a third group of clouds starting beyond 1.5 kpc (see Fig. 14). The Local arm and the split are both strong in extinction, but the Local arm is better traced with masers and star formation regions, especially at large distance (Reid et al. 2014), while the split is stronger in CO (Miville-Deschênes et al. 2016). According to the Cepheid study of Skowron et al. (2018), there seems to be an age difference between these structures and the split might be younger.

The new map reveals a particularly compact and well-delineated foreground part of Sagittarius-Carina that extends in the fourth quadrant and at 0 < l < 30° in the first quadrant. This structure, labeled inner in Fig. 14, is mainly visible within and below the Plane and has the same orientation as the split. A second and similarly compact outer part of Sagittarius–Carina is located at larger distance and predominantly 50–150 pc above the Plane. At variance with the inner part, it is oriented in the direction of rotation.

A large portion of the Perseus arm is perceptible in the second and third quadrants, but it appears to be significantly more fragmented than the other arm segments. In the second quadrant it is not detected above 120 < l < 145°, in agreement with the small number of molecular clouds in this area, but in contrast with HII regions and masers. Whether the absence of dust and CO is related to past star formation or if there is a bias in dust detection deserves further study.

The distribution of the dust in directions perpendicular to the Plane varies strongly from one region to the other, and there are striking structures that in some regions have a wavy pattern with periods on the order of 1–2 kpc, as revealed by images in vertical planes (Figs. 12 and 13). The largest amplitude of off-plane excursions is found in the vicinity of the Sun, within about 500 pc. It can hardly be an observing bias related to obscuring foreground clouds. However, additional data are required to preclude biases linked to the scarcity of targets at large distance from the mid-plane. Interestingly, stellar velocity patterns are also particularly complex within 500 pc (see, e.g., Kawata et al. 2018, Antoja et al. 2018, Gaia Collaboration 2018c). Finally, many more cloud structures and cavities devoid of dust are revealed by the 3D inversion, but describing them is beyond the scope of this paper.

In the future, improved parallax measurements and more numerous extinction estimates are expected, allowing inversions at better spatial resolution and the extension of the maps to larger distances. For these maps, bias estimates in parallax distance at the individual level will need to be included, and this is part of our current efforts. The models of Galactic evolution, which are currently being revolutionized by Gaia, will certainly benefit from detailed comparisons between the distributions and motions of stars and interstellar matter. 3D dust maps such as our inverted maps can be useful directly or indirectly as ingredients for new kinetic tomography techniques for the gaseous ISM (Tchernyshyov & Peek 2017) that use 3D cloud distributions and CO/HI radio spectral cubes and/or absorption data to assign velocities to interstellar clouds. This in turnwill allow detailed comparisons between the stellar and ISM motions and more constraints on star formation.

Fig. 14 Projections onto the Plane of the masers from Reid et al. (2014), Xu et al. (2016), and Honma et al. (2012) that possess parallax distances and are located within 50 pc from the Plane, superimposed on the dust distribution throughout the Galactic plane. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The dotted gray line represents the distance beyond which the final resolution of 25 pc cannot be achieved because targets are scarce (see text). We show the spiral arms derived by Reid et al. (2014) (dashed lines, inner and outer parts, together with the local spur of Xu et al. 2016) and Hou & Han (2014) (solid lines). In black we plot Perseus, in blue the Local arm, in magenta we show Sagittarius, and in cyan we plot Scutum.

Fig. 15 Same as Fig. 14 for star-forming regions from the compilation of Hou & Han (2014) and from Russeil (2003). Only objects with stellar distances (see text) and |z| < 50 pc are shown, and the size of the markers decreases for increasing |z|.

Fig. 16 Same as Fig. 14 for O stars from the recent work of Xu et al. (2018). Distances are from Gaia DR2. Only stars within 50 pc from the Plane are shown.

Fig. 17 Same as Fig. 14 for the CO clouds from Miville-Deschênes et al. (2016) (circles) and the nearby giant molecular clouds listed by Hou & Han (2014) (stars). The black-and-white scale indicates the CO equivalent width. Only objects with |z| < 50 pc are shown, and the size of the markers decreases for increasing |z|.

Acknowledgements

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 also 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 was supported by the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Reddening profiles and comparisons with profiles from Pan-STARRS 1 and 2MASS

Green et al. (2018) have made available reddening profiles derived from Pan-STARRS 1 (PS1) and 2MASS for directions spaced by about 10 arcmn in the whole sky region accessible to Pan-STARRS. It is particularly interesting to compare 3D maps based on different photometric systems and different techniques. In order to perform comparisons between our results, which are based on Gaia G, G_BP, and G_RP, and 2MASS, and the Green et al. (2018) reddening profiles, we have integrated our computed differential reddening in a series of directions. We uniformly converted the color excess E(B–V) from the Green et al. (2018) map into A₀ by multiplying by a factor 0.88 × 3.1, that is, assuming R = 3.1 and following the scaling factors recommended by the authors. For all the comparisons we used their maximum probability density estimate (the best-fit) distance-reddening curve. Figure A.1 displays comparisons within the Plane every 5° in the longitude interval l = 10–230°. In addition to the Green et al. (2018) central profile at the exact coordinates (l,b), we added their profiles for longitudes l ± 0.25° and b ± 0.25° as well as the average between the five profiles. Our maximum resolution is also indicated as a function of distance on the right scale. Both the adjacent directions and the resolution allow a better interpretation of the comparisons. The figure shows that at distances where our achieved resolution is 25 or 50 pc, our reddening profile compares well with either the central PS1 profile or the average profile. Marked differences appear where the resolution we achieved is poorer (from 100 pc) and in very complex regions with very many intervening clouds, for example, between longitudes l = 75° and 85° beyond 500–1000 pc. In these cases, clouds are similarly located, but our extinction is underestimated. On the other hand, the good agreement observed at short distances where the PS1 flag starts to be OFF suggests that our measurements within several hundred parsec continue PS1 profiles well.

Fig. A.1 Extinction profiles integrated through the 3D distribution of differential extinction for sightlines within the Plane (black) and reddening profiles from Green et al. (2018) scaled to A0 (see text, red solid line). We also show the Green et al. (2018) profiles for four surrounding directions at ±0.25° in longitude and latitude (dotted red line) and their average in non-flagged and all intervals (thick and thin dashed red line, respectively). Our achieved resolution (see Sect. 3) is also shown (pale blue line, right scale. Coordinates are indicated in the graph.

References

All Figures

	Fig. 1 Example of an intrinsic color-color relation fit (top) and its residuals (bottom) for G − GBP vs. G − K. In red we plot the outliers that were removed in the process.
In the text

	Fig. 2 A0 extinction distribution (mag) in low-extinction selections. In blue we show the local bubble selection (ϖ > 14 mas, corresponding to the first 70 pc), and in green we plot a high-latitude selection (|b| > 80°). The median residuals are smaller than 0.03 mag with a standard deviation of 0.2 mag.
In the text

	Fig. 3 Dereddened Hess diagram for a subset of stars with a relative parallax uncertainty smaller than 10%. The grayscale corresponds to the square root of the stellar density. The Red Clump is the most prominent feature. Its shape can be compared to Fig. 2 of Ruiz-Dern et al. (2018).
In the text

	Fig. 4 Comparison between the extinction A0 derived here and AV determined by Santiago et al. (2016) using APOGEE data.
In the text

Fig. 5 Left panel: 250 × 250 × 250 pc bins similar in volume to those used in the first step of the inversion. Shown are the cells centered on the Plane. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The colorcode corresponds to the number of target stars per bin, i.e., a quantity proportional to the targetspace density. Only bins with more than ten targets are used at this stage. Middle panel: 12.5 × 12.5 × 12.5 pc bins used in the last step of the inversion. Only bins with more than two targets are used. Right panel: resolution achieved in the Plane based on the hierarchical inversion. In regions of high target space density, essentially close to the Sun, the resolution is 25 pc and is achieved after the last step. In regions of scarce targets, mainly at large distances, the resolution is coarse, may reach 500 pc, and is achieved at the first steps. These regions are not improved during the following steps.

In the text

	Fig. 6 Inverted 3D distribution of differential extinction in the Galactic plane and as it is obtained after each of the first four iterative steps. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The achieved resolution is300, 200, 100, and 50 pc, respectively.
In the text

	Fig. 7 Biases on the parallax distance in parsec from a Monte Carlo simulation adapted to our dataset.
In the text

	Fig. 8 Top: Planck map of the Galactic dust optical thickness at τ353 GHz. Middle: estimated extinction toward all targets located beyond 1 kpc. Bottom: 3D differential extinction integrated from the Sun to the limits of the computational volume.
In the text

Fig. 9 Dust density along the Plane, based on the 3D distribution obtained after the last inversion step. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The color scale represents the differential extinction in units of mag per parsec. The dashed black line represents the distance beyond which the final resolution of 25 pc is not achievable due to target scarcity (see text). The dotted white contours correspond to a differential extinction of 0.003 mag.pc−1 and delimit the dense areas.

In the text

	Fig. 10 Differential extinction in horizontal planes 300, 250, 200, 150, 100, and 50 pc below the Galactic plane. The color coding and the significance of the black line are the same as in Fig. 9.
In the text

	Fig. 11 Same as Fig. 10 for planes above the Galactic plane. We note that the orientation of the Carina–Sagittarius region in the fourth quadrant, especially at Z = +100, +150 pc, differs by a large angle (~ 45°) from the orientation of the compact closer part shown in Fig. 9.
In the text

Fig. 12 Dust distribution in vertical planes containing the Sun and for longitudes spaced every 10°. The quantity represented is the square root of the differential extinction. The color scale is indicated in Fig. 9. The dotted gray line represents the distance beyond which the final resolution of 25 pc cannot be achieved because there are too few targets (see text). The dashed light gray contours correspond to a differential extinction of 0.003 mag.pc−1. Note the wavy structure throughout the Plane for longitudes 200, 220, and 250°.

In the text

	Fig. 13 Continuation of Fig 12. Note the wavy structure throughout the Plane, especially for longitudes 270 and 330°.
In the text

Fig. 14 Projections onto the Plane of the masers from Reid et al. (2014), Xu et al. (2016), and Honma et al. (2012) that possess parallax distances and are located within 50 pc from the Plane, superimposed on the dust distribution throughout the Galactic plane. The Sun is at the center of the figure at (0,0). The Galactic center is to the right. The dotted gray line represents the distance beyond which the final resolution of 25 pc cannot be achieved because targets are scarce (see text). We show the spiral arms derived by Reid et al. (2014) (dashed lines, inner and outer parts, together with the local spur of Xu et al. 2016) and Hou & Han (2014) (solid lines). In black we plot Perseus, in blue the Local arm, in magenta we show Sagittarius, and in cyan we plot Scutum.

In the text

	Fig. 15 Same as Fig. 14 for star-forming regions from the compilation of Hou & Han (2014) and from Russeil (2003). Only objects with stellar distances (see text) and |z| < 50 pc are shown, and the size of the markers decreases for increasing |z|.
In the text

	Fig. 16 Same as Fig. 14 for O stars from the recent work of Xu et al. (2018). Distances are from Gaia DR2. Only stars within 50 pc from the Plane are shown.
In the text

	Fig. 17 Same as Fig. 14 for the CO clouds from Miville-Deschênes et al. (2016) (circles) and the nearby giant molecular clouds listed by Hou & Han (2014) (stars). The black-and-white scale indicates the CO equivalent width. Only objects with |z| < 50 pc are shown, and the size of the markers decreases for increasing |z|.
In the text

Fig. A.1 Extinction profiles integrated through the 3D distribution of differential extinction for sightlines within the Plane (black) and reddening profiles from Green et al. (2018) scaled to A0 (see text, red solid line). We also show the Green et al. (2018) profiles for four surrounding directions at ±0.25° in longitude and latitude (dotted red line) and their average in non-flagged and all intervals (thick and thin dashed red line, respectively). Our achieved resolution (see Sect. 3) is also shown (pale blue line, right scale. Coordinates are indicated in the graph.

In the text