Deep Full-sky Coadds from Three Years of WISE and NEOWISE Observations

TL;DR

This paper presents the creation of the deepest full-sky maps at 3-5 microns from three years of WISE/NEOWISE data, enabling improved deep photometry for astronomical surveys and target selection.

Contribution

It introduces a novel reprocessing pipeline that combines all available WISE/NEOWISE exposures to produce high-resolution, deep full-sky coadds with artifact removal for enhanced astronomical analysis.

Findings

Coadds are ~3 times deeper than previous AllWISE stacks.

Forced photometry depths are increased by 0.56 and 0.46 magnitudes in W1 and W2.

The new maps facilitate uniform target selection and searches for rare objects.

Abstract

We have reprocessed over 100 terabytes of single-exposure WISE/NEOWISE images to create the deepest ever full-sky maps at 3-5 microns. We incorporate all publicly available W1 and W2 imaging - a total of ~8 million exposures in each band - from ~37 months of observations spanning 2010 January to 2015 December. Our coadds preserve the native WISE resolution and feature depth of coverage ~3 times greater than that of the AllWISE Atlas stacks. Our coadds are designed to enable deep forced photometry, in particular for the Dark Energy Camera Legacy Survey (DECaLS) and Mayall z-Band Legacy Survey (MzLS), both of which are being used to select targets for the Dark Energy Spectroscopic Instrument (DESI). We describe newly introduced processing steps aimed at leveraging added redundancy to remove artifacts, with the intent of facilitating uniform target selection and searches for rare/exotic objects (e.g. high-redshift quasars and distant galaxy clusters). Forced photometry depths achieved with these coadds extend 0.56 (0.46) magnitudes deeper in W1 (W2) than is possible with only pre-hibernation WISE imaging.

Tables

Table 1. Table 1: Description of bright star mask bits.

Bit	Description	Scan Direction
0	W1 bright star	south
1	W1 bright star	north
2	W2 bright star	south
3	W2 bright star	north

Full text

††pubyear: 2016††pagerange: Deep Full-sky Coadds from Three Years of WISE and NEOWISE Observations–References

1 Introduction

The Wide-field Infrared Survey Explorer (WISE; Wright et al., 2010) was designed to map the entire sky at mid-infrared wavelengths with sensitivity far exceeding that of its predecessors, IRAS (Wheelock et al., 1994) and DIRBE (Boggess et al., 1992). As a source of high-quality, full-sky imaging, its observations have an extremely wide range of applications, from near-Earth asteroids (e.g. Connors et al., 2011) to the most luminous galaxies in the universe (e.g. Tsai et al., 2015).

WISE is a 0.4 meter telescope onboard a satellite in low-Earth orbit, launched in late 2009. Always pointing near 90*∘* solar elongation while making successive scans in ecliptic latitude at fixed ecliptic longitude, WISE maps the entire sky once every six months when operational. From 2010 January to 2010 August, WISE carried out a full-sky mapping in each of four broad mid-infrared bandpasses centered at 3.4$\mu$m (W1), 4.6$\mu$m (W2), 12$\mu$m (W3) and 22$\mu$m (W4).

Due to the depletion of solid hydrogen cryogen, the W3 and W4 channels were rendered unusable as of 2010 September and 2010 August, respectively. However, WISE continued surveying in W1 and W2 through 2011 January as part of the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE; Mainzer et al., 2011) mission. In 2011 February WISE was placed in hibernation, temporarily ceasing data acquisition. In December 2013, WISE recommenced surveying the sky in W1 and W2, carrying out the NEOWISE-Reactivation (NEOWISER; Mainzer et al., 2014) mission. Despite the multi-year hiatus, NEOWISER images are of essentially the same high quality and sensitivity as exposures acquired pre-hibernation (Mainzer et al., 2014). In March 2015, NEOWISER published its first year’s worth of single-exposure images and frame-level source extractions. In March 2016, a second year of such data were made public.

Because NEOWISER is an asteroid hunting and characterization mission, it does not deliver any coadded data products. Nevertheless, the astrophysics research community has recognized the tremendous value of coadded data products incorporating NEOWISER images (e.g. Faherty et al., 2015).

