A THREE-DIMENSIONAL MAP OF MILKY WAY DUST

Green et al. (2014) describes how we use photometry of stars to determine the distribution of dust in the Galaxy. We begin by separating the sky into small sightlines. In a small region of the sky across which the dust column does not vary much angularly, more distant stars should be more heavily reddened than nearby stars, as all the stars lie along the same dust column. We therefore begin by calculating a probabilistic reddening and distance estimate for each star in the sightline. We then determine the amount of dust reddening in each distance bin along the sightline, under the constraint that the line of sight distance versus reddening profile has to be consistent with our distance and reddening inferences for all the stars in the sightline.

Green et al. (2014) states our map-making method generally, without choosing specific priors on the distribution of reddening in 3D. Here, we will briefly restate the formalism used in Green et al. (2014), and then define our 3D reddening priors.

Our goal of inferring the line of sight dust reddening on the basis of stellar photometry begins with inferring the distance and reddening of each star based on its photometry. Label each star in the line of sight by a number i, and denote the probability density for an individual star to lie at distance modulus, μ, and reddening, E, given photometry ${\boldsymbol{m}}$, by

$$\begin{eqnarray}&&p({\mu }_{i},{E}_{i}| {{\boldsymbol{m}}}_{i}).\end{eqnarray} \tag{ 1 }$$

We precompute $p({\mu }_{i},{E}_{i}| {{\boldsymbol{m}}}_{i})$ using a kernel density estimate of MCMC samples. This technique, similar to that used in Hogg et al. (2010), is not guaranteed to converge to exactly the target density. It is possible, however, to substitute a mathematically correct, but somewhat more computationally expensive algorithm, pseudo-marginal MCMC, (Beaumont 2003; Andrieu & Roberts 2009), in which a new noisy estimate of a marginalized likelihood term is computed at each Markov chain step. However, we expect any possible bias introduced by precomputing Equation (1) once to be small due to the relatively large number of stars in each sightline. Compared to possible inaccuracies in our stellar model, for example, we expect the inexactness of our sampling method to be a small effect.

As explained in Green et al. (2014), we place a flat prior on E_i for each individual star when precomputing Equation (1), as the prior on reddening will come in when we combine information from all of the stars in a sightline. We now parameterize the line of sight reddening by a set of parameters, ${\boldsymbol{\alpha }}$, so that the cumulative reddening out to a distance modulus μ can be written as:

