A FULL-SKY, HIGH-RESOLUTION ATLAS OF GALACTIC 12 μm DUST EMISSION WITH WISE

The WISE satellite followed a 95 minute, Sun-synchronous, polar orbit about Earth, always pointing close to zenith and imaging at nearly 90° solar elongation. Over a single orbit, WISE scanned all ecliptic latitudes β at a particular ecliptic longitude λ, and then all β at longitude (λ + 180°). Successive 8.8 s, 08 frames overlap by ∼10% in the scan direction, while corresponding frames in successive orbits have ∼90% overlap.

Because the All-Sky Release spans just under seven months, the vast majority of the sky at low ecliptic latitude has a resulting integer coverage of 10–12. For our purpose of mapping diffuse emission, this substantial redundancy is useful in compensating for single-frame artifacts such as cosmic rays and satellite streaks. However, for typical locations on the sky, all 10–12 WISE epochs occur within a single ∼20 hr time period. Unfortunately, the Moon and bright solar system objects do not move far enough on this timescale to be rejected.

All WISE data products utilized throughout this work, including single exposures and catalogs, are drawn from the All-Sky Data Release (Cutri et al. 2012). The All-Sky Release features two distinct options for accessing WISE images. The "Atlas images" are a set of 18,240 co-adds, 15 on a side. While it might seem convenient for these co-adds to form the basis of our 12 μm dust map, various considerations led us to conclude that we were best served by starting our analysis at the single-exposure level, then assembling our own custom stacks. First, each Atlas co-add is designed only to minimize the appearance of exposure boundaries within its footprint, but will not necessarily agree with its neighboring Atlas tiles on their overlaps. Additionally, the Atlas images contain all compact sources observed by WISE, which require extremely careful PSF subtraction in order to isolate the diffuse cirrus. The Atlas PSF must be defined in a complex way and would be difficult both to model and apply, given the multiple epochs combined at each Atlas pixel and the fact that the Atlas images have been smoothed.

Because of these Atlas co-add considerations, we have instead obtained WISE imaging data by downloading all ∼1.5 million W3 exposure sets, including the "intensity," "mask," and "uncertainty" files. The single exposures are referred to as "Level 1b" (L1b) images, are archived by the NEOWISE project (Mainzer et al. 2011), and are the lowest-level publicly available WISE imaging data. L1b images are square, 1016 pixels on a side, and have a platescale of 275 pixel⁻¹. The x (y) axis of each L1b exposure is oriented perpendicular (antiparallel) to the scan direction.

To remove compact sources, we also make use of the WISE All-Sky Source Catalog (Section 3) and WISE All-Sky Known Solar System Object Possible Association List (Section 5.3).

We discard a small subset of L1b exposures which, based on metadata alone, we know in advance will not be worthy of further processing. Our cuts generally follow the recommendations and cautions of Cutri et al. (2012). Specifically, we reject frames with saa_sep ⩽−5°, or moon_sep ⩽12°, or w3msknumsat > 500,000 pixels, or with dtanneal ⩽ 1000 s at high ecliptic latitude. These cuts remove 4.9% of the L1b images.

The WISE All-Sky Source Catalog was engineered to emphasize completeness over purity. Consequently, a striking feature of the W3-detected sample (w3snr > 2) is the tendency for multitudes of dubious, faint sources to trace any sufficiently bright, filamentary dust structure (see Figure 1). Leaving this issue unaddressed, the filaments in PSF-subtracted images take on a distinctly "chewed up" appearance, compromising our estimate of the dust cirrus.

Zoom In Zoom Out Reset image size

We were unable to identify any trivial cuts on e.g., signal-to-noise or source color that rectified this problem. However, we observed that the sources clustered on filaments tend to have a curve of growth (COG) more similar to that of a constant surface brightness profile than that of a PSF (see Figure 2). We measured the COG via a series of aperture magnitudes, w3mag_n, for n = 1, ..., 6, available in the All-Sky Catalog. The PSF COG was computed by integrating our model from Section 4 in the relevant apertures. We were able to satisfactorily reject nebula fragments by cutting those sources with COG more consistent with constant surface brightness than with the W3 PSF model. We also discarded sources with an undefined error w3sigm_1 on their innermost aperture magnitude, w3mag_1. The PSF versus constant surface brightness comparison was calculated in the sense of minimum absolute deviation, as the majority of All-Sky Catalog errorbars (w3sigm_n) were undefined for the faint sources of interest, precluding a formal likelihood analysis.

Zoom In Zoom Out Reset image size

The objects thus removed constituted a non-negligible fraction of All-Sky Catalog sources with w3snr > 2. Prior to the COG cut, the All-Sky Catalog contained ∼167 million sources detected in W3; applying the COG filter reduced the number of remaining W3 sources to ∼41 million.

By employing our detailed model of the PSF wings (see Section 4), we were able to derive custom bright source parameters that yielded cosmetically optimal PSF subtractions. Our procedure ignores the saturated core of each bright source, fitting only an annulus containing the wings.

