NEOSURVEY 1: INITIAL RESULTS FROM THE WARM SPITZER EXPLORATION SCIENCE SURVEY OF NEAR-EARTH OBJECT PROPERTIES

We selected our sample targets based on Spitzer Cycle 11 observability from the 2014 October 24 version of the Minor Planet Center NEO list¹¹ . Observability was checked using the NASA Horizons system¹² for each day in Cycle 11 (2015 February through 2016 September); Horizons was accessed with the Callhorizons ¹³ Python module, which retrieves ephemerides and orbital elements from the Horizons system based on user queries. An object was considered a potential target when its solar elongation with respect to Spitzer complied with the observatory constraints (825–120°), its positional uncertainty was less than 150'' at 3σ (so that the target would fall within the detector FOV), it had a proper motion of less than 1'' s⁻¹ (though in practice no objects are excluded by this constraint), and its distance from the Galactic plane was larger than 15°. These requirements are based on experience from our successful ExploreNEOs program (Trilling et al. 2010). From this list of potential observations and targets, we removed all targets that had already been observed by Spitzer (ExploreNEOs: Trilling et al. 2016), WISE (Mainzer et al. 2011a, 2012, 2014a), or Akari (Usui et al. 2011; Hasegawa et al. 2013).

Spitzer's Infrared Array Camera (IRAC; Fazio et al. 2004) is the only instrument available in Spitzer's Warm Mission, providing imaging at 3.6 and 4.5 μm (CH1 and CH2, respectively). The field of view is 5.2'. Fluxes from NEOs at CH2 are always dominated by thermal emission (Figure 1), but the ratio of reflected light to thermal emission in CH1 cannot be determined if the albedo at 3.6 μm is not known (Mainzer et al. 2011a). In our previous ExploreNEOs program (Trilling et al. 2016) we found that CH1 observations are generally not useful for thermal modeling. We therefore do not observe in CH1 in this program, saving more than 50% on our total time request¹⁴ .

Zoom In Zoom Out Reset image size

For each potential target we predicted the thermal-infrared flux density in IRAC CH2 as a function of time. Our flux predictions were based on the absolute optical magnitude H, a measure of D²p_V, as reported by Horizons (D is diameter, and p_V is albedo). H values for NEOs are of notoriously low quality and tend to be skewed to brighter magnitudes (Jurić et al. 2002; Romanishin & Tegler 2005). To account for this, we assumed an H offset (ΔH) of (0.6, 0.3, 0.0) mag for (low, nominal, high) fluxes, respectively. Optical fluxes were calculated from $H+{\rm{\Delta }}H$ together with the observing geometry. In our observation planning we assumed that asteroids are 1.4 times more reflective at 4.5 μm than in the V band (Trilling et al. 2008; Harris et al. 2011; Mainzer et al. 2011a), although in our analysis pipeline we use an updated value that is derived from NEOWISE data, as described below. Thermal fluxes also depend on p_V (D is determined from H and p_V) and η (see Section 3). We assumed p_V = (0.4, 0.2, 0.05) for (low, nominal, high) thermal fluxes. The nominal η value was determined from the solar phase angle α using the linear relation given by Wolters et al. (2008), which is generally in agreement with the newer results of Mainzer et al. (2011a) and the approach used in this work (see Section 3.3, below); 0.3 was added (subtracted) for low (high) fluxes to capture the scatter in the empirical relationship derived in Wolters et al. (2008). From these parameters, we calculated the predicted thermal fluxes using the Near-Earth Asteroid Thermal Model (Harris 1998; discussed further in Section 3). We calculated (low, nominal, high) predicted fluxes in the 4.5 μm IRAC band for each day of Cycle 11. The resulting NEATM fluxes were convolved with the CH2 passband to yield "color-corrected" in-band fluxes.

After removing all dates where an NEO's high predicted flux could saturate the detector (using saturation values from the IRAC Warm Mission Handbook), we identified a 5-day window centered on the peak brightness during Cycle 11. We used the minimum flux calculated for the low predicted values in this 5-day window to calculate our integration times and build our Astronomical Observation Requests (AORs). This 5-day window dictates the timing constraint that we place on the AOR. A window of this size allows good scheduling flexibility while allowing us to use the shortest possible integration times. Using the low prediction values here ensures that our signal-to-noise ratio (S/N) requirement will be met or exceeded. Our experience has shown that the ratio of observed to low predicted flux is in the range of 1–10, so our conservative prediction technique is appropriate.

