© R. H. Leike and T. A. Enßlin 2019

1 Introduction

Emission and extinction by Galactic dust is a prominent astronomical foreground at many wavelengths. Therefore, knowing its distribution on the 2D sky and in 3D is essential for many astronomical observations. However, dust is also interesting to be studied on its own, as it provides information on the physical conditions in the interstellar medium and informs us about star forming regions. Dust has been mapped out by surveys for a long time, the first notable contribution being Burstein & Heiles (1978). The dust reconstruction by these latter authors, as in most reconstructions of the dust distribution so far, was focused on mapping the dust in 2D on the sky. This can be done by looking at the sky in infrared (IR) wavelengths, where it is dominated by dust emission. However, when mapping out dust using IR emission one is biased by the radiation field, as dust emission is not only proportional to dust density but also to the amount of starlight absorbed by the dust. Furthermore, dust maps that were produced by mapping IR light might contain extended IR sources that are not within our galaxy, as was shown by Chiang & Ménard (2019). On the other hand, a hypothetical cold dust cloud cannot be seen in IR, leading to systematic errors in the analysis of distant targets, for example quasars or the cosmic microwave background (CMB).

For accurate analyses of objects in our Galactic vicinity, it is vital to have a 3D dust map as a foreground model, which informs us about regions that cannot be observed, or can only be observed with less fidelity, due to dust obscuration. The first nonparametric reconstruction of Galactic dust in 3D was published by Arenou et al. (1992). Since then there have been many attempts to chart the dust density in 3D in increasing resolution and accuracy and for an ever greater part of our galaxy (Vergely et al. 1997, 2001; Chen et al. 2018; Rezaei Kh et al. 2017, 2018a,b; Gontcharov 2012, 2017; Sale et al. 2014; Sale & Magorrian 2018). A large driving force for 3D dust reconstruction and astronomy in general comes from large surveys like 2MASS (Skrutskie et al. 2006), Pan-STARRS (Kaiser et al. 2002) and SDSS/APOGEE (Albareti et al. 2017) and WISE (Wright et al. 2010). These surveys provide photometric measurements and some of them provide spectra for thousands of stars, from which the calculation of photometric distances is possible. There are two 3D dust reconstructions based on these data sets that are closely related to the approach taken in this paper: in Lallement et al. (2018) reddening data from 71 000 sources was used to perform a 3D reconstruction using Gaussian process regression. The resulting dust map covers a 4 kpc square of the Galactic plane and 600 pc in perpendicular direction with a voxel size of (5 pc)³. In Green et al. (2018) a 3D dust map was produced by combining the star data of Pan-STARRS and 2MASS, binning it in angular and distance bins, and performing independent Bayesian reconstructions per angular bin. The result is a dust map that covers three quarters of the sky to a distance up to 2 kpc. This reconstruction shows artificial radial structures called the fingers of God effect, in analogy to the well-known phenomenon in cosmology (Jackson 1972; Tully & Fisher 1978). One way to mitigate this effect is to use higher-accuracy parallax information.

A prominent new survey is performed by the Gaia collaboration (Brown et al. 2016). In its second data release (DR2, Gaia Collaboration 2018) accurate parallaxes for roughly 2 billion stars are published. The provided catalog also contains estimates of extinction coefficients for a subset of about 88 million stars, using spectral information of the three energy bands of the Gaia satellite. Due to the limited spectral information, the accuracy of the extinction coefficients estimated for individual sightlines is quite low. For this reason it is recommended by Andrae et al. (2018) to not use the information of individual sightlines but only the joint information of several sightlines. Even though the data quality of individual sightlines is rather low, the sheer amount of data points and the accuracy of the parallaxes outweigh this limitation, as our work shows.

Three-dimensional dust reconstructions have never been performed using Gaia data only. Instead, Gaia data have only been used for accurate parallax measurements and the more accurate spectral information of other surveys was used (Lallement et al. 2018, for example).

In this paper we present a 3D dust reconstruction using Gaia DR2 data only. The results of the reconstruction are provided at the CDS and online¹, and can be used under the terms of the ODC-By 1.0 license. The inference of the unknown dust density from the extinction data is performed by a critical filter, a method for Gaussian process regression with nonparametric kernel learning, first published in Enßlin & Frommert (2011). While the statistical model used here is, apart from minor details, equivalent to the model introduced in that paper, the algorithm to arrive at approximate posterior summary statistics is quite different. The relevant numerical method used in this paper is outlined in Knollmüller & Enßlin (2019), which describes a general method to derive posterior summary statistics for high-dimensional Bayesian inference problems. For a theoretical discussion of the underlying inference framework of information field theory we refer to Enßlin (2019). The algorithm was implemented using the Python package NIFTy5, which is the newest version of the software package NIFTy (Selig 2013; Steininger et al. 2017; Arras et al. 2019a)². Even though mathematical theory, statistical motivations, and numerical details are distributed over the aforementioned papers, this paper is entirely self-contained by describing the whole method.

In Sect. 2 we discuss which part of the Gaia data we used. We introduce our statistical model of the interstellar dust density as well as of the measurement in Sect. 3. In Sect. 4 we present a test application of the algorithm using synthetic data, verifying the predictive power of the algorithm. The main results of the dust density reconstruction using Gaia data are presented in Sect. 5. This section also contains a brief recommendation on how to use our results. Our dust reconstruction is compared to other 3D dust density reconstructions in Sect. 6 and in Sect. 7 we summarize the findings of this paper.

2 Data

We used the data from the Gaia DR2 catalog by Gaia Collaboration (2018) to reconstruct the Galactic dust in the nearby interstellar medium. From the Gaia data archive we extract the parallaxes, the G-band extinction, the latitude and longitude, as well as their respective uncertainties. A plot of the full Gaia extinction data set can be seen in Fig. 1.

We select sources according to the following criteria:

• 1.

  the above mentioned data are available for the source;

• 2.

  the parallax $\tilde{\omega }$ is inside a 600 pc cube around the Sun;