Specifically, we fit custom fluxes for those ∼15k sources in the All-Sky Source Catalog brighter than w3mpro = 2 (>5 Jy). This threshold has no special significance with regard to the WISE instrumentation, and its choice is driven by pragmatic concerns relating to e.g., our available computational resources. To fit each bright source, we consider an annulus about its centroid with outer radius 21. The inner boundary of the annulus was determined via a flood-fill of pixels near saturation (>20,000 DN), starting from the centermost pixel. After dilating this saturation mask by 3 L1b pixels, we excluded the resulting interior region from the flux fit. We then perform a linear least squares fit to the annulus pixel values. The model is the sum of a rescaled PSF and a constant offset representing the background level. Thus, there are two parameters: the background offset and PSF rescaling factor, from which we infer the source flux. For each source, such a fit is performed for each L1b image in which the entire fitting annulus falls within the image boundaries. The L1b pixels of each fit are inverse variance weighted based on the corresponding L1b uncertainty image. The flux assigned is taken to be the median of the individual L1b fluxes. Typically 9 L1b fits per source are performed.

We also computed custom centroids for bright sources in order to achieve cosmetically optimal PSF subtractions. After computing our custom fluxes, we recompute the centroid by fixing the background and flux, but fitting x, y offsets via Levenberg–Marquardt chi-squared minimization.

Figure 3 shows a representative sample of 25 bright sources drawn from ∼300 deg² of the sky which we processed in turn with the raw All-Sky Catalog parameters and our custom parameters. On the whole, our custom parameters clearly provide substantial improvement in the monopole and dipole components of the residuals. For the purposes of this work, we only care about the cosmetic quality of our point source subtracted maps, and have therefore chosen bright source subtractions based on our custom parameters. We caution, however, that our custom bright source parameters have not been thoroughly vetted against other mid-infrared calibrators, whereas those of the official All-Sky Source Catalog have been (Cutri et al. 2012).

Zoom In Zoom Out Reset image size

We cannot rely upon one single sample of WISE sources to model all portions of the W3 PSF. For instance, the extremely bright sources that have high signal-to-noise in their outer wings will be saturated in their cores, and thus useless in modeling the innermost several FWHM.

W3 sources begin to saturate in their very central pixels at 3.5 < w3mpro < 4.0. In order to leverage the brightest non-saturated sources, we modeled the PSF's inner 115 (10.6 FWHM, the "core") using ∼4600 All-Sky Catalog sources with 4 < w3mpro < 5 and |b| > 40°.

To fully reap the benefits of PSF subtraction at the single-exposure level, our model must include spatial PSF variation across the chip, as a function of source centroid (x_c, y_c). We begin by gathering cutouts from individual L1b images in which each source in our sample appears. These cutouts include those in which the PSF core is only partially contained within the image. This permits us to model the PSF variation all the way to the chip's edges. Typically 10 such cutouts per source were collected. We sinc shift each cutout in both x and y directions by a fraction of a pixel such that all cutouts are precisely centered on the source centroid. Additionally, each L1b cutout is background subtracted and rescaled to an arbitrary common source flux. The (x_c, y_c) location of each cutout's centroid within its original L1b image is also stored.

For each pixel location with respect to the PSF model centroid, the set of rescaled pixel values available from these cutouts is fit with a third order polynomial surface in (x_c, y_c). We chose third order because our measurements of the PSF FWHM as a function of position in the chip revealed a dependence that required third order terms. The per-pixel fit is performed independent of other neighboring pixels and carried out via weighted linear least squares, including iterative outlier rejection. The weights are inverse variances based on propagation of errors from the L1b uncertainty masks. The result is a polynomial model of the PSF core wherein the zeroth order term represents the PSF at the center of the chip, and higher order terms correct for PSF variation across the chip.

Our procedure for modeling the PSF wings is quite similar to that employed in modeling the PSF core (Section 4.1). Our model for the PSF wings consists of a rectangular region 149 on a side, and thus characterizes the W3 PSF behavior over a diameter of ∼140 FWHM.

In assembling a sample of sources with which to model the PSF wings, we aimed to use approximately the ∼100th–1100th brightest high-latitude sources on the sky. Specifically, our sample was defined to be those sources with −2 < w3mag < 1, and with |b| >15°. This sample contains 1046 unique sources. We found that the handful of w3mag < −2 sources were predisposed toward anomalous defects such as bleeding, and thus would contribute systematic problems to the PSF model despite having nominally high signal-to-noise in the wings.

We again collect cutouts as in Section 4.1, but model the per-pixel spatial variation with only first order terms in x_c, y_c. Higher order terms did not appear necessary and tended to destabilize the per-pixel polynomial fit in outer regions of the wings where signal-to-noise is low. The resulting model of the PSF core and wings, tapered according to Section 4.5, is shown on a logarithmic stretch in Figure 5. The diffraction spikes are not treated specially relative to other pixels in the model of the PSF wings.

Zoom In Zoom Out Reset image size

The optical ghost location with respect to the PSF core depends strongly on the centroid coordinates of its parent source, particularly x_c. As x_c ranges from one edge of the image to the other, x_c ≈ 0 to x_c ≈ 1015, the ghost centroid shifts in position relative to x_c by 05 (∼11 L1b pixels). We initially attempted to model this translation of the ghost relative to x_c with an approach similar to that in Sections 4.1 and 4.2, employing high-order per-pixel polynomials. However, we found this model did not satisfactorily capture the ghost's translation and became unstable for x_c near the L1b edges.

