FIVE YEARS OF MID-INFRARED EVOLUTION OF THE REMNANT OF SN 1987A: THE ENCOUNTER BETWEEN THE BLAST WAVE AND THE DUSTY EQUATORIAL RING

About 10 years after its explosion on 1987 February 23, SN 1987A has evolved from a supernova (SN), when its radiative output was dominated by the release of radioactive decay energy in the ejecta, into a supernova remnant (SNR), when its radiative output became dominated by the interaction of its blast wave with the inner equatorial ring (ER). The ER is located at a distance of about 0.7 lt-yr from the center of the explosion and could have been produced by mass loss from a single rotating supergiant (Heger & Langer 1998) or by a merger event in a binary system that also formed the two outer rings (Morris & Podsiadlowski 2009). The transition from SN to SNR was observed at wavelengths ranging from radio to X-rays and is summarized in Figure 18 of Bouchet et al. (2006).

In this paper, we report on the continuing evolution of the ∼5–30 μm spectrum and the 3.6, 4.5, 5.8, 8.0, and 24 μm photometric fluxes from SN 1987A, spanning the ∼5 year period from day ∼6000 until day ∼8000 after the explosion. Initial reports and analysis of the IR evolution were presented by Bouchet et al. (2004, 2006) and Dwek et al. (2008). In Section 2, we present the IR data obtained by the Spitzer satellite. The evolution of the IR emission and the dust composition are presented in Section 3. In Section 4, we derive the plasma conditions from the IR observations. The comparison of the IR emission with the X-ray emission and the evolution of their flux ratio are discussed in Section 5. A brief summary of the paper is presented in Section 6.

Spitzer's observations using the Infrared Array Camera (IRAC), Infrared Spectrograph (IRS), and Multiband Imaging Photometer for Spitzer (MIPS) were conducted annually during the first two years (2004 and 2005) and then roughly every six months for the remainder of its cryogenic mission. SN 1987A was also observed incidentally with IRAC and MIPS by the SAGE survey (Meixner et al. 2006) and other projects during 2005.

Figures 1 and 2 illustrate the IRAC data at the final observations (∼day 8000). The IRAC data are all of good quality, i.e., the data are not adversely affected by any artifacts induced by bright sources in the field or in prior observations. Initial tests of generating mosaics from the basic calibrated data (BCD) showed no significant difference from the post-BCD mosaics. Therefore, aperture photometry was performed on the post-BCD mosaics in order to obtain the broadband evolution of SN 1987A at 3.6, 4.5, 5.8, and 8 μm. These measurements and their statistical uncertainties are shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 1 and 2 also illustrate the MIPS 24 and 70 μm images of SN 1987A at day 7983. There is no evident point source at 70 μm. As for IRAC, there seemed to be no advantage to reprocessing the BCD, thus aperture photometry was performed on the post-BCD mosaics. The 24 μm flux densities are shown in Figure 3. At day 7983, a 70 μm upper limit was determined by adding an artificial point source to the MIPS image at the location of the SN. We started with a source flux that was clearly visible above the background confusion and then reduced the flux until the source was no longer detected using the FIND routine of the IDLASTRO library. The results indicate the S(70μm) < 0.09 Jy.

A significant complication with the IRAC and MIPS photometry is that the SN is not resolved from the nearby Star 2 and Star 3. If the stars are assumed to have mid-IR magnitudes that are equal to the K magnitudes (15.06 and 15.80, respectively, from Walborn et al. 1993), then their combined IR flux densities should be 0.41, 0.26, 0.16, 0.09, and 0.01 mJy at 3.6, 4.5, 5.8, 8, and 24 μm, respectively. These assumptions may be accurate for Star 2, which is identified as a B2 III star (Walborn et al. 1993; Scuderi et al. 1996). However, Walborn et al. (1993) report both an IR excess and a variability by at least Δm = 1 mag for Star 3. The fluxes in Figure 3 are plotted after subtraction of the emission of Stars 2 and 3 under these assumptions.