• 3.

  the relative parallax error is sufficiently low, $\tilde{\omega }/{\sigma }_{\tilde{\omega }}>5;$

• 4.

  Priam flag 0100001 or 0100002.

The last two criteria are suggested by Andrae et al. (2018). There are about 3.7 million stars selected by these criteria. Figure 1 shows a sky average of the data points used in the reconstruction. In this data plot one can observe structures present also in other dust maps, for example the Planck dust map (Fig. 1).

Fig. 1 Top panel: natural logarithm of the thermal dust emission map produced by Planck Collaboration I (2019). The color-scheme was saturated to visually match that of the Gaia dust extinction data shown in the middle panel. Bottom panel: subset of Gaia extinction data within a 600 pc cube centered on the Sun; the data used in this paper. The scale is a natural logarithm of the extinction data in magnitudes. Data points in the same direction were averaged.

3 Model

3.1 Algorithm

The algorithm is derived from Bayesian reasoning. In Bayesian reasoning, information that some data d provides about a quantityof interest s is calculatedaccording to Bayes theorem: $$P(s|d)=\frac{P(d|s)P(s)}{P(d)}.$$(1)

We note that the quantity of interest can be a (possibly high-dimensional) vector; in our case it is the dust density for every point in space (256³ degrees of freedom after discretization). There are three main ingredients necessary for the inference of the quantity of interest s:

• 1.

  The likelihood P(d|s) of the data d given a realization of the quantity of interest s. We describe our likelihood in Sect. 3.2.

• 2.

  The prior P(s) describing the best available knowledge about the quantity of interest s in absence of data. We describe our prior in Sect. 3.3.

• 3.

  An inference algorithm that yields a statistical summary of P(s|d) given the joint distribution P(d, s) = P(d|s)P(s) of d and s. We use the inference algorithm described in Knollmüller & Enßlin (2019).

The main quantity of interest s is the logarithmic G-band dust extinction cross-section density s = ln(αρ∕pc), henceforth referred to as the logarithmic dust density. Hereby, ρ denotes the actual dust mass density and α the average G-band dust cross section per mass. The value of α is uncertain,which is why we report extinction densities, also referred to as dust pseudo-densities.

We approximate the posterior with a Gaussian $$\begin{array}{ll}Q(s)\hfill & =\mathcal{G}(s-m,D)\hfill \\ \hfill & =\frac{\text{exp}\left(-\frac{1}{2}{(s-m)}^{†}{D}^{-1}(s-m)\right)}{{\left|2\pi D\right|}^{\frac{1}{2}}},\hfill \end{array}$$

by adopting a suitable mean m and uncertainty dispersion D. The approximation is obtained by minimizing the Kullback–Leibler divergence (Kullback & Leibler 1951), $$\text{KL}(Q,P)={\displaystyle \int \text{d}}Q\text{\hspace{0.17em}ln}\frac{Q}{P},$$(4)

with respect to the parameters of Q. This approach is known as variational Bayes (Nasrabadi 2007) or Gibbs free energy approach (Enßlin & Weig 2010). The approach we take in finding the unknown approximate posterior mean m and covariance D of Eq. (3) is described in detail in Knollmüller & Enßlin (2019). It can be summarized as follows:

• 1.

  Calculate the negative log-probability −log(P(s, d)) for the problem, disregarding normalization terms like the evidence P(d).

• 2.

  Perform coordinate transformations of the unknown quantities until those are a priori Gaussian distributed with unit covariance (Kucukelbir et al. 2017; Knollmüller & Enßlin 2018).

• 3.

  Choose the class of approximating distributions to be Gaussian with variable mean m and covariance $D={(\mathbb{1}+M_m)}^{-1}$, where M_m is the Fisher information metric at the current m. Here 1 is the contribution of the prior which was transformed in step 2 to have unit covariance. This uncertainty dispersion is a lower bound to the uncertainty (Cramér 1946; Rao 1947) and has been shown to be an efficient technique to take cross-correlations between all degrees of freedom into account without having to parameterize them explicitly (Knollmüller & Enßlin 2019).

• 4.

  Minimize Eq. (4) with respect to m using the Newton Conjugate gradient as a second-order scheme with the covariance D of Q as curvature. The expectation value with respect to Q is hereby approximated through a set of samples drawn from the approximating distribution Q. Second-order minimization by preconditioning with the inverse Fisher metric is also called natural gradient descent (Amari 1997) in the literature.

3.2 Likelihood

The likelihood P(d|s) can be split into two parts: the first states how the true extinction depends on the dust density and the second states how the actual data is distributed given the true extinction on that line of sight. The first part, which we call the response R, states how the unknown dust extinction density ρ imprints itself on the data. The extinction of light on the ith line of sight L_i is given by the line integral $${(A_{\text{G}})}_i={\left[R(\rho )\right]}_i={\displaystyle {\int }_{L_i}\text{d}}l\alpha \rho (l).$$(5)

Here α is the average dust cross section per unit of mass and the line of sight L_i = L_i(ω) is dependent on the true parallax ω. As noted in Sect. 3.1, the value of α is uncertain and we reconstruct the dust extinction density $s=\text{\hspace{0.17em}ln}\left(\alpha \rho /\text{pc}\right)$ instead.

The extinction is additive because the extinction data are given in the magnitudes scale, which is logarithmic. The true parallax ω of the star is uncertain. We assume the true parallax ω to be Gaussian distributed around the published parallax $\tilde{\omega }$ with a standard deviation equal to the published parallax error σ_ω: $$P(\omega |\tilde{\omega },{\sigma }_{\omega })=\mathcal{G}(\omega -\tilde{\omega },{\sigma }_{\omega }^2).$$(6)

The parallaxes of Gaia DR2 were shown to be Gaussian distributed with incredible reliability by Luri et al. (2018). However, it was also noted by the same authors that there can be outliers. By restricting ourselves to close-by sources for which G-band extinction values are published, we expect to have cut out most of the outliers.