We thus resorted to a ghost model in which the ghost shape remains constant, but can translate with respect to x_c and can be amplified/suppressed by an overall rescaling factor f_g that varies with (x_c, y_c). To create such a model, we again started by extracting many L1b cutouts. The same 1046 sources used to model the PSF wings were also used to model the ghost. We began by computing the centroid of the ghost in each cutout, and fit a second-order polynomial model to the ghost's offset relative to (x_c, y_c). We then used this ghost translation model to rectify all the cutouts to a common fiducial ghost centroid with respect to (x_c, y_c). The rectified cutouts allowed us to compute the ghost profile via a weighted average with iterative outlier rejection. To determine f_g, we modeled each rectified cutout as a sum of an offset and f_g times the ghost profile. Finally, we fit a second-order polynomial in (x_c, y_c) to the set of f_g values we obtained. With f_g = 1 for (x_c, y_c) at the chip center, our resulting model always has 0.85 < f_g < 1.06.

4.4.1. First Latent

4.4.2. Second Latent

Our models for the PSF wings and ghost, as well as the latents, are defined on rectangular footprints, and may not necessarily asymptote identically to zero precisely at their very edges. Thus, subtracting our rectangular models would introduce extra power in modes oriented parallel to the x, y axes of the chip. For this reason we taper the aforementioned models to zero before performing subtractions. This is achieved by premultiplying the models by weights that decline smoothly from unity in each model's center to zero at each model's boundary. We denote the angular extent of a particular model in the x (y) direction θ_x (θ_y). In order to define our tapering weight we first compute a radial coordinate r relative to the model's central pixel:

$$\begin{equation} r = \sqrt{(\theta _y/\theta _x)^2 \Delta x^{2} + \Delta y^{2}}, \end{equation} \tag{ 1 }$$

where Δx and Δy are displacements relative to the model centroid. For r ⩽ r_in, the tapering weight is unity, and for r ⩾ r_out the tapering weight is zero. For r_in < r < r_out, the weight transitions from 1 to 0 linearly with respect to r. Table 1 lists θ_x, θ_y, r_in, r_out for the models that we taper prior to subtraction.

Table 1. Model Sizes and Tapering Parameters

Model	θx	θy	rin	rout
PSF wings (bright)	149	149	67	74
PSF wings (faint)	53	53	19	26
Ghost	62	40	17	20
First latent	149	149	67	74
Second latent	124	124	55	62

Download table as:  ASCIITypeset image

In order to relate magnitudes (w3mag) to L1b counts (DN), we must accurately determine the "zero-point magnitude" m₀, such that a source with integrated flux of 1 count corresponds to w3mpro = m₀.

We chose to define our m₀ with respect to an aperture of diameter 18. We then conducted aperture photometry on the same set of high-latitude, unsaturated, moderately bright sources used to model the PSF core. In order to determine m₀, we performed a least-squares fit with iterative outlier rejection to the trend of w3mpro versus 2.5log₁₀(F_DN), F_DN signifying the total counts within the 18 aperture. The result determined m₀ = 17.645. Reassuringly, the best-fit slope was indeed very nearly unity. More generally, we found no evidence for significant nonlinearity of w3mpro for unsaturated sources at any point throughout this work.

We found that m₀ was sufficiently stable with time throughout the All-Sky Release portion of the WISE mission to be considered a constant. Similarly, we did not find any evidence for temporal variation of the PSF shape.

Each raw L1b intensity image is accompanied by a corresponding *-w3-msk-1b.fits bit-mask provided by the WISE team. We will refer to such files as msk files, so as not to confuse them with our own custom L1b masks. The first step in our L1b processing is to create a mask of bad pixels in the raw L1b image which we must interpolate over, based on the msk file. We dilate the cosmic ray mask by a kernel of sidelength 3 pixels. Bad pixels are then taken to be those flagged in the dilated cosmic ray mask, or with any of the static bad pixel msk bits set, or with a NaN value. We immediately interpolate over this mask before proceeding.

Next we move on to correcting for static pointlike sources. From our custom catalog, we retrieve a list of all those sources whose PSF wings or ghost image, as represented by our models, have non-zero overlap with the L1b image at hand. This includes objects with centroids that fall outside the boundaries of the L1b image.

For each point source in the list, we compute its centroid coordinates (x_c, y_c) based on the L1b WCS and create a PSF model including the core, wings and ghost. The core and wings are computed from the per-pixel polynomials in (x_c, y_c) described in Sections 4.1 and 4.2. The ghost image is derived based on the model of Section 4.3, taking (x_c, y_c) into account through the ghost centroid translation and the amplification factor f_g. When computing these models, we always bound x_c, y_c to lie within the image, so as to not extrapolate beyond the limitations of our various polynomial corrections in x_c, y_c. We renormalize the PSF according to each source's w3mag value and the prescription of Section 4.6.

Before finally subtracting this PSF model from the L1b image, two further steps are taken: (1) we sinc shift appropriately such that the observed and model source centroids match at the sub-pixel level (2) we apply the tapering prescription of Section 4.5 so as not to introduce rectangular boundaries into the model-subtracted images. These two steps are also applied before subtracting the first and second latent models as described in Sections 5.2.1 and 5.2.2 below.

It was not necessary to apply the full 149 model of the PSF wings to arbitrarily faint sources. Rather, for sources deemed faint, we subtracted a 53 cutout of the PSF core/wings. This smaller model does not include the ghost. We set the threshold distinguishing bright versus faint at w3mag =4.2 for typical coverage, though the exact cutoff was scaled on a source-by-source basis to account for coverage w3m.

For sufficiently bright sources, the innermost pixels of the PSF subtraction tended to show strong residuals. Consequently, we interpolate over such regions, using a circular mask of size that varies with w3mag. Such interpolation only occurs for sources with w3mpro ⩽ 9. We flagged pixels interpolated over in this manner with bit 4 of our mask. These interpolations affect a very small fraction, <0.25%, of the sky.

Based on the appropriately shaped and scaled PSF model derived for each source, a corresponding bit-mask is created. Saturated pixels within the PSF core (>30,000 DN) are flagged with bit 0. Bright regions of the PSF profile (>250 DN) are flagged with bit 2. Pixels within the ghost model that have values greater than 10⁻⁴ times the peak value of the PSF core are flagged with bit 1. This mask bit is meant to trace the general outline of the entire ghost model. Especially bright pixels within the ghost model (>20 DN) are flagged with bit 5.

5.2.1. First Latent

5.2.2. Second Latent

5.2.3. Third, Fourth Latents

5.2.4. Diffraction Spike Mask

To some extent, the outlier rejection steps in our mosaicking process (Section 7) eliminate moving objects, without any special treatment of SSOs in particular. However, there are many cases in which moving objects proceed slowly enough so as to contaminate the same location on the celestial sphere at multiple epochs. Furthermore, the WISE team provides a convenient catalog of potential SSO detections—the WISE All-Sky Known Solar System Object Possible Association List (Cutri et al. 2012, henceforth PAL). For these reasons, we chose to incorporate steps for masking and flagging SSOs into our L1b processing pipeline.

The PAL uses the orbits of known SSOs to calculate the expected locations of these objects within each L1b exposure. If there exists a WISE detection of some source nearby to one of these predicted locations, the measured properties of that source will be recorded in the PAL. Importantly, each PAL entry contains an L1b frame identification, coordinates of the source within that L1b image, and the source's profile-fit magnitude, again designated w3mpro. In some cases, the associations may be inertial sources, or spurious detections, as PAL entries are based on single-epoch photometry. For our purposes, sample purity is not of much concern. Even if we interpolate unnecessarily in some cases, the redundant WISE coverage will tend to render such a mistake imperceptible.

For each SSO in a given L1b exposure, we produce a mask centered about the SSO centroid that scales in size with the detected object's brightness. These masks each contribute to bit 6 of our L1b bit-mask (see Table 2), which denotes the pixels we have interpolated over within the L1b image to remove the SSO core/wings. The mask about an SSO's centroid consists predominantly of a circular region, with diameter varying from 03 (7 L1b pixels) for w3mpro ⩾ 10 to 35 (83 L1b pixels) for w3mpro ⩽ 1.4. In addition, we apply our diffraction spike model of Section 5.2.4 to SSOs with w3mpro < 3. We interpolate over these SSO diffraction spikes at the L1b level, and the affected pixels are also flagged with bit 6. We considered PSF subtracting SSOs rather than simply interpolating over them, but found that the single-epoch PAL photometry was too unreliable and yielded a cosmetically inferior result.

In addition to addressing each SSO's core, wings and diffraction spikes, we take further actions with regard to artifacts of very bright SSOs. For SSOs with w3mpro < 2.5, we mask and interpolate over an elliptical region surrounding the object's ghost image. SSO ghost interpolation is registered in our bit-mask with bit 11. Similarly, for SSOs brighter than w3mpro = 1, we interpolate over the first latent and correspondingly set mask bit 12. Figure 6 provides an example of our SSO masking and interpolation results.

Zoom In Zoom Out Reset image size

Table 2 lists all mask bits included in the WSSA tiles of our data release. Throughout Section 5, bits 0–6, 8–13 have been explained. Bit 14 flags pixels that are saturated (>30,000 DN), but not marked as suspect by any of the msk file static mask bits. This can be useful downstream in identifying regions of extremely bright nebulosity that have reached saturation. Bits 7, 15–22 are only defined at the co-add level and are explained in Section 6.5.

Before writing our processed L1b files, we smooth to 15'' FWHM, trim the smoothed image to 1000 × 1000, and bin down by a factor of two in each dimension. We write the outputs of our L1b processing pipeline as multi-extension FITS files. For each raw L1b file in the All-Sky Release we write a corresponding processed file. The primary extension is the fully "cleaned" image, where we have removed point sources and their artifacts and SSOs. Extension 1 is the "dirty" image, which is merely a smoothed and rebinned version of the raw L1b image after interpolating as described in Section 5.1. Extension 2 is our bit-mask of Section 5.4, after dilation and rebinning that mimics the smoothing and rebinning of the first two extensions. These processed L1b images are not included in our data release, but will be made available upon request.

Overlapping L1b images may differ by more than a simple relative offset. This is especially true in regions that were imaged six months apart with opposite scan directions or near the ecliptic poles. In such cases, the zodiacal light can differ substantially between epochs. Even though the zodiacal light at a given epoch generally varies on spatial scales larger than a degree, this foreground is so strong that its gradient is readily apparent in single L1b exposures. As such, to ensure that our stacking can match neighboring exposures with a single offset, it is important to compute a correction for each exposure that removes time-dependent first order structure.

To arrive at such a gradient correction, we compare each L1b exposure to SFD 100 μm on the same footprint. For each cleaned L1b exposure, we interpolate generously over bright source artifacts using our masks, then smooth to 61 FWHM to match i100 resolution and bin down 5× in each dimension. We denote each pixel with coordinate x_i, y_i in this smoothed, rebinned W3 image as w_i. We then fit the following model to w_i:

$$\begin{equation} w_i = f_{wi}i_i + \nabla _x x_i + \nabla _y y_i + C, \end{equation} \tag{ 2 }$$

where f_wi, ∇_x, ∇_y, and C are the scalar parameters to be fit and i_i are the i100 values corresponding to w_i. The fit is performed via unweighted linear least squares, with iterative outlier rejection. ∇_x, ∇_y then specify x, y gradients that make the processed L1b exposure most consistent with Galactic emission, free of slowly varying time-dependent foregrounds. f_wi is a factor that allows small-scale dust structure to be nulled, accounting for the fact that conversion factor between W3 and i100 will be different for different L1b pointings. Similarly, C absorbs any constant offset arising from the WISE and i100 zero points. Our model for the time-dependent gradient in each exposure is computed from ∇_x, ∇_y, and subtracted before the mosaicking process begins. We only desire the gradient terms, and make no use of f_wi, C downstream.

Despite our metadata cuts (Section 2.3), some L1b images with severe systematic problems remain and could potentially corrupt the stacking process if allowed to contribute to our co-adds. To remove such images, we employ a pairwise comparison between overlapping L1b exposures.

Specifically, each gradient-corrected, processed L1b image is compared to up to 20 of its nearest neighbors at the locations of N_side = 2048 HEALPix pixel centers on which the frames overlap. For each pair, we compute the rms of the differences between the sampled values. Then, for each image, we compute the median of these rms values, $\mathcal {R}_{{\rm med}}$. We then reject the highest $\mathcal {R}_{{\rm med}}$ frame with $\mathcal {R}_{{\rm med}} > 500$ DN, updating $\mathcal {R}_{{\rm med}}$ for all other frames after each rejection, until no frames remain with $\mathcal {R}_{{\rm med}} > 500$ DN. ∼0.4% of images are rejected by the pairwise comparison, but these frames would have inflicted a disproportionate, negative influence on our stacks given the severity of their systematic problems.

For each WSSA tile, we co-add the constituent L1b exposures onto the relevant astrometric footprint, with pixel scale of 5625 pixel⁻¹. 125 × 125 co-adds are thus 8000 pixels on a side. Such tiles are well-sampled because we have smoothed our processed L1b images to 15'' FWHM.

We obtain a list of L1b exposures that will overlap the tile's astrometric footprint, sorted by proximity to the tile center. We apply the gradients of Section 6.1 to each cleaned L1b image, and then use the pairwise rejection step to reject the typically small number of problematic exposures. Starting with the exposure closest to the center, we add the cleaned L1b exposures into the stack one at a time, applying a single offset to make each exposure best agree with the existing stack. This offset is simply the median value of the difference between the exposure at hand and the existing stack on their overlap. As we stack, we also keep track of the minimum and maximum single-exposure pixel values at each pixel in the co-add footprint. We also store an integer coverage map, incrementing the relevant pixels each time an L1b image is added to the stack. For locations in the completed co-add that have coverage >2, we subtract out the average of the minimum and maximum images at those locations, and adjust the coverage map appropriately. We will refer to the average of the minimum and maximum images as the "transient artifact image." This outlier rejection step eliminates satellite streaks and any other unusual single-epoch features not specifically addressed by our L1b processing pipeline. We also apply the per-frame offsets to the dirty extension of our L1b images to create a corresponding stack with all point sources, bright source defects, and SSOs present.

With our initial WSSA tile stacks in hand, we can perform some further refinements, using the existing stack as a reference. There are two typically low-level problems we address by making use of our reference stack: (1) per-quadrant polynomial offsets in individual L1b exposures and (2) striping artifacts in L1b images.

6.4.1. Per-quadrant Polynomial Warps

We have already addressed scalar intensity offsets and gradient corrections on a per-exposure basis. But within a given exposure, the four quadrants can be offset relative to each other due to four-amplifier detector readout. These offsets can be simple scalars, but in some cases involve higher order polynomials (e.g., a gradient across one quadrant that is absent in the other quadrants). We refer to these polynomial offsets as per-quadrant "warps." For fixed position on the celestial sphere within our co-adds, such warps in the various contributing L1b quadrants are not correlated. As such, any warping will be suppressed by ≳90% in our initial stacks.

Thus, we may reasonably consider our initial stacks to be nearly free from the effect of such per-quadrant warping, and as a result use them as reference templates to correct individual L1b quadrants. Doing so can help eliminate any jagged L1b quadrant edges that persist at low levels in the initial stacks. For each tile, we gather the same sample of L1b exposures as in Section 6.3. We then rebuild the mosaic from scratch, one L1b quadrant at a time, first fitting a per-quadrant warp relative to the reference, and then subtracting the warp model before adding the corrected quadrant into the revised stack.

Our warp model is simply a polynomial in the relevant cleaned L1b quadrant's x, y coordinates. For each such quadrant, we fit a fourth order polynomial to the difference between the cleaned L1b image and the reference stack. The fit is performed via unweighted linear least squares. With the best-fit polynomial computed, we evaluate the χ² goodness of fit between the reference stack and the corrected L1b quadrant. We take the expected pixel variance in each pixel of the cleaned L1b image to be:

$$\begin{equation} \sigma _i^2 = f_{s}^2\left(R^2 + \frac{N_{i}}{g}\right) + (\delta g N_{i})^2 + (\nabla N|_{i} \delta _{\theta })^2, \end{equation} \tag{ 3 }$$

where N_i is the cleaned L1b pixel value in DN, R = 16.94 DN is the readnoise, g = 6.83 e⁻/DN is the gain, f_s is a constant that accounts for our smoothing/rebinning, δg is an estimate of the fractional gain variation, ∇N|_i is the magnitude of the reference image gradient evaluated at pixel i, and δ_θ is an estimate of the L1b astrometric error. N_i is floored at zero in the second term, since it is possible that e.g., a bright source residual could have negative values. ∇N|_i is calculated based on the reference stack, not the L1b image being compared to the stack. Based on studies of point sources presented in Cutri et al. (2012), we chose δg = 0.025 and δ_θ = 04. The first two terms account for statistical errors in the cleaned L1b image. The second term allows for the possibility that the per-quadrant L1b gain may vary slightly. The final term accounts for astrometric imperfections that could cause the L1b exposure to be out of alignment with the reference stack. We assume that statistical errors in the reference stack are negligible relative to those in the L1b quadrant of interest.

If the χ²/dof is greater than 7 for a particular L1b quadrant after applying the best-fit polynomial warp, that quadrant is excluded from the stack. 2% of L1b quadrants were thus rejected. Figure 7 shows a region in which L1b boundaries in the initial stack of Section 6.3 are corrected by our quadrant warping procedure.

Zoom In Zoom Out Reset image size

6.4.2. Destriping

L1b exposures display low-level striping. Complicating matters, this striping occurs with different amplitudes and on different wave numbers from exposure to exposure, and even from quadrant to quadrant within single exposures. Still, a particular set of relatively few modes tend to be enhanced, as the stripes within a given quadrant are typically oriented parallel to the scan direction and have wavelengths of ∼25–4'. To avoid high-pass filtering diffuse structure, our destriping never modifies modes larger than 46. We assume the striping to be an additive rather than multiplicative defect.

As such, our first step toward correcting L1b striping is to make a mask in Fourier space of the typically contaminated modes, using ∼20,000 cleaned L1b images at very high |b|, where there are few compact sources and diffuse structure is minimal. For each of these L1b exposures, we compute the isotropic Fourier power spectrum in each quadrant. We label any mode "contaminated" which has a modulus-squared that is a >5σ outlier relative to the appropriate isotropic value. We construct a mask of frequently contaminated modes by counting the number of occurrences of such outliers at each pixel within the quadrant fast Fourier transform (FFT).

For each quadrant that successfully passes the quadrant warping step of Section 6.4.1, we attempt a destriping correction before finally adding the quadrant into the reference-corrected stack. We take the FFT of the quadrant under consideration and identify any >5σ outliers relative to the isotropic power spectrum. We then modify the FFT values of those outliers which coincide with modes we previously flagged as frequently contaminated. We rescale such outliers to have a modulus that matches the appropriate isotropic value, leaving the phase unchanged. Transforming back to the spatial domain yields a destriped image.

For most L1b quadrants, this correction removes all apparent striping artifacts. However, there are cases in which the correction can go awry. For example, the residuals of bright source artifacts can be aligned with the scan direction, wreaking havoc with our analysis meant to deal with low-level striping of the same orientation. We therefore incorporated destriping into our reference comparison process, in order to robustly identify such problematic instances. Specifically, we perform a chi-squared comparison between each L1b quadrant and the reference stack, both before and after destriping. We use the pixel variances of Section 6.4.1. We then apply the destriping correction only if (1) the destriping correction improves the χ²/dof and (2) more individual pixels within the quadrant have their χ² improved than worsened by destriping. Since the striping is most apparent in regions of very low signal/background, and the destriping correction is very successful in such regions, we find this "do no harm" approach to be satisfactory. Figure 8 shows an example of our destriping correction.

Zoom In Zoom Out Reset image size

The following subsections describe mask bits defined at the mosaic level, rather than the L1b level.

6.5.1. Moon Mask

While our metadata cuts remove L1b frames pointed near the Moon (moon_sep ⩽ 12°), off-axis Moon-glow can adversely affect exposures with much larger moon_sep values. In fact, prominent off-axis Moon illumination features extend out to many tens of degrees away from the Moon's location. Accounting for the WISE survey strategy, this illumination detrimentally contaminates a non-negligible fraction, ∼5%–10%, of the sky. Adverse effects on our stacks arise because the extended Moon-glow violates our assumption that the differences between exposures can be accounted for by an offset and linear gradients. For example, Moon-glow diffraction spikes can imprint strong, higher-order structure within a single L1b frame (see Cutri et al. 2012, Section 2.4b).

We derived an empirical model of W3 off-axis Moon illumination in order to flag portions of our stacks which suffer serious Moon-glow corruption. In ecliptic coordinates, Moon-contaminated regions trace out a series of narrow vertical bands, a few degrees wide and generally separated by ∼30° in λ. We did not attempt to characterize sub-exposure Moon-glow structure, but instead created a model of the additive offset to the WISE zero point due to Moon illumination. Our model is not a function of (λ, β), but rather Moon angular separation (moon_sep) and Moon position angle relative to the WISE optics (termed ϕ_m = moonpa−w3pa). To arrive at such a model, we generated a full-sky W3 map, taking the L1b zero point at face value. We then subtracted 1.2× the Kelsall model to roughly cancel the zodiacal light. Next, we smoothed to 61, and subtracted scaled SFD i100 assuming i100/W3 = 17.5, to approximately cancel diffuse Galactic emission. We designated particular, narrow ranges of λ as Moon-affected, and created a Moon-free prediction by interpolating over such regions at constant β. Restricting to |b| > 15°, the difference between the Moon-free prediction and the corresponding L1b data constitutes our Moon-glow estimate, which yields the profile shown in Figure 9 when binned in (moon_sep, ϕ_m). We also created analogous profiles in Moon phase and Moon distance to enable corrections for these factors when predicting Moon contamination for individual L1b pointings.

Zoom In Zoom Out Reset image size

As we stack, at each pixel in the co-add, we keep an integer count of the number of Moon-affected frames contributing to that pixel. We define Moon-affected frames to be those for which our model predicts moon contamination of S_m > 100 DN, and $S_m/\mathcal {R}_{100}>450$. Pixels with Moon-contamination count >3 at the end of the stacking process are flagged with bit 15, as listed in Table 2.

In rare cases typically associated with moon contamination, all L1b quadrants at a particular sky location are rejected during the reference comparison step. Mask bit 19 labels such locations, at which we have tapered in the reference stack to avoid leaving regions of missing data in the final mosaic. 0.15% of the sky is affected.

6.5.2. Resolved Compact Sources

In regions with low integer frame coverage <5, compact sources are listed only in the WISE All-sky Reject Table, and are excluded from the All-Sky Source Catalog (Cutri et al. 2012). Rather than attempt to identify valid sources among Reject Table entries, we opted instead to flag such regions with mask bit 21. 1.2% of the sky is affected.

6.6.1. Big Objects

6.6.2. Low Ecliptic Latitude

Our mosaicking process only attempts to create a smooth image over each tile footprint by best matching each L1b exposure with its neighbors. However, this process does not ensure that neighboring tiles agree on their overlaps. Further, relics of the zodiacal light and Moon contamination that are not perfectly removed by our SFD-based gradient correction can persist, leading to wandering of the tile zero points on relatively large scales of several degrees.

We investigated the possibility of obtaining a zero point by taking the L1b zero point at face value, and then subtracting out models for the time-dependent foregrounds. Despite our best efforts to subtract the Moon-glow model of Section 6.5.1 and various customizations of the Kelsall zodiacal light model (Kelsall et al. 1998), such an approach ultimately yielded a cosmetically unappealing result. Based on the Planck comparison we describe below, we estimate that, at high |b|, the true W3 Galactic emission is ∼0.5% as bright as the |β| = 90° zodiacal light. While the Kelsall model appeared reasonable at the several percent level, this remained inadequate due to the overwhelming brightness of the zodiacal light at 12 μm.

We therefore resorted a zero point based on Planck (Planck Collaboration et al. 2013a). Henceforth, the term Planck will refer to Planck 857 GHz (350 μm), corrected for zodiacal emission according to Planck Collaboration et al. (2013c). To ensure this Planck map contained only diffuse Galactic emission, we also interpolated over the SFD compact source mask. The great advantage of Planck relative to other full-sky, mid-IR and far-IR datasets is the very faint zodiacal foreground. The basic strategy will be to replace modes on scales of several degrees in our WISE tiles with appropriately rescaled and smoothed Planck.

