Charting nearby dust clouds using Gaia data only

TL;DR

This paper presents a detailed 3D map of local Galactic dust within 600pc using Gaia data, revealing the spatial correlation spectrum and enabling better understanding of dust's impact on astronomical observations.

Contribution

It introduces a hierarchical Bayesian model to reconstruct 3D dust density and its correlation spectrum solely from Gaia data, advancing local dust mapping techniques.

Findings

Reconstructed 3D dust density within 600pc of the Sun.

Derived a spectral index of 3.1 for the dust density power spectrum.

Provided publicly available dust maps and spectral data.

Abstract

Aims: Highly resolved maps of the local Galactic dust are an important ingredient for sky emission models. In nearly the whole electromagnetic spectrum one can see imprints of dust, many of which originate from dust clouds within 300pc. Having a detailed 3D reconstruction of these local dust clouds enables detailed studies, helps to quantify the impact on other observables and is a milestone necessary to enable larger reconstructions, as every sightline for more distant objects will pass through the local dust. Methods: To infer the dust density we use parallax and absorption estimates published by the Gaia collaboration in their second data release. We model the dust as a log-normal process using a hierarchical Bayesian model. We also infer non-parametrically the kernel of the log-normal process, which corresponds to the physical spatial correlation power spectrum of the log-density. Results: Using only Gaia data of the second Gaia data release, we reconstruct the 3D dust density and its spatial correlation spectrum in a 600pc cube centered on the Sun. We report a spectral index of the logarithmic dust density of $3.1$ on Fourier scales with wavelengths between 2pc and 125pc. The resulting 3D dust map as well as the power spectrum and posterior samples are publicly available for download.

Equations

$$P(s|d)=\frac{P(d|s)P(s)}{P(d)}$$

$$Q(s)$$

$$=\frac{\text{exp}\left(-\frac{1}{2}(s-m{)}^{†}{D}^{-1}(s-m)\right)}{{\left|2\pi D\right|}^{\frac{1}{2}}},$$

$$\text{KL}(Q,P)=\int \text{d}Q\text{ln}\frac{Q}{P}$$

$$(A_G{)}_i={\left[R(\rho )\right]}_i={\int }_{L_i}\text{d}l\alpha \rho (l).$$

$$P(\omega |\tilde{\omega },{\sigma }_{\omega })=\mathcal{G}(\omega -\tilde{\omega },{\sigma }_{\omega }^2)$$

$$P((A_G{)}_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_1},\rho )=$$

$$=\int \text{d}{\omega }_{i}P((A_G{)}_i,{\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i},\rho )$$

$$=\int \text{d}{\omega }_{i}P((A_G{)}_i|{\omega }_i,{\tilde{\omega }}_i,{\sigma }_{{\omega }_i},\rho )P({\omega }_i|\widetilde{{\omega }_i},{\sigma }_{{\omega }_i},\rho )$$

$$=\int \text{d}{\omega }_{i}P((A_G{)}_i|{\omega }_i,\rho )P({\omega }_i|{\tilde{\omega }}_i,{\sigma }_{{\omega }_i})$$