The Dark Energy Spectroscopic Instrument (DESI, Levi et al., 2013; DESI Collaboration et al., 2016a, b) represents an important application which stands to benefit from access to deep WISE stacks, and which ultimately drives many of our design considerations in building coadded WISE/NEOWISE data products. DESI will select millions of luminous red galaxy (LRG) targets from $r$+$z$+W1 photometry, while its quasar targeting will additionally make use of $g$ and W2 fluxes. To obtain high-quality photometrically selected samples of these targets, DESI requires the deepest and cleanest possible W1/W2 coadds.

DESI targeting employs a ‘forced photometry’ approach which measures WISE fluxes for every optically detected source, fixing the centroid and morphology to those obtained from much higher resolution optical imaging. This technique has already been successfully applied to optical catalogs from SDSS (Lang et al., 2016) and DECam (Schlegel et al., 2015; Meisner et al., 2017a). The Lang et al. (2016) results have played a significant role in eBOSS LRG and quasar selection (Prakash et al., 2016; Prakash et al., 2015; Myers et al., 2015).

In the eBOSS/DESI targeting applications, WISE forced photometry has been performed on ‘unWISE’ coadds (Lang, 2014; Meisner et al., 2017a). The original unWISE coadds of Lang (2014) used the available WISE data from observations conducted in 2010 and early 2011. The W1/W2 coadds of Meisner et al. (2017a) are based on an adaptation of the Lang (2014) unWISE coaddition pipeline, folding in the first year of NEOWISER imaging. Meisner et al. (2017a) thereby doubled the depth of coverage while eliminating the dominant artifacts found in the Lang (2014) coadds.

Here we update the results of Meisner et al. (2017a), adding in the most recently published year of NEOWISER W1/W2 exposures. We highlight our recent processing improvements which aid in the elimination and flagging of remaining artifacts. Our latest coadds represent a 50% increase in depth of coverage relative to those of Meisner et al. (2017a), and a 200% increase relative to those of Lang (2014) in W1 and W2.

In $\S$2 we list the input data used for this work. In $\S$3 we provide an overview of unWISE coaddition and our image processing strategy/philosophy. In $\S$4 we highlight newly added improvements to our coaddition pipeline. In $\S$5 we display and validate key aspects of our new W1/W2 full-sky maps. We conclude in $\S$6.

2 Data

Our coaddition proceeds from the least-processed form of publicly available WISE imaging, namely the “Level 1b" (L1b) single-exposure framesets. We downloaded a local copy of every publicly available W1 and W2 frameset, including those from the All-Sky, 3band Cryo, NEOWISE, NEOWISER year one (NEO1) and NEOWISER year two (NEO2) releases (Cutri et al., 2012, 2013, 2015). For each frameset, a -int- FITS file gives the sky intensity, while a -unc- FITS file provides per-pixel uncertainty estimates, and a -msk- FITS file contains a bitmask flagging artifacts such as bad pixels and cosmic rays. In all, we downloaded $\sim$52 TB of L1b data products per band, totaling 49 terapixels of inputs.

To validate photometry derived from our coadds, we make use of the AllWISE Source Catalog (Cutri et al., 2013), as well as W1/W2 forced photometry from Lang et al. (2016) and Data Release 4 (DR4) of the DESI imaging Legacy Survey111http://legacysurvey.org. DR4 is currently internal to the DESI collaboration, but is expected to become public in the near future..

3 Processing Overview

3.1 Atlas versus unWISE Coadds

There are two main types of full-sky coadded WISE data products presently available. First, there are the Atlas stacks created by the WISE/NEOWISE team. These images are intentionally blurred by the WISE point spread function (PSF), since their primary purpose is to enable deep source extraction for e.g. the AllWISE Source Catalog222Once detections are in hand, however, the AllWISE catalog measurements for each source are derived via joint modeling of all relevant single-exposure images.. The most recent full-sky set of Atlas coadds is that of the AllWISE release and is based solely on pre-hibernation exposures.

Alternatively, the unWISE (Lang, 2014) line of coadds uses Lanczos interpolation to preserve the native WISE resolution. The unWISE stacks are designed for forced photometry, a use case for which any blurring of the native WISE PSF would be suboptimal. This work is part of an ongoing effort to upgrade the W1/W2 unWISE coadds with each new release of NEOWISER L1b exposures, enhancing the depths achieved and using added redundancy to remove artifacts.

3.2 unWISE Philosophy

unWISE coaddition is meant to be lightweight/minimalist, with the plurality of compute time consumed by reading in the large number of L1b frames which overlap a given coadd footprint, and subsequently writing the output coadd images. As such, no truly robust outlier rejection methods (e.g. median filtering the resampled L1b intensities at each coadd pixel location) are employed, as these would dramatically increase the total computational cost.

Given that we have already produced full-depth coadds including all exposures through the first-year NEOWISER release (Meisner et al., 2017a), one might imagine a scheme in which we update our coadds by only processing the “new” exposures, adding them onto the existing coadd. Such an approach would not fully take advantage of improved outlier rejection enabled by the two newly added sky passes from second-year NEOWISER observations. For instance, we would inherit artifacts present in previous iterations of the unWISE coadds, rather than leveraging the newly added redundancy of the latest set of exposures to improve artifact removal (see e.g. Figure 5). Such an approach would also likely require us to implement a series of contrived shortcuts/hacks.

Therefore, with every new release of NEOWISER exposures, we opt to rebuild all of our full-depth coadds from scratch. Even when jointly processing all exposures spanning 2010 January to 2015 December, as we have chosen to do, the total full-sky computational expense is only of order tens of thousands of CPU hours.

In combining observations which span the full WISE lifetime, we also strive to achieve the best possible relative calibration of all exposures. To accomplish this, we adopt a custom photometric calibration using repeat measurements of calibrator sources near the ecliptic poles, as described in $\S$4 of Meisner et al. (2017a). We have now performed this photometric calibration analysis for the second year of NEOWISER observations. The results are shown in Figure 1.

One final component of our philosophy in generating WISE coadds is to throw away as little data as possible. This motivated our work in Meisner et al. (2017a) to recover exposures corrupted by scattered moonlight. In this same vein, we retain all frames at/near the ecliptic poles, despite the added computational cost associated with exceptionally large depth of coverage in these regions. As a result, the tiles at the north/south ecliptic poles have peak integer coverage of $\sim$17,000 frames per band, for a total integration time of $\sim$1.5 days.

3.3 unWISE Implementation Details

unWISE coadds employ the same set of 18,240 tile centers as do the AllWISE Atlas stacks. We adopt a pixel scale matching that of the W1/W2 L1b exposures themselves, 2.75*′′* per pixel. Each coadd astrometric footprint has its $-x$ ($+y$) axis aligned with celestial east (north). We refer to a coadd astrometric footprint as a ‘tile’, and identify tiles by their coadd_id values, which are strings encoding their central (RA, Dec) coordinates, e.g. ‘0000p000’. Our coadds are 1.56*∘* (2048 pixels) on a side.

unWISE coaddition proceeds by identifying all exposures that overlap the tile of interest, then resampling every relevant exposure onto the coadd pixel grid using Lanczos interpolation. With all resampled single-exposure pixels in memory, a “first round” of coaddition computes a simple mean and standard deviation of the contributing L1b intensity values at each coadd pixel location. These first-round coadd mean and standard deviation images are then used to enable outlier rejection during a second and final round of coaddition. For full details refer to Lang (2014) and Meisner et al. (2017a).

4 Coaddition Updates/Improvements

The foremost improvement achieved in this work relative to the coadds of Meisner et al. (2017a) is the $\sim$50% enhancement in depth of coverage obtained by including second-year NEOWISER observations. However, because we uniformly reprocess all exposures following every NEOWISER release, each new, deeper set of coadds also represents an opportunity to build further processing improvements into the unWISE coaddition pipeline. The following subsections describe newly introduced facets of the unWISE processing and data products.

4.1 Min/Max Rejection

As mentioned previously, unWISE coaddition does not use robust statistics to perform outlier rejection, instead relying on naive first-round coadd mean and standard deviation images to flag problematic pixels in the L1b exposures. This computationally cheap approach is vulnerable to the appearance of extreme outliers in a small number of exposures. The most conspicuous example of this phenomenon arises due to large numbers of cosmic ray strikes in frames affected by the South Atlantic Anomaly, which mainly impacts observations at $-45^{\circ}<\delta<-10^{\circ}$. This leads to first-round coadd results which are strongly contaminated by cosmic rays, and therefore are relatively ineffective at flagging outliers (see left panel of Figure 2).

To combat such corruption of our first-round coadds, we have now implemented a “min/max” rejection step during first-round coaddition. Specifically, for each pixel in coadd space, we precompute which L1b exposures contribute the highest and lowest single-frame intensity values at that location. Then, when computing the first round coadd mean and standard deviation, we disregard the minimum and maximum single-exposure values at each coadd pixel location. Precomputing the minimum and maximum values at each coadd pixel location requires $\sim$0.01 seconds per exposure, which is negligible relative to the overall runtime of coaddition, and typically leads to $\sim$10 additional seconds of CPU time per tile at low ecliptic latitude.

The right panel of Figure 2 shows an example of dramatically reduced cosmic ray contamination in a first round coadd mean image achieved by incorporating our min/max rejection step. Obtaining a cleaner first-round estimate of the per-pixel mean and standard deviation improves our ability to reject outliers during the second round of coaddition. This is illustrated by the example region shown in Figure 3, where min/max rejection has eliminated the imprint of cosmic ray strikes from our final coadd image.

Note that our min/max rejection step does not necessarily result in a reduction of integer coverage for our final outputs at every coadd pixel, since it is only applied during the first round coaddition process. In future iterations of our full-depth coadds, as the typical integer coverage continues to increase, we may investigate more aggressive min/max rejection in which the highest/lowest $N$ single-exposure values at each coadd pixel location are discarded, with $N>1$.

4.2 Bitmasks

One feature which has been absent from previous unWISE data releases is a pixel-level bitmask corresponding to each coadd image. Such a bitmask data product could be very helpful for selection of spectroscopic samples (e.g. for DESI), in order to improve efficiency by flagging sources associated with artifacts. It could also be quite useful in aiding rare object searches, by indicating when unusual objects might in fact be of instrumental rather than astrophysical origin. For instance, a close (W1-only) variant of the bitmask procedure presented here was used to great effect in the Planet Nine search of Meisner et al. (2017b).

In this work, we address these needs by creating a pixel-level bright star mask output for each coadd tile, with suffix -msk.fits.gz. We construct one such bitmask image per tile, containing information about bright stars in both W1 and W2, rather than generating two separate bitmasks with one corresponding to each band.

The simplest bright star mask one might imagine would flag a circular region about every sufficiently bright star, with the radius scaled based on each star’s total flux according to some empirically calibrated prescription. Indeed, the eBOSS collaboration currently employs masks of bright WISE stars created in this way (Prakash et al., 2016; Myers et al., 2015). Our bitmasks are designed to go beyond this simplistic approach, taking into account both the WISE PSF and the WISE scan direction to create highly detailed masks which include features such as diffraction spikes and optical ghosts.

As the basis for our bitmasks, we must select samples of very bright WISE sources in each of W1 and W2. We do so by making use of the positions and fluxes in the AllWISE catalog. We create a full-sky list of W1-bright sources by selecting objects from the AllWISE catalog with w1mpro $<$ 9.5. For W2, we adopt a brightness threshold of w2mpro $<$ 8.3. Our masks are primarily designed with extragalactic/cosmology applications in mind, and therefore are optimized for high Galactic latitude sky regions. In locations far enough off the Galactic plane to be part of DESI’s footprint ($|b_{gal}|>18^{\circ}$), our bright star sample contains on average $\sim$29 ($\sim$11) sources per square degree in W1 (W2).

For each coadd_id astrometric footprint, we begin constructing its bitmask by identifying all bright stars which may have profiles at least partially falling inside of the tile boundaries. We properly account for all bright stars which overlap the tile footprint of interest with any portion of their PSF wings/diffraction spikes, even those for which the centroid falls outside of the tile boundaries. For each bright star, we then render a model of its appearance given its AllWISE centroid and flux. To render these models, we make use of the W1 and W2 PSFs of Meisner & Finkbeiner (2014). These PSFs, displayed in Figure 4 of Lang (2014), are $14.9^{\prime}$ on a side, and therefore extend quite far into the wings, capturing details such as diffraction spikes and the W2 optical ghost.

One subtlety involved in rendering our bright star profiles is that the Meisner & Finkbeiner (2014) PSFs provide models of the single-exposure PSF, not the effective PSF obtained after resampling onto the unWISE tile astrometric footprints and coadding. For instance, in L1b images, the diffraction spikes always emanate from bright sources at (45*∘, 135∘, 225∘, 315∘*) degrees from the +$x_{L1b}$ direction. However, this is not the case for unWISE tiles, which are oriented along the equatorial cardinal directions, whereas the L1b exposures are very nearly oriented exactly along the ecliptic cardinal directions.

A further subtlety is that the Meisner & Finkbeiner (2014) PSF models require specification of detector ($x_{L1b}$, $y_{L1b}$) coordinates, since they incorporate PSF variation across the single-exposure field of view (FOV). Therefore, before rendering any bright star models during bitmask construction, we precompute mean versions of the W1 and W2 Meisner & Finkbeiner (2014) PSFs, averaging over ($x_{L1b}$, $y_{L1b}$) detector location.

A final subtlety is that, except for a small area near $|\beta|=90^{\circ}$, there will be two discrete WISE scan directions contributing to each coadd, corresponding to the ecliptic north and south directions. The PSF is not perfectly symmetric with respect to swapping the scan direction. The most notable asymmetry in the bands of interest is the W2 ghost, which is offset by $\sim$5*′* from its parent bright source centroid, at a position angle that is fixed relative to the scan direction. In practice this means that, in our coadds, the ghost will appear on opposite sides of its parent bright star in exposures with opposite scan directions. This explains why the doughnut-shaped ghost in Figure 4 appears twice.

In order to create a coadd-level model rendering of a bright star, we compute the angle between the scan direction and celestial north, rotate the FOV-averaged PSF by this angle, and scale the rotated profile so that its total flux matches that quoted by the AllWISE catalog. For each bright star, we create two such renderings, one for each scan direction. Pixels in the bright star models above a threshold of 13.2 (38.3) Vega nanomaggies per square arcsecond in W1 (W2) are flagged according to the mask bits listed in Table 1. At very high latitude, $|b_{gal}|>80^{\circ}$, $\sim$0.4% of the total sky area is masked on average. This fraction ramps up toward lower $|b_{gal}|$, reaching 1.6% at $|b_{gal}|$ = 18*∘*, the closest DESI will observe to the Galactic plane.

Our pixel-level bitmasks have now been propagated into the DESI imaging Legacy Survey DR4 catalogs. An example bitmask image near an extremely bright star is shown in Figure 4. Much room remains to extend our present bitmasks in future unWISE coadd releases, making them even more intricately detailed and elaborate. The primary avenues for doing so are:

• •

Near the Galactic plane, our current bitmasks flag an excessively large fraction of pixels, rendering the masks unhelpful in these regions. In the future, we could counteract this issue by making our w?mpro thresholds dependent on e.g. Galactic latitude, so that a larger flux would be required to trigger masking of a star at low $b_{gal}$.

• •

As seen in Figure 4, the diffraction spikes of extremely bright stars can sometimes extend beyond the $14.9^{\prime}$ boundary of the Meisner & Finkbeiner (2014) PSF models. In the future, we could calibrate look-up tables of diffraction spike length versus magnitude in W1 and W2, and use these to create geometric masks composed of lines emanating from the centroids of extremely bright stars at angles of (45*∘, 135∘, 225∘, 315∘*) relative to ecliptic north. Such a procedure was employed in W3 by Meisner & Finkbeiner (2014), as described in their $\S$5.2.4.

• •

“Latents” represent another class of defects associated with stars bright enough to saturate in their cores. Latents are persistence artifacts which appear as diffuse positive blobs at the detector positions of saturated pixels, but in the frame immediately after imaging of the parent bright star. The locations of all latents can be predicted exactly given L1b metadata and a catalog of bright stars, so it would be possible to reliably flag these artifacts by adding new mask bits.

• •

Our assumption that there are exactly two discrete scan directions corresponding to ecliptic north and south is violated at the ecliptic poles, where continuous coverage yields a continuum of scan directions. As a result, very near $|\beta|$=90*∘*, bright star diffraction spikes become spread out into discs in the coadds and are correspondingly attenuated. One approach for taking these effects into account has been described in $\S$5.2.4 of Meisner & Finkbeiner (2014).

5 Results

5.1 Images

Relative to the coadds of Meisner et al. (2017a), we have in general increased the depth of coverage by $\sim$50% while further mitigating artifacts. The mean integer coverage over the entire sky is 108 (107) frames in W1 (W2), and every pixel has integer coverage of at least 33 (30) in W1 (W2). Figure 5 illustrates that folding in a third full year of observations has reduced the impact of scattered moonlight, the dominant systematic problem with the W1/W2 unWISE coadds on large angular scales. Because the Moon affects different ranges of ecliptic longitude during different sky passes, any remaining traces of scattered moonlight contamination become further suppressed as additional NEOWISER coverage is incorporated. Figure 6 shows a typical example of the reduction in statistical noise achieved for coadds at low ecliptic latitude.

5.2 Catalogs

DR4 of the DESI imaging Legacy Survey, hereafter referred to as simply ‘DR4’, contains the only currently available photometry based on the full-depth W1/W2 coadds presented in this work. This photometry is forced rather than WISE-selected, adopting source locations and morphologies derived from deep optical data. DR4 covers a $\sim$4,500 square degree extragalactic footprint limited to $\delta\gtrsim+30^{\circ}$, as shown in Figure 3.19 of DESI Collaboration et al. (2016a).

In validating our coadds, we first seek catalog-level confirmation of the decrease in statistical noise suggested by the images in Figure 6. We examine objects with W1/W2 forced photometry available from both Lang et al. (2016), based on the original unWISE coadds, and DR4 of the Legacy Survey. We select a comparison sample drawn from the region $170^{\circ}<\alpha<230^{\circ}$, $42.5^{\circ}<\delta<48^{\circ}$, using a 1*′′* match radius, and restricting to objects with DR4 forced photometry signal-to-noise in the range 10 $\pm$ 1. For these 127,000 (144,000) sources in W1 (W2), we find that forced photometry flux uncertainties from Legacy Survey DR4 are smaller than those of Lang et al. (2016) by median factors of 1.68$\times$ (1.53$\times$). These values indicate that forced photometry based on the present coadds is 0.56 (0.46) mags deeper in W1 (W2) than forced photometry based solely on pre-hibernation W1/W2 imaging. For both the W1 and W2 samples, the median increase in integer coverage is a factor of 2.94$\times$, leading us to expect flux uncertainties reduced by 1.71$\times$ under the assumption that all WISE observations have maintained the same sensitivity regardless of mission phase.

We additionally check that forced photometry fluxes derived from our new coadds are consistent with the AllWISE catalog. We select a comparison sample drawn from the region $170^{\circ}<\alpha<230^{\circ}$, $42.5^{\circ}<\delta<48^{\circ}$, using a 1*′′* match radius. We further require a one-to-one DR4-AllWISE match, DR4 morphological type PSF and w?cc_map = 0 in the band of interest. This yields 626,000 (612,000) sources in W1 (W2). Figure 7 summarizes our DR4-AllWISE comparison for this sample. There is good overall agreement across a wide range of fluxes, and the (DR4 $-$ AllWISE) offset asymptotes to 4.1 ($-$1.6) mmag in W1 (W2) toward the bright end. The faint end upturn has previously been noted in Lang et al. (2016) and Meisner et al. (2017a), and is due to Malmquist bias, as the AllWISE sources are WISE-selected while the DR4 objects are not.

5.3 Data Access

The full-depth W1/W2 coadds described in this work are publicly available via the unWISE web interface at http://unwise.me.

6 Conclusion

We have reprocessed all publicly available W1/W2 observations ever acquired to create the deepest full-sky maps at 3$-$5$\mu$m. New processing steps have been introduced, focused on enabling clean target selection and rare object searches. We have also validated our new coadds using forced photometry from DR4 of the DESI imaging Legacy Survey. It will be important to continue updating our full-depth W1/W2 stacks with each future release of additional NEOWISER exposures, as doing so represents a crucial step toward realizing the full potential of the entire WISE imaging data set.

Acknowledgements

This work has been supported by grant NNH17AE75I from the NASA Astrophysics Data Analysis Program.

This research 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 research 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 research has made use of the NASA/ IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

The National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, provided staff, computational resources, and data storage for this project.