$$\begin{eqnarray}&&E(\mu ;{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ 2 }$$

The probability density of ${\boldsymbol{\alpha }}$ is then given by a product over line integrals through the individual stellar probability density functions, following $E(\mu ;{\boldsymbol{\alpha }})$:

$$\begin{eqnarray}&&p({\boldsymbol{\alpha }}| \{{\boldsymbol{m}}\})\propto p({\boldsymbol{\alpha }})\;\displaystyle \prod _{i}\int d{\mu }_{i}\;p({\mu }_{i},E({\mu }_{i};{\boldsymbol{\alpha }})| \;{{\boldsymbol{m}}}_{i}).\end{eqnarray} \tag{ 3 }$$

We parameterize the reddening in each line of sight as a monotonically increasing, piecewise-linear function in distance modulus. Dividing up the line of sight into bins of equal width in distance modulus, we sample the increase in reddening in each bin, so that

$$\begin{eqnarray}&&{\boldsymbol{\alpha }}=\{{\rm{\Delta }}{E}_{k}| k=1,2,\ldots ,{N}_{\mathrm{bins}}\},\end{eqnarray} \tag{ 4 }$$

where k denotes the bin index, and ${N}_{\mathrm{bins}}$ is the number of distance bins. We place a log-normal prior on the reddening in each bin:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{k}={e}^{{\epsilon }_{k}},\ \mathrm{where}\ {\epsilon }_{k}\sim {\mathcal{N}}({\bar{\epsilon }}_{k},{\sigma }_{\epsilon }).\end{eqnarray} \tag{ 5 }$$

In each bin, we set ${\bar{\epsilon }}_{k}$ so that the mean reddening in each voxel matches what would be expected from the smooth disk component of Drimmel & Spergel (2001), up to a constant normalization that is the same for the entire map. In order to fully define the priors, there are therefore two global parameters that must be set:

• 1.  
  ${dE}(B-V)/{ds}{| }_{s=0}$, the local normalization of the dust reddening per unit distance, and
• 2.  
  ${\sigma }_{\epsilon }$, the width of the log-normal prior on dust column density in each voxel.

We fit these two numbers so that our priors, in projection, produce two-dimensional dust maps with the same overall normalization and standard deviation as the SFD dust map. We find that ${\sigma }_{\epsilon }=1.4$ and a local reddening per unit distance of ${dE}(B-V)/{ds}$ = 0.2 mag kpc⁻¹ roughly matches the mean and variance of the SFD dust map over large angular scales.

Figure 1 shows a random realization of our 3D dust priors, next to the SFD dust map. In each bin, we limit $\bar{\epsilon }$ to the range $-12\leqslant \bar{\epsilon }\leqslant -4$. On the lower end, the limit helps our fit to converge by preventing the prior from becoming too stringent. On the upper end, the limit prevents our fit from inferring large amounts of dust in the absence of data. Regions where $\bar{\epsilon }$ has been limited to $-4$ are shaded in Figure 1. Our priors do not impose spatial correlations across lines of sight, and thus in the absence of data, produce cloud-free maps. The detailed cloud structures that emerge from our 3D dust modeling therefore derive entirely from the data, rather than from the priors.

Zoom In Zoom Out Reset image size

It is worth noting that our priors diverge from what one should expect for the real distribution of dust in a number of ways. In order to render the problem of fitting the 3D distribution of dust tractable, we fit the dust column along each sightline separately. We know, however, that in real life, dust density has spatial correlations. Moreover, there is nothing special about the Sun's place in the Galaxy, but our dust map voxelizes the sky into pencil beams centered on the solar system. This of course makes practical sense, since angles are much easier to measure than distances in astronomy, but it is again an unreal feature of our voxelization. A more realistic model would treat the dust density as a continuous field, or voxelize the Galaxy in a way that treats the angular and radial directions equally, and would impose correlations between dust density in nearby points in space (see, e.g., Lallement et al. 2013; Sale & Magorrian 2014). This entails significantly more algorithmic and computational complexity than the method used here, and we defer such work to the future. Within the constraints of our present setup—independent sightlines with pencil-beam-like voxels—our priors attempt to reasonably trace the properties of the Galaxy, including the mean dust density in each voxel and the overall variance in dust column across the sky.

A basic assumption of our model, as described in Green et al. (2014), is that within a given pixel, the dust density varies only with distance, and not with angle. If the pixels are sufficiently small, this is a good assumption. We are limited, however, in how small we can make each pixel by the need to include enough stars in each pixel to probe the dust density at a range of distances. Increasing the angular resolution of the map decreases the number of stars in each pixel, effectively decreasing the distance resolution of the map. We have found that we obtain the best results for pixels containing a few hundred stars, and we vary the resolution of our pixels across the sky to obtain approximately the same number of stars in each pixel (see Figure 2).

Zoom In Zoom Out Reset image size

Note that varying the pixel size across the sky in this way technically violates the principles of Bayesian inference. From a forward-modeling point of view, the distribution of dust influences the number of stars that are observed in any given region of the sky. Sale (2015), for example, discusses how a catalog of stars observed in a survey can be described as a Poisson point process whose rate is determined by the distribution of stars in the Galaxy, the three-dimensional distribution of dust, and the survey selection function. Yet we are using the number of stars observed in each part of the sky to determine how to pixelize the sky. Our pixelization is therefore set, in part, by an observable (the density of stars across the sky) that is a consequence of the model. Nevertheless, as we are treating each pixel independently, and we would like to keep pixels as small as possible without reducing the number of stars per pixel below a few hundred, this violation is difficult to avoid. We expect the impact of this violation to be small in most regions of the sky. In regions with large sub-pixel variation in dust reddening, however, our map may be biased toward lower reddenings, as we preferentially detect stars in regions of the pixel with lower reddening.

The typical resolution of the map is $6\buildrel{\,\prime}\over{.} 8\times 6\buildrel{\,\prime}\over{.} 8$, corresponding to an ${\mathtt{nside}}=512$ HEALPix pixelization (Gorski et al. 2005). At this resolution, there can still be significant power in the dust density spectrum below the pixel scale. This can pose problems for our method, especially in the vicinity of dense clouds and filamentary structures, where the sub-pixel angular variation is largest. In order to deal with sub-pixel angular variation in the dust density, we relax our assumption that all stars lie along the same dust column by allowing each star to deviate from the local "average" dust column by a small amount. The reddening of star i is parameterized as

$$\begin{eqnarray}&&{E}_{i}=(1+{\delta }_{i})E({\mu }_{i};{\boldsymbol{\alpha }}),\end{eqnarray} \tag{ 6 }$$

where ${\delta }_{i}$ is the fractional offset of the star from the local dust column, $E({\mu }_{i};{\boldsymbol{\alpha }})$.

In effect, our model is therefore that within each HEALPix pixel, the reddening is a white noise process, with a mean that increases piecewise linearly with distance. Each star samples this white noise process at a particular distance and angular position within the pixel. The parameter ${\delta }_{i}$ is then understood as the fractional residual (from the mean reddening in the pixel at the given distance) of the reddening column at the angular location and distance of star i.

We put a Gaussian prior on ${\delta }_{i}$, with zero mean and standard deviation dependent to the scale of the pixel (allowing more variation in larger pixels) and the local dust column (allowing greater fractional variation in regions of greater reddening). In detail,

$$\begin{eqnarray}&&p({\delta }_{i}| {\mu }_{i},{\boldsymbol{\alpha }})={\mathcal{N}}({\delta }_{i}| 0,{\sigma }_{\delta }),\end{eqnarray} \tag{ 7 }$$

with

$$\begin{eqnarray}&&{\sigma }_{\delta }={aE}({\mu }_{i};{\boldsymbol{\alpha }})+b.\end{eqnarray} \tag{ 8 }$$

Here, a and b are parameters that we set in order to match the variation we see at the given pixel angular scale in the Planck radiance-based two-dimensional dust map. We compute the rms scatter within HEALPix pixels of difference scales, finding that the scatter is well described by setting the coefficients a and b to

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\;a=0.88{\mathrm{log}}_{10}(\displaystyle \frac{\phi }{1^{\prime} })-2.96,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\mathrm{log}}_{10}\;b=0.58{\mathrm{log}}_{10}(\displaystyle \frac{\phi }{1^{\prime} })-1.88,\end{eqnarray} \tag{ 10 }$$

where ϕ is the angular pixel scale, defined as the square-root of the pixel solid angle.

We have found through trial and error that it is preferable to impose a minimum scatter of $\sim 10\%$ on the in-pixel dust column. Additionally, if we allow ${\delta }_{i}$ to approach unity, we clearly risk the possibility of scattering a star to a negative dust column. We therefore never allow ${\sigma }_{\delta }\gt 0.25$.

Although the complication introduced in this section adds an additional parameter, ${\delta }_{i}$, for each star, it can be achieved with minimal additional computational resources. We are introducing a Gaussian scatter in the reddening of each star from the "average" reddening in the pixel, and an appropriate Gaussian smoothing of the individual stellar probability density surfaces, $p({\mu }_{i},{E}_{i}| {{\boldsymbol{m}}}_{i})$, achieves this effect. Appendix D details how the individual stellar probability surfaces are smoothed.

In summary, our model of each sightline contains the following elements:

• 1.  
  The increase in the "average" reddening in each distance bin: ${\boldsymbol{\alpha }}=\{{\rm{\Delta }}{E}_{k}| k=1,\ldots ,{n}_{\mathrm{bins}}\}$.