$$=\int \text{d}{\omega }_i\delta \left((A_G{)}_i-{\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_G{)}_i-\int \text{d}{\omega }_i{\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_G{)}_i-\int \text{d}{\omega }_i{\int }_{L_i(0)}\text{d}{l}_i\alpha \text{mathds}{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_G{)}_i-{\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).$$

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

$$\text{mathds}{1}_{[a,b]}(x)=\{\begin{array}{ll}1& \text{ for }x\in [a,b]\\ 0& \text{ otherwise}\end{array}$$

$$P(\widetilde{A_G}|\rho )$$

$$\text{for }d$$

$$\alpha \rho$$

$$\text{where }s$$

$$\sqrt{P_s(k)}=$$

$$\text{exp}\left(({\varphi }_m{\sigma }_m+\bar{m})\text{log}(k)+{\varphi }_y{\sigma }_y+\bar{y}+{\mathbb{F}}_{\text{log}(k)t}^{-1}\left(\frac{a}{1+\text{nicefrac}{t}^2{t}_0^2}\tau (t)\right)\right)$$

$$P(d,\xi ,\tau ,\varphi )$$

$$\mathcal{G}((\xi ,\varphi ,\tau {)}^T,\text{mathds}1)$$

$$\text{where }{\left[R(\alpha \rho )\right]}_i$$

$$\alpha \rho$$

$$\sqrt{P_s(k)}(\varphi ,\tau )$$

$$\text{exp}\left(({\varphi }_m{\sigma }_m+\bar{m})\text{log}(k)+{\varphi }_y{\sigma }_y+\bar{y}+{\mathbb{F}}_{\text{log}(k)t}^{-1}\left(\frac{a}{1+\text{nicefrac}{t}^2{t}_0^2}{\tau }_t\right)\right)$$

$${\text{TG}}_{x_{\text{min}},x_{\text{max}},\sigma ,\bar{x}}(x)$$

$$\text{for }x$$

$${\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}}.$$

$$r$$

Full text

*Methods. *To infer the dust density we use parallax and extinction estimates published by the Gaia collaboration in their second data release. We model the dust as a log-normal process using a hierarchical Bayesian model. We also infer non-parametrically the kernel of the log-normal process, which corresponds to the physical spatial correlation power spectrum of the log-density.

*Results. *Using only Gaia data of the second Gaia data release, we reconstruct the 3D dust density and its spatial correlation spectrum in a 600pc cube centered on the Sun. We report a spectral index of the logarithmic dust density of $3.1$ on Fourier scales with wavelengths between $2$pc and $125$pc. The resulting 3D dust map as well as the power spectrum and posterior samples are publicly available for download.

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 about 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). Their dust reconstruction, as 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 infra-red wavelengths, where it is dominated by dust emission. However, when mapping out dust using infra-red emission one is biased by the radiation field, as dust emission is not only proportional to the dust density but also to the amount of starlight absorbed by the dust. Furthermore, dust maps that were produced by mapping infrared light might contain extended infrared sources that are not within our galaxy, as was shown by Chiang & Ménard (2018). On the other hand, a hypothetical cold dust cloud cannot be seen in infra-red, 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 only be observed with less fidelity, due to dust obscuration. The first non-parametric reconstruction of galactic dust in 3D published is Arenou et al. (1992). Since then there have been many attempts to chart the dust density in 3D in increasing resolution, accuracy and for an ever greater part of our galaxy (Vergely et al., 1997, 2001; Chen et al., 2018; Kh et al., 2017, 2018a, 2018b; Gontcharov, 2012, 2017; Sale et al., 2014; Sale & Magorrian, 2018). A large driving force for 3D dust reconstruction and astronomy in general are 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 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 has been 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\,\text{pc})^{3}$.

In Green et al. (2018) a 3D dust map is 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 2kpc. 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 more accurate parallax information.

A prominent new survey is performed by the Gaia collaboration (Brown et al., 2016). In its second data release (DR2, Brown et al. (2018)) accurate parallaxes for roughly 2 billion stars are published. The provided catalogue also contains estimates of extinction coefficients for a subset of about 88 million stars, using spectral information of the Gaia satellite’s three energy bands. 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.

So far, 3D dust reconstructions have never been performed using solely Gaia data, instead Gaia data has been used for its accurate parallax measurements only 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 online on https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/dust.html or by its doi:10.5281/zenodo.2577337, 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 non-parametric kernel learning, first published in Enßlin & Frommert (2011). While the statistical model used here is up to minor details equivalent to the model introduced in that paper, the algorithm to arrive at an 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 (2018). 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)111The version of NIFTy used for this reconstruction is available on https://gitlab.mpcdf.mpg.de/ift/NIFTy .. 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 Sec. 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 Sec. 3. In Sec. 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 Sec. 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 Sec. 6. In Sec. 7 we summarize the findings of this paper.