In order to accomplish this goal, we must derive a map of the spatially varying coefficient f_pw that converts 350 μm emission to 12 μm emission. To derive f_pw we correlate cleaned WISE with Planck on a bandpass between 4' and 15'. We break each tile into many 05 × 05 rectangular patches and perform a linear least squares fit to the bandpass filtered WISE sub-image, where the model is bandpass Planck multiplied by f_pw, plus a constant offset. For each such fit, we also compute the Pearson correlation coefficient c of bandpass WISE versus bandpass Planck. We perform such fits over the full sky. For a variety of reasons, including poor signal to noise in Planck or WISE, or systematic issues in Planck or WISE, c is often closer to zero than to unity. In such cases, the fit is driven to f_pw → 0 for reasons that are instrumental rather than astrophysical. In order to obtain a map of f_pw on large scales, we make a cut on c > 0.7 and smooth the map of f_pw at those locations to 15°. With WISE in W3 DN and Planck in MJy sr⁻¹ f_pw, varies between 2.2 and 5.8. Converting WISE to MJy sr⁻¹ according to Cutri et al. (2012) so that both data sets have the same units, this translates to 0.036 <f_pw < 0.095.

To replace the large scale modes of WISE with those of Planck, we apply a 125 median filter to both data sets, and then smooth both to 2°. We inpainted over regions with strong molecular emission, defined as Planck Type 3 CO >15 K km s⁻¹ (Planck Collaboration et al. 2013b). We subtract 2° WISE from each tile, and then add in smoothed Planck, multiplied by f_pw. In order to derive tiles free of edge effects from the smoothing process, we initially built our WSSA stacks on 14° × 14° padded astrometric footprints.

Lastly, we convert the WISE tiles with Planck-based zero point from W3 DN to MJy sr⁻¹ according to the prescription of Cutri et al. (2012, Section 4.4h).

We release the 430 W3 WSSA tiles as multi-extension FITS files. Table 3 lists the extensions included. Figure 10 shows three of these file extensions for an example region with prominent diffuse emission. Refer to the tile headers for further details.

Zoom In Zoom Out Reset image size

The "dynamic calibration sky offset correction" (Cutri et al. 2012, Section 4.4a) introduced several problematic features into the L1b frames which we did not address for this initial 12 μm tile release. In particular, long-term latents due to bright sources are mirrored in acausal fashion to affect frames temporally preceding bright source detection. Examples of such features can be seen in Figure 8 as low-level, dark "divots" ∼1'–2' in size. These acausal latents are the dominant remaining imperfection in many regions of our co-adds with low dust signal. In a related problem, the dynamic sky offset procedure can often assign up to ∼4% of pixels in a given L1b image to exactly the median pixel value within the frame. Because these issues would be better addressed with lower-level data products than the L1b images, we have so far opted not to address them.

The relevant DIRBE data product to compare with our W3 map is the band 5 "Zodi-Subtracted Mission Average" (ZSMA5). Because DIRBE pointed at a variety of solar elongations, including some <90°, the DIRBE-observed zodiacal light is even more pronounced and complex than that in W3. Arguably, the most prominent feature in ZSMA5 is the residual zodiacal light left over after subtracting the DSZA model. The ecliptic plane is obvious in such residuals, and at high |b| there remains a clear gradient of average flux with β, with magnitude of ∼6 × 10⁻³ (MJy sr⁻¹) deg⁻¹, or ∼0.2 MJy sr⁻¹ between |β| ∼ 60° and |β| ∼ 90°. By comparison, our Planck-based far-IR zero point suggests that at |b| ≳ 80°, the entire non-IPD 12 μm signal is ∼0.03 MJy sr⁻¹. There is also a significant offset between ZSMA5 and our map in regions near the ecliptic poles, with amplitude ∼0.6 MJy sr⁻¹. While it might be interesting to investigate these differences, they are immaterial to our long-term goal of predicting dust extinction. Thanks to our far-IR zero point, our W3 map only shows traces of zodiacal light in rare circumstances, very near the ecliptic plane, where the zodi can vary strongly on scales smaller than 2°. Our W3 map is also superior to ZSMA5 in that point sources remain unsubtracted in the latter, impacting a substantial fraction of the sky. On the other hand, DIRBE band 5 is truly a full-sky map, whereas ∼5% of diffuse W3 emission is badly contaminated by Moon-glow and/or solar system planets.

IRAS 12 μm has a separate "reject" data set at low |β| that is not consistent with the remainder of their map. Thus, our W3 map improves upon IRAS 12 μm by virtue of (1) its consistent, zodi-free zero point (2) its factor of 16 enhancement in angular resolution and (3) our removal and masking of compact sources.

We thank Roc Cutri for his guidance in making use of the WISE data products. We also thank Bruce Draine for his careful reading of an early version of this manuscript.

We acknowledge support of NASA grant NNX12AE08G for this research. This research was conducted with Government support under and awarded by DoD, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a. This research made use of the NASA Astrophysics Data System (ADS) and the IDL Astronomy User's Library at Goddard.⁵

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication also makes use of data products from NEOWISE, which is a project of the Jet Propulsion Laboratory/California Institute of Technology, funded by the Planetary Science Division of the National Aeronautics and Space Administration. This publication makes use of data products from the Two Micron All Sky Survey, 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.