We do not reconstruct the actual positions of the stars in our reconstruction, and therefore we have to marginalize them out to obtain the response: $$\begin{array}{l}P({(A_{\text{G}})}_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_1},\rho )=\\ ={\displaystyle \int \text{d}}{\omega }_{i}P({(A_{\text{G}})}_i,{\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i},\rho )\\ ={\displaystyle \int \text{d}}{\omega }_{i}P({(A_{\text{G}})}_i|{\omega }_i,{\tilde{\omega }}_i,{\sigma }_{{\omega }_i},\rho )P({\omega }_i|\widetilde{{\omega }_i},{\sigma }_{{\omega }_i},\rho )\\ ={\displaystyle \int \text{d}}{\omega }_{i}P({(A_{\text{G}})}_i|{\omega }_i,\rho )P({\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i})\\ ={\displaystyle \int \text{d}}{\omega }_i\delta \left({(A_{\text{G}})}_i-{\displaystyle {\int }_{L_i({\omega }_i)}\text{d}}{l}_i\alpha \rho (l_i)\right)P({\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i})\\ \approx \delta \left({(A_{\text{G}})}_i-{\displaystyle \int \text{d}}{\omega }_i{\displaystyle {\int }_{L_i({\omega }_i)}\text{d}}{l}_i\alpha \rho (l_i)P({\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i})\right)\\ =\delta \left({(A_{\text{G}})}_i-{\displaystyle \int \text{d}}{\omega }_i{\displaystyle {\int }_{L_i(0)}\text{d}}{l}_i\alpha {\mathbb{1}}_{[0,\frac{1}{\omega }]}(l_i)\rho (l_i)\mathcal{G}({\omega }_i-{\tilde{\omega }}_i,{\sigma }_{{\omega }_i}^2)\right)\\ =\delta \left({(A_{\text{G}})}_i-{\displaystyle {\int }_{L_i(0)}\text{d}}{l}_i\alpha \rho (l_i){\text{ sf}}_{\mathcal{G}}\left(\frac{\frac{1}{l_i}-\widetilde{{\omega }_i}}{{\sigma }_{{\omega }_i}}\right)\right).\end{array}$$(7)

Here $${\text{sf}}_{\mathcal{G}}(x)=1-{\displaystyle {\int }_{-\infty }^x\text{d}}t\frac{1}{\sqrt{2\pi }}\text{exp}\left(\text{-}\frac{\text{1}}{\text{2}}{\text{t}}^{\text{2}}\right)$$(8)

denotes the survival function of a standard normal distribution and $${\mathbb{1}}_{[a,b]}(x)=\{\begin{array}{ll}1\hfill & \text{ for }x\in [a,b]\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}$$(9)

denotes the indicator function for the closed interval [a, b]. We note that we made an approximation where we replace the true extinction with the expected extinction. As a consequence the lines of sight are smeared out by the parallax uncertainty in our approximation. This smoothing can be regarded as a first-order correction for the uncertainty of the parallax and was already used by Vergely et al. (2001). A fully Bayesian analysis would treat the true parallax as unknown and infer these along the other unknowns, but this is beyond the scope of this work.

For the algorithm, the integral in Eq. (7) is discretized into a weighted sum, such that each voxel contributes to the line integral over L_i with weight equal to the length of the line segment of L_i within that voxel. That weight is further discounted by the probability $P(l|\tilde{\omega },{\sigma }_{\omega })$ of that voxel being on the line of sight. Applying the response R thus takes $\mathcal{O}(N_{\text{data}}{N}_{\text{side}})$ operations, where N_data = 3 661 286 is the number of data points used in the reconstruction and N_side = 256 is the number of voxels per axis. Due to the large number of data points, evaluating the response on a computer turns out to be numerically expensive. The inference algorithm (see Sect. 3.1) is a minimization for which this response has to be evaluated many times. To make the dust inference feasible in a reasonable amount of time we restricted our reconstruction to a 600 pc cube centered on the Sun.

The second part of the likelihood states how the published data $\widetilde{A_{\text{G}}}$ is distributed given the true extinctions ${(A_{\text{G}})}_i$. We use the data likelihood recommended by the Gaia collaboration in Andrae et al. (2018). This likelihood assumes the data $\widetilde{A_{\text{G}}}$ to be distributed according to a truncated Gaussian with a global variance $N={\left(0.46\text{\hspace{0.17em}mag}\right)}^2\mathbb{1}$. This leads to the likelihood $$\begin{array}{l}P(\widetilde{A_{\text{G}}}|\rho )={\displaystyle \prod _i\frac{\mathcal{G}({\widetilde{A_{\text{G}}}}_i-R{(\rho )}_i,N_{ii})}{{\text{cdf}}_{\mathcal{G}(R{(\rho )}_i,N_{ii})}(A_{\text{G}}^{\text{max}})-{\text{cdf}}_{\mathcal{G}(R{(\rho )}_i,N_{ii})}(A_{\text{G}}^{\text{min}})}}\\ \text{for }d\in [A_{\text{G}}^{\text{min}},A_{\text{G}}^{\text{max}}].\end{array}$$

Here ${\text{cdf}}_{\mathcal{G}(R{(\rho )}_i,N_{ii})}$ denotes the cumulative density function of a normal distribution with mean $R{(s)}_i$ and variance N_ii. We took the boundaries of the truncated Gaussian to be $A_{\text{G}}^{\text{min}}=0$ and $A_{\text{G}}^{\text{max}}=3.609\text{\hspace{0.17em}mag}$ as recommended in Andrae et al. (2018).

3.3 Prior

We assume the dust density to be a positive quantity that can vary over orders of magnitude. The dust is assumed to be spatially correlated, statistically homogeneous, and isotropic. The statistical model is constructed to be as general as possible with these two properties in mind.

