Resolving nearby dust clouds

Aims: Mapping the interstellar medium in 3D provides a wealth of insights into its inner working. The Milky Way is the only galaxy, for which a detailed 3D mapping can be achieved in principle. In this paper, we reconstruct the dust density in and around the local super-bubble. Methods: The combined data from surveys such as Gaia, 2MASS, PANSTARRS, and ALLWISE provide the necessary information to make detailed maps of the interstellar medium in our surrounding. To this end, we use variational inference and Gaussian processes to model the dust extinction density, exploiting its intrinsic correlations. Results: We reconstruct a highly resolved dust map, showing the nearest dust clouds at a distance of up to 400\,pc with a resolution of 1\,pc. Conclusions: Our reconstruction provides insights into the structure of the interstellar medium. We compute summary statistics of the spectral index and the 1-point function of the logarithmic dust extinction density, which may constrain simulations of the interstellar medium that achieve similar resolution.


Introduction
Although dust contributes only a small fraction in terms of mass, it is an important constituent of the interstellar medium (ISM) that is observable in many wavebands of the electromagnetic spectrum. Dust efficiently absorbs and scatters ultra-violet and visible range photons, obscuring large parts of the Galaxy and hiding star forming regions at these wavelengths. The dust absorbed energy is re-emitted in the infrared to microwave bands, offering a diagnostic for physical conditions of the ISM. The microwave emission of dust is a significant foreground to the Cosmic Microwave Background (CMB).
Dust plays a role in many processes that drive galactic evolution. Grain surfaces can adsorb material from interstellar gas and act as catalytic sites for chemical reactions. Stars, including the most massive ones, are observed to form from dusty molecular clouds. Thermal emission from dust grains can be an important cooling channel for these clouds and grains can drive their chemistry, suggesting that dust plays an important role in regulating the star formation process. Photons absorbed by dust can convey radiation pressure to interstellar matter or, if they are energetic enough, eject electrons, contributing to the heating of interstellar gas.
Finally, the distribution of dust can be used as a tracer of other quantities. A significant portion of the observed Galactic gamma rays in the GeV-range originates in dense clouds, where it is produced by hadronic interactions of cosmic rays with gas. This can be seen e.g. in the morphology of cosmic rays with hadronic spectrum from FERMI (Selig et al. 2015). Dust can be used to trace these dense clouds and identify gamma ray production sites. Another example is the magnetic field structure of the Galaxy, which is imprinted in the dust density, as dust filaments tend to be aligned to the line of sight magnetic field (Panopoulou et al. 2016). Dust also reveals the large scale dynamics and structure of the Galaxy, as the gravitational and differential rotation imprints on the filaments of dust.
Studying how dust is distributed in the Galaxy can not only provide understanding of its contents and structure but also into its inner workings and aid in the interpretation of observations in dust affected wavebands. Most 3D mapping efforts so far have aimed at reconstructing the distribution of dust in our galaxy on large scales. This is interesting as it reveals the structure of our galaxy such as spiral arms. Some notable recent contribution in this direction was provided by Green et al. (2019), who map three quarters of the sky using Gaia, 2MASS and PANSTARRS data using importance sampling on a gridded parameter space and by assuming a Gaussian process prior. Lallement et al. (2019) reconstructed a map extending out to 3 kpc with a 25 pc resolution based on Gaia and 2MASS data with Gaussian process regression. Chen et al. (2018) reconstructed a map extending out to 6 kpc with a 0.2 kpc radial resolution based on Gaia, 2MASS and WISE data with random forest regression. This paper can be regarded as a follow-up to . Some derivations are kept short here, and we advise Leike & Enßlin (2019) as a co-read for the statistically inclined reader. We focus on reconstructing only the nearby dust clouds, within ∼ 400 pc. While this prohibits revealing spiral arms, it enables us to achieve higher resolution. This way, we hope to be able to constrain simulations of the ISM, which achieve similar resolution. The map might also prove relevant for foreground corrections to the CMB, especially for CMB polarization studies. It was shown that most of the Galactic infrared polarization at high latitudes (|b| > 60) comes from close by regions around 200-300 pc (Skalidis & Pelgrims 2019). Corrections maps so far were based on infrared observations, and could be biased