• 2.  
  The distance modulus, ${\mu }_{i}$, stellar type, ${{\rm{\Theta }}}_{i}$ and fractional offset, ${\delta }_{i}$ from the "average" reddening of each star i.

The reddening of star i is therefore determined both by ${\boldsymbol{\alpha }}$ and ${\delta }_{i}$. Together with the distance and type of the star, one can obtain model apparent magnitudes for the star, and thus a likelihood:

$$\begin{eqnarray}&&p({{\boldsymbol{m}}}_{i}| {\mu }_{i},{{\boldsymbol{\Theta }}}_{i},{\delta }_{i},{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ 11 }$$

We also have per-star priors on distance and stellar type, given by a smooth model of the distribution of stars throughout the Galaxy, and a prior on the offset of each star from the local reddening column:

$$\begin{eqnarray}&&p({\mu }_{i},{{\boldsymbol{\Theta }}}_{i},{\delta }_{i}| {\boldsymbol{\alpha }})=p({\mu }_{i},{{\boldsymbol{\Theta }}}_{i})p({\delta }_{i}| {\mu }_{i},{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ 12 }$$

We finally have a log-normal prior on the increase in reddening in each distance bin along each sightline, $p({\boldsymbol{\alpha }})$, whose mean is chosen to match a smooth model of the distribution of dust throughout the Galaxy.

The posterior on the line of sight reddening along one sightline is then given by

$$\begin{eqnarray}p({\boldsymbol{\alpha }}\;| \{{\boldsymbol{m}}\}) & = & p({\boldsymbol{\alpha }})\displaystyle \prod _{i=1}^{{n}_{\mathrm{stars}}}p({{\boldsymbol{m}}}_{i}| {\mu }_{i},{{\boldsymbol{\Theta }}}_{i},{\delta }_{i},{\boldsymbol{\alpha }})\\ & & \times p({\mu }_{i},{{\boldsymbol{\Theta }}}_{i})p({\delta }_{i}| {\mu }_{i},{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ 13 }$$

We pre-compute the likelihood and prior terms for the individual stars, and then sample in ${\boldsymbol{\alpha }}$. The details of how this is done are given in Green et al. (2014), and in Appendix D of this paper.

In order to create the 3D dust map, we first probabilistically infer the distance and reddening to each star individually, with no information about the 3D structure of the dust. After creating the 3D dust map, however, we have a very strong constraint on how reddening should increase with distance, and this should impact our inferences about individual stellar distances and reddenings. In order to improve our stellar parameter inferences, we should replace our initial assumption about stellar reddening (a flat prior) with a new one that favors stellar reddenings close to the measured distance-reddening relation along the sightline. This can be done, in practice, by reweighting the probability density functions we initially calculated for each star. In this section, we therefore define a reweighting of the Markov Chain samples for the individual stars which takes into account the line of sight reddening.

Let us first consider the case in which we fix the line of sight reddening profile. In the formalism used here, that means that we fix the parameters ${\boldsymbol{\alpha }}$, expressing the line of sight reddening profile as $E(\mu ;{\boldsymbol{\alpha }})$. The reddening of an individual star is then determined by the distance modulus of the star, μ, as well as the fractional deviation, δ, of the stellar reddening from the local dust column. As before, there are also parameters representing the stellar type, ${\boldsymbol{\Theta }}$. In this formulation, the posterior density of the stellar parameters, given its photometry ${\boldsymbol{m}}$ and the line of sight reddening profile, is determined by

$$\begin{eqnarray}&&p(\mu ,{\boldsymbol{\Theta }},\delta | {\boldsymbol{m}},{\boldsymbol{\alpha }})\propto p({\boldsymbol{m}}| \mu ,{\boldsymbol{\Theta }},\delta ,{\boldsymbol{\alpha }})\\ &&\quad p(\mu ,{\boldsymbol{\Theta }})p(\delta | \mu ,{\boldsymbol{\alpha }}).\end{eqnarray} \tag{ 14 }$$

We already have Markov chain samples in distance, reddening, and stellar type for each star, for a model which does not take the line of sight reddening into account. We would like to apply weights to these samples, so that they correspond to the model sketched out directly above. As shown in detail in Appendix E, the correct reweighting of the stellar Markov chain samples, assuming a particular line of sight reddening profile ${\boldsymbol{\alpha }}$, is

$$\begin{eqnarray}&&{w}_{k}\propto \displaystyle \frac{{\mathcal{N}}({\delta }_{k}| 0,{\sigma }_{\delta })}{E({\mu }_{k};{\boldsymbol{\alpha }})},\end{eqnarray} \tag{ 15 }$$

where k indexes the sample, and ${\delta }_{k}$ is the fractional offset of the sample reddening, E_k, from the line of sight reddening, $E({\mu }_{k};{\boldsymbol{\alpha }})$, at distance ${\mu }_{k}$. The numerator in the above weight corresponds to a prior on the offset of the stellar reddening from the local dust column, while the denominator comes from the Jacobian transformation from reddening to fractional offset from the local reddening.

Reweighting each of the stellar parameter samples by the above factor, and normalizing the sum of the weights to unity, we obtain an inference for the stellar parameters, conditioned on a particular line of sight reddening profile, $E(\mu ;{\boldsymbol{\alpha }})$. We marginalize over the line of sight reddening profile by repeating this procedure for each sampled value of ${\boldsymbol{\alpha }}$, summing the weight applied to each stellar sample.

It might be objected that since we infer ${\boldsymbol{\alpha }}$ using the photometry of all the stars in the given pixel, the samples we drew for ${\boldsymbol{\alpha }}$ in our initial processing are already dependent on the photometry for the star whose samples we are reweighting. Put more formally, the prior we use here when marginalizing over line of sight reddenings is the inference we obtained earlier, $p({\boldsymbol{\alpha }}\;| \{{\boldsymbol{m}}\})$, which is conditional on all the photometry in the pixel, $\{{\boldsymbol{m}}\}$. We are using the photometry of a given star to infer the line of sight reddening, and then again to infer the stellar parameters, conditional on that line of sight reddening.

This would be a problem if the photometry of any given star significantly affected the line of sight reddening inference. However, if photometry from a large number of stars informs the line of sight reddening, then no one star should have a significant impact on the inferred line of sight reddening profile. In a more formally correct formulation of the problem, we would first infer the line of sight reddening using all the stars in the pixel except the star whose parameters we wish to infer, and we would then infer the parameters for that star, conditional on the inferred line of sight reddening profile. As each pixel contains hundreds of stars, however, we expect our procedure to approximate this formally correct procedure closely.

Figure 3 shows the result of reweighting the samples for one star. Depending on the line of sight reddening profile and the distribution of the unweighted stellar parameter samples, individual stellar inferences can be tightened dramatically by taking the line of sight reddening into account. Our knowledge of the parameters describing one star is therefore dependent not only on its photometry, but also on the photometry of its neighbors in the sky. Because neighboring stars lie along the same dust column, inferences for nearby stars are coupled through the requirement that the dust column increase with distance.

Zoom In Zoom Out Reset image size

PS1 is a 1.8 m optical and NIR telescope located on Mount Haleakala, Hawaii (Hodapp et al. 2004; Kaiser et al. 2010). The telescope is equipped with the GigaPixel Camera 1, consisting of an array of 60 CCD detectors, each 4800 pixels on a side (Onaka et al. 2008; Tonry & Onaka 2009). From 2010 May to 2014 April, the majority of the observing time was dedicated to a multi-epoch $3\pi$ steradian survey of the sky north of $\delta =-30^\circ$ (K. C. Chambers 2015, in preparation). The $3\pi$ survey observes in five passbands ${g}_{{\rm{P}}1}$, ${r}_{{\rm{P}}1}$, ${i}_{{\rm{P}}1}$, ${z}_{{\rm{P}}1}$, and ${y}_{{\rm{P}}1}$, similar to the SDSS (York et al. 2000), with the most significant difference being the replacement of the Sloan u band with a NIR band, ${y}_{{\rm{P}}1}$. The PS1 filter set spans 400–1000 nm (Stubbs et al. 2010). The images are processed by the PS1 Image Processing Pipeline (Magnier 2006), which performs automatic astrometry (Magnier et al. 2008) and photometry (Magnier 2007). The data is photometrically calibrated to better than 1% accuracy (Schlafly et al. 2012; Tonry et al. 2012). The $3\pi$ survey reaches typical single-exposure depths of 22 mag (AB) in ${g}_{{\rm{P}}1}$, 21.5 mag in ${r}_{{\rm{P}}1}$ and ${i}_{{\rm{P}}1}$, 20.8 mag in ${z}_{{\rm{P}}1}$, and 20 mag in ${y}_{{\rm{P}}1}$. The resulting homogeneous optical and NIR coverage of three quarters of the sky makes the Pan-STARRS1 data ideal for studies of the distribution of the Galaxy's dust.

2MASS is a uniform all-sky survey in three NIR bandpasses, J, H, and K_s (Skrutskie et al. 2006). The survey derives its name from the wavelength range covered by the longest-wavelength band, K_s, which lies in the longest-wavelength atmospheric window not severely affected by background thermal emission (Skrutskie et al. 2006). The survey was conducted from two 1.3 m telescopes, located at Mount Hopkins, Arizona and Cerro Tololo, Chile, in order to provide coverage for both the northern and southern skies, respectively. The focal plane of each telescope was equipped with three $256\times 256$ pixel arrays, with a pixel scale of $2^{\prime\prime} \times 2^{\prime\prime}$. Each field on the sky was covered six times, with dual 51 ms and 1.3 s exposures, achieving a $10\sigma$ point-source depth of approximately 15.8, 15.1 and 14.3 mag (Vega) in J, H and K_s, respectively. Calibration of the survey is considered accurate at the 0.02 mag level, with photometric uncertainties for bright sources below 0.03 mag (Skrutskie et al. 2006).

The addition of 2MASS data requires an expansion of the stellar model described in Green et al. (2014) to cover the 2MASS J, H, and K_s passbands, as well as an estimate of the survey selection function for 2MASS. We address these additions to our model in Appendices A and B.

We match each PS1 source to the nearest source in the 2MASS point-source catalog, rejecting stars at greater than 2'' separation. We require detection in at least four passbands, two of which must be PS1 passbands. In order to reject extended sources, we require that ${m}_{\mathrm{psf}}-{m}_{\mathrm{aperture}}\lt 0.1\;\mathrm{mag}$ in at least two PS1 passbands. We additionally reject sources flagged as extended sources in 2MASS.

Individual PS1 passbands are rejected if they have a photometric uncertainty greater than $0.2\;\mathrm{mag}$, or if the sources are beyond or close to the saturation limit for PS1 (here considered 14.5, 14.5, 14.5, 14, and $13\;\mathrm{mag}$ in ${{grizy}}_{{\rm{P}}1}$, respectively). In 2MASS passbands, we make the recommended "high-reliability catalog" selection cuts,⁸ in addition to the following requirements:

• 1.  
  ${\mathtt{contamination}}/{\mathtt{confusion}}\ {\mathtt{flag}}=0$,
• 2.  
  ${\mathtt{galaxy}}\ {\mathtt{contamination}}\ {\mathtt{flag}}=0$.

Our final catalog contains 798,611,689 sources, of which 32% are detected in four passbands, 49% in five passbands, and 19% in six or more passbands.

We divide the sky into HEALPix pixels (Gorski et al. 2005), adjusting the pixel scale in order to keep the number of stars per pixel roughly constant (see Figure 2). Our procedure is to begin with ${\mathtt{nside}}=64$ pixels, and then subdivide each pixel recursively, as long as the number of stars exceeds some threshold, dependent on the pixel scale. We use thresholds given in Table 1, chosen to allow us to reach a resolution of ${\mathtt{nside}}=512$ with a relatively small number of stars, but to avoid going to higher resolutions unless the stellar density is much higher. We reject pixels with fewer than ten stars. Such pixels comprise a negligible fraction of the sky. In all, we assign just under 800 million stars, covering just over three quarters of the sky, to 2.4 million pixels, with an average of 327 stars per pixel.

Table 1.  Pixelization

	Max.	Solid Angle	# Of
nside
	Stars Pixel−1	At Resolution (${\mathrm{deg}}^{2}$)	Pixels
64	200	77	91
128	250	90	430
256	300	11,980	228,373
512	800	16,071	1,225,471
1024	1200	2957	901,971
2048	⋯	66	80,956
Total	⋯	31,240	2,437,292

Download table as:  ASCIITypeset image

In Figures 4 and 5, we present the differential reddening in four spherical shells, centered on the Sun. Each panel shows the median dust reddening in a different range of solar-centric distances. Due to perspective, dust at high Galactic latitudes resides nearby, as Galactic dust lies in a thin disk. We recover the wealth of structure seen in the ISM across a wide range of scales, from thin filaments to large cloud complexes. Large, coherent cloud complexes appear at consistent distances.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 gives a closer view of the anticentral region (${\ell }\sim 180^\circ$). Different features appear clearly in each rendered distance slice. The Perseus, Taurus, and Auriga cloud complexes dominate the anticentral region in the closest distance slice, while the Orion molecular cloud complex (${\ell }\sim 210^\circ$, $b\sim -15^\circ$) and the California nebula (${\ell }\sim 160^\circ$, $b\sim -8^\circ$) appear very strikingly in the second distance slice. Of particular interest is the ring-like structure that Orion A and B appear to be embedded within. This ring-like structure is only apparent when the background dust is removed. In particular, the northeast portion of the ring is confused with the plane of the Galaxy in projection. Schlafly et al. (2014a) discusses the "Orion ring," including possible formation scenarios for the ring, in greater depth.

Zoom In Zoom Out Reset image size

Figure 7 shows the Galactic plane from ${\ell }=60^\circ$ to $155^\circ$ in more detail. The Cepheus flare, which lies at the center of the frame, at $95^\circ \lesssim {\ell }\lesssim 110^\circ$, $b\approx 15^\circ$, separates into two components at different distances, visible in the first and third panels. Using a modified version of the method used in this paper along individual sightlines, Schlafly et al. (2014b) places the two components of the Cepheus flare at distances of $360\pm 35\;\mathrm{pc}$ and $900\pm 90\;\mathrm{pc}$. The dust associated with Cygnus X (${\ell }\approx 80^\circ$, $b\approx 0^\circ$) appears in the third panel, along with a wealth of fine structure along the Galactic plane.

Zoom In Zoom Out Reset image size

We note that each pixel is fit independently, and our only prior assumption about the spatial structure of the dust is that if forms a diffuse disk, as shown in Figure 1. The detailed spatial structure in the interstellar medium that our analysis reveals indicates that the PS1 and 2MASS photometry dominates over our priors out to a distance of several kiloparsecs and reddening of $E(B-V)\approx 1.5\;\mathrm{mag}$. With the assumption of spatial correlations between neighboring pixels, we expect that one would be able to significantly reduce the uncertainty in the map, and achieve better distance resolution (see, e.g., Sale & Magorrian 2014).

Our 3D dust map is based on measurements of stellar distances and reddenings. Beyond the most distant stars, we have no sensitivity to dust, and in front of the nearest stars, we have no information about the distance to the dust. Therefore, we estimate the minimum and maximum distance to which our map is reliable by locating the nearest and farthest stars in each pixel. Outside of this distance range, our line of sight reddening inferences are dominated by our priors. Using the improved stellar parameter inferences (described in Section 2.4), we define the minimum reliable distance in each line of sight as the distance out to which there are ${N}_{\mathrm{closer}}$ observed main-sequence stars, and the maximum reliable distance as the distance beyond which there are ${N}_{\mathrm{farther}}$ observed main-sequence stars.

For this calculation, we count each Markov Chain sample in stellar distance as a fraction of a star. We exclude stars that fail to converge, for which the model is a very poor match to the data (as determined by the Bayesian evidence for the star; see Green et al. 2014), or which are inconsistent with the inferred line of sight reddening profile. We consider a star inconsistent with the line of sight reddening inference if none of the 100 stored Markov Chain samples of the stellar distance and reddening is within a fractional distance ${\sigma }_{\delta }$ (the modeled intra-pixel scatter in the dust column; see Section 2.4) of the line of sight reddening profile. Such objects are likely not well fit by any of our stellar templates, or alternatively signal that there is more variation in the dust column at fine angular scales than we allow.

In determining the minimum and maximum reliable distances in the 3D dust map, we use only main-sequence stars. This is because we consider our inferences for giants to be less reliable than our inferences for dwarfs. The colors and luminosities of giants depend more strongly on metallicity and age, the latter of which we do not model.

In the left two panels of Figure 8, we show the results obtained by requiring ${N}_{\mathrm{closer}}=2$ and ${N}_{\mathrm{farther}}=10$. The results are qualitatively similar for other choices of these parameters, as long as they are small compared to the typical number of stars in a pixel. The closest reliable distance is set almost entirely by the angular pixel scale. The nearby density of stars is, to a very rough estimation, uniform, meaning that the distance to the closest star in a pixel is a function primarily of the solid angle covered by the pixel. Accordingly, the top left panel of Figure 8 is essentially a map of pixel solid angle, with boundaries in distance following pixel scale boundaries.

Zoom In Zoom Out Reset image size

The farthest reliable distance of the 3D dust map is strongly influenced not only by the pixel scale, but also by the distribution of stars throughout the Galaxy, the 3D distribution of dust and the survey depth. The survey depth plays a larger role in the far limit because the most distant observed stars lie preferentially near the limiting survey magnitude. The closest stars, by contrast, are distributed more evenly in apparent magnitude. Accordingly, the high-latitude patches with particularly shallow maximum reliable distances (e.g., the shallow patch near ${\ell }=205^\circ$, $b=-55^\circ$) correspond to areas which are missing one or more PS1 bands, or which have fewer observation epochs in the PS1 $3\pi$ survey, and therefore worse coverage.

As a check on these estimated nearest and farthest reliable distance maps, we construct equivalent maps using simulated stellar catalogs. For each pixel in our map, we use our Galactic priors and the estimated limiting magnitudes for ${{grizy}}_{{\rm{P}}1}$ and 2MASS ${{JHK}}_{s}$ to generate an equal number of stars as actually observed. We then determine the distance in front of which there are ${N}_{\mathrm{closer}}$ and beyond which there are ${N}_{\mathrm{farther}}$ main-sequence stars. In order to remove the complicating factor of the line of sight dust inference, we restrict this test to $| b| \gt 15^\circ$ and to pixels for which $E{(B-V)}_{{\tau }_{353}}\lt 0.05\;\mathrm{mag}$, and assume for our simulated catalogs that there is zero dust extinction. In these comparisons, we find that a halo number density close to the value given in Jurić et al. (2008) is required to reproduce inferred map depths of Figure 8. The results for a halo strength of ${n}_{{\rm{h}}}=0.004$, binned down to ${\mathtt{nside}}=64$, are shown in the right two panels of Figure 8.

These results indicate that while we dialed down the halo strength in our priors, the data, as reflected in our photometric stellar inferences, nonetheless prefers a stronger halo. In future work, we will investigate the implications of our stellar and line of sight dust inferences for a global Galactic structure model.

As described in Section 2.4, we reweight the naively inferred parameters for each star in order to take the line of sight reddening profile into account. In order to demonstrate the improvement in stellar inferences after this re-weigthing procedure, we compare our reddening inferences to a set of independently determined stellar reddening standards, as in Green et al. (2014). For this, we use the set of stellar reddening standards developed by Schlafly & Finkbeiner (2011). The Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), part of SDSS-II, used moderate-resolution spectra to classify 240,000 stars. Empirically adjusted SDSS colors based on model atmospheres can then be compared with observed SDSS photometry to obtain reliable reddening estimates. We select the same sample of 200,000 SEGUE target stars as Schlafly & Finkbeiner (2011), which excludes white dwarfs and M dwarfs (the former because they are not contained in our model, and the latter because of their unreliable spectral classification).

In Figure 9, we compare reddening samples drawn from our PS1+2MASS-based stellar inferences to samples drawn from the SEGUE-based estimates. which have Gaussian uncertainties. The agreement between the two reddening estimates improves significantly after reweighting our PS1+2MASS-based stellar inferences. The improvement is most pronounced at low reddenings. The improvement is negligible for $E(B-V)\gtrsim 1\;\mathrm{mag}$, due to the fact that we allow the individual stars to deviate from the local line of sight reddening profile by about 10%. Unlike the SEGUE-based reddenings, the PS1+2MASS-based reddening inferences are constrained to be non-negative, leading to a negative slope in the residuals near zero reddening.

Zoom In Zoom Out Reset image size

It is possible to obtain a 2D dust map from our 3D map by projecting out distance, i.e., by taking the cumulative reddening out to some large distance. Such a map is of particular interest for extragalactic astronomy, where Milky Way dust is essentially a foreground screen to be removed. Several all-sky 2D maps of dust reddening already exist, among them Burstein & Heiles (1982), based on H i emission and galaxy number counts, SFD and two recent reddening maps derived from Planck data (Planck Collaboration et al. 2014), which model dust optical depth and temperature from FIR emission, and then calibrate a dust optical depth to reddening relation. NIR stellar colors have also been used to determine dust reddening more directly, as in the NICE/NICER/NICEST family of algorithms (Lada et al. 1994; Lombardi & Alves 2001; Lombardi 2009, respectively), and Rowles & Froebrich (2009).

Schlafly et al. (2014c) compares a 2D projection of an earlier version of our 3D dust map with the widely used SFD reddening map, as well as the newer Planck reddening maps. We repeat a number of the same tests for our new 3D dust map.

We begin, however, with a comparison between the map presented in this paper and the map presented in Schlafly et al. (2014c). The most important differences between the 3D maps used here and in Schlafly et al. (2014c) are the addition of NIR 2MASS photometry in our newer map, and that we sample here from the full posterior distribution on line of sight reddening, rather than finding the maximum-likelihood line of sight reddening. Because Schlafly et al. (2014c) requires that each star be detected in ${g}_{{\rm{P}}1}$, and because we incorporate NIR 2MASS photometry alongside PS1 photometry, our map reaches to deeper dust extinctions. This is apparent in the right panel of Figure 10, where our new map predicts more reddening both in the inner Galactic plane and in dense dust clouds off the plane, such as Orion, Taurus and Perseus.

Zoom In Zoom Out Reset image size

At reddenings below $E(B-V)\lesssim 0.08\;\mathrm{mag}$, our new map predicts significantly less reddening than Schlafly et al. (2014c), as can be seen in the left panel of Figure 10. Beyond $E(B-V)\approx 0.1\;\mathrm{mag}$, our reddening scales agree very closely, with our new map predicting $\sim 4\%$ more reddening. The scatter between the two reddening measures comes to $\sim 0.04\;\mathrm{mag}$, with a maximum median offset of $0.04\;\mathrm{mag}$ at $E(B-V)\approx 0.08\;\mathrm{mag}$, decreasing gradually to an offset of less than $0.01\;\mathrm{mag}$ between the two measures at $E(B-V)=1\;\mathrm{mag}$. The behavior of the residuals suggests that at very low reddenings, our new 3D dust map prefers essentially zero reddenings too strongly. This may be due to the stronger priors on dust reddening used in the present work, or due to slight differences in our new combined PS1-2MASS stellar locus, versus the PS1 stellar locus used in Schlafly et al. (2014c).

Next, we compare our inferred cumulative reddening out to 5 $\mathrm{kpc}$ with the SFD map, the Planck ${\tau }_{353\;\mathrm{GHz}}$-based reddening (hereafter, denoted simply as ${\tau }_{353}$, with units of magnitudes of $E(B-V)$), and the Planck radiance-based map (hereafter denoted ${\mathcal{R}}$, likewise with units of magnitudes of $E(B-V)$). All three of these maps are derived by modeling dust emission, and should therefore have different types of systematic errors than our stellar reddening-based dust map.

The upper panel of Figure 11 shows the SFD reddening across the footprint of our map, while the lower panel shows the residuals after subtracting off our 5 $\mathrm{kpc}$ map. Off the Galactic plane, the residuals are small, while larger residuals are found in the plane of the Galaxy, where there is significant dust past 5 $\mathrm{kpc}$, and where PS1 does not necessarily detect stars beyond all of the dust. Near the Galactic center, one noticeable feature is a blue halo, where we infer more dust than SFD. The residuals are correlated with dust reddening in this area, indicating a scale offset between the two maps, rather than a constant offset.

Zoom In Zoom Out Reset image size

The "blue halo" is again apparent if we compare our map to the Planck ${\tau }_{353}$ map. In Figure 12, we compare our inferred reddening at 5 $\mathrm{kpc}$ with SFD, ${\tau }_{353}$, and ${\mathcal{R}}$. The left panels map the residuals across the PS1 survey footprint on a square-root stretch, emphasizing small-amplitude differences. In the left panels, the Planck-based maps have been scaled by a constant factor to match our inferred PS1 reddening for $| b| \gt 20^\circ$.

Zoom In Zoom Out Reset image size

As the "blue halo" occurs in the direction of the nearby Aquila Rift, it is possible that the cloud has anomalous dust properties. For example, the grain size distribution or composition may vary, so that ${R}_{V}=3.1$ is not a good assumption in the Aquila Rift. Another possibile explanation for the "blue halo" is that our stellar inferences may be systematically biased toward greater reddenings in this direction due to limitations in our Galactic model. Our priors do not, for example, include a radial metallicity gradient in the disk components of the Galaxy, which could lead to biased metallicity estimates for stars toward the Galactic center. We leave this question for future investigation.

In the right panels of Figure 12, we plot the residuals of our inferred reddening to 5 $\mathrm{kpc}$ with SFD and the Planck emission-based maps as a function of reddening. Here, we do not scale the Planck maps by any factor to bring them into alignment with our map. To conduct this comparison, we compare the emission-based maps to multiple random realizations drawn from the posterior on reddening to 5 $\mathrm{kpc}$ from our 3D dust map. We restrict our comparison to high-Galactic-latitude regions ($| b| \gt 30^\circ$), and cut out the ecliptic plane ($| \beta | \lt 20^\circ$), where imperfectly subtracted Zodiacal light might contaminate the emission-based reddening maps. When comparing our PS1-based reddening with SFD, we place the Planck ${\tau }_{353\;\mathrm{GHz}}$-based reddening on the x-axis, as its errors are uncorrelated with both maps on the y-axis. When displaying the residuals between our PS1-based reddening and the Planck reddening maps, we place SFD along the x-axis for the same reason.

For reddenings above $E(B-V)\gtrsim 0.05\;\mathrm{mag}$, we see broadly similar residuals as found in Schlafly et al. (2014c), with different behaviors above and below $E(B-V)\approx 0.15\;\mathrm{mag}$. Above $E(B-V)\approx 0.15\;\mathrm{mag}$, our map predicts about 10% less reddening than SFD, but is in good agreement with ${\tau }_{353}$. As in Schlafly et al. (2014c), we find an overall difference in scale between our PS1-based map and the Planck ${\mathcal{R}}$-based map. For $E(B-V)\lesssim 0.05\;\mathrm{mag}$, we see the same residual between our map and the emission-based maps as found between our map and Schlafly et al. (2014c). For these small reddenings, our map favors essentially zero reddening too heavily.

The exact behavior of the residuals depends on which regions of the sky are masked in the analysis, indicating that there are spatially dependent systematic differences in the residuals between our reddening map and emission-based maps. However, the essential features remain the same, with different slopes in the residuals below and above $E(B-V)\approx 0.15\;\mathrm{mag}$, and our map favoring lower reddening below $E(B-V)\approx 0.05\;\mathrm{mag}$.

Marshall et al. (2006) developed a method to determine the reddening-distance relation along individual lines of sight by comparing 2MASS $J-{K}_{s}$ stellar colors to those of simulated catalogs based on the Besançon model of the Galaxy (Robin et al. 2003). Marshall et al. (2006) then applied this method to a regular grid of sightlines separated, by 15' covering the region $| {\ell }| \lt 100^\circ$, $| b| \lt 10^\circ$. The result is a 3D map of reddening in the inner Galaxy, extending to a maximum extinction of $\sim 1.4-3.75\;\mathrm{mag}$ in the 2MASS ${K}_{{\rm{s}}}$ band (equivalent to $\sim 4.5-12\;\mathrm{mag}$ in $E(B-V)$), and a maximum distance of $\sim 7\;\mathrm{kpc}$. Because Marshall et al. (2006) use only giants in their analysis, the dust map they produce has little information in the nearest kiloparsec. We will refer to this map as the "Marshall map."

Our dust map overlaps with the Marshall map in the approximate region $0^\circ \lt {\ell }\lt 100^\circ$, $| b| \lt 10^\circ$. Figure 13 shows the cumulative reddening at increasing distances in both the Marshall map and our 3D dust map. When converting from extinction to reddening, we assume ${A}_{K{\rm{s}}}=0.320\times E(B-V)$, as calculated by Yuan et al. (2013) for a 7000 K source spectrum at $E(B-V)=0.4\;\mathrm{mag}$, using the Cardelli et al. (1989) reddening law and assuming ${R}_{V}=3.1$. We mask our map beyond our predicted maximum reliable distance, as determined in Section 4.2. In the large-distance limit, our maps show good qualitative agreement outside of the masked areas. Figure 14 shows the differential reddening in the two maps in bins of increasing distance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the nearest two to three kiloparsecs, our map shows clearly differentiated structures at discrete distances that are spread over several distance bins in the Marshall dust map. For example, the Cygnus rift (located at ${\ell }\sim 80^\circ$, $b\sim 0^\circ$) appears clearly in our map in the distance bin spanning $0.5-1\;\mathrm{kpc}$, while in the Marshall map, it is spread over all the distance bins closer than $\sim 3\;\mathrm{kpc}$, and is therefore not cleanly separated from superimposed dust structures at greater distances. The greater distance resolution of our map at nearby distances is most likely due to the fact that we use both main-sequence stars and giants, while Marshall et al. (2006) relies solely on giants, which are saturated nearby and form a larger fraction of the observable stellar population at greater distances.

In addition to better distance resolution in the first few kiloparsecs, the greater source density of PS1 relative to 2MASS allows us to achieve better angular resolution than the Marshall map. Figure 15 demonstrates the difference in angular resolution between the two maps. In most of the region of overlap between our reddening map and the Marshall map, we achieve an angular resolution of $3\buildrel{\,\prime}\over{.} 4$, as opposed to the constant 15' resolution of the Marshall map. This allows us to resolve detailed filamentary structure not seen in the latter map. As dust reddening can vary significantly on small angular scales, this increased angular resolution will be important in practice for correctly de-reddening extinguished sources in regions with complex dust structure.

Zoom In Zoom Out Reset image size

Deep in the plane of the Galaxy, where high extinction reduces source counts, our finer angular resolution limits the distance to which our map can trace dust, compared to the Marshall map. Although the PS1 $3\pi$ survey is deeper than 2MASS, the 2MASS passbands are less affected by dust extinction, and the advantage of PS1 decreases in regions of high extinction. In such regions, the greater depth of PS1 does not fully compensate for our smaller pixels, limiting the depth to which we trace dust deep in the Galactic plane. The fact that our map derives its most accurate reddening information from main-sequence stars also limits the maximum extinction to which it is reliable.

Combining distance and reddening estimates for $\sim 23,000$ stars with the assumption of spatial correlation in dust density, Lallement et al. (2013) infer reddening in 3D out to a distance of $\sim 800\;\mathrm{pc}$ from the Sun. We find close morphological agreement between our 3D dust map and that of Lallement et al. (2013), with some differences which are worth taking note of.

In Figure 16, we show the distribution of dust and stars in a slice $25\;\mathrm{pc}$ above the Galactic plane, level with the Sun. We show the median dust density. In order to generate the stellar locations, we draw a sample at random from the improved distance posterior of each star (see Section 2.4), and select only those stars that lie within $5\;\mathrm{pc}$ of the chosen plane. For display purposes, we only display one out of every thousand stars.

Zoom In Zoom Out Reset image size

Just like Lallement et al. (2013), we find cavities in the dust density in the directions of ${\ell }=70^\circ$ and $225^\circ$. The overall morphology of the dust structure in the right planel of Figure 16 matches that of Figure 1 in Lallement et al. (2013).

The most obvious difference between our maps and those of Lallement et al. (2013) is the different voxel shapes employed in our work and theirs. Lallement et al. (2013) use small cubic voxels, and assume a spatial correlation function that favors similar dust densities in nearby voxels. This allows them to densely sample the reddening distribution, with voxels that are not directly constrained by stellar reddening measurements being constrained by neighboring voxels. In contrast, we infer the dust distribution in each line of sight separately, without assuming spatial correlations in the dust density, as laid out in Green et al. (2014). While our map has excellent angular resolution, it has distance bins with a width of about 25%, giving our voxels their pencil-beam shape.

Although we see roughly the same structures, such as cavities in the reddening distribution centered on ${\ell }=70^\circ$ and $225^\circ$, the distances we derive for a number of clouds is greater than the distances Lallement et al. (2013) find. In particular, while Lallement et al. (2013) place the Cygnus rift between 500 and $600\;\mathrm{pc}$, we place it at a distance of 800–$1000\;\mathrm{pc}$.

The basic data product our map contains is samples of the differential reddening in 31 distance bins, in each of 2.4 million pixels. The data structure is thus

$$\begin{eqnarray*}&&(2.4\times {10}^{6}\;\mathrm{pixels})\times (500\ \mathrm{samples})\times (31\ \mathrm{distance}\ \mathrm{bins}).\end{eqnarray*}$$

This allows us to determine the probability density of the cumulative reddening to any distance (within a few kiloparsecs) in the $3\pi$ steradians covered by the PS1 survey.

We encourage users to use the Markov chain samples of reddening versus distance provided by our interface in their analyses, as the samples contain the full statistical information generated by our method. These samples can be queried using our web API, downloaded as an ASCII table for individual lines of sight, or accessed directly in the complete data cube provided for download.

An additional, and larger, data set that we produce while generating the 3D dust map is a library of photometrically determined stellar parameters. As described in Green et al. (2014), we infer the probability distribution of the distance modulus, reddening, metallicity, and absolute ${r}_{{\rm{P}}1}$-magnitude for each star. We thus have a second data cube with shape

$$\begin{eqnarray*}&&(8\times {10}^{8}\ \mathrm{stars})\times (100\ \mathrm{samples})\times (4\ \mathrm{parameters}).\end{eqnarray*}$$

In addition, we store quality assurance information for each star, including whether the Markov Chain converged during the fitting procedure, and the Bayesian evidence for the stellar model, which is similar to the ${\chi }^{2}$ statistic in maximum-likelihood fitting. Point sources with poor evidence are likely either of stellar types not contained in our model, such as very young stars or binary systems, or are not stars (e.g., white dwarfs, quasars, unresolved galaxies).