To reflect the positivity and to allow variations of the dust density by orders of magnitude we assume the dust density ρ to be a priori of log-normal distribution with $$\begin{array}{l}\alpha \rho ={\rho }_0\text{\hspace{0.17em}exp}(\text{s}),\\ \text{where }s\curvearrowleft \mathcal{G}(s,S)\end{array}$$

is assumed to be Gaussian distributed with Gaussian process kernel S. Here, ${\rho }_0=\frac{1}{1000}$ pc is a constant introduced to give ρ the correct unit and to bring it to roughly the right order of magnitude. By using an exponentiated Gaussian process we allow the dust density to vary by orders of magnitude while simultaneously ensuring that it is a positive quantity. In Eq. (13), S is the prior covariance. If we assume no point or direction to be special a priori, then according to the Wiener-Khinchin theorem S can be fully characterized by its spatial power spectrum $S_{k{k}^{\prime }}=2\pi \delta (k-k^{\prime })P_s(k)$. We nonparametrically infer this power spectrum P_s(k) as well. There are two main motivations to reconstruct the power spectrum. From a physical perspective the power spectrum provides valuable insight into the underlying processes. From a signal-processing point of view, many linear filters can be identified with a Bayesian filter that assumes a certain prior power spectrum. The optimal linear filter is obtained when the power spectrum used for the filter is exactly equal to the power spectrum of the unknown quantity (Enßlin & Frommert 2011). However, the power spectrum of the unknown quantity is usually also unknown, and therefore this must be reconstructed well. While this argumentation holds for linear filters, certainly many aspects of it carry over to nonlinear filters such as the reconstruction performed in this paper.

Figure 2 depicts the hierarchical Bayesian model for the extinction data $\widetilde{A_{\text{G}}}$ resulting from the logarithmic dust density ln(ρ), which itself is shaped by the power spectrum P_s(k).

Our statistical model for the power spectrum P_s(k) is a falling power law with Gaussian distributed slope and offset modified by differentiable nonparametric deviations. It is, apart from minor details³, an integrated Wiener process (Doob 1953) on log–log-scale.

This is realized by the following formula, $$\begin{array}{l}\sqrt{P_s(k)}=\\ \text{exp}\left(({\varphi }_{\text{m}}{\sigma }_{\text{m}}+\overline{m})\text{log}(\text{k})\text{+}{\varphi }_{\text{y}}{\sigma }_{\text{y}}+\overline{y}+{\mathbb{F}}_{\text{log}(\text{k})\text{t}}^{\text{-1}}\left(\frac{\text{a}}{\text{1+}\frac{{\text{t}}^{\text{2}}}{{\text{t}}_{\text{0}}^{\text{2}}}}\tau (\text{t})\right)\right),\end{array}$$(14)

where ϕ_m, ϕ_y, and τ(t) are the parameters to be reconstructed, ${\sigma }_m=1,\overline{m}=-4,{\sigma }_y=2,\overline{y}=-16,{\sigma }_y=3,a=11,t_0=0.2$ are fixed hyperparameters, ${\mathbb{F}}_{\text{log}(\text{k})\text{t}}^{-1}$ denotes the inverse Fourier transform on log-scale, and $V={(600\text{\hspace{0.17em}pc})}^3$ is the total volume of the reconstruction. These hyperparameters settings were determined by trial and error such that data measured from a prior sample have roughly the same order of magnitude as the actual data and such that the dust density varies by more than one order of magnitude in prior samples.

In our reconstruction, the parameter τ for the smooth deviations of the log–log power spectrum was discretized using 128 pixels. The mathematical motivation to take Eq. (14) as a generative prior for power spectra is discussed in Arras et al. (2019b). As a rule of thumb, k-modes for which the data constrains the power spectrum very well will be recovered in great detail due to the nonparametric nature of the model. For k-modes on which the data provide little information, the power spectrum will be complemented by the prior which forces it into a falling power law whenever the data are not informative. If the actual physical process deviates strongly from a falling power law for the unobserved k-modes, the prior might artificially suppress or amplify the posterior uncertainty of the result, possibly biasing the uncertainty quantification. Figure 3 shows a few examples of prior samples of power spectra using our choice of hyperparameters. While the individual samples might not appear qualitatively different, it should on the one hand be noted that any kind of power spectrum is representable with our model given enough data and on the other hand that the figure depicts the power spectrum of the log-density on log–log scale. A small deviation in this figure can have a huge impact on the actual statistics. Reconstructing the power spectrum is equivalent to reconstructing the correlation kernel. We show our reconstructed normalized kernel as well as the one assumed by Lallement et al. (2018) in Fig. 4. Certain biases can appear when using a fixed kernel, for example introducing a characteristic length scale of the order of the FWHM of the kernel.

Puttingtogether likelihood and prior, the overall joint information Hamiltonian for our parameters ξ, τ, and ϕ is $$\begin{array}{l}P(d,\xi ,\tau ,\varphi )={\text{TG}}_{A_{\text{G}}^{\text{min}},A_{\text{G}}^{\text{min}},0.46,d}(R(\alpha \rho ))\mathcal{G}({(\xi ,\varphi ,\tau )}^T,\mathbb{1}),\\ \text{where }{\left[R(\alpha \rho )\right]}_i={\displaystyle {\int }_{L_i}\text{d}}l\alpha \rho (l),\\ \alpha \rho =\frac{1}{1000}\text{exp}\left({\mathbb{F}}_{\text{xk}}^{\text{-1}}\sqrt{{\text{P}}_{\text{s}}(\text{k})\text{V}}(\varphi ,\tau ){\xi }_{\text{k}}\right),\end{array}$$ $$\begin{array}{l}\sqrt{P_s(k)}(\varphi ,\tau )=\text{exp}(({\varphi }_{\text{m}}{\sigma }_{\text{m}}+\overline{m})\text{log}(\text{k})\text{+}{\varphi }_{\text{y}}{\sigma }_{\text{y}}+\overline{y}\\ +{\mathbb{F}}_{\text{log}(\text{k})\text{t}}^{-1}\left(\frac{a}{1+\frac{t^2}{t_0^2}}{\tau }_t\right)),\\ {\text{TG}}_{x_{\text{min}},x_{\text{max}},\sigma ,\overline{x}}(x)={\displaystyle \prod _i\frac{\mathcal{G}({\overline{x}}_i-x_i,{\sigma }^2)}{{\text{cdf}}_{\mathcal{G}(x_i,{\sigma }^2)}(x_{\text{max}})-{\text{cdf}}_{\mathcal{G}(x_i,{\sigma }^2)}(x_{\text{min}})}}\\ \text{for }x\in [x_{\text{min}},x_{\text{max}}],\end{array}$$