Data
For our 3D reconstruction, we use combined observational data of Gaia DR2, ALLWISE, PANSTARRS and 2MASS. These data-sets were combined and processed to yield one consistent catalogue with stellar parameters by Anders et al. (2019). We use these high level preprocessed data for our reconstruction. Table 1 contains a summary of the columns we extract from this data set. We further only select sources that are inside an 800 pc × 800 pc × 600 pc cube centered on the Sun. To determine whether a source is inside this cube, we use their 84% distance quantile dist 84 . We assume a Gaussian error on the parallax, with mean m ω and standard deviation computed from the distance quantiles as Furthermore, we apply the following selection criteria: (4) ph ∈ Table 2 (5) σ ω/m ω < 0.3 (6) av 05 av 16 .
In words, we selected only stars, which do have clean starhorse pipeline flags, a clean Gaia flag, a specific photo-flag, and sufficiently small parallax error. We require the constraint on the photo-flag, because we only derived the noise statistic for stars with this flag. For details see Sec. 3.2. Additionally, we excluded stars for which the 5% V-band extinction quantile is equal to the 16% quantile, as this suggests that the pipeline had difficulties for these sources. These criteria result in the selection of a total of 5 096 642 sources. Fig. 1 shows an inverse-noise weighted average of our data projected onto the sky. To consistently combine the information of many data points, it is crucial to know the likelihood of a data point given the true amount of extinction for that source. We call this likelihood of one data point given its true extinction the noise statistic to distinguish it from the likelihood of the whole data set given the true 3-dimensional dust extinction distribution, which contains additional operations (see Sec. 3 for details). Unfortunately, Anders et al. (2019) did not publish a noise statistic for their data-set, and a noise statistic is not readily derivable from posterior quantiles. This is because posterior quantiles a give very limited information on the distribution P(a * |a) of the true extinction a * , while a full noise statistic would be given by P(a|a * ). In particular, there is no natural way to derive an analytic form of P(a * |a), inhibiting the calculation of P(a|a * ). We will thus choose a different approach to infer the noise statistic which we describe in Sec 3.2. Fig. 1: A Mollweide projection of the G-band extinction optical depth a to all sources in the used dataset. For this healpix nside 128 plot, we average the data sources that are in the same pixel, using the inverse noise dispersion as weights. Pixels with no data appear in white.

Response
If given the true 3D extinction density s(x), we can compute the extinction a * i for each source i by computing the line integral R i where θ i is the position of the i-th source projected onto the unit sphere and ω * i is the true parallax of the source. The true parallax ω * i is assumed to be Gaussian distributed with error and mean computed from the 16% and 84% percentiles of the starhorse dataset according to Eqs. (1), (2): Given this uncertainty of the true source distance, we can compute the expected extinction density for the source i as a weighted line integral R i where cdf denotes the cumulative density function of Eq. 9 with r = (ω * ) −1 . We compute the line integral of Eq. (10) on the fly for every step, using a parallelized fortran code 1 . The uncertainty of the true position of the source introduces a source dependent supplementary noise contribution σ 2 i . This uncertainty arises due to the uncertainty of the true source distance,   (2019) which introduces uncertainty on the line of sight extinction even when given the true extinction density s. The standard deviation of this supplementary noise contribution can be computed as The last inequality holds as P(ω * i |ω 16,i , ω 84,i ) has strictly more variance than P(ω * i |ω 16,i , ω 84,i , s). We sample this additional noise contribution before every step of our algorithm. We do this by drawing M = 20 samples j of parallaxes ω j i according to the statistic given by Eq. (9). We then compute as the sample variance of the extinction estimate using the samples j and the current reconstructed dust extinction density s. This error correction was not done in . However, for this paper the smaller data uncertainty and slightly higher parallax error of the sources raises the importance of computing this error correction, while the use of the new code for the response enables its calculation.

