A 3D Dust Map Based on Gaia, PanSTARRS 1 and 2MASS
Abstract
We present a new threedimensional map of dust reddening, based on Gaia parallaxes and stellar photometry from PanSTARRS 1 and 2MASS. This map covers the sky north of a declination of , out to a distance of several kiloparsecs. This new map contains three major improvements over our previous work. First, the inclusion of Gaia parallaxes dramatically improves distance estimates to nearby stars. Second, we incorporate a spatial prior that correlates the dust density across nearby sightlines. This produces a smoother map, with more isotropic clouds and smaller distance uncertainties, particularly to clouds within the nearest kiloparsec. Third, we infer the dust density with a distance resolution that is four times finer than in our previous work, to accommodate the improvements in signaltonoise enabled by the other improvements. As part of this work, we infer the distances, reddenings and types of 799 million stars. We obtain typical reddening uncertainties that are 30% smaller than those reported in the Gaia DR2 catalog, reflecting the greater number of photometric passbands that enter into our analysis. Our 3D dust map can be accessed at doi:10.7910/DVN/2EJ9TX or through the Python package dustmaps. Our catalog of stellar parameters can be accessed at doi:10.7910/DVN/AV9GXO.
0000000154172260]Gregory M. Green
\move@AU\move@AF\@affiliationKavli Institute for Particle Astrophysics and Cosmology,
Stanford University
452 Lomita Mall,
Stanford, CA 943054060, USA
0000000235697421]Edward Schlafly
\move@AU\move@AF\@affiliationLawrence Berkeley National Laboratory
One Cyclotron Road
Berkeley, CA 94720, USA
\move@AU\move@AF\@affiliationHubble Fellow
000000022250730X]Catherine Zucker
\move@AU\move@AF\@affiliationHarvard Astronomy,
HarvardSmithsonian Center for Astrophysics
60 Garden St.,
Cambridge, MA 02138, USA
0000000250659896]Joshua S. Speagle
\move@AU\move@AF\@affiliationHarvard Astronomy,
HarvardSmithsonian Center for Astrophysics
60 Garden St.,
Cambridge, MA 02138, USA
000000032808275X]Douglas Finkbeiner
\move@AU\move@AF\@affiliationHarvard Astronomy,
HarvardSmithsonian Center for Astrophysics
60 Garden St.,
Cambridge, MA 02138, USA
1 Introduction
Dust is both a critical foreground for many astrophysical measurements and a tracer of Galactic structure and starforming regions. In the ultraviolet, optical and nearinfrared, dust causes both extinction and reddening. It is necessary to correct for these effects in order to measure intrinsic luminosities or colors of obscured objects. For extragalactic astronomy, a twodimensional map of integrated dust extinction and reddening is sufficient, while for sources embedded in the Milky Way, a distancedependent correction for dust obscuration is necessary. Dust emits strongly in the mid and farinfrared, and thus provides an important foreground for the Cosmic Microwave Background, further motivating accurate maps of the distribution of Galactic dust.
Farinfrared emission of dust can be used to construct twodimensional maps of dust reddening. SFD98 models farinfrared dust emission as a modified blackbody, deriving maps of dust optical depth and temperature. Assuming a linear relation between dust optical depth at 100 m and reddening at optical wavelengths, SFD calibrates their optical depth map against measurements of from elliptical galaxies, obtaining a fullsky map of Galactic dust reddening at optical wavelengths. Because it is not possible to determine the distance to the emitting dust from farinfrared intensity alone, this method does not recover the threedimensional distribution of dust. In addition, the relation between farinfrared optical depth and optical reddening depends on the dust composition and grainsize distribution, and thus varies throughout the Galaxy, introducing systematic errors into the resulting reddening maps.
One method of overcoming these limitations is to use optical stellar photometry to trace dust. Stellar distances can be determined (both from spectral models and from geometric parallax measurements), allowing one to trace reddening as a function of both position on the sky and distance. In addition, stellar colors are a more direct measurement of optical reddening, obviating the need to convert from farinfrared optical depth to reddening at optical wavelengths. Several works have now used optical and nearinfrared stellar photometry to trace dust in 3D. In lieu of an extensive summary, a few works will be mentioned here. Marshall2006 uses nearinfrared colors of postmainsequence stars in the Two Micron All Sky Survey (Skrutskie2006, hereafter “2MASS”) to map dust in 3D, while Berry2012 uses the optical and nearinfrared colors of stars in the Sloan Digital Sky Survey (York2000).
These methods have matured with the introduction of probabilistic models that infer both the distribution of dust and the types and distances of stars. Sale2014IPHAS applies a hierarchical Bayesian model of the distribution of dust, as well as stellar distances and types, to the INT Photometric H Survey of the Northern Galactic Plane (Drew2005IPHAS, “IPHAS”). Green2015 models PanSTARRS 1 and 2MASS photometry of 800 million stars over threequarters of the sky to create a map that extends to several kiloparsecs. Green2018 improves upon this map, with both updated PanSTARRS 1 photometry and methodological improvements. Lallement2014 uses catalog distances and reddenings of stars to trace dust, but applies a Gaussian process prior to the spatial distribution of dust, enforcing smoothness on small spatial scales.
Recently, the availability of geometric parallax measurements from Gaia (GaiaCollaboration2016Mission) has given a significant boost to the field. Whereas most previous studies of the 3D distribution of dust have relied on stellar spectral energy distribution models to determine distances to stars, Gaia more directly constrains the distances to the hundreds of millions of stars for which it precisely measures parallaxes. In addition, Gaia measures stellar optical photometry, which Andrae2018DR2Apsis uses to estimate stellar reddenings and extinctions. Using a method similar to that of Andrae2018DR2Apsis, Chen2018 estimates stellar reddenings from Gaia photometry and parallaxes, producing an allsky map of dust reddening in 3D. LeikeEnsslin2019 combines the extinction estimates from Andrae2018DR2Apsis with a Gaussian process model for the spatial distribution of dust to map reddening within a few hundred parsecs of the Sun. Lallement2019 combines Gaia parallaxes with photometry from Gaia and 2MASS to estimate stellar distances and reddenings, and models the spatial distribution of dust as a Gaussian process, obtaining a map that extends to several kiloparsecs in the midplane of the Galaxy.
In the following, we present a new 3D map of dust reddening, based on Gaia parallaxes and broadband optical and nearinfrared photometry. We make several major improvements over Green2018, incorporating Gaia parallaxes into our model, improving the distance resolution of the map by a factor of four, and imposing a Gaussian process prior on the distribution of dust. For stars with precisely measured Gaia parallaxes, we obtain dramatically more precise distance determinations than is possible with broadband photometry alone. The imposition of a Gaussian process prior on the distribution of dust not only regularizes the resulting map, but also allows us to more precisely determine distances to dust clouds within 1 kpc.
Our method produces not only a map of the 3D distribution of dust reddening, but also a catalog of distances, reddenings and types for 799 million stars. Our 3D dust map may be accessed at doi:10.7910/DVN/2EJ9TX, while our catalog of stellar parameters may be accessed at doi:10.7910/DVN/AV9GXO. In addition, our dust map can be accessed through the Python package dustmaps (Green2018dustmaps).
The paper is organized as follows. Section 2 reviews our method for inferring stellar properties (Section 2.1) and the distribution of dust along each sightline (Section 2.2), and then describes a new method for imposing a correlated prior on the dust, linking nearby sightlines (Section 2.3). This latter process yields a map in which dust densities in nearby voxels, spanning different sightlines, are correlated, dependent on the distance between the voxels. Section 3 describes the survey data we use, from Gaia, PanSTARRS 1 and 2MASS. Section 4 presents our new 3D dust map, compares it to previous dust maps, assesses evidence in our map for spiral structure in the Milky Way, and presents our catalog of stellar parameters. Section 5 discusses ways forward for dust mapping. Finally, we conclude in Section 7.
2 Method
Here, we sketch out our method for determining the threedimensional distribution of dust, as well as stellar types and distances. In brief, we use stars as tracers of the dust column, modeling the apparent magnitudes and parallaxes of stars as a function of their type, distance and foreground reddening. We then group the stars into discrete sightlines, using the stellar distances and reddenings to constrain the dust as a function of distance.
The method for inferring types and distances for single stars is much the same as in Green2015 and Green2018, with the addition of Gaia parallaxes. However, we have significantly updated our method for sampling the dust reddening along individual sightlines, allowing us to increase the number of distance bins by a factor of four while decreasing the computational runtime, and we have implemented a new iterative scheme for imposing a prior with correlations between nearby sightlines.
2.1 SingleStar Inference
For an individual star with measured photometry, , and parallax, , we wish to determine the star’s type, , distance modulus, , and reddening, . We can write down the posterior probability density of these parameters, given the data, as
(1) 
Since the photometry and parallax are observed, fixed quantities, we can treat as a constant, so
(2) 
If we are given the type, distance modulus and reddening of a star, we can look up its model magnitudes, , using the stellar templates described in Green2015, and predict its parallax, which is simply a function of distance modulus, . Comparison of these model observables with the observed photometry and parallax yields the likelihood term, , which is a multivariate Gaussian:
(3)  
(4) 
where and are the uncertainty in the observed magnitudes and parallax, and is the extinction in the observed passbands, assumed to be a linear function of reddening, . We assume that for each star, extinction is given by
(5) 
where is the “extinction vector,” relating a scalar reddening to the extinction in each passband. We use the empirical reddening vector from Schlafly2016ExtinctionVector, converting it to an extinction vector by requiring that (Indebetouw2005) and that when . The latter choice puts our measure of reddening on a similar scale as SFD (SchlaflyFinkbeiner2011). The resulting extinction vector is given in Table 2.1. This extinction vector is slightly different from the one we use in Green2018, as we normalize the vector in that work by setting extinction in the WISE W2 band to zero. The difference between the two extinction vectors is most pronounced in the 2MASS passbands.
We place priors on the distribution of stars of different types throughout the Galaxy, yielding a joint prior on distance modulus and type, . These priors are described in Green2014, and are summarized in Table 2.1. For reasons that will become clear in Section 2.2, we place a flat prior on the reddening of each star: .
We evaluate the posterior density of each star, , on a grid of distance modulus in steps of 0.125 mag and in steps of 0.01 mag. Note that stellar type, , has been integrated out here. Our procedure for evaluating this posterior density quickly over a regular grid is explained in Appendix A. This new method is both faster and more accurate than the evaluation technique used in our previous work, which involved kernel density estimation on samples drawn from an MCMC chain.
As in Green2015 and Green2018, we smooth the stellar probability density functions along the reddening axis with a Gaussian kernel, which is equivalent to treating reddening within each pixel as a white noise process at any given distance. This takes into account the smallscale power in the dust density spectrum, allowing reddening to vary across a single angular pixel. The Gaussian smoothing applied is a fraction of the reddening in the pixel, as explained in Section 2.2 of Green2015. We make one change to the smoothing method developed in Green2015, in that we increase the minimum reddening smoothing percentage to 15% and the maximum smoothing percentage to 50%.
Not all of the point sources that we model are indeed stars – our sample is inevitably contaminated by compact galaxies and quasars. We attempt to filter out galaxies by enforcing a compactness criterion on point sources (see Section 3), but at the faint end of our source catalog, this criterion becomes less effective, both because morphology of sources becomes more difficult to determine and because galaxies become more numerous in relation to stars. There are additionally certain stellar types, such as O and B stars, which are not included in our stellar template library. In order to filter out such sources, we make a goodnessoffit cut after modeling the stellar parameters. In detail, we filter out stars for which the minimum (i.e., for which no stellar template accurately models the observed photometry). For the purposes of this cut, we count stellar parallax as an additional photometric band. In previous work, we cut on the Bayesian evidence of the model for each individual star. Our new cut is similar in spirit, but does not depend on the choice of priors, relying instead only on the likelihood. Approximately 1.4% of point sources fail this goodnessoffit cut, and are therefore not used to determine lineofsight reddening.
2.2 SingleSightline Inference
In this section, we will describe how to use the stars within one individual pixel to infer the dust density as a function of distance. The way in which we do this for one isolated sightline is described in detail in Green2018, but we will recap it here.
We split the sightline up into discrete distance bins, and parameterize the logarithm of the increase in reddening in the distance bins by . We want to know what distribution of dust along the line of sight and what choice of stellar types and distances is consistent with the observed stellar magnitudes and parallaxes. The posterior on the lineofsight dust distribution is given by
(6) 
where denotes the photometry and parallaxes for all the stars in the sightline, and is a constant for any given observed data. Assuming that the stars are independent, we get
(7) 
The term encodes our prior expectations about the amount of dust in each distance bin along the sightline. As in Green2018 and Green2015, we put a lognormal prior on the jump in reddening in each distance bin, based on a simple, smooth model of the distribution of dust in the Milky Way. In Section 2.3, we describe how to impose a prior that links the dust in nearby sightlines.
We have a term for each star . In what follows, we’ll drop the subscript , and it is assumed we are discussing a single star. We can expand this term as
(8)  
(9) 
Note the similarity between the term in the integrand and the righthand side of Eq. (2). In fact, these two terms are exactly identical. The observed stellar photometry and parallax are dependent only on the stellar type, distance and reddening. Given a lineofsight distribution of reddening, , and stellar distance modulus, , we can calculate the stellar reddening, . Thus,
(10) 
where . Recall from Section 2.1 that we placed a uniform prior on the stellar reddening, . As a consequence,
(11) 
We define the function
(12) 
which is precomputed for each star on a grid in and , according to the method laid out in Appendix A. In terms of this function, Eq. (9) becomes
(13) 
Plugging this all into Eq. (7), we obtain
(14) 
Conceptually, the above means that the posterior on the lineofsight reddening is a product of the prior and a line integral tracing the distancereddening curve through each individual stellar posterior. This expression is generally quick to compute, once the precomputed stellar posteriors, , are available.
We model the lineofsight reddening as a step function, with encoding the increase in reddening in each distance bin. The jumps in reddening are discretized, as integer multiples of 0.01 mag. The distance bins are spaced evenly in distance modulus, from to 19 mag, with a bin spacing of 0.125 mag. Relative to our previous work (Green2014; Green2015; Green2018), we find that our discretized parameterization of reddening significantly speeds up the evaluation of the line integrals in Eq. (14), as well as the convergence of our MCMC sampler.
2.3 Correlated prior
In the previous section, we showed how to determine the distancereddening relation along one isolated sightline. However, 3D dust maps created sightlinebysightline display a number of undesirable and unphysical features.
Because distance is fundamentally more difficult to determine than angular position on the sky, dust clouds in 3D dust maps constructed from individual sightlines are typically highly elongated in the radial direction. This is an unphysical feature, which violates the Copernican principle. We might reasonably expect, a priori, for dust clouds to be oriented randomly with respect to the Sun (or at least, to not be oriented in some special configuration with respect to the Sun).
Our technique for measuring the distribution of dust along the line of sight depends on bracketing the dust with foreground and background stars. Because of the stochastic nature of the distribution of stars throughout the Galaxy, a particular dust cloud may be tightly constrained in distance by foreground and background stars in one sightline, but poorly constrained in distance in a neighboring sightline. This leads to 3D dust maps in which the precision of our distance estimate for the same cloud can vary significantly across different sightlines. A priori, we expect the dust density field to be smooth on small scales, so that if a dust cloud appears at some distance in one sightline, it is likely to appear at a similar distance in neighboring sightlines.
In order to solve these problems, several works have suggested placing Gaussian process priors on either the dust reddening density (i.e., the derivative of reddening with distance) (Rezaei2017Method), or on the logarithm of the dust reddening density (Lallement2014; LeikeEnsslin2019). There are advantages and disadvantages to each approach. If dust reddening density is a Gaussian process, then line integrals through the reddening density are themselves Gaussian, meaning that the integrated reddenings to observed stars are a Gaussian process. This feature is exploited by Rezaei2018Spirals to avoid having to explicitly model the dust reddening density throughout space, but rather only the integrated reddenings to observed stars. However, dust reddening density cannot be negative and varies over several orders of magnitude, from the diffuse interstellar medium to dense molecular cores. For these reasons, modeling the logarithm of the dust density as a Gaussian process is more physically natural, and is the approach taken in this paper. Both approaches encode the expected smoothness of the dust density field, and will tend to produce maps with round clouds, not stretched along the line of sight. The true interstellar medium is not a Gaussian process, as it contains filamentary structure that is not fully described by a twopoint correlation function, yet one may reasonably expect that a Gaussian process prior that correlates dust density in nearby sightlines will produce more physical results than a completely uncorrelated prior in which neighboring sightlines are independent of one another. A third possibility which is worth mentioning is to model the logarithm of the integrated extinction to each point in space as a Gaussian process (SaleMagorrian2014GP). This has the advantage of being less computationally taxing than modeling the logarithm of the dust extinction density as a Gaussian process, but is somewhat less physical, as it does not encode our knowledge that integrated extinction increases with distance, and even allows negative extinction density.
In Green2015 and Green2018, we modeled the logarithm of the reddening density in each voxel (one distance bin in one sightline) as a Gaussian, with a mean given by the smooth disk component of the model of the Galaxy from DrimmelSpergel2001, and a fixed variance in each bin. In this paper, we alter this prior by introducing covariances between the voxels, including covariances between voxels in different sightlines:
(15) 
where is the logarithm of the reddening density in a single voxel, is the total number of voxels (in multiple sightlines), is the expectation of the logarithm of the reddening density in voxel , given by a smooth Galactic model (the same as used in Green2018), and is a covariance matrix. The covariance between any two voxels is a function of the distance between their centers – this function is called the “kernel.” The kernel contains information about the power spectrum of the interstellar medium, and in principle could be derived from physical considerations (see, e.g., SaleMagorrian2014GP, which assumes Kolmogorov turbulence in the interstellar medium), or parameterized and inferred itself (as in LeikeEnsslin2019). In this work, we choose instead to fix the kernel, which reduces the computational complexity of the problem.
In Section 2.2, we worked out the posterior density on reddening in a single sightline. With our new Gaussian process prior, the joint posterior density of reddening in many sightlines is similar to Eq. (14):
(16) 
The joint prior, , is the Gaussian process prior over the logarithm of the reddening density in all the voxels in the entire volume under consideration, which spans multiple sightlines.
2.3.1 Covariance kernel
In this section, we describe the kernel, which sets the covariance between the logarithm of the reddening density in any two voxels. We choose a modified exponential kernel, for which covariance falls exponentially at large distances, but which smoothly asymptotes to a maximum value within distances . This effective inner radius limits the maximum covariance between any two voxels, which helps to make the problem more computationally tractable. Let and be the logarithm of the dust density in two different voxels, whose centers lie at and , respectively, and be the physical distance between the centers of the voxels. Our kernel is given by
(17) 
where is a distance scale over which the correlation between voxels decays, is the asymptotic value of the correlation coefficient for small physical separations, is the inner radius at which the correlation begins to smoothly asymptotes to , and defines how smoothly the cap is imposed. We set the variance in the logarithmic reddening density (i.e., the diagonal of the covariance matrix) to . This parameter thus sets the overall scale of variation in reddening from the smooth Galactic model. For (zero distance, representing the variance in a single voxel), the kernel is defined to be . We set , which yields a prior with approximately the same mean and variance in the expected reddening per unit distance as the prior chosen in Green2018, taking into account the increased distance resolution of our new map.
For the remaining kernel parameters, we choose , , and . The scale length, , is chosen to be similar to the transverse distance between neighboring sightlines, with a typical angular separation of , at a radial distance of 1 kpc, which comes to . The inner radius, , is set so that the prior is not dominated by correlations between the nearest distance bins, which are separated by smaller transverse distances. Although in theory, large correlations between nearby voxels are not a problem, they can pose difficulties for the iterative approach we develop in Section 2.3.3.
We give a summary of our chosen kernel parameters in Table 2.3.4.
2.3.2 Importance sampling
Our goal is to impose a correlated prior that links dust reddening densities across different sightlines. In order to do this, one naively has to infer the dust reddening density in all voxels simultaneously. A naive approach would involve calculating the covariance matrix, , between all of the voxels in the map, inverting it, and then using MCMC to sample the posterior density (Eq. 16) of the logarithmic reddening densities in all the voxels. With millions of sightlines, each with 100 distance bins, this would require us to invert an enormous covariance matrix, and then to sample from a space with hundreds of millions of parameters, a daunting task.
One way of reducing this problem to a feasible one, which we employ in this work, is to infer the dust reddening density only on restricted patches of sky, rather than over the entire sky at once. This limits the size of the covariance matrix which must be calculated and inverted, and also decreases the dimensionality of the parameter space which must be sampled. Rather than sampling from the reddening density over all voxels in the entire map simultaneously, we instead analyze small patches of sky independently. A second technique one might try is importance sampling, which we will flesh out here.
In importance sampling, one wishes to sample from a distribution , which might be expensive to compute. One begins instead by drawing points from a distribution which is easier to sample, which we will call . One then reweights each sample by the ratio
(18) 
termed the “importance weight.” The weighted samples, , are then representative of the distribution . Importance sampling works well when the distribution is a good approximation to the true distribution, . In such cases, the weights, , will be close to unity. If is a poor approximation to , then the samples, , will not be concentrated in the regions of high probability density and will therefore receive nearzero weights, while the few samples which happen to lie in highprobabilitydensity regions will receive large weights.
In our case, we have a Gaussian process prior on the logarithm of the reddening density in each voxel , :
(19) 
where is the number of voxels. If the correlations between voxels are small, we can approximate this distribution as a product of independent distributions:
(20) 
We can sample from this uncorrelated, approximate distribution, and then reweight the resulting samples using the ratio
(21) 
This method is only effective if the correlations between voxels are small, making the uncorrelated case a good approximation to the correlated case. However, in our case, the uncorrelated prior is sufficiently different from our desired prior that importance sampling would be inefficient.
2.3.3 An iterative scheme for imposing the correlated prior
Our central idea here will be to first sample from the uncorrelated case, and then over the course of several iterations, to resample the individual pixels in a way that progressively yields a better approximation to our desired Gaussian process prior.
Our iterative scheme begins by sampling from each sightline individually. For each sightline, , we obtain samples of , the logarithm of the dust reddening density in all the distance bins along the line of sight. We then iteratively update each sightline.
In order to update one sightline, we select a small patch of sky surrounding it. We will call the reddening density in the central sightline , and the reddening density in the neighboring sightlines . We resample the entire patch of sky, but treat the central sightline differently from the neighbors. The reddening profile in each neighbor is represented by a sample from the previous iteration, while the reddenings along central sightline are allowed to take on any value. The neighboring sightlines need only be resampled in order to impose the correct prior on the central sightline, and we only store samples of the reddening density in the central sightline.
We begin by initializing the reddening density in each neighboring sightline to a random sample from the previous iteration, and set the reddening density in the central sightline to an initial guess. We then alternate between updating the central and the neighboring sightlines.
To update one sightline, holding all the others fixed, we need to sample from the conditional distribution
(22) 
where is the sightline we are updating, and denotes the set of sightlines we are holding fixed. One useful property of a multivariate Gaussian is that the conditional probability densities are also Gaussian. If we fix the sightlines , the conditional prior on sightline , , is itself a Gaussian.
To update the central sightline, we hold all the neighboring sightlines fixed, and sample from the conditional distribution on the reddening in the central sightline. We use Markov Chain Monte Carlo sampling to explore the space of reddening densities along the central sightline, holding the neighbors fixed.
To update each neighboring sightline, we hold all the other neighbors and the central sightline fixed. We choose a new sample for the neighboring sightline, from the set of samples stored in the previous iteration. We wish the sample we draw to be distributed according to Eq. (22), with representing the neighbor to be updated, and representing all the other sightlines. However, the stored samples of from the previous iteration are drawn from the posterior on in the previous iteration, which we will denote . We therefore have to weight each stored sample according to the ratio of the conditional probability density of that sample, Eq. (22), to the posterior of that sample from the previous iteration:
(23) 
This is essentially an importancesampling step, with the conditional probability density (Eq. 22) playing the role of , and the stored posterior density of the sample from the previous iteration playing the role of in Eq. (18). The details of how we calculate the ratio of these probability densities is given in Appendix B.
The importancesampling step for the neighboring sightlines works well as long as the samples stored from the previous iteration are a decent representation of the conditional probability density of each neighbor, holding all the other neighbors and the central sightline fixed. If we adiabatically change the correlated prior from iteration to the next, then this is generally the case.
In each iteration, we jointly sample all the sightlines in a small patch surrounding the central sightline, alternating between updating the central sightline and updating each neighbor. We assess convergence by monitoring the autocorrelation of the reddening in the central sightline, as detailed in Appendix D. In order to speed up convergence, we use a technique that is similar to parallel tempering, with an ensemble of samplers that explore modified versions of the joint posterior on reddening density (see Appendix C).
2.3.4 Implementation details
In order to impose the Gaussian process prior, we choose the 32 nearest neighbors of each sightline. Each sightline in our map is thus sampled from a slightly different model, in which it is correlated with a patch of sky surrounding it. In each iteration, we sample the patch, as described above, and then discard the information from the neighboring sightlines in the patch, retaining only the Markov chain samples from the central sightline. After an initial pass over the sky, in which we sample the reddening along each sightline independently, we conduct four iterations with a correlated prior, increasing the correlation length with each iteration (see Table 2.3.4).
As described above, as we sample the reddening field in a patch of sky, we alternate between updating the central sightline and updating the neighbors. We call one such set of updates a “round.” In each round, we make 24000 MetropolisHastings proposals in the central sightline (equivalent to 20 proposals per distance bin), and five Gibbs steps in each of the neighboring sightlines. In the initial, uncorrelated iteration, we sample for 2500 rounds, discarding the first 500 as burnin. In subsequent, correlated iterations, we sample for 2250 rounds, discarding the first 450 as burnin. We save only the results from the central sightline, storing 100 samples.
The typical runtime per sightline is a function both of the number of stars in the sightline, as well as the number of neighboring sightlines. The runtime contains a term that is linear in the number of stars, due to the need to calculate one line integral per star in the lineofsight reddening posterior, Eq. (14). For reasons that may have to do with the need to store gridded stellar posterior densities in memory (specifically, with the number that can be stored in cache), there is also a small quadratic dependence of runtime on the number of stars. The runtime has a term that scales quadratically with the number of sightlines, due to the need to sample each neighboring sightline, and that the prior on each sightline depends on all of the neighbors. With 32 neighbors, running on a single core of an Intel Xeon ES2683v4, with a clock speed of 2.1 GHz and a cache size of 40 MB, we achieve a typical persightline runtime of
(24) 
where is the number of stars in the sightline. In the initial iteration, in which no neighboring pixels must be sampled, the typical runtime has the same dependence on the number of stars, but lacks the constant term. This runtime is dominated by the lineofsight sampler, as the time required for the grid evaluation of the individual stellar posterior densities is negligible, in comparison. The median number of stars per sightline is 130, yielding a typical runtime of for the initial, uncorrelated iteration, and for the correlated iterations.
Running on a mix of Intel and AMD cores on Harvard’s Odyssey cluster, the initial, uncorrelated iteration consumed approximately 700 million CPU seconds, while each successive iteration consumed approximately 1.5 billion CPU seconds, for a total runtime of 7 billion CPU seconds to generate the full map.
The method presented here makes a number of approximations to the Gaussian process model. Firstly, it only takes into account correlations between voxels which are at close angular separations on the sky, as seen from the Solar system, thereby giving the Sun a somewhat special position in our priors. Secondly, our voxels are highly elongated along the radial direction, and we calculate the distances between voxels using the distances between their centers. For this reason, distances between voxels are much smaller in the transverse directions than in the radial direction, and the correlations are correspondingly stronger along the transverse directions. We treat the correlations along the radial direction heuristically, through a modification of the inverse covariance matrix, as explained in Appendix C. Nevertheless, this method provides a relatively computationally inexpensive way to approximately apply a Gaussian Process prior to our map, and as will be seen in the following, provides a significantly improved result.
3 Data
3.1 Gaia
The European Space Agency’s Gaia spacecraft aims to map the 3D positions and kinematics – as well as spectral types – of a substantial fraction of Milky Way stars (GaiaCollaboration2016Mission). The mission will eventually provide astrometry (positions, proper motions and parallaxes) and optical spectrophotometry for over a billion stars (Liu2012Spectrophotometry), as well as radial velocity measurements of more than 100 million stars (Katz2004RVS; Wilkinson2005RVS). Gaia DR2 contains fiveparameter astrometric determinations for 1.3 billion sources (Lindegren2018DR2Astrometry).
Gaia’s spectrophotometry is an unusual middle ground between manyband photometry and lowresolution spectroscopy. Gaia takes dispersed images of sources in two bands, RP (330–680 nm) and BP (630–1050 nm). Gaia Data Release 2 (GaiaCollaboration2018DR2Contents) provides integrated RP and BP photometry of 1.4 billion sources is provided, which Andrae2018DR2Apsis combines with parallax and band photometry to estimate effective temperatures for 161 million stars and reddenings for 88 million stars. In future data releases, lowresolution RP and BP spectra will allow better characterization of the parameters of bright stars.
For our work here, the most valuable dataset provided by Gaia DR2 consists of the stellar parallax measurements. Parallaxes provide an independent means of estimating stellar distances from template fitting of spectral energy distributions. Particularly for bright and for nearby stars, Gaia parallaxes can much better constrain stellar distance than photometry alone. Stellar distancereddening posterior densities based on photometry are often bimodal, as mainsequence and giant stars of the same temperature can lie at vastly different distances. Even relatively uncertain parallaxes can be sufficient to break this degeneracy (see Fig. 1 of Zucker2019MolecularCloudsDR2).
We make two quality cuts on the Gaia DR2 catalog recommended by Arenou2018CatalogValidation:
(25)  
(26) 
where is related to the Gaia band magnitude, , by
(27) 
The first of these cuts removes sources with an insufficient number of Gaia observations to safely make a parallax measurement, while the second cut removes sources which are poorly fit by the astrometric model (and is essentially a cut on reduced ).
3.2 PanSTARRS 1
The bulk of our stellar photometry comes from the Panoramic Survey Telescope and Rapid Response System 1 (PanSTARRS 1, hereafter abbreviated as “PS1”), a 1.8 m optical and nearinfrared telescope located on Mount Haleakala, Hawaii (Chambers2016PS1). The telescope is equipped with the Gigapixel Camera #1 (GPC1), consisting of an array of 60 CCD detectors, each 4800 pixels on a side (Tonry2006GPC1; Onaka2008GPC1; Chambers2016PS1). From May 2010 to April 2014, the majority of the observing time was dedicated to a multiepoch survey of the sky north of declination (Chambers2016PS1, named the “ survey” for its footprint in steradians). The survey observed in four SDSSlike passbands (York2000), , , and , as well as an additional passband in the nearinfrared, . The entire filter set spans the range 400–1000 nm (Stubbs2010PS1lasercal).
Astrometry and photometry was extracted by the PS1 Image Processing Pipeline (Magnier2016PS1dataprocessing; Magnier2016PS1pixelanalysis; Magnier2016PS1photoastrocalib). PS1 photometry features a very uniform flux calibration, achieving better than 1% precision over the sky (Schlafly2012; Chambers2016PS1). In singleepoch photometry, the survey reaches typical depths of 22.0 mag (AB) in , 21.8 mag in , 21.5 mag in , 20.9 mag in , and 19.7 mag in (Chambers2016PS1). The PS1 Data Release 1 (DR1) occurred in December 2016, and provided a staticsky catalog, stacked images from the survey, along with other data products (Flewelling2016PS1database).
Because of its large footprint, homogeneous depth and excellent internal calibration, PS1 photometry provides an ideal dataset for mapping Galactic dust through stellar colors. As in Green2018, we use stellar photometry derived from the singleepoch photometry, which is very similar to PS1 DR1, with the primary difference being in the treatment of observations taken in nonphotometric conditions.
3.3 TwoMicron AllSky Survey
The Two Micron All Sky Survey (2MASS) is an allsky survey in three nearinfrared passbands, , and (Skrutskie2006), conducted from two 1.3 m telescopes at Mount Hopkins, Arizona and Cerro Tololo, Chile. The name of the survey derives from the wavelength range covered by the reddest passband, . Groundbased surveys beyond the atmospheric window at are hampered by severe background thermal emission (Skrutskie2006). The focal plane of each telescope was equipped with three pixel arrays, with a pixel scale. The entire sky was covered six times with dual 51millisecond and 1.3second exposures, in order to cover a wide range in apparent magnitudes. The survey achieved a 10 pointsource depth of 15.8, 15.1 and 14.3 mag (Vega) in , and , respectively. 2MASS photometry was calibrated to 0.02 mag accuracy, with persource photometric uncertainties for bright sources below 0.03 mag (Skrutskie2006).
In this work, as in Green2018, we make use of the “highreliability catalog,”^{3}^{3}3See the Explanatory Supplement to the 2MASS All Sky Data Release: https://old.ipac.caltech.edu/2mass/releases/allsky/doc/sec1_6b.html#composite and exclude sources that are marked as possibly having contamination from nearby point sources or galaxies.
3.4 Input Catalog
We generate a matched input catalog combining PS1 photometry, 2MASS photometry and Gaia parallaxes, using the Large Survey Database framework (Juric2012LSD). We use a matching radius of to match 2MASS sources to PS1, and a radius of to match Gaia DR2 sources to PS1. We require that each source be observed in at least four photometric passbands, at least two of which must be from PS1. We do not require that sources have measured parallaxes. To exclude extended sources, we require that the in at least two PS1 bands, where are magnitudes determined using PSF photometry, and are aperture magnitudes. Our resulting catalog contains 799 million sources, of which 19.5% are observed in at least one 2MASS passband, 62.4% have a Gaia parallax measurement, and 8.6% have a Gaia parallax measurement with . It should be noted, however, that even parallax measurements that do not confidently rule out zero parallax can put useful lower limits on distance. In our input catalog, 32.4% of sources are constrained at by their Gaia parallaxes to be more distant than 1 kpc.
We divide up the sky using a multiresolution HEALPix pixelization (Gorski2005). We assign the sources to angular pixels in the same manner as in Green2018, beginning with a coarse angular pixelization of the sky and recursively subdividing pixels with more than a given number of sources. Given that our spatially correlated priors allow information to be shared between nearby voxels, we more aggressively subdivide pixels than in Green2018, as reflected in the subdivision thresholds in Table 3.4 (“max. stars / pixel”), resulting in 4.21 million pixels (versus 3.42 million over the same footprint in our previous work). The properties of our input catalog are summarized in Tables 3.4 and 3.4.
4 Results
We obtain a map of dust reddening density, with angular sightlines of a typical scale of to , and 120 distance bins spaced logarithmically in distance from 63 pc to 63 kpc. Figs. 4–4 show differential reddening in four distance ranges, out to 5 kpc, while Fig. 4 shows the cumulative reddening out to the maximum distance in the map. Our map is probabilistic, and we obtain samples of the reddening density in each voxel. In order to transform this into a single value for plotting purposes, we first construct a map of median integrated reddening out to each voxel. When we wish to display reddening densities, we use the lineofsight derivative of this “median” map.
Figs. 4–4 reveal a wealth of structure at different scales and at different distances. The nearest dust structures appear off of the Galactic plane in Fig. 4 – Orion, Taurus, Perseus, California, the Aquila Rift and Ophiuchus. In Fig. 4, Monoceros R2 and Cygnus X are clearly visible. At increasing distances, in Figs. 4–4, dust features are increasingly confined to low Galactic latitude, as should be expected for the Galactic disk.
4.1 Dependence on correlation length
The Gaussian process prior, which couples the reddening density in nearby voxels, has a significant impact on the recovered reddening map. We can see this by comparing the results from the first iteration, in which the sightlines are independent, and the final iteration, in which the sightlines are coupled, as shown in Fig. 4.1. The correlated prior favors dust clouds which are isotropic, rather than stretched along the line of sight. As the transverse distance between neighboring sightlines grows with distance, the correlations between them are stronger in the nearest distance bins than in the farthest distance bins. The effect of our Gaussian process prior is thus greatest for the nearest clouds, particularly within a distance of 1 kpc.
The effect of the Gaussian process prior can be seen in a striking way in closeup views of individual clouds. Fig. 4.1 shows the reddening density in our correlated and uncorrelated maps (corresponding to the final and initial iterations, respectively) in the vicinity of the Orion B molecular cloud complex. The distance to the cloud complex is much better localized in the correlated map, whereas in the uncorrelated map, the dust is spread over a wider range of distances. Slices through the dust density field in the uncorrelated map show holes, caused by the model placing the cloud complex at different distances in neighboring sightlines. This noise is much reduced in the correlated map. The results for Orion B are typical of our results for clouds within kpc.
4.2 Comparison with 2D dust map: Planck Collaboration (2014)
Although we infer dust reddening in 3D voxels, we can compare our results with those of 2D dust maps by integrating out the distance direction in our map. In the following, we use the median of our integrated reddening at each location on the sky.
PlanckCollaboration2013 fits a modified blackbody model of dust emission to the Planck 857, 545 and 353 GHz bands, as well as the IRAS band. Planck14 then constructs two different maps of dust reddening by separately calibrating optical depth at 353 GHz and integrated radiance against reddening measurements of SDSS quasars. While farinfrared optical depth should be a better tracer of dust reddening than integrated radiance is, at high Galactic latitude, optical depth estimates can be noisy, due to covariances between the dust optical depth, temperature and (which describes the shape of the modified blackbody spectrum). In these comparisons, We therefore use the reddenings derived from integrated radiance, which we denote as . We transform these reddenings from to using the reddening curve of SchlaflyFinkbeiner2011.
In Fig. 4.2, we compare our map to the Planck radiancebased dust map out of the Galactic midplane (). In the midplane of the Galaxy, we would expect Planck to detect dust out to greater extinctions, as dust remains optically thin in farinfrared emission. By contrast, our method relies on observing stars behind the dust, and therefore cannot trace reddening to as great a depth. We compare the Planck14 map with both the initial iteration of our dust map, without spatially correlated priors, and the final iteration of our dust map, which has spatially correlated priors. For both versions of the map, we find good agreement with Planck14 out to a reddening of , with median residuals between the maps of .
At reddenings of , the final iteration of our map (with correlated spatial priors) matches Planck14 better than the initial iteration does (with completely independent sightlines), in that the residuals have a flat trend with increasing reddening. At these low reddenings, the trends in the residuals between the two maps are driven by possible zeropoint offsets in both maps, as well as the details of the priors we impose on the dust reddening. The introduction of the correlated spatial prior on dust reddening density should act to increase the signaltonoise ratio in our reddening measurements, as the reddening in each sightline is influenced by stars in nearby sightlines. If our prior favors too much reddening at high Galactic latitudes, then increased signaltonoise should drive our reddening measurements down, and closer to the correct answer.
4.3 Comparison with 3D dust maps
4.3.1 Green et al. (2018)
In Fig. 4.3.1, we compare our integrated reddenings at large distance (for both the initial and final iterations of our map) with those of Green2018. At reddenings greater than , our results are largely consistent, with median residuals of out to a reddening of . The major differences between the initial iteration and Bayes17 are the use of Gaia parallaxes in the former, as well as slightly different reddening vectors (particularly in the 2MASS passbands) and additional distance bins in the former, which has the effect of slightly altering our priors. In the final iteration of our dust map, regions of the sky with low reddening differ most from Bayes17 (and from the initial iteration), with more of the sky being placed at the lowest reddenings. The spatially correlated priors effectively cause information to be shared between nearby voxels. If the stellar photometry in one voxel constrains the reddening to be close to zero, this information propagates to the neighboring voxels, pulling down the inferred reddening. If the priors on individual voxels at high Galactic latitude favor more dust than is typically found in these regions, then the introduction of spatially correlated priors will tend to pull these regions closer to the true reddening.
4.3.2 Chen et al. (2019)
Chen2018 uses stellar parallax measurements from Gaia DR2, as well as optical and nearinfrared photometry from Gaia, WISE and 2MASS to trace dust reddening in the Galactic plane (). In Fig. 4.3.2, we compare our map with that of Chen19. Due to our use of deeper PanSTARRS 1 photometry, we trace dust density to greater distances. The greater distance resolution of our map also reveals much finer features. This greater distance resolution is enabled in part by our larger input catalog (799 million versus 56 million stars). The gross features revealed by both maps are similar, increasing our confidence that both maps are recovering real features in the interstellar medium. Neither map shows clear spiral structure – we discuss the possible relation of the features in our map to spiral structure in Section 4.4.
4.3.3 Leike & Enßlin (2019)
LeikeEnsslin2019 reconstructs the threedimensional structure of nearby dust using Gaia DR2 parallaxes and reddening estimates (the latter from Andrae2018DR2Apsis), treating the logarithm of the dust density as a Gaussian process. LE2019 applies their model to 3.7 million stars with welldetermined Gaia parallaxes in a cube centered on the Sun, obtaining the threedimensional distribution of dust, as well as its spatial correlation spectrum.
There are a number of important differences between our present work and that of LE2019. Unlike in our work, LE2019 uses a hierarchical Bayesian model in which the kernel parameters are inferred, rather than fixed beforehand. Their resulting correlation kernel has correlations on much larger scales than the kernel we assume, on the order of 30 pc. LE2019 uses a constant mean in their Gaussian process prior over all of their volume. The mean of our prior, by contrast, is set to follow a smooth Galactic disk, encoding our prior expectations about the overall structure of the Galaxy. A key difference between our methods is that in the present work, we approximate the Gaussian process prior in small angular patches surrounding each sightline, disregarding correlations between voxels that are separated by more than in angle on the sky. At a distance of 300 pc, this means that our prior ties voxels that are separated by a maximum physical distance of . This produces more spherical clouds than one would get with an uncorrelated prior (see Fig. 4.1), but does not fully exploit the information that is available in spatial correlations between nearby voxels. One consequence of this treatment is that some of the nearest clouds, such as the Ophiuchi and Taurus cloud complexes, are placed somewhat farther away in our map than they are generally taken to be. Zucker2019MolecularCloudsDR2, which uses a method which is more carefully tailored towards estimating the distance to individual clouds, obtains closer distances to these nearby clouds than does the map presented here, due to the greater number of stars per sightline and different priors used by Zucker2019MolecularCloudsDR2.
In Fig. 4.3.3, we compare orthographic projections of our map and LE2019. In Fig. 4.3.3, we compare our map to that of LE2019 in a Suncentered stereographic projection of the Galactic anticentral region. These comparisons highlight the differences between our methods. As is visible in the orthographic projections, LE2019 produces more spherical dust structures, due partially to the fact that it uses Cartesian voxels, while we use voxels that are elongated in the radial direction. As is visible in the Suncentered projection, however, our map has much finer angular resolution, due to the much larger number of stars that enter into our analysis. Because we do not require well constrained (or indeed any) Gaia parallaxes in our method, we are able to take into account a far larger number of stars, particularly at larger distances, where Gaia parallaxes become unavailable. The voxels in LE2019 are approximately 2.3 pc on a side, limiting the angular resolution that is achieved by their method, at present. Our map tends to capture more lowreddening dust structure than LE2019, as can be seen, for example, in the bottom right of Fig. 4.3.3. This may be due to the greater number of stars that go into our map, to the quality of our input stellar reddening estimates, or to differences in the choice of priors.
4.4 Spiral structure
Widely differing models of the spiral structure of the Milky Way have been proposed, with between two and four arms (GeorgelinGeorgelin1976; Drimmel2000; Vallee2008; ReidMentenZheng2009; HouHan2014), and with or without a “molecular ring” at 4 kpc from the Galactic center (Stecker1975; CohenThaddeus1977; DobbsBurkert2012). Spiral structure should be apparent not only in 3D maps of the Galaxy, but also in the distribution of interstellar gas and starforming regions in positionvelocity space (see, e.g., Fig. 3 of Dame2001). Though spiral structure is not readily visible in bird’seye views of our map (see Figs. 4.1 and 4.3.2), we can test whether the dust in our map preferentially falls along models of spiral arms developed using both position and velocity information.
XuReid2016 compiles kinematic information for masers in highmass starforming regions (HMSFRs), based on verylongbaseline interferometry in the radio, and develops a model of the Milky Way’s spiral arms. In Fig. 4.4, we compare parallax distances to these HMSFRs with our dust map. Visually, the HMSFRs preferentially lie in regions of higher dust density, as expected. In order to verify this numerically, we compare two Poisson point process models for the distribution of HMSFRs from XuReid2016:

The masers are independently and uniformly distributed in space.

The masers are independently distributed proportionally to the dust reddening density.
We carry out this analysis in a fixed box, defined by , , in a Cartesian representation of Solarcentered Galactic coordinates, excluding regions South of a declination of , which are not covered by our map. The ratio of the maximum likelihoods of the two models is given by
(28) 
where is the dust reddening density (in ) at the location of maser , and is the average reddening density over the entire volume. We find a likelihood ratio of , favoring the model in which the masers are distributed according to the reddening density.
We test that this correlation between dust density and maser position is not due to the largescale distribution of dust and masers (e.g., dust density generally being larger towards the inner Galaxy and HMSFRs also lying preferentially in this direction), but rather correlations between smallerscale features and maser positions. We do so by adding a small random displacement to the position of each maser ( pc in and , and pc in ), which should reduce any correlations between smallscale features, but keep correlations between largescale features intact. Comparing the likelihood ratio in Eq. (28) for many realizations of scattered maser positions to the likelihood ratio with the measured maser positions, we find that the measured maser positions lead to a higher ratio in greater than 99% of cases, indicating that the maser positions are correlated with smallscale features in our map.
Finally, we fit a Poisson point process model to the maser locations in which the masers are distributed according to the density , with again being the inferred dust reddening density in our “median” 3D dust map, and being a free parameter that sets a floor on the maser abundance. The likelihood of this model is maximized when satisfies
(29) 
With the XuReid2016 maser positions in the same volume as used above, we find that . The uniform floor is thus low compared to the average dust density, again indicating that the HMSFRs largely follow the dust density inferred in our 3D dust map.
Despite the correlation between our inferred dust density and the locations of HMSFRs, visual inspection of our dust map does not, in itself, reveal obvious evidence of spiral structure. This could be due to a number of factors, among them:

Our dust map does not extend beyond a few kiloparsecs. The overall spiral structure of our Galaxy might not be apparent on this scale.

The spiral structure may be less apparent in dust than in the distribution of starforming regions, particularly on small scales. This is the case in many external spiral galaxies, where the dust density can show complicated structure, with spurs and filaments connecting spiral arms.
Linear features in our map – for example, the feature running largely along the Local Arm in Fig. 4.4 – may correspond to spiral arms and spurs, but a dust map that extends to larger distances will be necessary to see any possible curvature in these features, which one would expect if they correspond to spiral arms. In particular, pushing deeper into the inner Galaxy may allow us to gain a better view of the overall structure of the Milky Way. Doing so will require deeper photometry in the nearinfrared than is provided by 2MASS.
4.5 Stellar Parameter Catalog
Inferring stellar distances, reddenings and types (absolute band magnitude and ) is a necessary first step towards generating a 3D map of reddening, as described in Section 2.1. We save samples of these stellar parameters, and make them available along with our 3D reddening map. In order to validate the accuracy of our stellar reddening estimates, we compare them to those estimated by Queiroz2018StarHorse (hereafter, “StarHorse”), which are based on APOGEE DR14 spectroscopy (Abolfathi2018APOGEEDR14), 2MASS photometry, APASS DR9 photometry (Henden2014APASS) and Gaia DR1 parallaxes and band photometry. We also compare our stellar reddenings to those of Andrae2018DR2Apsis (hereafter, “Andrae2018”), which are based on Gaia DR2 parallaxes and band, BP and RP photometry. StarHorse estimates , while Andrae2018 estimates both and . Here, we compare our reddenings to .
We first match our stellar reddening catalog to that of Andrae2018 using Gaia DR2 source_id, excluding stars for which the bestfit for our stellar model. This yields a catalog of 41.3 million sources. We convert our reddening estimates, which are in the arbitrary unit , to using the relation
(30) 
as implied by the extinction coefficients in Table 2.1. We then match our combined catalog to the StarHorse APOGEE DR14 catalog using a matching radius of , obtaining a catalog of 5350 stars.
In order to compare these catalogs on a uniform scale, we have to convert them to the same unit of reddening or extinction. In order to do this, we fit a linear relationship between the StarHorse estimates and our estimates, and a separate linear relationship between StarHorse and the Andrae2018 estimates of . Similarly to the method described in HoggBovyLang2010, we also fit nuisance parameters describing whether each star is an outlier and the overall distribution of outliers. We additionally fit a systematic error floor for our estimates and for Andrae2018’s estimates, assuming that the uncertainties reported by StarHorse are accurately estimated. The details of this model are described in Appendix E.
We obtain the relations
(31)  
(32) 
with a systematic error floor of 0.08 mag on our estimates, and a systematic error floor of 0.09 mag on the Andrae2018 estimates. With these relations, we are able to transform between the three reddening measurements, and compare them on a common scale. The left panel of Fig. 4.5 shows the residuals stellar recovered by our method and by StarHorse, as a function of , while the right panel shows the comparison between Andrae2018 and StarHorse (likewise in ). As StarHorse makes use of spectroscopically determined stellar atmospheric parameters, its extinction estimates are generally much more precise than those produced by our method or by Andrae2018. The residuals in Fig. 4.5 are thus dominated by errors in our measurements (in the left panel) or in Andrae2018’s measurements (in the right panel). As can be seen from these panels, the errors in our stellar extinction estimates are generally smaller than those of Andrae2018, due to the greater number of photometric passbands used in our study (up to eight, versus two). In Fig. 4.5, we directly compare the uncertainties in stellar reddening inferred by our method and Andrae2018, for the 41.3 sources in our matched catalog, with the recovered systematic error floors added in quadrature. In , the percentile of the reddening uncertainties is 0.11 mag for our estimates, and 0.27 mag for the Andrae2018 estimates. In this matched catalog, our median uncertainty is 0.086 mag, 30% less than the median uncertainty of 0.124 mag in Andrae2018. These represent the typical uncertainties that can be expected for relatively bright stars with wellmeasured Gaia parallaxes. However, we obtain individual stellar reddening estimates for all 799 million stars in our input catalog, including those without Gaia parallax measurements. We make this entire catalog available to the community, as described in Section 6.
Figs. 4.5 shows the spatial distribution of stars in our catalog, excluding those with poor goodnessoffit (maximumlikelihood ). This distribution is determined by the true number density of stars throughout the Galaxy, the sensitivities of the PS1, 2MASS and Gaia surveys, our selection function, and the distribution of interstellar dust. With a good understanding of the latter confounding effects, this catalog can be used as a basis for determining the true distribution of stars throughout the Galaxy. Fig. 4.5 shows the distribution of stars with absolute magnitudes and distances that are determined to better than 20%. Although a detailed treatment of selection functions is necessary to determine the structure of the Galaxy, the Galactic bulge is visible as in increase in stellar density near the Galactic center in the projection to the plane. Because the boundary of the PS1 survey passes essentially through the Galactic center, the bulge appears offcenter from the Galactic center. The relative lack of stars directly in the Galactic midplane () is due to dust obscuration, while the lack of bright stars close to the Sun is due to saturation in photometric surveys. This catalog will serve as a basis for future work on the structure of our Galaxy.
5 Discussion
The work presented here can be extended in a number of ways, both through improvements of the model and incorporation of more data. One methodological improvement would be to optimize the kernel parameters in the spatial prior. In this work, we have demonstrated that the incorporation of a Gaussian process prior on the logarithm of the dust reddening density yields significant improvements in the effective distance resolution of the map, with cloud distances becoming much more precise. However, the kernel that we impose is determined beforehand, rather than inferred from the data. An optimized kernel would allow us to extract more information about the dust reddening density from our data. As LeikeEnsslin2019 shows, with an optimized kernel, it is possible to extract a significant amount of information about the spatial distribution of dust from even relatively small stellar catalogs (3.7 million stars, albeit with well determined parallaxes, versus 800 million in this work). An additional improvement to the method here would be to take into account larger patches of sky surrounding each sightline, in order to better approximate the target Gaussian process prior. This is particularly important in the near regime, in which our angular patches correspond to smaller physical extent. Finally, although we increase the number of distance bins in this work by a factor of four over Green2018, our voxels are still elongated in the radial direction. This corresponds to the fundamental fact that we have much better knowledge of the angular distribution of dust than of its distribution with distance, but it makes implementation of spatially correlated priors – particularly those that obey the Copernican principle – more difficult.
A separate area in which significant improvement is achievable is in the dust reddening curve that we assume. At present, our work (like all other 3D dust maps that we are aware of) assumes a single, universal dust extinction curve for the entire Galaxy. However, the shape of the dust extinction spectrum varies significantly throughout the Galaxy (FitzpatrickMassa1986ExtinctionI; FitzpatrickMassa1988ExtinctionII; FitzpatrickMassa1990ExtinctionIII; FitzpatrickMassa2005ExtinctionIV; FitzpatrickMassa2007ExtinctionV; FitzpatrickMassa2009ExtinctionVI). Although there are multiple features in the dust extinction spectrum which can vary independently (such as the width and strength of the 2175 Å bump, the slope of the wavelengthextinction relation in the optical, and the shape of the rise in extinction in the far ultraviolet), the dust extinction spectrum is often parameterized as a singleparameter family of curves (Cardelli1989RV), using the ratio of total to selective extinction, , defined as:
(33) 
In the diffuse interstellar medium of the Milky Way, is typically taken to be approximately 3.1 (Cardelli1989RV). However, the value of can vary significantly throughout the interstellar medium. Schlafly2017MappingExtinctionCurve combines determinations for individual stars (Schlafly2016ExtinctionVector) with the Bayestar2015 3D dust map (Green2015) to map the spatial variation of the extinction curve throughout the Milky Way, finding variations on kiloparsec scales. As most dust maps that use stellar photometry rely primarily on changes in stellar colors, and are less sensitive to changes in absolute magnitude, they are essentially maps of dust reddening, rather than extinction. Extinction maps derived from these reddening measurements (e.g., using Eq. 33) without knowledge of the spatial dependence of will have largescale, spatially dependent systematic errors. Tackling this source of systematics is of particular importance to precision cosmological measurements (Nataf2016RVCosmology; Huterer2013).
One way to address variation in , without unduly complicating our model, is to treat the dust reddening density as a pair of numbers at each point in space. Schlafly2016ExtinctionVector decomposes the reddening vector into a mean component, and an orthogonal vector along which most of the variation in the reddening vector occurs. The virtue of this formulation, from the perspective of 3D dust mapping, is that the two components of reddening add linearly, simplifying their inclusion in our model. Zucker2019MolecularCloudsDR2, which measures distances to a catalog of molecular clouds using stellar photometry and Gaia parallaxes, takes advantage of this linearity to determine the distance and of the dust.
More generally, any systematics in reddening maps at high Galactic latitude are of importance to cosmology. Both this work and Green2018 have identified systematic trends in reddening when comparing with the Planck14 and SFD dust maps, particularly for reddenings of . These systematics could be due to the uncertain zeropoint of the dust reddening (e.g., what is the absolute reddening of some reference point on the sky?) and variation in across the highGalacticlatitude sky. An additional worry for cosmology is contamination from largescale extraGalactic structure. For farinfrared emissionbased maps, such as Planck14 and SFD, this contamination comes from thermal dust emission in external galaxies residing at a wide range in redshift (ChiangMenard2019). For maps based on stellar colors, such as ours, this contamination comes from quasars and unresolved galaxies misconstrued as stars (ChiangMenard2019). More effective rejection of these classes of objects (for example, by treating them explicitly in our model) will be important to providing reddening maps for use in precision cosmological measurements.
There are also additional datasets which can be applied to the dustmapping problem. Deeper nearinfrared photometry, in particular, will allow us to trace dust to greater extinctions. Dust maps which extend to greater distances, particularly in the Galactic midplane, will allow us to better uncover the overall structure of the Milky Way. As discussed in Section 4.4, with our present map, extending only a few kiloparsecs, the spiral structure of the Galaxy is not readily apparent. Incorporation of deep infrared photometry – e.g., from the Spitzer GLIMPSE surveys (Churchwell2009GLIMPSE), the ongoing UKIDSS Galactic Plane Survey (Lucas2008UKIDSSGPS) in the North and the VISTA Variables in the Via Lactea survey in the South (Minniti2010VVV), and from the newly published unWISE catalog (Schlafly2019unWISE) – would extend the range of our dust map, particularly in the inner Galaxy. This would allow us to directly see the overall structure of the spiral arms, rather than drawing curves through short segments of them, as we are currently limited to doing with available 3D dust maps.
Finally, the work presented in this paper is limited to the Galaxy north of a declination of . The ongoing Dark Energy Camera Plane Survey (Schlafly2018DECaPS, hereafter “DECaPS”) has already filled in the Southern Galactic plane in the region , and is in the process of being extended to . DECaPS observes in five PS1like passbands on the Dark Energy Camera, mounted on the Blanco 4m telescope in Cerro Tololo, Chile. Incorporation of the resulting photometry will allow us to map all of the Galactic midplane.
6 Accessing the map
The dust map presented here can be downloaded at doi:10.7910/DVN/2EJ9TX. The easiest way to interact with the map, however, is through the Python package dustmaps (Green2018dustmaps), which provides a uniform interface to a range of maps of dust reddening and extinction, including all of the dust maps discussed in this paper. The map can additionally be queried interactively at argonaut.skymaps.info.
In addition to the dust map, we provide our individual stellar inferences (described in Section 2.1). Samples of reddening, distance modulus, absolute magnitude and metallicity for 799 million stars can be downloaded at doi:10.7910/DVN/AV9GXO. This is an order of magnitude more stellar reddenings than are provided by the Gaia DR2 catalog (Andrae2018DR2Apsis). Our technique leverages a greater number of photometric passbands (between four and eight, compared to two independent passbands for Gaia DR2), and delivers typical reddening uncertainties that are 30% lower. For stars which are present in Gaia DR2, we provide a crossmatch to the corresponding Gaia DR2 source_id.
7 Conclusion
We have presented a new 3D map of dust reddening, covering three quarters of the sky out to a distance of several kiloparsecs. This map is based on stellar distances and reddenings, inferred from Gaia DR2 parallaxes and optical and nearinfrared photometry from PS1 and 2MASS. This dust map has four times the distance resolution as Green2018. Another improvement over our previous method is that we impose a Gaussianprocess prior on the logarithm of the dust reddening density, using an approximation that requires far less computational time than a naive approach would require. This prior encodes our knowledge that the spatial distribution of dust should be smooth on some scale, and has the effect of sharpening our distance determinations to dust clouds, particularly in the nearest kiloparsec, and reducing noise in our reddening estimates.
We have made both the 3D dust map, and the parameter inferences (distances, reddenings and types) for 799 million stars that underlie the map, available at doi:10.7910/DVN/2EJ9TX and doi:10.7910/DVN/AV9GXO, respectively. The dust map is also available through the Python package dustmaps (Green2018dustmaps). For the stellar parameter inferences, we provide matches to the Gaia DR2 catalog.
At high Galactic latitude, our integrated reddening estimates agree well with those of Planck14 and Bayestar17. We find median residuals between our map and those of the farinfrared radiancebased Planck14 reddening map of , out to a reddening of . Our new map matches Planck14 better at low column densities, due to the superior signaltonoise ratio of its reddening determinations.
Comparisons with other 3D dust maps show similar overall structures, though the maps have different resolutions and ranges of validity. We see the same basic features as Chen19 within a distance of , although we trace the dust distribution with finer distance resolution and to greater distances. Much of this difference is likely attributable to our use of deeper optical photometry from PS1, as well as the use of a Gaussianprocess prior. In the nearby regime, in a box centered on the Sun, we compare our map with that of LE2019. There are many differences between our methods, with LE2019 employing a Cartesian voxelization of the dust, in contrast to the spherical voxelization that we choose. LE2019 likewise imposes a Gaussianprocess prior on the logarithm of the dust reddening density, but infers the parameters of the dust reddening correlation function as part of a hierarchical model. The spatial resolution of the resulting LE2019 map is greater than ours, while the angular resolution of our map is greater.
The spiral structure of the Milky Way galaxy is uncertain, with a wide range of models having been developed. XuReid2016 uses HMSFRs with kinematic and distance measurements from embedded masers to trace spiral arms. We have shown that the masers compiled by XuReid2016 are preferentially distributed in regions of high reddening density in our map. However, in order to use the 3D distribution of dust reddening to assess models of spiral arms, it is necessary to extend our map beyond the current limits of a few kiloparsecs, particularly in the inner Galaxy. One way of achieving this will be the incorporation of deeper nearinfrared photometry, allowing us to see stars through far greater dust column densities.
8 Acknowledgements
The computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.
This work took part under the program MilkyWayGaia of the PSI2 project funded by the IDEX ParisSaclay, ANR11IDEX000302.
E.S. acknowledges support for this work by NASA through ADAP grant NNH17AE75I and Hubble Fellowship grant HSTHF251367.001A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 526555. D.F. and C.Z. acknowledge support by NSF grant AST1614941, “Exploring the Galaxy: 3Dimensional Structure and Stellar Streams.”
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.
The PanSTARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the PanSTARRS Project Office, the MaxPlanck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg, and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the HarvardSmithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under grant No. AST1238877, the University of Maryland, and Eötvös Loránd University (ELTE).
This publication makes use of data products from the 2MASS, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.
Software: Astropy (astropy2013; astropy2018), dustmaps (Green2018dustmaps), emcee (ForemanMackey2013emcee), Matplotlib (Hunter2007Matplotlib), NumPy (NumPy), SciPy (SciPy)
APPENDIX
A Gridevaluation of stellar posterior densities
In contrast to our approach in Green2014; Green2015 and Green2018, here we bruteforce evaluate the singlestar posterior probability densities on a regular grid, rather than approximating them using a kernel density estimate on Markov Chain Monte Carlo samples. This bruteforce approach is feasible – and indeed extremely fast – with a few approximations. For any given stellar type, the observed magnitudes depend only linearly on the reddening and distance modulus. As the uncertainties on our observed magnitudes are Gaussian, the likelihood of any parameters that the magnitudes depend on linearly is also Gaussian. In other words,
(A1) 
is Gaussian in and for fixed . Moreover, the covariance of this Gaussian depends only on the observational uncertainties and the reddening curve.
Our scheme for calculating the stellar posteriors is thus as follows:

For each point in a uniform grid in , calculate the and that maximize the likelihood . Call these values and .

Add a Gaussian centered on each and , weighted by .
The second step makes the approximation that over the Gaussian centered on a given , the likelihood in parallax and the prior on are constant. In our testing, this approximation typically only had a small effect.
Using this method, we evaluate over a grid in distance modulus and reddening, with a bin scale of 0.125 mag in distance modulus and 0.01 mag in reddening. This information is carried over into the lineofsight inference. However, in order to save disk space, we ultimately store samples from the posterior, discarding the gridevaluated posterior density after we have completed the lineofsight inference.
This scheme is fast. On a single core, with a library of 41000 stellar templates, we are typically able to evaluate the gridded posterior densities of 50 stars per second.
B Roundrobin sampling of neighboring pixels
To update each neighboring sightline, we choose a new sample from the previous iteration, conditioned on the reddening densities in the other pixels in the patch. This is essentially an importancesampling step, and we must choose the new sample according to the ratio of the new posterior density to the posterior density of the sample in the previous iteration. During each iteration, we must therefore store the posterior density for each sample of the central pixel, to be used when resampling pixels in the following iterations.
How does one calculate this posterior density? For a given sightline ,
(B2)  
(B3)  
(B4)  
(B5) 
where is a normalizing constant dependent only on the data, and can therefore be ignored. For simplicity, here denotes both photometry and parallaxes. The second term, , is the likelihood of , which does not change from iteration to iteration. The third term, , is the “effective prior” on generated by the neighboring sightlines, averaged over the range of values that the reddening in those sightlines can take.
In the first iteration, the calculation of this “effective prior,” , is trivial. As the pixels are independent, it is simply equal to the prior on . However, in subsequent iterations, we must estimate for each sample we store of , by calculating the average of over samples of drawn from the chain. This is done as a postprocessing step after we are done sampling the patch of sky. The computational complexity of this step is in the number of samples stored, because for each of the samples of we store, we have to average over for samples of .
In order to select a new sample for one of the neighboring pixels , we therefore select one of the stored samples from the previous iteration, with the probability of selecting sample being proportional to
(B6) 
where and are the stored likelihood and effective prior, respectively, for sample of sightline from the previous iteration.
If the stored samples for the neighboring pixels are a decent approximation to the new correlated posterior on neighboring sightlines, then this sampling procedure will work well. Over succeeding iterations, we slowly change the covariance kernel, from completely diagonal (i.e., independent sightlines) to the target kernel. If we change the covariance kernel too drastically from one iteration to the next, then the stored samples will be a poor representation of the posterior density, and will not be able to represent the state of the neighboring sightlines in the current iteration. As in importance sampling, this would manifest itself as the entropy of the weight distribution approaching zero.
C Improving Markov Chain mixing with the “deformation ladder”
The sampling procedure for a patch of sky described above can suffer from poor Markov Chain mixing. Our distance bins are significantly longer in the radial than in the angular direction (by a typical factor of 30), leading to the transverse correlations between neighboring voxels at the same distance being much greater than correlations between pixels at different distances. One possible failure mode of our sampling procedure may occur when there is a dust cloud at a particular distance, spanning several neighboring sightlines. The correlated prior will tend to favor solutions in which the dust is placed at the same distance in all the sightlines. The prior imposes a penalty if the dust cloud is moved forward or backward in distance in just one sightline. With the roundrobin sampling procedure we employ for each patch of sky, this means that it is difficult for the cloud to transition from one distance to another. The sampler becomes locked in a state with the cloud at one distance, and has to traverse a high potential barrier in order to reach solutions with the cloud at different distances.
In order to deal with this problem, we employ a method which is conceptually similar to parallel tempering. Parallel tempering works by creating a ladder of samplers, with the bottom rung sampling from the target distribution, and successively higher rungs sampling from progressively flattened versions of the distribution (using a sampling temperature parameter to control the flattening of the distribution). The highertemperature samplers are able to explore a larger region of parameter space, allowing the sampler to transition between different modes, while swaps between the states of different rungs on the ladder allow the lowertemperature samplers to benefit from the increased exploration at higher temperatures.
In our problem, parallel tempering performs badly beyond the initial iteration (without coupling between nearby sightlines), as the samples we store from previous iterations are a poor representation of the highertemperature distributions. We build a ladder of samplers, but employ a different parameterization to modify the target distribution. Our approach is to modify the inverse covariance matrix, in order to introduce correlations between voxels in neighboring distance bins. Instead of a temperature, we introduce a distancemixing parameter, , which allows us to smoothly interpolate between the approximately isotropic target distribution and distributions in which correlations in the radial direction are stronger than in the transverse direction.
The form of the deformation of the inverse covariance matrix is sketched out in Fig. C. The entries linking neighboring distance bins of different angular pixels increase with the distancemixing parameter, , according to the formula given in Fig. C. The entries linking neighboring distance bins of the same angular pixel (i.e., sightline) are fixed to zero. This ensures that when we sample a single sightline, keeping the neighbors fixed, the different distance bins are independent, and thus speeds up the sampling procedure significantly.
A the distancemixing parameter of indicates no change to the inverse covariance matrix. As increases, the penalty for moving a dust cloud forward or backward in a single distance bin typically decreases.
Just as in parallel tempering, we can treat the ladder of samplers as if they are sampling a single, augmented parameter space. If the original parameter space is denoted by , then the augmented parameter space is the Cartesian product of , with corresponding to distancemixing parameter , and being the number of rungs in the ladder. The probability density that we assign to the augmented space is given by
(C7) 
where uses the inverse covariance matrix deformed with parameter . Because the probability density of the augmented space is factorizable in this manner, we can sample in each subspace independently. However, we introduce swap steps between neighboring rungs of the ladder, in which we propose to exchange the values of and . The MetropolisHastings acceptance probability for this step is
(C8) 
For parallel tempering, this acceptance probability simplifies down to an elegant form, but in the more general case of a ladder of modified probability densities, the above acceptance ratio does not necessarily simplify.
Notice that the crossdistance coupling introduced in the higher rungs of the “deformation ladder” only affects the Gaussian process prior that couples nearby sightlines. In the initial iteration, however, we treat each sightline independently, so the “deformation ladder” would have no effect. In the initial iteration, we therefore use standard parallel tempering in order to improve convergence of our Markov Chain Monte Carlo sampler. In correlated iterations, we use a ladder with three rungs, with the distancemixing parameter increasing linearly from in the bottum rung to in the top rung. Note that the distancemixing parameter is nonzero in the bottom rung, meaning that our target posterior density is slightly modified, and contains small crossdistance correlations.
D Assessing convergence of the lineofsight reddening sampler
We test for convergence of our lineofsight sampling by measuring the autocorrelation time of certain functions of the lineofsight reddening density parameters, . There are many different functions that one could choose to try to assess chain mixing, but we choose three: the prior, the likelihood, and the amplitude of the first principal component of in the central pixel. This third function measures how quickly the sampler explores the direction in parameter space that accounts for most of the variance in the Markov Chain. Typically, this direction in parameter space corresponds to changing the distance to the dominant cloud along the line of sight. We require that the number of samples in the Markov chain be at least 20 times the autocorrelation time of each of these three functions. Fig. D illustrates these convergence criteria for both a converged and nonconverged case.
For pixels that fail this convergence criterion, we sample again with more robust settings (i.e., more steps and additional rungs in the deformation ladder). If a pixel fails to converge within three attempts, we flag the pixel and continue. Fewer than 0.25% of pixels fail this convergence test in any of the iterations.
E Fitting a line to data
Here, we describe our method of fitting a linear relation between two quantities with measurement uncertainties. In Section 4.5, we use this model to determine the relation between , as measured by StarHorse, and either our inferred or the inferred from Andrae2018. The model we describe here is very similar to that described by HoggBovyLang2010, with the addition of a systematic error floor in the measured values of .
Our observed data consists of measured values of and , with reported Gaussian uncertainties and . We assume that each data point is either “good,” in which case the true values of and (without measurement error) follow a linear relation,
(E9) 
or an “outlier,” in which case the true values of and are independently distributed Gaussian random variates, with means and , and standard deviations and . The prior probability of any given data point being an outlier is . We additionally assume that the reported values of are systematically underestimated, and that the true uncertainty in the measured value of is given by
(E10) 
In order to determine the calculate the likelihood of a single data point,
(E11) 
it is necessary to specify a prior on the true value of . We impose a flat prior, which yields a closedform solution for the likelihood, as long as the width of the prior is large compared to and contains the measured value, . This unfortunately introduces one more arbitrary parameter into the problem, , the width of the flat prior in . We choose to be approximately the same as the range of observed , and the recovered parameters do not depend sensitively on the choice of . It can be shown that the likelihood of a single data point in this model is given by
(E12) 
We assume that the data points have independently distributed errors, so that the likelihood of the entire data set is the product of the individual likelihoods.
In Section 4.5, we treat the StarHorse estimates of as , and our estimates or Andrae2018’s estimates as . The priors on our model parameters are given in Table E. We use the emcee Python package (ForemanMackey2013emcee) to sample from our model, using 90 walkers, 2500 steps per walker for burnin, and subsequently 5000 steps per walker. The resulting median inferred values of , and are given in Section 4.5.
bayestar2019