which is a truncated Gaussian, and $V={(600\text{\hspace{0.17em}pc})}^3$. The application, calculation of the gradient, and the application of the Fisher metric of Eq. (15) scales almost linearly with the number of voxels $N_{\text{side}}^3$. More specifically, it takes $\mathcal{O}(N_{\text{data}}{N}_{\text{side}}+N_{\text{side}}^3{\text{logN}}_{\text{side}})$ operations to evaluate Eq. (15), where N_data = 3 661 286 is the number of data points used in the reconstruction and N_side = 256 is the number of voxels per axis.

Fig. 2 Graphical representation of the data model for our reconstruction. The logarithmic dust density log (ρ) is a Gaussian process with a smooth Gaussian process power spectrum $S_{k{k}^{\prime }}=2\pi \delta (k-k^{\prime })P_s(k)$. The true extinctions AG are directly dependent on the dust density ρ on each line of sight. The measured extinctions $\widetilde{A_{\text{G}}}$ are assumed to be distributed around the true extinctions AG following a truncated Gaussian distribution as described in Sect. 3.2.

Fig. 3 Several prior samples of the logarithmic spatial correlation power spectrum in units of pc3.

Fig. 4 Log-normal process normalized two-point correlation reconstructed by our method (solid line) and imposed in the reconstruction by Lallement et al. (2018) (dashed line). One can see that the dust is assumed to be strongly correlated at a distance scale of up to about 30 pc. This plot shows normalized 1D cuts through the 3D Fourier transform of the log-normal spatial correlation power spectrum shown in Fig. 7.

4 Simulated data test

4.1 Datageneration

In this section, a test on simulated data is presented. This test enables us to compare the results of the reconstruction to a known ground truth. As ground-truth for the dust-density, a cube of interstellar medium simulated by the SILCC collaboration (Walch et al. 2015) was used, more specifically from the magneto-hydrodynamic simulation of the interstellar medium B6-1pc at 50 Myr published by Girichidis et al. (2018). This simulation result spans a cube of ${(512\text{\hspace{0.17em}pc})}^3$ in size. We computed our synthetic ground truth differential absorption ρ_mock from the gas density of the simulation ρ_sim via $${\rho }_{\text{mock}}\left(x-{(150,150,0)}^T\right)=\sqrt[3]{{\rho }_{\text{sim}}\left(\frac{512x}{600}\right){10}^{17}\frac{{\text{cm}}^3}{\text{g}}}\frac{1}{\text{pc}}.$$(20)

Thus, we stretch the 512 pc simulated cube to the 600 pc of our reconstruction, scale it with a constant factor, and shift it by 150 pc. The shift is performed in order to have an underdense region at the center. We also take the third root of the gas density in Eq. (20). There are two reasons for this. A practical motivation for taking the third root is that it reduces the dynamic range. If one does not do this, the sky will be dominated by one very small, but very strongly absorbing blob. A more physical motivation is that very dense regions lead to star formation. These forming stars again reduce the density by blowing the material out of these regions. This feedback mechanism was not included in the simulation by Girichidis et al. (2018) but was shown to have a strong impact on the gas density in a followup simulation by Haid et al. (2018). The third root can be seen as a very crude way of reducing the density in these overdense regions. To obtain the synthetic data from the ground truth differential extinction cube ρ_mock, the following operations were performed:

• 1.

  Sampling ground truth parallaxes ${\omega }_i\curvearrowleft \mathcal{G}({\omega }_i-{\tilde{\omega }}_i,{\sigma }_i^2)$ according to the parallax likelihood published by the Gaia collaboration.

• 2.

  Integrating the dust density from the center of the cube to the location of the sampled star location $\frac{1}{\omega }$ using the full resolution of 512³ voxels⁴.

• 3.

  Sampling an observed extinction magnitude according to the truncated Gaussian likelihood described in Sect. 3.2.

Fig. 5 Results of our test using simulated data. The rows show integrated dust extinction for sightlines parallel to the z-, x-, and y-axis, respectively.First column: test reconstruction and second column: ground truth synthetic extinction.

4.2 Results

We were able to recover a slightly smeared out version of the original synthetic extinction cube. Figure 5 shows integrations with respect to the x-, y-, and z-axis of the synthetic extinction cube and the reconstructed extinction cube are shown. This visually confirms the reliability of the reconstruction.

For a more quantitative analysis, we compared the reconstructed differential extinctions with the ground-truth voxel-wise. More specifically, we computed the uncertainty-weighted residual r $$\begin{array}{ll}r\hfill & =\frac{{\rho }_{\text{reconstruction}}-{\rho }_{\text{ground truth}}}{{\sigma }_{\rho }},\hfill \end{array}$$(21)

where ρ_{reconstruction} and σ_ρ are the posterior mean and standard deviation computed from the approximate posterior samples. The ground truth differential extinction was recovered well within the recovered approximate posterior uncertainty, apart from outliers which make up about 0.15% of the voxels. See Fig. 6 for a histogram of the uncertainty-weighted residual. Overall, the reconstruction seems very reliable on aqualitative and quantitatively level within the uncertainty for most of the voxels.

Fig. 6 Normalized histogram of the deviation of the reconstruction from the true solution, in sigmas (gray curve). The black curve is the probability density function of a standard normal distribution, which is plotted as a reference. We note that the values were clipped to the range from − 10 to 10, i.e., the bump in the gray curve at − 10 corresponds to outliers that can be up to 250 sigma. These outliers correspond to about 0.15% of all voxels.