Noise Statistic
The noise statistic specifies how probable an observed G-band extinction value is, given one would know the true amount of Gband extinction for that source. Since there is no detailed noise statistic published for the dataset we use, we have to construct it.
To do this, we look at regions of the sky where there is no significant amount of dust expected. These regions were identified by using the Planck dust map (The Planck Collaboration et al. 2018), more specifically the dust map from the COMMANDER pipeline of the 2014 Planck data release. Here, regions with less than exp(2) µK /rJ were taken to be dustless. This criterion selects 606 pixels of the healpix nside 256 dust map, corresponding to 0.077% of the sky. For every SH_PHOTOFLAG for which we have more than 100 values in these dustless regions, we calculate the mean m ph and standard deviation σ ph of all G-band extinctions. Using these values, we define the probability to measure an extinction a given the true extinction a * as P(a|a * , SH_PHOTOFLAG = ph) = G (a|a * + m ph , σ 2 ph ) .   Table 2 show our used means m ph and standard deviations σ ph for all used photoflags ph. As can be seen by investigating Table  2, the mean values deviate strongly from zero, and correcting the zero-point is vital to our reconstruction. Note that because we fix the noise statistic for an actual extinction value of zero, the reconstruction might be biased for high extinction values. We discuss some biases that could be attributed to this effect in Sec. 6.2.

Prior
We fold our physical knowledge into the prior of the dust extinction density. We choose the exact same prior model as in . We assume the extinction density s to be positive and spatially correlated. This can be enforced by assuming a log-normal Gaussian process prior where ρ 0 is the prior median extinction density and T is the correlation kernel of the Gaussian process τ. The prior median extinction density is a hyper-parameter of our model and we choose ρ 0 = 1 /1000pc −1 . We infer the kernel T during our reconstruction. This can be achieved by rewriting s in terms of a generative model where all ξ are a-priori standard normal distributed and T k (ξ T ) is a non-parametric model for the Fourier transformed correlation kernel T k , also called the spatial correlation power spectrum. One should note that this model is degenerate, any change in T k can be absorbed into ξ k instead, as only the product of these two fields enters the overall dust extinction density s. Because of this property, the reconstructed power spectrum T k does not have to be the empirical power spectrum of s x , that can be calculated by Fourier transforming and binning. To avoid misunderstandings and artifacts from the degenerate model, we mainly report the empirical power spectrum in this paper, which is computed from posterior samples of s x . We now focus on our model for the power spectrum T k . This model assumes the spatial correlation power spectrum to be a preferentially falling power-law, but allows for arbitrary deviations. It can be written as where the first part describes a linear function on log-log scale, i.e. a power law; and the second part describes the nonparametric deviations which are assumed to be differentiable on log-log-scale. The operation Exp * denotes the exponentiation of the coordinate system. More explicitly, F ln(k)t is Fourier transformation on log-log scale, and the function sym is defined as where

Algorithm
We combine the prior and the likelihood into one generative model of the data. We compute approximate posterior samples using Metric Gaussian Variational Inference (MGVI) (Knollmüller & Enßlin 2019). This variational approach alternates between drawing samples around the current estimate for the latent parameters and optimizing the current estimate using the average gradient of the samples. The final set of samples is used to derive an uncertainty estimate on all our maps as well as on all derived quantities. For further parallelization, we split the problem into the 8 octants. Each octant has size 410pc × 410pc × 310pc, such that they overlap for 20pc.
We hereby use a threefold parallelization scheme, parallelizing by octants, parallelizing by samples and a parallelized response. The latter two parallelizations are enabled by our new fortran implementation, which computes the arising line integrals (Eq. (10)) on the fly. This is in contrast to our previous paper , where we computed the line integral using sparse matrices. Computing the response on the fly takes approximately the same time, but does not have any additional memory requirements and therefore allows for parallelization and a larger reconstruction.
The total number of degrees of freedom is ≈ 417 Million, exceeding those of our previous map by a factor of 30. The total computation time was about 2 weeks of wall clock time, or about 0.5 Million CPUh on 1920 cores.
The final samples of the independently reconstructed octants are combined into the full reconstruction using a differentiable variance-preserving interpolation scheme. The details are described in Appendix A.
One noteworthy point is that we cut away the outer 30pc due to artifacts from periodic boundary conditions, resulting in a final map volume of 740 pc × 740 pc × 540 pc.

Results
We were able to reconstruct the nearby dust clouds. Fig. 2 shows various maps produced from our result and their relative uncertainty. The maps show tendrils and filaments of dust on scales as small as 2pc up to scales of several hundred parsecs, at which they become disconnected.
All octants inferred similar logarithmic convolution kernels, as can be seen in Fig.3. These correlation kernels were computed by taking a slice out of the reconstructed Fourier transformed square root power spectrum.
A comparison of the empirical power spectra of the 8 different octants can be found in Fig. 4. Most octants have very similar power spectra, only octant 3 deviates strongly. This octant, located at 180 < l ≤ 270 and b > 0 (disregarding the overlap), is strongly devoid of dust, explaining the significantly lower power spectrum.
For the power of the full-volume extinction density, we find a power law with spectral index of 2.52 ± 0.015 at scales from 2pc to 100pc. For the logarithmic power, we report a spectral index of 2.82 ± 0.022 at scales from 2.3pc to 125pc.
Using our reconstruction, we can determine distances to nearby dust clouds. We derive two distance maps. Fig. 5 shows the distance to the nearest dust clouds in all directions, as well as an uncertainty on that distance estimate. Note that we compute the distance by checking for the first voxel that exceeds a the threshold of 0.005 e-folds per pc of extinction density. Some of our samples do not exclude the existence of nearby dense clouds, which raises the uncertainty in the corresponding directions tremendously. Fig. 6 shows the distance to the densest dust clouds in all directions, as well as an uncertainty on that distance estimate. Note that the uncertainty estimate is quite high on the boundaries of dust clouds, as the reconstruction is uncertain which voxel is densest along these lines of sight.

Comparison
An implicit assumption of the algorithm is that the voxels are smaller than the achievable resolution. Phrased in physical terms, an increase in pixel resolution can be regarded as a renormalization and we need to reach the continuous limit, i.e the limit of negligible discretization effects, for the algorithm to work. This is a byproduct of the inference of the power spectrum, if the achieved posterior resolution is of the order of the imposed voxel resolution, then the reconstruction changes drastically from one voxel to another and the extinction of sources behind an affected voxel also changes dramatically at the boundary. This sudden change in extinction is not compatible with a falling spatial correlation power law in Fourier space, thus the reconstructed power will fall less steep than the real one. This would significantly hamper the ability of the algorithm to extrapolate between measurements. We avoid this behavior by significantly increasing resolution compared to our previous reconstruction . However, it is conceivable that the reconstruction would still benefit from increasing the amount of voxels. We recommend distrusting the smallest scales of our reconstruction, only at scales of 2 pc or larger can the result be considered to be stable. This resolution limit was deduced from the reconstructed logarithmic correlation kernels as seen in Fig.3. At this limit, the reconstructed correlation kernels of the logarithmic dust extinction density have fallen to 10%.
A comparison of our results to  can be found in Fig. 7; see Fig. 8 for a logarithmic version.
We compare our results to the map of Green et al. (2019). Fig. 9 shows column density comparisons of the two reconstructions. Fig. 10 shows the same column densities, but on a logarithmic scale. A more detailed comparison to Green et al. (2019) in angular coordinates can be seen in Fig. 11.
In contrast to our old map, we use the dataset of Anders et al. (2019), which provides more sources and tighter constraints on the parallax and G-band extinction than the previously used Gaia data. The new reconstruction has a volume of 800 pc×800 pc× 600pc, compared to the (600 pc) 3 cube in . Furthermore, using a designated fortran routine for the computation of the line of sight integrals lead to the necessary speedup to handle the additional data constraints and massively more degrees of freedom. Finally, in the new reconstruction the parallax error is propagated into the measurement error, causing extinction values with stars of high parallax error to be less informative. In Fig. 7 one can see dust column densities along Galactic x, y, and z coordinates. Both dust maps agree on the morphology of large dust clouds on large scales. However, the current dust map contains significantly more dust. Part of the reason is that the data we use in the reconstruction of this paper has higher resolution and lower noise, allowing more dust to be reconstructed. We also believe the data used in  to be slightly biased to underestimating the amount of dust, an effect that accumulates in a reconstruction that uses many data points. In contrast, the data used in this reconstruction might have a tendency to overestimate the amount of dust, despite our effort to calibrate the zero point (see Sec. 3.2).
We furthermore reconstruct our correlation kernel nonparametrically, which should lead to an unbiased estimate of the power spectrum. Fig. 12 shows power spectra of , the reconstruction of this paper, and of the reconstruction of Green et al. (2019). Our new reconstruction and Leike & Enßlin (2019) seem to have quite consistent power spectra. The general tendency of the falling power law is also remarkably consistent with Green et al. (2019), however at scales of a few parsec the power spectrum of Green et al. (2019) flattens which we believe to be an artifact of how we put their reconstruction on a cartesian grid, i.e. the boundaries of the reconstructions intrinsic voxels introduce steep cuts which flatten the resulting power spectrum. However, none of the power spectra are consistent within the uncertainties estimates. While this seems problematic, one has to bear in mind that all reconstructions focus on dust in differing regions, potentially explaining the difference in the power spectrum. In Fig. 13 we show the power spectra of the logarithmic reconstructions. These seem to be less consistent in general, however one has to bear in mind that the logarithmic power spectrum is dominated by regions of low dust content, as these occur more frequently. Using Gaia data our method was found to underestimate low dust regions, and we anticipate that with the starhorse data we tend to overestimate low-dust regions. Nonetheless we find that the spectral index of 2.82 ± 0.022 at scales from 2.3 pc to 125 pc is compatible with the empirical spectral index of  within a 2σ joint uncertainty margin. The spectral index of the empirical power spectrum of  is 3.2 ± 0.14. 2 The logarithmic power spectrum of Green et al. (2019) seems to be inconsistent with our measurements. However, this effect is probably due to how we treat the missing values in that map, where a quarter of the sky was not measured. We have to set these values and every possible choice will impact the derived power spectra. We chose to set them to 10 −7 , which has minimal impact on the power spectrum on linear scale, but biases the power spectrum of the logarithmic dust extinction density and could potentially explain the difference. Fig. 14 shows a histogram of dust extinction density per voxel. One can see a good agreement between the histogram of our old and our current reconstruction in the region between 10 −3 pc −1 and 10 −1 pc −1 . A dust extinction density of 10 −4 pc −1 integrated to the boundary of our simulation cube yields an integrated extinction of 0.046, which is below our noise level even when pooling the information of many stars. For this reason, we do not show the histogram below 10 −4 pc −1 as its shape is mostly dependent on how the reconstruction extrapolates into dustless region. From the histogram it can be seen that the dust density is well described by a log-normal distribution. Note that since we show only the part of the histogram that has high signal to noise, this result should be relatively unbiased by our choice of prior. The fitted log-normal model has a standard deviation of σ = 1.906 ± 0.009 and a mean of m = −9.79 ± 0.04.

Using the reconstruction
One should note that the reconstruction shows a non-negligible amount of dust in the local bubble. We believe that the level found is an artifact of our noise statistic. As described in Sec. 3, our data model involves some heuristics which might systematically affect the reconstruction. This causes estimates made with this data to be biased and we were not able to fully correct for this bias. When integrating the reconstructed dust density to 70 pc, we find that the nearby dust looks like a smeared out version of farther dust clouds, indicating that it is indeed an artifact related to systematic data biases.
The posterior samples of the extinction density are available for download under https://doi.org/10.5281/zenodo. 3750926 or by its DOI 10.5281/zenodo.3750926 . When using the reconstruction we advise to beware of systematic overestimations of dust, especially in the local bubble. When deriving numeric quantities, we advise to do so for every sample and then estimate the mean and standard deviation of the results in order to get an error estimation.

Implications
Our map can be used to constrain simulations of the ISM. For example, in simulations of radiatively cooling dust clouds in hot winds, it has been shown that dust density power spectra are flatter than was previously thought (Sparre et al. 2019). Our maps show power spectra compatible with these simulations, and morphologically similar structures. Our reconstructed spectral index of 2.82 ± 0.022 at scales from 2.3 pc to 125 pc could be used to constrain parameters of sub-grid models of simulations of the ISM. Furthermore we find the density histogram of the logarithm of the G-band dust extinction density in e-folds per parsec shown in Fig. 14 is well described by a log-normal distribution with standard deviation σ = 1.906±0.009 and mean m = −9.79±0.04 on extinction density scales from 10 −4 pc −1 to 1 pc −1 .
(e) (f) Fig. 2: Result of our 3D dust reconstruction. The first column shows dust extinction, the second shows the relative error. The first row shows the integrated extinction in e-folds in a Mollweide projection of the whole reconstructed box of 740 pc × 740 pc × 540 pc. The second row also shows integrated extinction in e-folds in the same box, but integrated normal to the Galactic plane instead of radially. The third row shows differential extinction in e-folds per parsec in a slice along the Galactic plane.      . We overplot a log-normal model that was fitted to our reconstructed logarithmic extinction density (dash-dotted line). The curve of the fit is described by f (x) ∝ exp −0.5σ −2 (ln(x) − m) with σ = 1.906 and m = −9.79 and follows our empirical distribution function closely.