Our photometric requirement is S/N ≥ 15, at which point, our detailed thermal modeling shows, the flux density uncertainty becomes small compared to overall model uncertainties. We have empirically measured the Warm Spitzer CH2 sensitivity (5σ) for faint moving objects to be 2.2 μJy in 5 hr and 1.8 μJy in 10 hr, scaling roughly as the square root of integration time for shorter integrations (Mommert et al. 2014b). This is less sensitive than predicted by SENS-PET due to contamination from background sources that cannot be completely subtracted out in the moving frame.

Our threshold is set to integration times of ≤10,000 s (2.7 hr), which produces an AOR clock time including overheads of 3.2 hr. All sample targets have a 4.5 μm thermal-infrared flux density of 7.2 μJy or greater. In 710.1 hr of total observation time, including all overheads, we will observe 597 NEOs. Our final sample is shown in Figures 2 and 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We use a standard moving single target AOR, in which the telescope tracks at the appropriate moving object rates for the object, and a medium cycling dither, with only the 4.5 μm field integrating on the NEO (as we did recently for 2011 MD; Mommert et al. 2014a). As described below, the 4.5 μm band provides the most information on the thermal emission from the NEO, whereas the 3.6 μm band can have a significant (and in general unknown) contribution from reflected light and does not significantly improve the derivation of an albedo and size. Not observing in CH1 yields significant time savings as well; since most NEOs are fainter in CH1 than in CH2 (Figure 1), observing in CH1 would more than double the needed time. Each target was assigned to one of six uniform observing templates with total frame times of 120, 480, 1200, 2400, 6000, and 10,000 s. The individual frame times are 12 s and 30 s for the first two bins, respectively, and 100 s for the remainder; this allows many dithers for each observation.

The frame times used for each object were chosen based on our estimate of the source flux in order to obtain the required sensitivity using the minimum total observation time. However, we have seen that occasionally the NEO is brighter than expected, possibly due to errors in the previously determined H magnitude, or, as in the case of Don Quixote, the result of cometary activity (Mommert et al. 2014c). We therefore use HDR mode for the 12 and 30 s frame times, which minimizes the chance of the object saturating. However, if saturation does occur, we can use a PSF-fitting technique to recover the NEO flux, as we did for Don Quixote¹⁵ . (Using HDR mode for the brightest targets increases the total overhead only a small amount.) For the 100 s frame times, where the objects are very likely much fainter than the saturation limit, we do not use HDR mode.

We have also checked for objects that move only slightly during the minimum observation time required to reach S/N ≥ 15. If an object moves only a small distance during the AOR, it is difficult to subtract the background field and perform accurate photometry, especially if it falls in close proximity to a bright source. We have therefore found all objects that would move less than ∼3 FWHMs during the observation, and we added additional integrations for these sources so that the object moves sufficiently to be able to measure and subtract the background field. This has been done for approximately 50 sources, which are all relatively bright targets and have short-duration AORs.

The objects in our target pool are distributed all over the sky and have an essentially uniform distribution of ecliptic longitudes (except at the Galactic plane crossings). No AOR is longer than 3.2 hr, including all overheads, and more than half the targets have AORs that are 30 minutes or less. As such, this program maximizes both total NEO yield and efficiency.

The data processing is performed in a manner similar to our previous reductions for ExploreNEOs (Trilling et al. 2010) and our other observations of NEOs (Mommert et al. 2014a, 2014b). The BCDs are downloaded from the Spitzer archive, and we use IRACproc/mopex in the moving object mode to mosaic the frames relative to the NEO position. For objects that are relatively bright compared to the background, and for low background fields, the outlier rejection process in the mosaicking software is sufficient to remove background objects. For fainter objects or regions of relatively high background, we often must subtract the background field before constructing the NEO mosaic. We first mosaic the data in the standard (nonmoving object) mode to construct an image of the background field (the NEO is masked from the individual frames before making the background mosaic). Then we subtract the background image from each BCD and mask any remaining strong residuals from bright background objects. We then construct the NEO mosaic from the background-subtracted BCDs. Once the NEO mosaic has been constructed, we perform aperture photometry to measure the NEO flux. We use an aperture correction derived from observations of standard stars during the Warm Mission to determine the flux.

Each mosaic is examined during the photometry process to ensure that there are no image artifacts present that will affect the photometry, and that the NEO is found in the field and not confused with background objects. The Spitzer Science Center produces a data delivery approximately every 2 weeks that contains ∼14 new NEO observations from the campaign that ended 2 weeks prior. In most cases, we can extract fluxes and perform thermal modeling within a few days of data delivery.

The measured flux density includes both thermal emission from the asteroid's surface and reflected solar light. The latter is generally small, but in some cases may be significant. We subtract this contribution from the measured flux densites using the method described by Mueller et al. (2011), but with an updated reflectance ratio between the optical and IR wavelengths of 1.6 ± 0.3, which is based on measurements from WISE by Mainzer et al. (2011a), who assumed, as do we, that NEO albedos are equal at 3.6 and 4.5 μm. This albedo ratio can differ throughout the NEO population, and our error bar on this ratio captures this variance. We subtract the (estimated) reflected light contribution before calculating the color correction because reflected sunlight has a different spectral shape than asteroid thermal emission; the reflected light spectral shape is very close to the one used in the IRAC flux calibration, so the corresponding color corrections for the reflected sunlight are negligible.

In our previous program, ExploreNEOs, we corrected solar flux densities by converting from IRAC in-band flux densities to monochromatic flux densities by integrating over the model surface temperature distribution and the product of the model spectral energy distribution and the IRAC CH2 bandpass (Mueller et al. 2011), requiring an iterative approach. By contrast, in the current program, we use a different approach in which we color-correct locally emitted thermal flux density from monochromatic to IRAC CH2 in-band fluxes using the exact model temperature of that surface element at the time the Planck equation is integrated over the surface temperature distribution. In other words, we derive separate color corrections for each spot on the surface, similar to the approach used in Mainzer et al. (2011b). Using this approach, the color correction is properly iterated as part of the fitting process. The derived color corrections are in the range 5%–20%; the only uncertainty associated comes from the ratio of reflectivity in the infrared compared to the optical (Mueller et al. 2011), and our error bars (described above) already capture our lack of knowledge of this reflectance ratio (Mueller et al. 2011). In the measured 4.5 μm flux, the reflected light fraction is between 0.01% and 29%, with a mean contribution of 1% and a median value of 0.7%. This new method produces color corrections that agree at the 1%–2% level with our previous approach (Mueller et al. 2011).

We obtain our optical H magnitudes from the ASTDyS¹⁶ database and improve them using statistical corrections provided by Pravec et al. (2012) over the given range of H magnitudes; uncertainties are also adopted from that work. For the range of H > 20.0, we apply a generic offset of 0.3 mag (we reduce the brightness relative to the catalog) and adopt an uncertainty of 0.3 mag, which is supported by Pravec et al. (2012), Jurić et al. (2002), and Hagen et al. (2016). In the future we will replace ASTDyS optical data with more accurate H magnitudes from ground-based surveys, where available.

Because we model single-band IRAC data, we are unable to use η as a free-floating fitting parameter. Instead, we have to make use of existing NEO thermal-infrared data to derive a robust estimate for η. As a result of asteroid surface roughness, thermal emission from the surface can be preferentially directed toward the Sun, which is referred to as infrared beaming. Energy conservation therefore requires that an observer at a high phase angle observes less thermal radiation than would be received from a smooth object, leading to a trend of increasing η with solar phase angle α. This trend has been quantified in the literature (Delbo' et al. 2003; Wolters et al. 2008; Mainzer et al. 2011a) using linear relations between the parameters. In general, all of ExploreNEOs, NEOSurvey, and Warm NEOWISE must use such a derived (not fitted) η, though NEOWISE was able to fit η directly for data obtained during their cryogenic mission, where some NEOs were detected at multiple thermal wavelengths, and both ExploreNEOs and Warm NEOWISE can fit both CH1 and CH2 thermal fluxes for very low albedo objects.

We also use a linear relation to derive the best-fit η based on its phase angle α, but here we take a more sophisticated approach to sample the structure and full range in the η distribution of NEOs. Based on values of η from NEOWISE (Mainzer et al. 2011a), which provides the largest published sample of fitted η values, we derive a new relation. Phase angles of the observed objects are derived by querying the Minor Planet Center database for WISE observations and represent average values. Figure 4 (left) shows the η data, the literature η–α relations, and our newly derived relation (see next paragraph). The NEOWISE fractional uncertainties in η are uniformly distributed throughout this space. The high degree of scatter in η across the whole range of phase angle is an indication of the heterogeneity in the thermal properties of NEOs and is a reminder that a simple linear relation is not sufficient to describe the variety in η. Instead, the whole range of η values must be taken into account in a statistical approach to properly interpret thermal-infrared data. Here we derive a formalism that provides the most likely η as a function of the phase angle, as well as an uncertainty distribution that samples the whole range of possible η values.

Zoom In Zoom Out Reset image size

From the left panel of Figure 4, we observe that there are no NEOs with η < 0.6 or η > 3.1 (0.6 is the smallest value that η can ever take on; Mommert et al. 2012). Furthermore, there is a clear grouping of low-η NEOs but a lack of low-η NEOs at high phase angles. A comparison of the η distributions derived by Mainzer et al. (2011a) and Wolters et al. (2008) shows that the large population of NEOs with low η is not a result of the thermal-infrared nature of the NEOWISE survey, as it is also present in the optically selected sample by Wolters et al. (2008). For the η relation used in this work, we adopt a slope of 0.01 deg⁻¹, which is in agreement with previous estimates and the slope of the lowest η values (we show below that a slightly different slope will not affect the results). We choose the intercept of our linear relation in such a way that it divides the sample of η values in Figure 4 (left) into two equally sized subsamples; the probability to overestimate η is the same as to underestimate it. Our new relation is $\eta (\alpha )=(0.01\,{\deg }^{-1})\alpha +0.87$. We define the residual Δη as the difference between the measured values of η and a specific η(α). The right panel of Figure 4 shows that Δη for the individual η(α) relations agree well with one another; despite the slightly different slopes, the shapes of the Δη histograms align. We fit a lognormal distribution to the residuals using a least-squares method. The fitted distribution has parameters μ = 0.655 and σ = 0.628, using the notation of Limpert et al. (2001). The median of the Δη distribution is 0.675. The upper and lower 1σ confidence interval of η uncertainties spans the range from −0.3 to +0.55, relative to the nominal value of η as derived by our η(α) relation (3σ interval: −0.55 to +1.7). Overall, the η–α parameter space shown in Figure 4 is poorly sampled, and we assume that the averaged lognormal distribution shown in the right panel of Figure 4 applies across all phase angles. The data in hand do not warrant or allow any more detailed treatment of these uncertainties.

Uncertainties in diameter and albedo are derived in a Monte Carlo approach very similar to the one described by Mueller et al. (2011). We perform 10,000 simulations with randomized parameters, including the measured flux densities that are varied based on the measured uncertainties; the absolute magnitude H, with Gaussian uncertainties from Pravec et al. (2012) or, for H > 20, as a Gaussian with 1σ uncertainties of 0.3; and the optical/IR reflectance ratio, also varied as a Gaussian. The beaming parameter η is varied within the lognormal distribution described above. Note that occasionally values η > 3.1 or η < 0.6 are drawn from the distribution; these extreme values are rejected. Clipping the η distribution affects the final results in a negligible way and is in agreement with the observations obtained by Mainzer et al. (2011a).

We adopt the median of the diameter and albedo distributions as our nominal solutions. We derive 1σ and 3σ confidence intervals as those ranges of values that bracket 68.3% and 99.7% of the values around the median, respectively. In Figure 5 we show representative cases of the diameter and albedo probability distributions that result from our thermal modeling.

Zoom In Zoom Out Reset image size

The beaming parameter η is not well constrained (Figure 4) and drives the overall uncertainties in our model results. For the results presented here, more than 90% of the total uncertainty in diameter results from uncertainty in η. For albedo, 80%–90% of the total uncertainty is due to η uncertainties, with uncertainties in H responsible for most of the remaining uncertainty. Uncertainties in other aspects (such as emissivity) are about 10 times smaller and are therefore negligible.

Finally, we make no attempt to correct for intrinsic NEO light curves for the targets in this program. The elapsed time of our observations spans the range from 10 minutes to 3.2 hr. For the longest exposures, we might integrate over an entire asteroid rotation period, as periods of seconds to a few hours are common for small asteroids, in which case we derive the mean diameter from our observations. For short exposures, we could easily sample just one part of the light curve, in which case our solution may represent the effective spherical diameter that corresponds to a light-curve minimum or maximum, thereby under- or overestimating the true mean diameter.

We estimate how many of our observations are affected by light-curve effects with a Monte Carlo approach. Taking into account the measured light-curve amplitude distribution of known NEOs (Warner et al. 2009), we randomly sample sine curves with different amplitudes to simulate the offset from the light-curve-averaged thermal flux of the target in our observations and produce diameter and albedo solutions that correspond to this instantaneous sampling. We find that 3%–6% of all trials have diameter discrepancies greater than our quoted diameter uncertainty of 40%. Hence, we expect the same fraction of our targets to have diameters affected by light-curve-induced discrepancies greater than 1σ. For our albedo solutions, 5%–9% of trials have solutions greater than our quoted albedo uncertainty of 70%. If we smooth the light curve over 20% of its period (which is a more realistic approximation of the situation), we find that fewer than 5% of our diameters and fewer than 9% of our albedos are affected by this light-curve sampling.

Derived albedos that are (much) greater than around 0.5 are suspicious, as such values are rarely seen for asteroids for which albedos have been reliably measured (for radar or spacecraft measurements, for example, as opposed to radiometric approaches), and are at odds with expected surface compositions. The most likely explanations for such high albedos are H magnitudes that are significantly incorrect and/or (optical or thermal) light-curve amplitudes and periods that produce a high albedo solution. In both cases, additional ground-based photometry can allow us to correct our albedo solution. We have initiated such a program using a range of ground-based telescopes, and in a future paper we will report revised optical magnitudes and/or light curves for these high-albedo objects, along with recalculated albedos. (Note that poor H values and/or light curves could also produce anomalously low albedos, which may be harder to identify. We will carry out ground-based observations of the objects for which we find the lowest albedos to determine whether a correction is also needed for those objects.)

A small number of our targets will be observed by NEOWISE—it is impossible to calculate which targets will overlap in the future, since the NEOWISE SPICE kernel (orbit) is only published a few months ahead—and we will use those targets as cross-checks. Furthermore, there are more than 100 NEOs that have been observed by NEOWISE that were also observed in our ExploreNEOs program (Trilling et al. 2016). For these common objects, we use our updated thermal model pipeline and used only CH2 data (as if the ExploreNEOs data were obtained in this NEOSurvey program) to rederive the diameter and albedo from our ExploreNEOs measurements, and we compare these to the published NEOWISE diameters and albedos (Mainzer et al. 2011a). The mean ratio of ExploreNEOs values to NEOWISE values for (diameter, albedo) is (0.9, 1.28). The median value in this ratio of ExploreNEOs values to NEOWISE values for (diameter, albedo) is (0.9, 1.02). This overall agreement, based on independent data measurements and modeling, is quite good, although we caution that for any individual object significantly different solutions may exist. For diameter, 70% of the newly modeled ExploreNEOs targets have 1σ error bars that overlap with NEOWISE values; 98% of the targets agree at 3σ. For albedo, 66% of the newly modeled ExploreNEOs targets have 1σ error bars that overlap with NEOWISE values; 98% of the targets agree at 3σ. This suggests both that the errors are largely random and that our error bar estimations are appropriate. The largest deviation is in albedo, where uncertainties in H magnitude typically introduce large errors, and where the two approaches frequently use somewhat different H magnitudes for their modeling. In the future we will carry out a more comprehensive comparison, considering targets that overlap with both ExploreNEOs and NEOSurvey, and using common H values for the thermal modeling of the different thermal measurements.