5 Results from Gaia data

We reconstruct the dust density in a 600 pc cube using 256³ voxels, resulting in a resolution of ${(2.34\text{\hspace{0.17em}pc})}^3$ per voxel. For our reconstructed volume we also infer the spatial correlation power spectrum of the log-density; see Fig. 7.

In Fig. 8b one can see a projection of the reconstructed dust onto the sky in Galactic coordinates. Figure 9b shows the corresponding expected logarithmic dust density.

In Fig. 8a the projection of the dust reconstruction on the Galactic plane is shown. This view is especially interesting to study the dust morphology as this projection introduces no perspective-based distortion. It is especially suited to identifying underdense regions such as the local bubble in high resolution. A logarithmic plot of the projection on the Galactic plane can be seen in Fig. 9a. We show integrated dust density for sightlines parallel to the x-, y-, and z-axis in Fig. 13.

We provide posterior uncertainty estimation via samples. One should note that these uncertainties might be underestimated due to the variational approach taken in this paper. One can see a map of the expected posterior variance of the sky projection in Fig. 10a and in the plane projection in Fig. 11a. A slice through the galactic plane and corresponding uncertainties can be seen in Fig. 14.

The results of the reconstruction are provided online⁵ and at the CDS, and can be used under the terms of the ODC-By 1.0 license. Proper attribution should be given to this paper as well as to the Gaia collaboration (Gaia Collaboration 2016).

Now, we give an overview of known systematic effects and advice on how to use the provided dust map. We do not recommend using the outer 15 pc of the reconstruction. Periodic boundary conditions were assumed for algorithmic reasons, which leads to correlations leaking around the border of the cube. The inferred prior correlation kernel (Fig. 4) suggests that correlations are vanishing after 30 pc.

Furthermore, we provide posterior samples. When performing a further analysis of our reconstruction we recommend doing so for every sample in order to propagate errors.

It was observed in Sect. 4.2 that there can be a small number of outliers, that is, differential dust extinction values that are much larger than the reconstructed value, by amounts that cannot be explained by the reconstructed uncertainty.

Finally, we anticipate a perception threshold that leads to the absence of extremely low-density dust clouds. The two main reasons for this are that the truncated Gaussian likelihood provides less evidence in the regime where the extinction is close to zero and that variational Bayesian schemes are known to underestimate errors. Studying a larger volume will shed further lighton this subject, as sightlines for more distant stars still provide information about nearby dust clouds. One should note that the overall Gaia extinction data provide 20 times more sightlines than were used in this reconstruction.

Fig. 7 Log-normal process spatial correlation power spectrum inferred in our reconstruction (solid line) as well as the imposed power spectrum of Lallement et al. (2018) (dashed line). The shaded area around the solid line indicates 1σ error bounds. The y-axis is in units of pc3. The functions can be interpreted as the a priori expected value of ${\left|\mathbb{F}\text{ln}(\rho )\right|}^2/V$, where V is the volume on which the density ρ is defined, and F is the Fourier transform. The region between 0.0008/pc and 0.426/pc is almost power-law-like with a slope of 3.1, the spectral index of the power law.

6 Discussion

Here we discuss qualitative, quantitative, and methodological differences to other dust-mapping efforts. Table 1 shows a detailed break down of methodological differences to other papers. There are three notable differences between our method and other methods that we would like to stress. On the one hand, the data set used here is one of the largest of its kind used so far. Additionally, we use a large amount of data while still taking 3D correlations into account. Furthermore, we reconstruct the spatial correlation power spectrum. The motivation and impact of this was already briefly discussed in Sect. 3.3

We compare our dust map to other maps. Comparisons to 2D dust maps are only possible on a qualitative level, since it is not clear whether structures visible in the 2D maps that are not present in the 3D map are simply further away or are too noisy in the data for the algorithm to pick them up. On a qualitative level it is possible to see several morphological similarities between our reconstruction in Fig. 8b and the Planck dust map (Planck Collaboration I 2019) in Fig. 1. These figures also show that many dust structures that are not inside the Galactic plane are local features.

The two 3D dust maps mentioned in Sect. 1 permit a more thorough analysis. Figure 8 shows a compilation of projected dust densities for our reconstruction as well as the reconstruction of Lallement et al. (2018) and Green et al. (2018)⁶, restricted to the same volume as the reconstruction discussed in this paper. A logarithmic version of this figure is provided in Fig. 9.

While our map seems to agree on large scales with the other maps, there seems to be a pronounced tension in the predictions of the position of some dust clouds compared to the reconstruction of Lallement et al. (2018). Compared to the map of Green et al. (2018) our map recovers the small scales significantly better and suffers far less from radial smearing. It should be noted that Green et al. (2018) mapped a significantly larger part of our galaxy, and that the region that overlaps with our map was declared as relatively unreliable by the authors themselves. The differences are probably due to the different nature of the data sets used. Gaia DR2 which was used in our reconstruction contains a much greater number of data points than the data set used for the other reconstructions. These data points, taken from Gaia DR2, have a very small parallax error. Additionally, our reconstruction takes the full 3D correlation structure into account.

Our reconstruction and that of Green et al. (2018) allow uncertainties to be quantified using samples. A plot of uncertainties of the dust density reconstructions projected into the Galactic plane can be seen in Fig. 11. Uncertainties of the dust density reconstructions in the sky projection can be seen in Fig. 10.

To quantify the dynamic range of the reconstruction and as a prediction on the variability of the logarithmic dust density we calculated histograms of dust density which show the number of voxels that can be attributed to each dust density. These histograms can be seen in Fig. 12a. One can see that the histogram of our reconstruction extends slightly more towards high dust densities and substantially towards low dust densities. This is possibly because our reconstruction is more sharply resolved, thus regions of high dust density get captured better and bleed less into the regions where dust is absent.

We characterize the agreement between pairs of those reconstructions using the heat maps of their voxel-wise value pairs. These heat maps can be seen in Fig. 12. For two perfectly agreeing reconstructions, the heat map shows a line with a slope of 1. Again it can be seen that the dust density in our reconstruction varies significantly more than in the two other reconstructions. While all maps agree more or less for high dust densities, our dust map exhibits a vastly greater volume with low dust density.

The reconstruction of Lallement et al. (2018) is performed using Gaussian process regression, as is ours. Therefore, one can compute theprior Gaussian process correlation power spectrum used in their reconstruction. Figure 7 shows both ours and their inferred power spectrum. These two power spectra agree more or less for the larger modes (low k), where the data arevery constraining.

One canempirically compute power spectra of the dust density using a Fourier transformation. A comparison plot with all the three mentioned reconstructions can be found in Fig. 15. This shows a white noise floor in the reconstruction of Green et al. (2018), which can visually also be seen as small scale structures in the plane projections shown in Figs. 8 and 9.

Fig. 8 Left column: integrated dust extinction from − 300 to 300 pc for sightlines perpendicular to the Galactic plane. The image covers a 600 pc cube centered around the Sun. The units are e-folds of extinction. The coordinates are Galactic Cartesian coordinates. The Sun is at coordinate (0, 0), the Galactic center is located towards the left of the plot, and the Galactic West is located towards the top. Right column: all-sky integrated dust extinction maps of the same region, but for sightlines towards the location of the Sun. First row: result of the reconstruction discussed in this paper, second row: reconstruction performed by Lallement et al. (2018), last row: reconstruction by Green et al. (2018).

Fig. 9 Natural logarithmic version of Fig. 8.

Fig. 10 Uncertainty of the reconstruction of this paper derived from posterior samples (first row) and of the reconstruction of Green et al. (2018) (second row), both in the sky projection. The uncertainties are in the same units as the corresponding maps in Fig. 8, while the logarithmic uncertainties are dimensionless. First column: variance for the dust extinction and second column: variance of the logarithmic projected dust density on a natural log-scale, which can be interpreted as a relative error.

Fig. 11 Posterior uncertainty of the reconstruction of this paper derived from samples (first row) and of the reconstruction of Green et al. (2018) (second row) in the plane projection. The uncertainties are in the same units as the corresponding maps in Fig. 8, while the logarithmic uncertainties are dimensionless. First column: variance for the dust extinction and second column: variance of the logarithmic projected dust density on natural log-scale which can be interpreted as a relative error.

7 Conclusions

We were able to obtain highly resolved map of the local dust density using only Gaia data. This map agrees on large scales with previously published maps of Lallement et al. (2018) and Green et al. (2018), but shows significant differences on small scales. These differences might to a large degree stem from the different data used. Our map shows many structures that are visible in the Planck dust map (Planck Collaboration I 2019).

In comparison to previous maps, we were able to improve on 3D resolution while still being mostly consistent on the large scales. A comparison to 2D maps like the Planck dust map seems to confirm the features present in our map.

We find that the logarithmic density of dust exhibits a power-law power spectrum with a 3D spectral index of 3.1, corresponding to a 1D index of 1.1. This is a significantly harder spectrum than that expected for a passive tracer in Kolmogorov turbulence, which would be a 1D index of 5∕3. The harder spectrum is probably caused by the sharp edges of the local bubble and other ionization or dust evaporation fronts.

In contrast to other dust reconstructions, we predict very low dust densities inside the local bubble. This discrepancy is possibly an artifact of our reconstruction as there are known dust clouds in our vicinity; for example the northern high-latitude shells (Puspitarini & Lallement 2012) and the local Leo cold cloud (Peek et al. 2011). The Leo cold cloud is however considerably smaller than a voxel of our simulation. The possibility that Gaia extinction estimates are biased for small distances can also not be excluded.

We hope that by providing accurate reconstructions of the nearby dust clouds, further studies of dust morphology will be possible, as well as the construction of more accurate extinction models for photon observations in a large range of frequency bands.

Fig. 12 Normalized histograms of dust densities (panel a) and heat maps of voxel-wise correlations (panels b–d). The solid line corresponds to our reconstruction, while the dashed line is the reconstruction of Lallement et al. (2018) and the dash-dotted line is the reconstruction of Green et al. (2018). The other three plots are heat maps of voxel-wise correlations between reconstructed logarithmic dust densities, where the colors show bin counts. The black line in the heat maps is the identity function, corresponding to perfect correlation.

Fig. 13 Reconstructed dust density in different projections. The rows show integrated dust extinction for sightlines parallel to the z-, x-, and y-axis, respectively. In the first row, the Galactic center is located towards the left of the plot, in the other two rows the Galactic north is located towards the top of the plot. The cube is in Galactic coordinates, thus the x-axis is oriented towards the Galactic center and the z-axis is perpendicular to the Galactic plane. First column: integrated G-band extinction ine-folds of extinction, second column: logarithmic version of the first column.

Fig. 14 Reconstructed dust pseudo-density in a slice of the Galactic plane. First plot: differential dust extinction in the plane containing the Sun. Second plot of first row: logarithmic version of the first plot. Second row: corresponding uncertainty maps. The unit of the dust is G-band extinction ine-folds per parsec. The coordinates are Galactic Cartesian coordinates.

Fig. 15 Empirical spatial correlation power spectra of the reconstructed mean dust density in parsecs. The black line was computed from our reconstruction, the dark-gray line was computed from the reconstruction of Lallement et al. (2018), and the light-gray line was computed from the reconstruction of Green et al. (2018). For the reconstruction of Green et al. (2018), unspecified voxel values on sightlines that lacked data were replaced with a zero.

Acknowledgements

We acknowledge helpful commentaries and suggestion from our anonymous referee as well as fruitful discussions with S. Hutschenreuter, J. Knollmüller, P. Arras and others from the information field theory group at the MPI for astrophysics, Garching. 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.

References

All Tables

All Figures

Fig. 1 Top panel: natural logarithm of the thermal dust emission map produced by Planck Collaboration I (2019). The color-scheme was saturated to visually match that of the Gaia dust extinction data shown in the middle panel. Bottom panel: subset of Gaia extinction data within a 600 pc cube centered on the Sun; the data used in this paper. The scale is a natural logarithm of the extinction data in magnitudes. Data points in the same direction were averaged.

In the text

Fig. 2 Graphical representation of the data model for our reconstruction. The logarithmic dust density log (ρ) is a Gaussian process with a smooth Gaussian process power spectrum $S_{k{k}^{\prime }}=2\pi \delta (k-k^{\prime })P_s(k)$. The true extinctions AG are directly dependent on the dust density ρ on each line of sight. The measured extinctions $\widetilde{A_{\text{G}}}$ are assumed to be distributed around the true extinctions AG following a truncated Gaussian distribution as described in Sect. 3.2.

In the text

	Fig. 3 Several prior samples of the logarithmic spatial correlation power spectrum in units of pc3.
In the text

Fig. 4 Log-normal process normalized two-point correlation reconstructed by our method (solid line) and imposed in the reconstruction by Lallement et al. (2018) (dashed line). One can see that the dust is assumed to be strongly correlated at a distance scale of up to about 30 pc. This plot shows normalized 1D cuts through the 3D Fourier transform of the log-normal spatial correlation power spectrum shown in Fig. 7.

In the text

	Fig. 5 Results of our test using simulated data. The rows show integrated dust extinction for sightlines parallel to the z-, x-, and y-axis, respectively.First column: test reconstruction and second column: ground truth synthetic extinction.
In the text

Fig. 6 Normalized histogram of the deviation of the reconstruction from the true solution, in sigmas (gray curve). The black curve is the probability density function of a standard normal distribution, which is plotted as a reference. We note that the values were clipped to the range from − 10 to 10, i.e., the bump in the gray curve at − 10 corresponds to outliers that can be up to 250 sigma. These outliers correspond to about 0.15% of all voxels.

In the text

Fig. 7 Log-normal process spatial correlation power spectrum inferred in our reconstruction (solid line) as well as the imposed power spectrum of Lallement et al. (2018) (dashed line). The shaded area around the solid line indicates 1σ error bounds. The y-axis is in units of pc3. The functions can be interpreted as the a priori expected value of ${\left|\mathbb{F}\text{ln}(\rho )\right|}^2/V$, where V is the volume on which the density ρ is defined, and F is the Fourier transform. The region between 0.0008/pc and 0.426/pc is almost power-law-like with a slope of 3.1, the spectral index of the power law.

In the text

Fig. 8 Left column: integrated dust extinction from − 300 to 300 pc for sightlines perpendicular to the Galactic plane. The image covers a 600 pc cube centered around the Sun. The units are e-folds of extinction. The coordinates are Galactic Cartesian coordinates. The Sun is at coordinate (0, 0), the Galactic center is located towards the left of the plot, and the Galactic West is located towards the top. Right column: all-sky integrated dust extinction maps of the same region, but for sightlines towards the location of the Sun. First row: result of the reconstruction discussed in this paper, second row: reconstruction performed by Lallement et al. (2018), last row: reconstruction by Green et al. (2018).

In the text

	Fig. 9 Natural logarithmic version of Fig. 8.
In the text

Fig. 10 Uncertainty of the reconstruction of this paper derived from posterior samples (first row) and of the reconstruction of Green et al. (2018) (second row), both in the sky projection. The uncertainties are in the same units as the corresponding maps in Fig. 8, while the logarithmic uncertainties are dimensionless. First column: variance for the dust extinction and second column: variance of the logarithmic projected dust density on a natural log-scale, which can be interpreted as a relative error.

In the text

Fig. 11 Posterior uncertainty of the reconstruction of this paper derived from samples (first row) and of the reconstruction of Green et al. (2018) (second row) in the plane projection. The uncertainties are in the same units as the corresponding maps in Fig. 8, while the logarithmic uncertainties are dimensionless. First column: variance for the dust extinction and second column: variance of the logarithmic projected dust density on natural log-scale which can be interpreted as a relative error.

In the text

Fig. 12 Normalized histograms of dust densities (panel a) and heat maps of voxel-wise correlations (panels b–d). The solid line corresponds to our reconstruction, while the dashed line is the reconstruction of Lallement et al. (2018) and the dash-dotted line is the reconstruction of Green et al. (2018). The other three plots are heat maps of voxel-wise correlations between reconstructed logarithmic dust densities, where the colors show bin counts. The black line in the heat maps is the identity function, corresponding to perfect correlation.

In the text

Fig. 13 Reconstructed dust density in different projections. The rows show integrated dust extinction for sightlines parallel to the z-, x-, and y-axis, respectively. In the first row, the Galactic center is located towards the left of the plot, in the other two rows the Galactic north is located towards the top of the plot. The cube is in Galactic coordinates, thus the x-axis is oriented towards the Galactic center and the z-axis is perpendicular to the Galactic plane. First column: integrated G-band extinction ine-folds of extinction, second column: logarithmic version of the first column.

In the text

	Fig. 14 Reconstructed dust pseudo-density in a slice of the Galactic plane. First plot: differential dust extinction in the plane containing the Sun. Second plot of first row: logarithmic version of the first plot. Second row: corresponding uncertainty maps. The unit of the dust is G-band extinction ine-folds per parsec. The coordinates are Galactic Cartesian coordinates.
In the text

Fig. 15 Empirical spatial correlation power spectra of the reconstructed mean dust density in parsecs. The black line was computed from our reconstruction, the dark-gray line was computed from the reconstruction of Lallement et al. (2018), and the light-gray line was computed from the reconstruction of Green et al. (2018). For the reconstruction of Green et al. (2018), unspecified voxel values on sightlines that lacked data were replaced with a zero.

In the text