IRS observations were reduced starting with the BCD and using the basic steps outlined in the IRS Data Handbook. For the low-resolution (R ∼ 60–120) data, we constructed a two-dimensional background frame using the observations of a specific off-source target position. IRSCLEAN_MASK was used to identify and clean rouge pixels in the background-subtracted BCD frames. The cleaned BCD data were averaged for each nod (positions of the source at 1/3 or 2/3 the length of the slit), and spectra were extracted from these super BCDs using the Spitzer IRS Custom Extraction (SPICE) software.¹⁰ After extraction, the spectra orders were merged to produce the full spectra shown in Figure 3 (top right panel). The high-resolution (R ∼ 600) data were processed in a similar fashion. The major difference was that instead of using a background from a dedicated off-source target, the background is derived from the average of two pointings that flank the SN on opposite sides (15'' and 25'' away for the short high-resolution (SH) and long high-resolution (LH) spectra, respectively). The position angle of these pointings varied as dictated by the scheduling of the observations, but the distance of the pointings from the SN remained constant. For both the high- and low-resolution data, we used the Spectroscopic Modeling Analysis and Reduction Tool¹¹ (SMART) to fit and extract line fluxes.

Figure 3 depicts the evolution of the IRAC and MIPS broadband photometry (top left panel), and the low-resolution IRS spectrum of SNR 1987A (top right panel), as a function of time. The relatively strong background may affect the continuum at λ>20 μm and the [Ne iii] 15.6 μm and [S iii] 18.7, 33.5 μm emission lines. Line emission from SN 1987A and its surrounding medium will be presented elsewhere (R. G. Arendt et al. 2011, in preparation). The figures show that the IR flux is rising smoothly in all bands, increasing by a factor of ∼3 between days 6190 and 7980. The bottom left panel of the figure shows the spectra normalized by least-squares fit to the one obtained on day 7554. The spectrum maintained a nearly constant shape throughout the 1800 days of observations. This constancy enables the use of any IRAC or 24 μm MIPS band as a proxy of the total integrated IR flux from the ER. The bottom right panel in the figure shows the ratio between the integrated IR flux and the 24 μm flux from the ER. The ratio is constant with a value of 1.64 to within <10%. The constancy in the spectral shape also places strong constraints on the parameters of the hot X-ray gas in which the grains are embedded.

Figure 4 (left panel) shows that the mean spectrum of the ER can be well fitted by astronomical silicate grains radiating at a single temperature of ∼180 K. The relative ratio between the two peaks of the silicate features provides strong constraints on the temperature. The mass of the silicate dust is 1.2 × 10⁻⁶ M_☉ on day 7554 and should scale linearly with the IR flux for other epochs (see Table 1). The figure also shows the presence of an IR excess at wavelengths between ∼5 and 8 μm that cannot be attributed to emission from the silicate grains. This excess emission was first detected by Bouchet et al. (2006), but was not further pursued because of the shorter integration time and larger noise of that early Spitzer spectrum. An unidentified K-band (2.2 μm) continuum at a level of ∼0.4 mJy, detected by Kjær et al. (2007), may be part of this emission component. This emission component, now clearly present in all epochs, can be approximated in the 5–8 μm wavelength region by a ν^α power law with a spectral index of α ≈ −1.8. This rules out thermal bremsstrahlung, which has a well-defined spectral index of α ≈ −0.1 in this frequency regime, as a possible source. Synchrotron radiation can also be ruled as a source of this emission component. The observed radio synchrotron emission (Zanardo et al. 2010) has a spectral index of ≈−0.8, and its extrapolation from a value of 300 mJy at 1.4 GHz on day 7000 will make a negligible contribution to the 8 μm emission around that epoch. We therefore conclude that this emission arises from an additional hot dust component that is tightly correlated with the silicate emission. Its spectral shape should be featureless at longer wavelengths as not to reduce the goodness of the fit to the observed 9.7 and 18 μm silicate emission features. The right panel of the figure shows the fit of a blackbody spectrum, modified by a λ⁻² emissivity law and a temperature of ∼350 K, to this component. A λ⁻² emissivity is consistent with metallic or carbon grain composition, rather than dielectric or silicate composition (Bohren & Huffman 1983).

Zoom In Zoom Out Reset image size

