Resolving nearby dust clouds
AAstronomy & Astrophysics manuscript no. main c (cid:13)
ESO 2020April 16, 2020
Resolving nearby dust clouds
R. H. Leike , M. Glatzle , and T. A. Enßlin Max Planck Institute for Astrophysics, Karl-Schwarzschildstraße 1, 85748 Garching, Germany Ludwig-Maximilians-Universität, Geschwister-Scholl Platz 1, 80539 Munich, Germany Physik-Department, Technische Universität München, James-Franck-Str. 1, 85748 Garching, GermanyReceived XXXX, accepted XXXX
ABSTRACT
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 localsuper-bubble.
Methods.
The combined data from surveys such as Gaia, 2MASS, PANSTARRS, and ALLWISE provide the necessary informationto make detailed maps of the interstellar medium in our surrounding. To this end, we use variational inference and Gaussian processesto 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 of1 pc.
Conclusions.
Our reconstruction provides insights into the structure of the interstellar medium. We compute summary statistics of thespectral index and the 1-point function of the logarithmic dust extinction density, which may constrain simulations of the interstellarmedium that achieve similar resolution.
Key words.
ISM: dust, extinction – Galaxy: local interstellar matter – methods: data analysis
1. 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 electromagneticspectrum. Dust e ffi ciently absorbs and scatters ultra-violet andvisible range photons, obscuring large parts of the Galaxy andhiding star forming regions at these wavelengths. The dust ab-sorbed energy is re-emitted in the infrared to microwave bands,o ff ering a diagnostic for physical conditions of the ISM. Themicrowave emission of dust is a significant foreground to theCosmic Microwave Background (CMB).Dust plays a role in many processes that drive galactic evolu-tion. Grain surfaces can adsorb material from interstellar gas andact as catalytic sites for chemical reactions. Stars, including themost massive ones, are observed to form from dusty molecularclouds. Thermal emission from dust grains can be an importantcooling channel for these clouds and grains can drive their chem-istry, suggesting that dust plays an important role in regulatingthe star formation process. Photons absorbed by dust can conveyradiation pressure to interstellar matter or, if they are energeticenough, eject electrons, contributing to the heating of interstellargas.Finally, the distribution of dust can be used as a tracer ofother quantities. A significant portion of the observed Galacticgamma rays in the GeV-range originates in dense clouds, whereit is produced by hadronic interactions of cosmic rays with gas.This can be seen e.g. in the morphology of cosmic rays withhadronic spectrum from FERMI (Selig et al. 2015). Dust can beused to trace these dense clouds and identify gamma ray produc-tion sites. Another example is the magnetic field structure of theGalaxy, which is imprinted in the dust density, as dust filamentstend to be aligned to the line of sight magnetic field (Panopoulou et al. 2016). Dust also reveals the large scale dynamics and struc-ture of the Galaxy, as the gravitational and di ff erential rotationimprints on the filaments of dust.Studying how dust is distributed in the Galaxy can not onlyprovide understanding of its contents and structure but also intoits inner workings and aid in the interpretation of observations indust a ff ected wavebands. Most 3D mapping e ff orts so far haveaimed at reconstructing the distribution of dust in our galaxyon large scales. This is interesting as it reveals the structureof our galaxy such as spiral arms. Some notable recent contri-bution in this direction was provided by Green et al. (2019),who map three quarters of the sky using Gaia, 2MASS andPANSTARRS data using importance sampling on a gridded pa-rameter space and by assuming a Gaussian process prior. Lalle-ment et al. (2019) reconstructed a map extending out to 3 kpcwith a 25 pc resolution based on Gaia and 2MASS data withGaussian process regression. Chen et al. (2018) reconstructed amap extending out to 6 kpc with a 0 . ∼
400 pc. While this prohibits revealing spiral arms, itenables us to achieve higher resolution. This way, we hope tobe able to constrain simulations of the ISM, which achieve simi-lar resolution. The map might also prove relevant for foregroundcorrections to the CMB, especially for CMB polarization stud-ies. It was shown that most of the Galactic infrared polarizationat high latitudes ( | b | >
60) comes from close by regions around200-300 pc (Skalidis & Pelgrims 2019). Corrections maps sofar were based on infrared observations, and could be biased
Article number, page 1 of 16 a r X i v : . [ a s t r o - ph . GA ] A p r & A proofs: manuscript no. main through di ff erent starlight illumination or di ff ering dust temper-atures.
2. Data
For our 3D reconstruction, we use combined observational dataof Gaia DR2, ALLWISE, PANSTARRS and 2MASS. Thesedata-sets were combined and processed to yield one consistentcatalogue with stellar parameters by Anders et al. (2019). Weuse these high level preprocessed data for our reconstruction.Table 1 contains a summary of the columns we extract fromthis data set. We further only select sources that are inside an800 pc ×
800 pc ×
600 pc cube centered on the Sun. To deter-mine whether a source is inside this cube, we use their 84% dis-tance quantile dist . We assume a Gaussian error on the paral-lax, with mean m ω and standard deviation computed from thedistance quantiles as m ω =
12 (1 / dist + / dist ) (1) σ ω =
12 (1 / dist − / dist ) . (2)Furthermore, we apply the following selection criteria:SH_OUTFLAG = =
000 (4)ph ∈ Table 2 (5) σ ω / m ω < . (cid:44) av . (7)In words, we selected only stars, which do have clean starhorsepipeline flags, a clean Gaia flag, a specific photo-flag, and suf-ficiently small parallax error. We require the constraint on thephoto-flag, because we only derived the noise statistic for starswith this flag. For details see Sec. 3.2. Additionally, we excludedstars for which the 5% V-band extinction quantile is equal to the16% quantile, as this suggests that the pipeline had di ffi cultiesfor these sources.These criteria result in the selection of a total of 5 096 642sources. Fig. 1 shows an inverse-noise weighted average of ourdata projected onto the sky. To consistently combine the infor-mation of many data points, it is crucial to know the likelihood ofa data point given the true amount of extinction for that source.We call this likelihood of one data point given its true extinc-tion the noise statistic to distinguish it from the likelihood ofthe whole data set given the true 3-dimensional dust extinctiondistribution, which contains additional operations (see Sec. 3 fordetails). Unfortunately, Anders et al. (2019) did not publish anoise statistic for their data-set, and a noise statistic is not read-ily derivable from posterior quantiles. This is because poste-rior quantiles a give very limited information on the distribu-tion P ( a ∗ | a ) of the true extinction a ∗ , while a full noise statisticwould be given by P ( a | a ∗ ). In particular, there is no natural wayto derive an analytic form of P ( a ∗ | a ), inhibiting the calculationof P ( a | a ∗ ). We will thus choose a di ff erent approach to infer thenoise statistic which we describe in Sec 3.2. Fig. 1: A Mollweide projection of the G -band extinction opticaldepth a to all sources in the used dataset. For this healpix nside128 plot, we average the data sources that are in the same pixel,using the inverse noise dispersion as weights. Pixels with no dataappear in white.
3. Likelihood
If given the true 3D extinction density s ( x ), we can compute theextinction a ∗ i for each source i by computing the line integral R i a ∗ i = R ω ∗ i ( s ) = (cid:90) ω ∗ i s ( r θ i )d r , (8)where θ i is the position of the i -th source projected onto the unitsphere and ω ∗ i is the true parallax of the source. The true parallax ω ∗ i is assumed to be Gaussian distributed with error and meancomputed from the 16% and 84% percentiles of the starhorsedataset according to Eqs. (1), (2): ω ∗ i (cid:120) G ( ω ∗ i |
12 ( ω , i + ω , i ) ,
14 ( ω , i − ω , i ) ) . (9)Given this uncertainty of the true source distance, we cancompute the expected extinction density for the source i as aweighted line integral R i (cid:10) a ∗ i (cid:11) P ( ω ∗ i | ω , i ,ω , i ) = R i ( s ) = (cid:68) R ω ∗ i ( s ) (cid:69) P ( ω ∗ i | ω , i ,ω , i ) = (cid:90) ω ∗ i s ( r θ i )(1 − cdf( r | ω , i , ω , i ))d r , (10)where cdf denotes the cumulative density function of Eq. 9 with r = ( ω ∗ ) − . We compute the line integral of Eq. (10) on the flyfor every step, using a parallelized fortran code .The uncertainty of the true position of the source introduces asource dependent supplementary noise contribution (cid:98) σ i . This un-certainty arises due to the uncertainty of the true source distance, https://gitlab.mpcdf.mpg.de/mglatzle/gda_futils Article number, page 2 of 16. H. Leike et al.: Resolving nearby dust cloudsname in Anders et al. (2019) our notation explanationdist16 dist
16% distance quantiledist50 dist
50% distance quantiledist84 dist
84% distance quantileag50 a G -band extinction quantileSH_PHOTOFLAG ph photo-bands used for data pointSH_GAIAFLAG SH_GAIAFLAG output flag of GaiaSH_OUTFLAG SH_OUTFLAG output flag of the starhorse pipeline Table 1: data columns extracted from Anders et al. (2019)which introduces uncertainty on the line of sight extinction evenwhen given the true extinction density s . The standard deviationof this supplementary noise contribution can be computed as (cid:98) σ i = Var (cid:2) P (cid:0) a ∗ i | ω , i , ω , i , s (cid:1)(cid:3) = Var (cid:34)(cid:90) d ω ∗ i P (cid:0) a ∗ i , ω ∗ i | ω , i , ω , i , s (cid:1)(cid:35) = Var (cid:34)(cid:90) d ω ∗ i P (cid:0) a ∗ i | s , ω ∗ i (cid:1) P (cid:0) ω ∗ i | ω , i , ω , i , s (cid:1)(cid:35) ≤ Var (cid:34)(cid:90) d ω ∗ i P (cid:0) a ∗ i | s , ω ∗ i (cid:1) P (cid:0) ω ∗ i | ω , i , ω , i (cid:1)(cid:35) = Var (cid:34)(cid:90) d ω ∗ i δ (cid:16) a ∗ i − R ω ∗ i ( s ) (cid:17) P (cid:0) ω ∗ i | ω , i , ω , i (cid:1)(cid:35) . (11)The last inequality holds as P ( ω ∗ i | ω , i , ω , i ) has strictly morevariance than P ( ω ∗ i | ω , i , ω , i , s ). We sample this additionalnoise contribution before every step of our algorithm. We do thisby drawing M =
20 samples j of parallaxes ω ji according to thestatistic given by Eq. (9). We then compute (cid:98) σ i = M (cid:88) j R ω ji ( s ) (12)as the sample variance of the extinction estimate using the sam-ples j and the current reconstructed dust extinction density s .This error correction was not done in Leike & Enßlin (2019).However, for this paper the smaller data uncertainty and slightlyhigher parallax error of the sources raises the importance of com-puting this error correction, while the use of the new code for theresponse enables its calculation. The noise statistic specifies how probable an observed G -bandextinction value is, given one would know the true amount of G -band extinction for that source. Since there is no detailed noisestatistic 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 sig-nificant amount of dust expected. These regions were identifiedby using the Planck dust map (The Planck Collaboration et al.2018), more specifically the dust map from the COMMANDERpipeline of the 2014 Planck data release. Here, regions with lessthan exp(2) µ K / rJ were taken to be dustless. This criterion selects606 pixels of the healpix nside 256 dust map, corresponding to0 . m ph and standard deviation σ ph of all G -band extinc-tions. Using these values, we define the probability to measurean extinction a given the true extinction a ∗ as P ( a | a ∗ , SH_PHOTOFLAG = ph) = G ( a | a ∗ + m ph , σ ) . (13) ph = SH_PHOTOFLAG mean m ph standard deviation σ ph GBPRP 0.493 0.439GBPRPJHKs 0.131 0.259GBPRPJHKs
Table 2: SH_PHOTOFLAG values and the corresponding meanand standard deviations for sources in dustless regions. Regionsare considered as dustless if the Planck dust map shows weakeremission than exp(2) µ K / rJ .Table 2 show our used means m ph and standard deviations σ ph for all used photoflags ph. As can be seen by investigating Table2, the mean values deviate strongly from zero, and correctingthe zero-point is vital to our reconstruction. Note that becausewe fix the noise statistic for an actual extinction value of zero,the reconstruction might be biased for high extinction values. Wediscuss some biases that could be attributed to this e ff ect in Sec.6.2.
4. Prior
We fold our physical knowledge into the prior of the dust extinc-tion density. We choose the exact same prior model as in Leike& Enßlin (2019). We assume the extinction density s to be pos-itive and spatially correlated. This can be enforced by assuminga log-normal Gaussian process prior s x = ρ exp( τ x ) , (14) τ (cid:120) G ( τ | , T ) , (15)where ρ is the prior median extinction density and T is the cor-relation kernel of the Gaussian process τ . The prior median ex-tinction density is a hyper-parameter of our model and we choose ρ = / pc − . We infer the kernel T during our reconstruc-tion. This can be achieved by rewriting s in terms of a generativemodel s x = ρ exp( F (cid:112) T k ( ξ T ) ξ k ) , (16)where all ξ are a-priori standard normal distributed and T k ( ξ T )is a non-parametric model for the Fourier transformed correla-tion 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 twofields enters the overall dust extinction density s . Because of thisproperty, the reconstructed power spectrum T k does not have tobe the empirical power spectrum of s x , that can be calculatedby Fourier transforming and binning. To avoid misunderstand-ings and artifacts from the degenerate model, we mainly report Article number, page 3 of 16 & A proofs: manuscript no. main the empirical power spectrum in this paper, which is computedfrom posterior samples of s x . We now focus on our model forthe power spectrum T k . This model assumes the spatial correla-tion power spectrum to be a preferentially falling power-law, butallows for arbitrary deviations. It can be written as (cid:112) T k ( ξ T ) = Exp ∗ Exp (cid:2) ( m s + σ s ξ s )ln( k ) + m + σ ξ + F ln( k ) t sym( A / (cid:16) + ( t / t ) (cid:17) ξ φ ( t )) (cid:3) , (17)where the first part describes a linear function on log-logscale, i.e. a power law; and the second part describes the non-parametric deviations which are assumed to be di ff erentiable onlog-log-scale. The operation Exp ∗ denotes the exponentiation ofthe coordinate system. More explicitly, F ln( k ) t is Fourier transfor-mation on log-log scale, and the function sym is defined as f : [0 , b ] → R (18)sym( f )( x ) = ( f ( x ) − f (2 b − x )) (cid:12)(cid:12)(cid:12) [0 , b ] , (19)where f (cid:12)(cid:12)(cid:12) M denotes the restriction of the domain of the function f to M . The function sym is required to deal with the periodicboundary conditions introduced by the Fourier transform. De-tails can be found in the appendix of Arras et al. (2019). Thehyper-parameters of the model are ( A , t , m s , σ s , m , σ ) whichwe chose to be (11 , . , − , , − ,
3) in complete analogy toLeike & Enßlin (2019).
5. Algorithm
We combine the prior and the likelihood into one genera-tive model of the data. We compute approximate posteriorsamples using Metric Gaussian Variational Inference (MGVI)(Knollmüller & Enßlin 2019). This variational approach alter-nates between drawing samples around the current estimate forthe latent parameters and optimizing the current estimate usingthe average gradient of the samples. The final set of samples isused to derive an uncertainty estimate on all our maps as well ason all derived quantities.For further parallelization, we split the problem into the 8octants. Each octant has size 410pc × × ≈
417 Million,exceeding those of our previous map by a factor of 30. The totalcomputation time was about 2 weeks of wall clock time, or about0 . ff erentiablevariance-preserving interpolation scheme. The details are de-scribed in Appendix A.One noteworthy point is that we cut away the outer 30pc dueto artifacts from periodic boundary conditions, resulting in a fi-nal map volume of 740 pc ×
740 pc ×
540 pc.
6. Results and Discussion
We were able to reconstruct the nearby dust clouds. Fig. 2 showsvarious maps produced from our result and their relative uncer-tainty. The maps show tendrils and filaments of dust on scalesas small as 2pc up to scales of several hundred parsecs, at whichthey become disconnected.All octants inferred similar logarithmic convolution kernels,as can be seen in Fig.3. These correlation kernels were computedby taking a slice out of the reconstructed Fourier transformedsquare root power spectrum.A comparison of the empirical power spectra of the 8 di ff er-ent octants can be found in Fig. 4. Most octants have very sim-ilar power spectra, only octant 3 deviates strongly. This octant,located at 180 < l ≤
270 and b > . ± .
015 at scales from 2pcto 100pc. For the logarithmic power, we report a spectral indexof 2 . ± .
022 at scales from 2 . .
005 e-folds per pc of extinction density. Someof our samples do not exclude the existence of nearby denseclouds, which raises the uncertainty in the corresponding direc-tions tremendously. Fig. 6 shows the distance to the densest dustclouds in all directions, as well as an uncertainty on that dis-tance estimate. Note that the uncertainty estimate is quite highon the boundaries of dust clouds, as the reconstruction is uncer-tain which voxel is densest along these lines of sight.
An implicit assumption of the algorithm is that the voxelsare smaller than the achievable resolution. Phrased in physi-cal terms, an increase in pixel resolution can be regarded as arenormalization and we need to reach the continuous limit, i.ethe limit of negligible discretization e ff ects, for the algorithm towork. This is a byproduct of the inference of the power spec-trum, if the achieved posterior resolution is of the order of theimposed voxel resolution, then the reconstruction changes dras-tically from one voxel to another and the extinction of sourcesbehind an a ff ected voxel also changes dramatically at the bound-ary. This sudden change in extinction is not compatible with afalling spatial correlation power law in Fourier space, thus thereconstructed power will fall less steep than the real one. Thiswould significantly hamper the ability of the algorithm to ex-trapolate between measurements. We avoid this behavior by sig-nificantly increasing resolution compared to our previous recon-struction Leike & Enßlin (2019). However, it is conceivable thatthe reconstruction would still benefit from increasing the amountof voxels. We recommend distrusting the smallest scales of ourreconstruction, only at scales of 2 pc or larger can the resultbe considered to be stable. This resolution limit was deducedfrom the reconstructed logarithmic correlation kernels as seen inFig.3. At this limit, the reconstructed correlation kernels of thelogarithmic dust extinction density have fallen to 10%. Article number, page 4 of 16. H. Leike et al.: Resolving nearby dust clouds
A comparison of our results to Leike & Enßlin (2019) can befound 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 reconstruc-tions. Fig. 10 shows the same column densities, but on a loga-rithmic 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 Anderset al. (2019), which provides more sources and tighter con-straints on the parallax and G -band extinction than the previ-ously used Gaia data. The new reconstruction has a volume of800 pc ×
800 pc × cube in Leike& Enßlin (2019). Furthermore, using a designated fortran rou-tine for the computation of the line of sight integrals lead to thenecessary speedup to handle the additional data constraints andmassively more degrees of freedom. Finally, in the new recon-struction the parallax error is propagated into the measurementerror, causing extinction values with stars of high parallax errorto be less informative. In Fig. 7 one can see dust column densitiesalong Galactic x , y , and z coordinates. Both dust maps agree onthe morphology of large dust clouds on large scales. However,the current dust map contains significantly more dust. Part of thereason is that the data we use in the reconstruction of this paperhas higher resolution and lower noise, allowing more dust to bereconstructed. We also believe the data used in Leike & Enßlin(2019) to be slightly biased to underestimating the amount ofdust, an e ff ect that accumulates in a reconstruction that usesmany data points. In contrast, the data used in this reconstruc-tion might have a tendency to overestimate the amount of dust,despite our e ff ort to calibrate the zero point (see Sec. 3.2).We furthermore reconstruct our correlation kernel nonpara-metrically, which should lead to an unbiased estimate of thepower spectrum. Fig. 12 shows power spectra of Leike & Enßlin(2019), the reconstruction of this paper, and of the reconstruc-tion 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 remark-ably consistent with Green et al. (2019), however at scales ofa few parsec the power spectrum of Green et al. (2019) flattenswhich we believe to be an artifact of how we put their reconstruc-tion on a cartesian grid, i.e. the boundaries of the reconstructionsintrinsic voxels introduce steep cuts which flatten the resultingpower spectrum. However, none of the power spectra are consis-tent within the uncertainties estimates. While this seems prob-lematic, one has to bear in mind that all reconstructions focuson dust in di ff ering regions, potentially explaining the di ff erencein the power spectrum. In Fig. 13 we show the power spectra ofthe logarithmic reconstructions. These seem to be less consistentin general, however one has to bear in mind that the logarithmicpower spectrum is dominated by regions of low dust content, asthese occur more frequently. Using Gaia data our method wasfound to underestimate low dust regions, and we anticipate thatwith the starhorse data we tend to overestimate low-dust regions.Nonetheless we find that the spectral index of 2 . ± .
022 atscales from 2 . σ joint uncer-tainty margin. The spectral index of the empirical power spec-trum of Leike & Enßlin (2019) is 3 . ± . The logarithmic Note that Leike & Enßlin (2019) reports a spectral index of 3 . ff erent values but enables us to derive uncertainty estimates for allcompared maps in the same way. power spectrum of Green et al. (2019) seems to be inconsistentwith our measurements. However, this e ff ect is probably due tohow we treat the missing values in that map, where a quarter ofthe sky was not measured. We have to set these values and everypossible choice will impact the derived power spectra. We choseto set them to 10 − , which has minimal impact on the powerspectrum on linear scale, but biases the power spectrum of thelogarithmic dust extinction density and could potentially explainthe di ff erence.Fig. 14 shows a histogram of dust extinction density pervoxel. One can see a good agreement between the histogramof our old and our current reconstruction in the region between10 − pc − and 10 − pc − . A dust extinction density of 10 − pc − integrated to the boundary of our simulation cube yields an inte-grated extinction of 0 . − pc − as its shape is mostlydependent on how the reconstruction extrapolates into dustlessregion. From the histogram it can be seen that the dust densityis well described by a log-normal distribution. Note that sincewe show only the part of the histogram that has high signal tonoise, this result should be relatively unbiased by our choice ofprior. The fitted log-normal model has a standard deviation of σ = . ± .
009 and a mean of m = − . ± . One should note that the reconstruction shows a non-negligibleamount of dust in the local bubble. We believe that the levelfound is an artifact of our noise statistic. As described in Sec. 3,our data model involves some heuristics which might systemati-cally a ff ect the reconstruction. This causes estimates made withthis data to be biased and we were not able to fully correct for thisbias. When integrating the reconstructed dust density to 70 pc,we find that the nearby dust looks like a smeared out version offarther dust clouds, indicating that it is indeed an artifact relatedto systematic data biases.The posterior samples of the extinction density are availablefor download under https://doi.org/10.5281/zenodo.3750926 or by its DOI 10.5281 / zenodo.3750926 . When usingthe reconstruction we advise to beware of systematic overesti-mations of dust, especially in the local bubble. When derivingnumeric quantities, we advise to do so for every sample and thenestimate the mean and standard deviation of the results in orderto get an error estimation. Our map can be used to constrain simulations of the ISM. Forexample, in simulations of radiatively cooling dust clouds in hotwinds, it has been shown that dust density power spectra are flat-ter than was previously thought (Sparre et al. 2019). Our mapsshow power spectra compatible with these simulations, and mor-phologically similar structures. Our reconstructed spectral indexof 2 . ± .
022 at scales from 2 . G -band dust extinction density in e-folds per parsec shownin Fig. 14 is well described by a log-normal distribution withstandard deviation σ = . ± .
009 and mean m = − . ± . − pc − to 1 pc − . Article number, page 5 of 16 & A proofs: manuscript no. main
7. Conclusion
We were able to reconstruct the dust clouds within ∼
400 pc ofthe Sun down to a resolution of 2 pc, improving in resolution andvolume on our previous reconstruction (Leike & Enßlin 2019).The resulting map is public and can be downloaded; see Sec. 6.3for details. Distances to and densities of all dust clouds largerthan 2 pc are expected to be well constrained by the reconstruc-tion. We report our estimate on the power spectrum of the dustextinction density as well as the logarithmic density. Further-more, we provide a histogram of dust densities in the interstel-lar medium and find them to be well described by a log-normalmodel. We hope that our diverse summary statistics allow simu-lations of the ISM to be constrained.
Acknowledgements.
We acknowledge fruitful discussions with S. Hutschen-reuter, J. Knollmüller, P. Arras, A. Kostic, and others from the information fieldtheory group at the MPI for Astrophysics, Garching. We acknowledge the sup-port by the DFG Cluster of Excellence "Origin and Structure of the Universe".The prototypes for the reconstructions were carried out on the computing facili-ties of the Computational Center for Particle and Astrophysics (C2PAP). We ac-knowledge support by the Max-Planck Computing and Data Facility (MPCDF).The main computation for the reconstructions were carried on compute clus-ter FREYA. This work has made use of data from the European Space Agency(ESA) mission
Gaia ( ), processed bythe Gaia
Data Processing and Analysis Consortium (DPAC, ). Funding for the DPAC hasbeen provided by national institutions, in particular the institutions participatingin the
Gaia
Multilateral Agreement.
References
Anders, F., Khalatyan, A., Chiappini, C., et al. 2019, Astronomy & Astrophysics,628, A94Arras, P., Frank, P., Leike, R., Westermann, R., & Enßlin, T. 2019, arXiv preprintarXiv:1903.11169Chen, B., Huang, Y., Yuan, H., et al. 2018, Monthly Notices of the Royal Astro-nomical Society, 483, 4277Green, G. M., Schlafly, E. F., Zucker, C., Speagle, J. S., & Finkbeiner, D. P.2019, arXiv e-prints, arXiv:1905.02734Knollmüller, J. & Enßlin, T. A. 2019, arXiv preprint arXiv:1901.11033Lallement, R., Babusiaux, C., Vergely, J., et al. 2019, Astronomy & Astro-physics, 625, A135Leike, R. & Enßlin, T. 2019, Astronomy & Astrophysics, 631, A32Panopoulou, G., Psaradaki, I., & Tassis, K. 2016, Monthly Notices of the RoyalAstronomical Society, 462, 1517Selig, M., Vacca, V., Oppermann, N., & Enßlin, T. A. 2015, Astronomy & As-trophysics, 581, A126Skalidis, R. & Pelgrims, V. 2019, Astronomy & Astrophysics, 631, L11Sparre, M., Pfrommer, C., & Vogelsberger, M. 2019, Monthly Notices of theRoyal Astronomical Society, 482, 5401The Planck Collaboration et al. 2018, arXiv preprint arXiv:1807.06205
Article number, page 6 of 16. H. Leike et al.: Resolving nearby dust clouds(a) (b)(c) (d)(e) (f)
Fig. 2: Result of our 3D dust reconstruction. The first column shows dust extinction, the second shows the relative error. The firstrow 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 ofradially. The third row shows di ff erential extinction in e-folds per parsec in a slice along the Galactic plane. Article number, page 7 of 16 & A proofs: manuscript no. main
Fig. 3: Reconstructed correlation kernels for the di ff erent oc-tants. Note that the logarithmic dust extinction in our model isthe result of an a-priori normal distributed field that is foldedwith these kernels, dependent on the octant. The octants are ar-ranged such that octant i = b + b + b (for b i ∈ { , } ) ex-tends in positive x -direction if and only if b =
0, in positive y -direction if and only if b = z -direction if andonly if b =
0. Note that all kernels fall to about 10% in the first2 pc.Fig. 4: Empirical power spectra of the dust extinction densityof the eight octants. The octants are arranged such that octant i = b + b + b (for b i ∈ { , } ) extends in positive x -directionif and only if b =
0, in positive y -direction if and only if b = z -direction if and only if b = Article number, page 8 of 16. H. Leike et al.: Resolving nearby dust clouds
Fig. 5: A Mollweide projection showing the distance to the first voxel of our reconstruction that exceeds an extinction estimateof 0 .
005 e-folds per parsec (left side) and corresponding uncertainty map (right side). Directions for which the threshold is neverreached appear in white.Fig. 6: A Mollweide projection showing the distance to the voxel with the highest extinction estimate in that direction (left side)and corresponding uncertainty map (right side).
Article number, page 9 of 16 & A proofs: manuscript no. main(a) (b)(c) (d)(e) (f)
Fig. 7: Comparison of column densities of our current reconstruction (left column) and Leike & Enßlin (2019) (right column). Therows show integrated dust extinction for sightlines parallel to the z - x - and y -axis respectively. Article number, page 10 of 16. H. Leike et al.: Resolving nearby dust clouds(a) (b)(c) (d)(e) (f)
Fig. 8: As Fig. 7 but on logarithmic scale.
Article number, page 11 of 16 & A proofs: manuscript no. main(a) (b)(c) (d)(e) (f)
Fig. 9: Column density comparison of our current reconstruction (left column) and that of Green et al. (2019) (right column). Therows show integrated dust extinction for sightlines parallel to the z - x - and y -axis respectively. Note that for Green et al. (2019) weshow the integrated extinction only if more than 50% of the projected voxels exist in the reconstruction. Article number, page 12 of 16. H. Leike et al.: Resolving nearby dust clouds(a) (b)(c) (d)(e) (f)
Fig. 10: As Fig. 9 but on logarithmic scale.
Article number, page 13 of 16 & A proofs: manuscript no. main(a) (b)(c) (d)(e) (f)
Fig. 11: Comparison of integrated extinction of our reconstruction (left column) and that of Green et al. (2019) (right column) insky projection. The rows show integrated dust extinction out to the boundary of our 740 pc ×
740 pc ×
540 pc box in an all sky view(first row), as well as two selected directions towards the Galactic anticenter (middle row) and center (last row).
Article number, page 14 of 16. H. Leike et al.: Resolving nearby dust clouds
Fig. 12: Empirical power spectra of the dust extinction densityof this paper (solid line), Leike & Enßlin (2019) (dashed line)and the reconstruction of Green et al. (2019) (dotted line).Fig. 13: Empirical power spectra of the logarithmic dust ex-tinction density of this paper (solid line), Leike & Enßlin (2019)(dashed line) and the reconstruction of Green et al. (2019) (dot-ted line). Fig. 14: Histogram of dust extinction density per voxel of thispaper (solid line), Leike & Enßlin (2019) (dashed line) and thereconstruction of Green et al. (2019) (dotted line). We overplota log-normal model that was fitted to our reconstructed logarith-mic extinction density (dash-dotted line). The curve of the fit isdescribed by f ( x ) ∝ exp (cid:16) − . σ − (ln( x ) − m ) (cid:17) with σ = . m = − .
79 and follows our empirical distribution functionclosely.
Article number, page 15 of 16 & A proofs: manuscript no. main
Appendix A: interpolation scheme
To parallelize the reconstruction, we reconstructed the eight oc-tants of the coordinate system independently, with a 20 pc over-lap region. To get one final reconstruction, we have to glue thesereconstructions together and specify how we deal with the over-lap region. We do so using a di ff erentiable, variance preservinginterpolation scheme, i.e. if the octants are di ff erentiable thenthe result is di ff erentiable; and the final samples have at least thevariance that the individual reconstructions imposed. We com-pute the uncorrected interpolated logarithmic extinction samples τ (cid:48) ( x ) j from the samples of the eight octants o i ( x ) j as τ (cid:48) ( x ) j = (cid:88) i w i ( x ) o i ( x ) j . (A.1)Hereby the weights w i ( x ) can be computed as w i ( x ) = (cid:89) k = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b k ( i ) − f (cid:32) x k − (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (A.2)where f ( x ) = x ∈ ( −∞ , x (3 − x ) x ∈ (0 , b k ( i ) denotes the k -th digit of i in binary format. Noteworthyproperties of this scheme are that the weights sum to one ∀ x (cid:88) i w i ( x ) = , (A.4)and the polynomial f is the unique polynomial of degree 3 suchthat f (0) = , (A.5) f (1) = , (A.6) ∂ f ∂ x (0) = , (A.7) ∂ f ∂ x (1) = . (A.8)By using this interpolation scheme only voxels that have a co-ordinate x k of which the absolute value is at most 8 pc get non-zero contributions from more than one octant, or put in a dif-ferent way: We cut away the outermost 2 pc of all reconstruc-tions, which mitigates artifacts from periodic boundary condi-tions. From the preliminary logarithmic extinction density τ (cid:48) ( x ) j we can compute the logarithmic sample mean ¯ τ ( x ) as¯ τ = N (cid:88) j τ (cid:48) j . (A.9)Here N denotes the number of samples. The variance of τ (cid:48) j isartificially low at overlapping regions, as independent sampleswere averaged. We correct for this e ff ect and compute the overalllogarithmic extinction density samples τ ( x ) j as τ ( x ) j = τ (cid:48) ( x ) j − ¯ τ ( x ) (cid:112)(cid:80) i w i ( x ) + ¯ τ ( x ) . (A.10)(A.10)