2 Data

We used the data from the Gaia DR2 catalogue by Brown et al. (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:

the above mentioned data are available for the source 2. 2.

the parallax $\widetilde{\omega}$ is inside a 600pc cube around the Sun 3. 3.

the relative parallax error is sufficiently low, $\widetilde{\omega}/\sigma_{\widetilde{\omega}}>5$ 4. 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. Fig. 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).

3 Model

3.1 Algorithm

The algorithm is derived from Bayesian reasoning. In Bayesian reasoning information some data $d$ provides about a quantity of interest $s$ is calculated according to Bayes theorem:

[TABLE]

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^{3}$ degrees of freedom after discretization). There are three main ingredients necessary for the inference of the quantity of interest $s$:

The likelihood $P(d|s)$ of the data $d$ given a realization of the quantity of interest $s$. We describe our likelihood in Subsec. 3.2. 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 Subsec. 3.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=\text{ln}(\alpha\rho/\text{pc})$, henceforth called the logarithmic dust density. Hereby $\rho$ denotes the actual dust mass density and $\alpha$ the average G-band dust cross section per mass. The value of $\alpha$ is uncertain, which is why we report extinction densities, also called dust pseudo-densities, instead.

We approximate the posterior with a Gaussian

[TABLE]

by adopting a suitable mean $m$ and uncertainty dispersion $D$. The approximation is obtained by minimizing the Kullback-Leibler divergence (Kullback & Leibler, 1951)

[TABLE]

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:

Calculate the negative log-probability $-\text{log}(P(s,d))$ for the problem, disregarding normalization terms like the evidence $P(d)$. 2. 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. 3.

Choose the class of approximating distributions to be Gaussian with variable mean $m$ and covariance $D=(\mathds{1}+M_{m})^{-1}$, where $M_{m}$ is the Fisher information metric at the current $m$. Here $\mathds{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. 4.

Minimize Eq. (4) with respect to $m$ using Newton Conjugate gradient as 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.

A description of the used likelihood and the prior follows.

3.2 Likelihood

The likelihood $P(d|s)$ can be split into two parts, one part states how the true extinction depends on the dust density and one part that 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 $\rho$ imprints itself on the data. The extinction of light on the $i$-th line of sight $L_{i}$ is given by the line integral

[TABLE]

Here $\alpha$ is the average dust cross section per unit of mass and the line of sight $L_{i}=L_{i}(\omega)$ is dependent on the true parallax $\omega$. As noted in subsection 3.1, the value of $\alpha$ is uncertain and we reconstruct the dust extinction density $s=\text{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 $\omega$ of the star is uncertain. We assume the true parallax $\omega$ to be Gaussian distributed around the published parallax $\widetilde{\omega}$ with a standard deviation equal to the published parallax error $\sigma_{\omega}$:

[TABLE]

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, thus we have to marginalize them out to obtain the response:

[TABLE]

Here

[TABLE]

denotes the survival function of a standard normal distribution and

[TABLE]

denotes the indicator function for the closed interval $[a,b]$. Note that we did an approximation where we replace the true extinction by 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}$ exactly equal to the length of the line segment of $L_{i}$ within that voxel while being discounted by the probability $P(l|\widetilde{\omega},\sigma_{\omega})$ of that voxel being on the line of sight. Applying the response $R$ thus takes $\mathscr{O}(N_{\text{data}}N_{\text{side}})$ operations, where $N_{\text{data}}=3\,661\,286$ is the number of data points used in the reconstruction and $N_{\text{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 Sec. 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 600pc cube centered on the Sun.

The second part of the likelihood states how the published data $\widetilde{A_{G}}$ is distributed given the true extinctions $(A_{G})_{i}$. We use the data likelihood recommended by the Gaia collaboration in Andrae et al. (2018). This likelihood assumes the data $\widetilde{A_{G}}$ to be distributed according to a truncated Gaussian with a global variance $N=\left(0.46\,\text{mag}\right)^{2}\mathds{1}$. This leads to the likelihood

[TABLE]

Here $\text{cdf}_{\mathscr{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_{G}^{\text{min}}=0$ and $A_{G}^{\text{max}}=3.609\,\text{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 and 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 $\rho$ to be a-priori log-normal distributed with

[TABLE]

is assumed to be Gaussian distributed with Gaussian process kernel $S$. Here $\rho_{0}=\nicefrac{{1}}{{1000\,\text{pc}}}$ is a constant introduced to give $\rho$ 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_{kk^{\prime}}=2\pi\delta(k-k^{\prime})P_{s}(k)$. We non-parametrically 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 insights 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, thus one has to reconstruct it as 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.

Fig. 2 depicts the hierarchical Bayesian model for the extinction data $\widetilde{A_{G}}$ resulting from the logarithmic dust density $\text{ln}(\rho)$, 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 non-parametric deviations. It is up to minor details222The amplitude model given by Eq. 14 is not exactly equivalent to an integrated Wiener process, but shown by (Arras et al., 2019b) to be equivalent to it in a certain limit while still allowing a numerically stable transformation of the prior to a white Gaussian. an integrated Wiener process (Doob, 1953) on log-log-scale.

This is realized by the following formula:

[TABLE]

Here $\phi_{m}$, $\phi_{y}$ and $\tau(t)$ are the parameters to be reconstructed, $\sigma_{m}=1,\bar{m}=-4,\sigma_{y}=2.,\bar{y}=-16,\sigma_{y}=3.,a=11,t_{0}=0.2$ are fixed hyperparameters, $\mathbb{F}^{-1}_{\text{log}(k)t}$ denotes the inverse Fourier transform on log-scale, and $V=(600\,\text{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 has 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 $\tau$ 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 non-parametric 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 is 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. Fig. 3 shows a few examples of prior samples of power spectra using our choice of hyperparameters. While the individual samples might not look too different qualitatively, it should be noted on the one hand 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.

Putting together likelihood and prior, the overall joint information Hamiltonian for our parameters $\xi$, $\tau$, and $\phi$ is

[TABLE]

[TABLE]

is a truncated Gaussian, and $V=(600\,\text{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 $\mathscr{O}(N_{\text{data}}N_{\text{side}}+N_{\text{side}}^{3}\text{log}{N_{\text{side}}})$ operations to evaluate Eq. 15, where $N_{\text{data}}=3\,661\,286$ is the number of data points used in the reconstruction and $N_{\text{side}}=256$ is the number of voxels per axis.

4 Simulated Data Test

4.1 Data generation

In this section, a test on simulated data is presented. This test enables comparing the results of the reconstruction to a known ground truth. As ground truth dust density public data from 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 with size $(512\,\text{pc})^{3}$. We computed our synthetic ground truth differential absorption $\rho_{\text{mock}}$ from the gas density of the simulation $\rho_{\text{sim}}$ via

[TABLE]

Thus we stretch the $512\,\text{pc}$ simulated cube to the $600\,\text{pc}$ of our reconstruction, scale it with a constant factor, and shift it by $150\,\text{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 into 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 $\rho_{\text{mock}}$, the following operations were performed:

Sampling ground truth parallaxes $\omega_{i}\curvearrowleft\mathscr{G}(\omega_{i}-\widetilde{\omega}_{i},\sigma_{i}^{2})$ according to the parallax likelihood published by the Gaia collaboration. 2. 2.

Integrating the dust density from the center of the cube to the location of the sampled star location $\nicefrac{{1}}{{\omega}}$ using the full resolution of $512^{3}$ voxels333Note that the simulation of which the data is used was performed on an adaptive grid. The full resolution of $512^{3}$ is only realized in the high density regions.. 3. 3.

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

4.2 Results

We were able to recover a slightly smeared out version of the original synthetic extinction cube. In Fig. 5 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$

[TABLE]

where $\rho_{\text{reconstruction}}$ and $\sigma_{\rho}$ 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 a qualitative and quantitatively level within the uncertainty for most of the voxels.

5 Results from Gaia Data

We reconstruct the dust density in a 600pc cube using $256^{3}$ voxels, resulting in a resolution of $(2.34\,\text{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. 8(b) one can see a projection of the reconstructed dust onto the sky in galactic coordinates. Fig. 9(b) shows the corresponding expected logarithmic dust density.

In Fig. 8(a) 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 spot 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. 9(a). 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. 10(a) and in the plane projection in Fig. 11(a).

5.1 Using the reconstruction

The results of the reconstruction are provided online on https://wwwmpa.mpa-garching.mpg.de/~ensslin/research/data/dust.html, or by its doi:10.5281/zenodo.2577337, 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 et al, 2016).

We give an overview of known systematic effects and advice on how to use the provided dust map.

• •

We do not recommend to use the outer 15pc 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 30pc.

• •

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

• •

It was observed in Sec.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.

• •

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 light on this subject, as sightlines for more distant stars still provide information about nearby dust clouds. One should note that the overall Gaia extinction data provides 20 times more sightlines than were used in this reconstruction.

6 Discussion

Here we discuss qualitative, quantitative, and methodological differences to other dust mapping efforts. In table 1 a detailed break down of methodological differences to other papers are shown. There are three notable differences of our method to other methods that we would like to stress.

The here used dataset is one of largest one used so far. 2. 2.

We use a high amount of data while still taking 3D correlations into account. 3. 3.

We reconstruct the spatial correlation power spectrum. The motivation and impact of this already briefly discussed in Sec. 3.3

[TABLE]

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 of our reconstruction in Fig. 8(b) to the Planck dust map (The Planck Collaboration et al., 2018) 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 Sec. 1 permit a more thorough analysis. Fig. 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)444It should be noted that there is a new version (Green et al., 2019) that appeared during the revision of this paper., restricted to the same volume as the reconstruction discussed in this paper. A logarithmic version of this figure is provided by 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) we recover the small scales significantly better and suffer 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 to be not that reliable by the authors themselves. The differences are probably due to the different nature of the used datasets. The Gaia DR2 data used in our reconstruction has a vastly higher amount of data points than those 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 as well as the reconstruction of Green et al. (2018) permit quantifying uncertainties 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 how many voxels have which dust density. These histograms can be seen in Fig. 12(a). 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 how much pairs of those reconstructions agree by the heatmaps of their voxel-wise value pairs. These heatmaps can be seen in Fig. 12. For two perfectly agreeing reconstructions the heatmap would show a line with slope 1. Again it can be seen that the dust density in our reconstruction varies significantly more than in the two other reconstruction. While all maps agree more or less for high dust densities, our dust map exhibits vastly more volume with low dust density.

The reconstruction of Lallement et al. (2018) is performed using Gaussian process regression, as is ours. Thus one can compute the prior Gaussian process correlation power spectrum used in their reconstruction. Fig. 7 shows both our inferred power spectrum as well as their assumed power spectrum. These two power spectra agree more or less for the larger modes (low $k$), where the data is very constraining.

One can empirically compute power spectra of the dust density using a Fourier transformation. A comparison plot with all the three mentioned reconstruction 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.

7 Conclusions

We provide a 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 also shows significant differences on small scales. These differences might to a large degree stem from the different data used. Our map shows many structures visible in the Planck dust map (The Planck Collaboration et al., 2018). 2. 2.

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. 3. 3.

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 as 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. 4. 4.

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. 5. 5.

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.

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.