The modified blackbody is of course an idealized representation of this secondary dust component. Figure 5 shows a two-component fit to the IR spectrum when the secondary component was taken to consist of either carbon, pure iron, magnetite (Fe₃O₄), or iron sulfide (FeS) grains. All components, combined with the main silicate dust component, provided acceptable fits to the overall IR spectrum of the ER. The temperature of the secondary dust component is significantly higher than that of the silicate grains, ranging from ∼370 K for pure iron grains to ∼460 K for carbon grains. Some other possible secondary dust components, in particular iron and aluminum oxides, were ruled out because these grain species produce strong emission features at wavelengths longward of ∼8 μm that are not present in the Spitzer IRS spectra.

Zoom In Zoom Out Reset image size

The problem of accounting for the widely different temperatures of the two dust components, even though they may be collisionally heated by the same gas, will be addressed in Section 4.

Figure 6 shows the plasma densities required to collisionally heat the different grain types to their respective temperatures given in Figure 5 as a function of grain radius. The grains are assumed to be embedded in an ionized plasma with a temperature of 5 × 10⁶ K, reflecting the average temperature of the soft X-ray component (see Section 5). If we require both dust components to arise from the same gas, then only the carbon grains support a range of plasma densities (enclosed between the two horizontal dashed lines) that overlap with that of the silicate grains. So the IR spectrum of the ER is consistent with the emission from a population of small (a ≲ 0.05 μm) carbon grains intermixed with a population of large (a ≳ 0.20 μm) silicate grains embedded in a plasma with a density n_e ≈ n_H ≈ (2–4) × 10⁴ cm⁻³. Requiring a common environment for the two dust components rules out the other dust compositions (Fe, Fe₃O₄, and FeS) as viable dust candidates. Alternatively, the emission could arise from these non-carbonaceous grains, provided they reside in a significantly denser phase (n_e ≳ 2 × 10⁵ cm⁻³) of the X-ray emitting gas. This phase then has to be tightly correlated with the lower density one in order to explain the tight correlation between the IR light curves of the two dust emission components.

Zoom In Zoom Out Reset image size

In principle, the second dust component could consist of grains that have formed in the SN ejecta. There is ample evidence for the formation of dust in SN 1987A (Moseley et al. 1989; Lucy et al. 1991; Wooden et al. 1993). In particular, the disappearance of IR fine structure line emission from Fe has been interpreted as evidence for the formation of iron dust in the ejecta (Dwek et al. 1992; Wooden 1997). However, it is unlikely that this second dust component is ejecta dust. The total IR luminosity of this second dust component is about 10³⁶ erg s⁻¹, calculated for a distance of 50 kpc to the SN. This luminosity is higher by a factor of ∼10² from that expected from the radioactive decay of ⁴⁴Ti, assuming that ∼10⁻⁴ M_☉ of ⁴⁴Ca was produced in the explosion (Woosley et al. 1989). So the second dust component can only be heated by a reverse shock expanding through the ejecta. As shown below, the IR light curve depends on the rate at which the mass of the dust is swept up and its rate of destruction. These rates depend on the medium morphology and density, and the shock velocity. The tight correlation between the light curves is unlikely to arise from separate shocks propagating into vastly different environments such as the ejecta and the ER. It is therefore very unlikely that the emission from the second dust component arises from SN-condensed dust in the ejecta.

The potential presence of carbon with the silicate dust in the ER is surprising, defying the common paradigm that either carbonaceous or silicate dust can form in the wind of the progenitor, depending on the C/O abundance ratio in the outflow. The observations therefore suggest either that CO formation did not exhaust all the carbon in the outflow or that the co-existence of the two dust components is a manifestation of a binary origin of the ER. The progenitor of SN 1987A could have shared a common envelope with a less massive carbon-rich star. The interaction between them could have led to the formation of the ER as well as the two outer rings observed by the Hubble to be in the vicinity of the SN (Morris & Podsiadlowski 2009).

A rise in the IR emission is expected as the shock sweeps up more dust during its propagation into the dusty circumstellar medium. The rate of increase of the IR emission depends on the morphology of the medium into which the shock is expanding, the expansion rate of the shock, and on the possibility that dust may be destroyed in the postshock gas. The expected increase of the IR emission with time is schematically illustrated in Figure 7. The left column represents the case in which the grain lifetime is sufficiently long that it can be ignored, that is, τ_sput > t', where t' is the residence time of the dust in the hot gas. The right column, for which τ_sput < t', represents the case in which the grain lifetime is sufficiently short that an equilibrium is established between the rate of grain destruction and the rate at which the dust is replenished by the expanding shock. For both cases, the figure depicts a shock expanding into a one-dimensional protrusion (top), a two-dimensional disk, or ring (middle), and a homogeneous sphere (bottom). The dependence of the volume of the swept-up gas is presented for two different phases in the evolution of a blast wave: (1) the free-expansion phase, during which the velocity, v_f, is constant and (2) the Sedov–Taylor phase, which commences after most of the kinetic energy of the explosion has thermalized, during which the velocity v_s ∼ t'^−3/5 (Zel'Dovich & Raizer 1967).

Zoom In Zoom Out Reset image size

When τ_sput > t', following the left column of Figure 7, the intensity of the IR emission scales as the mass of swept-up dust, which for a constant dust-to-gas mass ratio simply scales as the volume of the shocked gas. For a freely expanding blast wave, the volume of the shocked gas, V_f, will increase as t', t'², and t'³, for the three morphological scenarios, respectively. During the Sedov phase, the volume of the shocked gas, V_s, will be increasing more slowly as t'^2/5, t'^4/5, and t'^6/5, respectively.

When τ_sput < t', following the right column of Figure 7, the dust column density through the shock will eventually reach a constant value, roughly equal to v_s × τ_sput, where τ_sput is the grain destruction time in the postshock gas. The IR intensity behind a freely expanding shock will therefore reach a constant value if the shock is expanding into a one-dimensional protrusion and will increase as t' or t'² when it expands into a ring or sphere, respectively. During the Sedov phase, the volume of the shocked gas, V_s, will actually be decreasing as t'^−3/5, t'^−1/5, or increasing more slowly as t'^1/5, for the three morphological scenarios, respectively.

Figure 8 depicts the evolution of the IRAC and 24 μm MIPS fluxes. The fluxes were fitted by a power-law fit of the form (t − t₀)^α, where t is the time since the explosion and α and t₀ are free parameters. The value of t₀ can be interpreted as the time when the shock first encountered the ring. The values of {α, t₀} were found to be {0.6, 5700 days} and {0.87, 5500 days} for the IRAC and MIPS fluxes, respectively.

Zoom In Zoom Out Reset image size

The dust flux scales as the volume, and the value of α derived from the light curve is consistent with V_s ∼ t'^4/5 (left column of Figure 7), indicating that grain destruction is unimportant. More specifically, the 24 μm flux increases at a rate that is consistent with a Sedov–Taylor blast wave expanding into a circumstellar disk/ring. As shown below, the inferred dust lifetime is consistent with currently used values for the sputtering rates of dust grains in hot plasmas.

The value of t₀ is consistent with the time inferred from X-ray observations that showed a deceleration of the radial expansion of the SN around day 6000 (Racusin et al. 2009) and the upturn in the soft X-ray luminosity that occurred around day 6200 (Park et al. 2006). The IR and X-ray observations therefore offer a consistent scenario in which both fluxes originate from the dust and gas shocked by the blast wave that has started to penetrate the main body of the ER around day 6000 after the explosion.

5.1. The Evolution of IRX

Analysis of the Chandra X-ray observations of the ER shows that the spectrum is generated by two distinct emission components: a soft component with temperatures ∼0.3–0.6 keV and a harder component characterized by temperatures between ∼2and5 keV (Zhekov et al. 2009; Park et al. 2006). The emission from the soft component arises from the gas heated by the shock that is transmitted into the ER. The hard component arises from the lower density gas interior to the ER that is shocked twice: once by the advancing SN blast wave and second time by the shock reflected from the ER. We assume here that the IR emission arises from the soft X-ray component generated by the shock that is penetrating the ER.

Figure 9 (left panel) depicts the evolution of the soft X-rays and IR fluxes and their ratio, IRX, as a function of time. The IRS observations only cover the ∼6804–8000 days epoch, so we used the MIPS 24 μm observations that cover the additional ∼6190–6800 days epoch as a proxy for the integrated 5–30 μm emission. The figure also depicts the day 6067 datum point, for which we used the scaled mid-IR photometry from Gemini T-ReCS 10.4 μm observations (Bouchet et al. 2006) as a proxy for the integrated IR flux.

Zoom In Zoom Out Reset image size

The radial expansion of the soft X-ray component slowed down around day ∼6000, which was interpreted as evidence for the transition of the blast wave from the free expansion to the Sedov–Taylor phase as it enters the main body of the ER (Park et al. 2005; Racusin et al. 2009). The T-ReCS datum suggests that this transition has also manifested itself in the IR emission, re-enforcing the consistency between the X-ray and IR evolution of the remnant.

The evolution of IRX is shown in more detail in the right panel of the figure. It has a roughly constant value of ∼2.5 throughout the entire ∼6000–8000 days epoch. An exception is the 6190 datum point. We believe that the anomalously large value of IRX at this point is caused by a relative dip in the X-ray flux starting near day 5800 and ending with the rise in the soft X-rays near day 6200. The large value of IRX also suggests that a significant fraction of the refractory elements in the ER are locked up in dust.

A constant value of IRX indicates that the IR and X-ray fluxes evolve at the same rate. This requires the dust-to-gas mass ratio to be uniform throughout the medium into which the shock is expanding and that neither grain destruction nor gas cooling affect the IR and X-ray emission from the shocked gas. Alternatively, the constancy in the ratio could be the result of both processes, grain destruction and gas cooling, being important. In this case, the IR and X-ray emission have attained a steady state in which the losses to grain destruction and cooling are balanced by the influx of new material such that the volume of radiating gas and dust remains constant in time. However, the rate of increase in the IR fluxes suggests that the volume of the shock-heated dust is increasing, so that the second scenario is a priori not likely to be the cause for the constancy in IRX. In the following, we examine in more detail the physical condition of the shocked gas that can lead to a constant IRX.

5.2. Grain Destruction in the Hot Gas

The lifetime of a dust grain of radius a moving at a high velocity through a hot gas with temperatures between ∼10⁶and10⁷ K is given by Nozawa et al. (2006) and Dwek et al. (2008):

$$\begin{equation} \tau _{{\rm sput}}({\rm yr}) = A_{\rm sput}\, {a(\mu \mbox{m})\over n_{\rm H}(\rm \, cm ^{-3})}, \end{equation} \tag{ 1 }$$

where A_sput = 3.3 × 10⁵, 9.0 × 10⁵,  5.6 × 10⁵,  2.8 × 10⁵, and 3.9 × 10⁵, respectively, for silicate, carbon, iron, FeS, and Fe₃O₄ dust, calculated for grains moving at velocities ≳500 km s⁻¹ through a 5 × 10⁶ K gas. The values of the velocity and temperature characterize those of the shock giving rise to the soft X-ray component. Figure 10 depicts the sputtering lifetime of the dust as a function of grain radii and gas density. The figure shows that the lifetime of the silicate grains is between ∼4 and ∼15 years, depending on the combination of grain size and gas density. In contrast, for the secondary dust component, this same range of gas temperatures and densities yields significantly shorter grain lifetimes, between ∼0.4 and ∼1 year, which are too short to explain the observed temporal evolution of the secondary dust component.

Zoom In Zoom Out Reset image size

5.3. Cooling of a Shock-heated Dusty Plasma

Figure 11 depicts the equilibrium atomic cooling rate, Λ(T_e), of a gas with ER composition (thick line, left panel) as a function of electron temperature, T_e. The ER abundances were taken from Mattila et al. (2010) and were derived from Anglo-Australian Telescope and Very Large Telescope observations of the SN between days ∼1400 and ∼5000 after the explosion. The thin lines in the figure show the cooling of the plasma via gas–grain collisions, assuming a pure silicate dust composition. The silicate dust-to-gas mass ratio was assumed to be 0.0029, half the solar value. The different curves assume that all the dust is represented by grains with the indicated radii. Dust cooling dominates the cooling via atomic transitions at temperatures above ∼10⁶ K, depending on grain size. The value of IRX is the ratio of these two cooling rates. From the fact that the initial value of IRX is only ∼2.5, we can conclude that the grain sizes in the ER must be relatively large, with typical radii of ≳0.2 μm, consistent with the sizes required to heat the dust to a temperature of ∼180 K (see Figure 6). Using this lower limit on the grain radii in Equation (1) gives a postshock density of about (2–4) × 10⁴ cm⁻³, where we assumed an upper limit of 0.5 μm on the grain radii.

Zoom In Zoom Out Reset image size

The equilibrium cooling may not represent the actual cooling of a shocked gas, which may be far from ionization equilibrium. The right panel in Figure 11 depicts the temperature profile behind shocks of different velocities expanding into a homogeneous dusty gas with a preshock hydrogen density of 10⁴ cm⁻³ and ER composition. The shock models include the effects on non-equilibrium ionization (NEI) and the effect of dust cooling. The ∼600–700 km s⁻¹ shocks well represent the 3–7 keV temperature of the soft X-ray-emitting plasma (Park et al. 2006; Zhekov et al. 2009). The figure shows that the cooling time of this shocked gas is about 12–20 years for a preshock gas density of n_H = 10⁴ cm⁻³. The IR observations suggest that the postshock density is ∼(2–4) × 10⁴ cm⁻³, giving a gas cooling time of ∼20–40 years.

5.4. Scenarios for the Constancy of IRX

5.5. The Mystery of the Secondary Grain Component

In this paper, we presented the mid-IR evolution of SN 1987A over a five year period spanning the epochs between days ∼6000 and 8000 since the explosion. Its radiative output during this epoch is dominated by the interaction of the SN blast wave with the pre-existing ER. The main results of this paper can be briefly summarized as follows.

• 1.  
  The ∼8–30 μm mid-IR spectrum is dominated by emission from ∼180 K silicate dust, collisionally heated by the hot X-ray emitting gas with a temperature and density of ∼5 × 10⁶ K and ∼(2–4) × 10⁴ cm⁻³, respectively. The mass of the radiating dust is ∼1.2 × 10⁻⁶ M_☉ on day 7554 and scales linearly with IR flux.
• 2.  
  A secondary emission component dominates the spectrum in the ∼5–8 μm region. Its intensity and spectral shape rule out any possible gas or synchrotron emission mechanism as the source of this emission. It must therefore be attributed to a secondary dust component radiating at temperatures above ∼350 K.
• 3.  
  The overall shape of the ∼5–40 μm dust spectrum has not changed during the observations, suggesting that the density and temperature of the soft X-ray emitting gas have not significantly changed during more than five years of IR observations. The constancy in the spectral shape of the IR emission also suggests that the mass ratio of the silicate to the secondary dust component remained roughly constant during this period.
• 4.  
  The evolution of the IRAC and MIPS fluxes can be described by a power law in time since the first shock-ER encounter. The silicate emission increases as t^0.87, consistent with X-ray observations, suggesting that the blast wave has transitioned from a free expansion to the Sedov phase and is now expanding into the main body of the ER.
• 5.  
  The infrared-to-X-ray flux ratio, IRX, is constant with a value of ∼2.5 throughout this epoch. The magnitude of IRX shows that the cooling of the shocked gas is dominated by IR emission from the collisionally heated dust with radii ≳0.2 μmand that a significant fraction of the refractory elements in the ER should be depleted onto dust.
• 6.  
  The constancy of IRX is consistent with the premise that neither grain destruction by sputtering nor cooling of the shocked gas has played a significant role in the evolution of the IR and X-ray emission. This scenario is consistent with the sputtering rates currently used in the literature as expressed in Equation (1).
• 7.  
  The presence of a secondary dust component, radiating at significantly higher temperatures than the silicate dust, suggests that the grain radii or IR emissivities of this component must be significantly smaller than those of the silicates. Their sputtering lifetime could therefore be significantly shorter than that of the silicate grains and their evolution distinctly different from that of the silicates, especially around the ∼6100 days time interval. However, a significantly shorter grain lifetime contradicts the nearly constant or slightly decreasing flux ratio between the hotter secondary component and the silicate grain component. Hence, the nature of the secondary dust component remains a mystery.

Five years of continued monitoring of SN 1987A have led us to conclude that there is no current evidence for grain destruction. This conclusion stands in contrast to the results of our earlier analysis (Dwek et al. 2008), which were based only on two epochs of data at early times, and which we now see to be somewhat anomalous. Continuing observations will reveal evidence for grain destruction in the ER and may elucidate the nature of the mysterious secondary dust component.

We benefited from useful comments by Dick McCray and Svetozar Zhekov. We also thank the referee Diane Wooden for her careful reading of the manuscript and many constructive suggestions. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA.