THE WIRED SURVEY. III. AN INFRARED EXCESS AROUND THE ECLIPSING POST-COMMON ENVELOPE BINARY SDSS J030308.35+005443.7

The WISE InfraRed Excesses around Degenerates (WIRED) survey is designed to detect infrared (IR) excesses around white dwarfs (WDs) using photometry from the WISE mission (for a complete description of the WISE mission, see Wright et al. 2010). Dust, low-mass companions, and cyclotron radiation from accreting magnetic WDs all emit at mid-IR wavelengths, providing a rich variety of sources to be discovered. Because of the all-sky coverage of WISE, WIRED provides a more systematic and less biased search for IR excesses around WDs than those performed with targeted Spitzer observations.

In Debes et al. (2011a, Paper I), we characterized the IR excess discovered around GALEX 1931+0117 using WISE photometry. In Debes et al. (2011b, Paper II), we cross-correlated visible, near-IR, and WISE photometry of WDs discovered in the Sloan Digital Sky Survey (SDSS) DR7 Preliminary WD catalog (Kleinman 2010) to discover excesses around WDs due to dust disks as well as low-mass stellar and substellar companions. In the process, we have discovered a WD+dM system, SDSS J030308.35+005444.1 (hereafter SDSS J0303+0054), that shows an infrared excess above the photosphere of the M star companion. SDSS J0303+0054 has a W1 − W3 color of 1.8 ± 0.1, much redder than the average W1 − W3 = 0.5 ± 0.2 of M4–M6 dwarfs detected by WISE with W3 magnitudes having signal-to-noise (S/N) > 5 (Kirkpatrick et al. 2011).

SDSS J0303+0054 is an eclipsing post-common envelope binary first discovered as a WD+dM system in the SDSS DR4 WD catalog (Eisenstein et al. 2006; Silvestri et al. 2006). Pyrzas et al. (2009, hereafter P09) detected significant radial velocity variations in the SDSS J0303+0054 system, as well as eclipses of the WD from optical photometric monitoring. They determined that SDSS J0303+0054 is a DC+dM4 binary with a 3.2 hr orbit. From the combination of radial velocities and light-curve modeling, they constrained the mass of the WD, assumed to have a helium-dominated atmosphere, to 0.88–0.95 M_☉ and the companion's mass to between 0.22 and 0.28 M_☉. The T_eff of the WD is not well constrained due to a lack of absorption features that might otherwise be detected if not for the overwhelming flux of the M dwarf companion; however, P09 estimated that it must be <8000 K based on spectral decomposition of the two stellar components.

Observations at u, g, and i were also conducted by Parsons et al. (2010), showing almost complete eclipses of the WD at short wavelengths, a brightening of the system during eclipse of the WD at u and g, and ellipsoidal variations in i. Tappert et al. (2011a) also found evidence for multiple Hα components, which could be caused by a combination of activity/irradiation of the M dwarf and the accretion of material on the WD surface.

Based on current prescriptions of orbital migration after common envelope (CE) evolution (e.g., Schreiber & Gänsicke 2003), SDSS J0303+0054 should have left the CE with a period of ≈0.2 days and slowly evolved through magnetic braking to its current orbit (Zorotovic et al. 2010), where the M dwarf does not overfill its Roche lobe and hence no mass transfer through the L1 point should occur in this system. Once the companion fills its Roche lobe and mass transfer through the L1 point begins, the binary becomes a cataclysmic variable (CV) and circumbinary dust is a distinct possibility. Some CVs have been found to show evidence for cold dust in a circumbinary disk (Hoard et al. 2009). This dust could originate from condensation of gaseous or dusty material that is lost from the inner binary, for example during classical nova outbursts, in a wind from the accretion disk and/or secondary star, or as "spillage" from the mass transfer process between the two stars. Numerous single WDs show IR excesses due to dust (e.g., Farihi et al. 2010; Xu & Jura 2012; Debes et al. 2011b, and references therein), likely produced by the tidal destruction of asteroids perturbed into the WD Roche limit by remnant planetary systems (Debes & Sigurdsson 2002; Jura 2003; Bonsor et al. 2011; Debes et al. 2012). Given the large mass and relatively large cooling age of the WD in SDSS J0303+0054 (∼3 Gyr; P09), the presence of an infrared excess is an intriguing development in the study of this binary system.

We present follow-up observations of the SDSS J0303+0054 system in Section 2. In Section 3, we construct a spectral energy distribution (SED) of the WD+dM in order to accurately determine the SED of the IR excess. We construct some simple models of the excess in Section 4. In Section 5, we present our follow-up near-IR observations of SDSS J0303+0054 as well as time-series WISE photometry to determine the possible photometric variability in this binary and discuss the implications of our work in Section 6.

In addition to GALEX, SDSS DR7, Two Micron All Sky Survey (2MASS), UKIDSS, and WISE photometry, follow-up near-IR observations of SDSS J0303-0054 were obtained in clear conditions on 2011 October 20 and 21 with the 4.2 m William Herschel Telescope using the Long-Slit Intermediate Resolution Infrared Spectrograph (LIRIS; Manchado et al. 1998). The follow-up data were obtained to search for photometric variability on orbital timescales and to obtain additional contemporaneous photometry of the system. Additionally, we obtained the All Sky Single Exposure photometry⁷ and WISE 3-Band Cryo Single Exposure photometry⁸ in W1 and W2 of SDSS J0303+0054 in the WISE All-Sky and 3-Band Cryo catalogs in order to track possible short-term variability.

On the first night of the LIRIS run, a single nine-point dither pattern of images was obtained in the J, H, and K_s bands, with total exposure times of 18 s in each filter. This was followed immediately by 10 identical and consecutive sequences at H, then again at K_s. On the second night, 50 such consecutive sequences were executed at K_s only, spanning approximately 1.2 hr in total elapsed time. On both nights, four standard star fields (ARNICA, 11 stars total; Hunt et al. 1998) were observed in a similar manner for photometric zero-point calibration. The data were reduced in the standard manner, by subtracting a median sky from each image in the dithered stack, flat fielding (using sky flats), then averaging and recombining frames. All data taken in the K-band filter were flux-calibrated using the ARNICA K-band standard star photometry, and the error introduced by this should be significantly smaller than the 5% absolute calibration uncertainty.

LIRIS suffers from what is known as a detector reset anomaly, which appears in certain frames as a discontinuous jump (in dark current) between the upper and the lower two quadrants. To remove this unwanted signal, after flat fielding and sky subtraction, the detector rows were collapsed into a median column (with real sources rejected), and subsequently subtracted from the entire two-dimensional image. The resulting fully reduced frames exhibit smooth backgrounds, free of the anomalous gradient.

Figure 1 shows LIRIS K_s, WISE W1, W2, and W3 images of SDSS J0303+0054. It is detected with an S/N of >100, 43, 42, and 11, respectively. Investigation of the point-spread function (PSF) in each WISE band does not reveal any extension: the science target in each band appears to be consistent with the expected PSF FWHM. Furthermore, inspection of higher spatial resolution UKIDSS and LIRIS images shows no indication of any neighboring sources, limiting contamination from comparably bright sources to angular separations <035. Table 1 lists the GALEX, SDSS, near-IR, and WISE photometry for this system, the Modified Julian Date (MJD) of each observation after a heliocentric correction and its implied orbital phase for the binary, as well as our best-fit model described below.

Zoom In Zoom Out Reset image size

3.1. Empirical M Dwarf Photometry

3.2. Fitting Models to the SED

In order to determine the significance of the possible nature of the IR excess, one must accurately subtract off the emission from the two stellar components. We first converted the photometry of SDSS J0303+0054 into units of flux density as in Debes et al. (2011b), and then simultaneously fitted empirical SEDs of M dwarfs and model SEDs of WDs to the observed photometry. WD SEDs were constructed from He-rich cooling models kindly provided by P. Bergeron that include the GALEX and WISE bands (Bergeron et al. 1995; Holberg & Bergeron 2006). We compared our models with GALEX FUV through J photometry of SDSS J0303+0054, determined a median scaling between the observed and model photometry, and minimized the χ²_ν value to determine the best fit in log g, T_eff, d, and spectral type. The resulting best-fit parameters and their 98% confidence intervals are shown in Table 2. Our best-fit model requires a WD with log g = 8.5, T_eff = 8940 K, and an M4.5 companion at a distance of 134 pc—the expected photometry is listed in Table 1. The inferred cooling age of the WD from our models is 2 Gyr. The resulting gravity and distance to the WD are dependent on our assumption that the companion has the same luminosity as a typical field M4.5 dwarf and are thus relatively uncertain given the inherent uncertainties in our sample of field dwarf absolute magnitudes and colors. However, the observations of P09 quite accurately determined both WD and M dwarf radii and masses, which we use in Section 4 and list in Tables 2 and 3. Figure 2 shows a comparison of the observed photometry and our best-fitting model. The photometry shows excess in the H band, with all three sets of near-IR data showing varying levels of excess above the photosphere. However, the three near-IR epochs are marginally inconsistent with each other: relative to 2MASS, the LIRIS H photometry is inconsistent at the 4σ level, and the UKIDSS K photometry is inconsistent at the 3σ level. These are small enough to conceivably be caused by underestimating the photometric errors, or could signal photometric variability in the system. We address that further in Section 5. Subtracting off the model stellar SED, we find that the W1, W2, and W3 bands have excesses of 22σ, 24σ, and 10σ above the expected combined stellar photospheres at wavelengths >3 μm.

Zoom In Zoom Out Reset image size

Table 2. SDSS J0303+0054 Fitted Fundamental Parameters

Fundamental Parameter	Value	P09 Constraints
WD log g	8.5+0.5− 1.0	8.47–8.6
WD mass (M☉)	0.9+0.3− 0.6	0.88–0.95
WD Teff (K)	8940+820− 400	<8000
WD cooling age (Gyr)	2.0+0.8− 1.6	...
Companion spectral type	M4.5+1.0− 1.5	M4–M5.5
χ2ν	0.7	...
Inferred distance (pc)	134+206− 63	...

Download table as:  ASCIITypeset image

Table 3. SDSS J0303+0054 Adopted Fundamental Parameters

Fundamental Parameter	White Dwarf	M dwarf
Mass (M☉)	0.89	0.25
Radius (R☉)	0.009	0.26
Inferred Teff	8960	3020
Inferred distance (pc)	135	...
Luminosity (L☉)	4.6 × 10−4	5.0 × 10−3

Download table as:  ASCIITypeset image

We now investigate the origin of the IR excess around this system. IR excesses around presumed single WD systems arise from (1) cool, low-luminosity companions, (2) dust disks, or (3) cyclotron emission. Integrating over the excess and the fitted secondary SED (assuming a Rayleigh–Jeans approximation to the flux beyond 12 μm) results in L_IR/L_⋆ = 0.3. The lack of any evidence for strong magnetic fields from the WD suggests that it is unlikely to be due to cyclotron emission (e.g., as in the magnetic CV EF Eri, see Hoard et al. 2007). The fact that the excess is brighter than the M companion at long wavelengths rules out a tertiary low-luminosity companion as the cause of the excess. We instead examine if a circumbinary, optically thin, dust disk can qualitatively match the observed excess. The large emitting surface area implied by the relatively large L_IR/L_⋆ requires a circumbinary, rather than a circumprimary or circumsecondary disk—the implied radius of such a disk is larger than the orbital separation of the binary and thus would not be stable around either component (Artymowicz & Lubow 1994). In order to model the SED of the two stellar components and a circumbinary dust disk, we use the procedure developed to model IR emission from circumbinary dust observed around the CV V592 Cas (Hoard et al. 2009). The model is designed to predict emission from several components in a CV, including an optically thick accretion disk and an optically thin circumbinary disk (Hoard et al. 2007; Brinkworth et al. 2007). In the case of SDSS J0303+0054, there is no accretion disk and so the only constituents are the WD, M companion, and the dust.

The adopted fundamental parameters of the WD and M companion that we use in our model are given in Table 3. These values are consistent, to within our fitting uncertainties, to the mass and gravity we derived independently for the WD from the SED alone and are consistent with the constraints placed by P09 from spectral decomposition and light-curve analysis. This combination of parameters gives a statistically indistinguishable χ²_ν from our best-fitting model.

Using these assumptions, our model for the circumbinary disk is shown in Figure 3 compared only to our roughly contemporaneous LIRIS and WISE data. Because of the small number of photometric points, we do not constrain a formal best fit—given the significant degeneracies present in SED modeling, we seek only to obtain a qualitative match to the data to determine whether the circumbinary disk hypothesis is plausible. In general, our simple model reproduces the observed excess, starting at H and extending out to 12 μm. The model underpredicts the W1 point, which could be due to structure in the disk that deviates from the basic structure we assumed. The inner edge of a circumbinary disk is defined by the innermost stable orbit before tidal forces from the binary dominate a dust grain's orbit and ejects or accretes the dust grain (Artymowicz & Lubow 1994). Therefore, we assume a tidal truncation radius of 218 R_WD (=1.96 R_☉, 1.7 times the binary separation as described in Hoard et al. 2009), with an inner radius temperature of 1400 K, and an outer temperature of 50 K (at 18,500 R_WD = 0.77 AU) and a disk coplanar with the binary system. The inferred mass of the dust is 3 × 10²⁰ g assuming 1 μm silicate grains. If micron-sized silicate grains were present, however, one might also expect the W3 band photometry to be brighter due to a strong silicate feature at ∼10 μm. This could be accounted for by assuming a smaller outer radius to the disk, but at this time a specific composition for the dust is not constrained by the data, but is used to infer a rough order of magnitude on the mass of the dust. Furthermore, 10 μm silicate features are not observed in dust around CVs (Brinkworth et al. 2007), implying that the dust is larger than 1 μm. Assuming the dust all originated from a single body composed solely of silicates, the radius would be 21 km.

Zoom In Zoom Out Reset image size

Aperture photometry of all science targets and standard stars in our LIRIS observations was performed using an r = 375 aperture radius and sky annuli of 5''–75 in size, including extinction corrections. The derived JHK zero points for both nights agree to within 0.01 mag, but the scatter in the 11 standard stars each night was 0.05 mag. All photometry was calibrated using the J-, H-, and K-band system of the standard star data. Absolute fluxes for SDSS J0303+0054 were derived using the derived zero points for each night, and also by using 2MASS photometry for five relatively bright comparison stars in the field. The photometry thus derived for SDSS J0303+0054 using both methods agrees very well, within 0.02 mag.

The comparison stars were also used to search for variability in time-series relative photometry in the K band. The flux ratio of SDSS J0303+0054 relative to the comparison stars was calculated for all photometry resulting from each imaging sequence, with a variance-weighted average of each set of resulting magnitudes for SDSS J0303+0054 to create a final light curve of the system. The Heliocentric Julian Date was calculated from the midpoint of each photometric observation and converted to a phase using the ephemeris published in P09. Based on the uncertainty in the period, the phases we calculate should be accurate to within <10⁻² or about 83 s.

WISE observed SDSS J0303+0054 over the course of 1 day on 2010 January 29 and another day on 2010 August 6 for a total of 23 visits. Each visit had a detection of the system, with S/N of 30 and 10 in the W1 and W2 bands, respectively. We phased the data in the same manner as our K-band data.

Figure 4 shows the resulting K-, W1-, and W2-band light curves for SDSS J0303+0054. In both plots, Δm is negative when the system is brighter and positive when it is dimmer. We find that the light curve in the W1 and W2 photometry shows tentative evidence of variability, especially in a minor linear trend at phase ∼0.2, but the scatter in the data is large compared to this possible trend. The K-band light curve is clearly variable: the system varies at the level of 8% peak-to-trough. Based on our models, we expect the M dwarf to contribute close to 82% of the total flux at H and K and 50% (41%) of the flux at W1 (W2).

Zoom In Zoom Out Reset image size

Parsons et al. (2010) found similar variations in i-band light curves with an amplitude of 10% peak-to-trough and interpreted them to be due to ellipsoidal variations due to tidal deformation of the M dwarf by the WD. In i band, the M dwarf should contribute 80% of the flux in the system, comparable to H and K. If both H and K vary by ≈10% during an orbital period, this could account for some of the variation observed in the near-IR photometry. The phase of our H-band observations overlaps with the minimum we observe at K, potentially explaining the dimmer value compared to 2MASS and UKIDSS. We have also overplotted the calculated phases for the 2MASS and UKIDSS K photometric data. The 2MASS K_s photometry falls roughly at a phase in the orbit where the M dwarf is expected to be close to an inflection point of the light curve and close to its median brightness. If we extrapolate the behavior of the LIRIS light curve to the orbital phase of the UKIDSS photometry, we would expect the UKIDSS K magnitude to be ≈13.31, almost 0.2 mag brighter than what was observed in 2005. If the UKIDSS photometry is truly dimmer, this would imply variability in the amount of dust in the system (or of the M dwarf) on a timescale much longer than the orbital period. Observations of a full orbital phase in the K band, as well as long-term monitoring of this system, are necessary to confirm these unresolved issues.

The variability at K is in phase with the ellipsoidal variations observed at i and consistent with what one would expect based on the tidal distortion of the M star from the WD. Based on Roche geometry calculations neglecting limb darkening and using the mass and radius values we adopted in our circumbinary dust model, we would expect ellipsoidal variations peak-to-trough of 12% (which scales to 10% if the secondary contributes 84% of the flux at K). The shape of the light curve in Figure 4 deviates slightly from a pure sine function; this is also seen at shorter wavelengths (P09; Parsons et al. 2010). We further discuss our results in the context of all the light curves for this system in Section 6.

The WISE observations of SDSS J0303+0054 represent the first likely detection of dust in orbit around a non-interacting post-common envelope binary. If the dust disk is aligned with the orbit of the binary, SDSS J0303+0054 could be a fruitful target for high-cadence near-IR transit searches. Our photometry was sensitive to transits of depths >50 mmag and durations >3 minutes. With the same photometric accuracy and sampling at <30 s, planets with radii of ∼5 R_⊕ could be detected if they orbited at the inner edge of the disk and transited with durations of ∼1 minute.

The presence of ellipsoidal variations deviating from pure sinusoidal behavior, now seen at K with our new observations, also opens up the possibility that the structure of the circumbinary disk can be directly probed by detailed analysis of long-wavelength time-series photometry. The circumbinary disk is not symmetrically heated throughout an orbital period, and if there exist significant azimuthal density or compositional variations, these could be detected as asymmetries in the light curve. Since the disk contributes 15% of the light at K, such variations could be detected at the few percent level.

One possibility is that SDSS J0303+0054 is the binary equivalent of isolated WD dusty disks, which most likely form from the tidal disruption of rocky bodies a few tens of kilometers in radius. Given the relatively long cooling age of this system compared to other dusty WDs and the inferred T_eff of the WD, this would be comparable to the coldest recorded dusty single WD observed, G166-58 (Farihi et al. 2010). The amount and location of dust present at SDSS J0303+0054 would not be detectable around a single cool WD. The presence of the M companion's additional luminosity, which far exceeds that of the WD, enhances the detectability of the circumbinary dust.

If the dust is caused by a tidally disrupted asteroid, then SDSS J0303+0054 is reminiscent of Kepler-16b, which contains a 0.6 M_☉ primary, a 0.2 M_☉ companion, and a Saturn-mass companion in a circumbinary orbit (Doyle et al. 2011). Such a binary, with a 41 day orbit, would eventually evolve into a CE binary. The remains of the circumbinary planetary system could feed a similar dust disk to that seen around SDSS J0303+0054. However, several potential problems exist for such a system. First, the Poynting–Robertson (P-R) drag timescale might be prohibitively short for such a disk. The P-R drag timescale can be approximated by (Hansen et al. 2006)

$$\begin{eqnarray} T_{\rm P\hbox{-}R} &=& 100 \left(\frac{s}{1\,\mbox{$\mu $m}}\right)\left(\frac{\rho _s}{3\,\rm {g\,cm}^{-3}}\right)\left(\frac{r}{10^{11}\,\rm {cm}}\right)^2 \nonumber\\ &&\times \left(\frac{L_\star }{10^{-3}\,L_\odot }\right)^{-1} \rm {yr}, \end{eqnarray} \tag{ 1 }$$

where L_⋆ is the luminosity of the star, s and ρ_s are the average grain size and density respectively, and r is the distance from the central star. As noted previously, the M dwarf dominates the luminosity of the system. Assuming a bolometric luminosity of 5 × 10⁻³ L_☉ (Rebassa-Mansergas et al. 2007), and assuming a minimum distance of the disk from the M dwarf of 0.76 R_☉,  the P-R timescale can have a range of values since the dust orbits the binary's center of mass. At the closest distance to the M dwarf, the P-R timescale is 6 yr, while at the furthest distance away from the dust disk inner edge the P-R timescale can be as long as 100 yr. The relatively short P-R timescale implies that dust must be replenished at a rate of 5 × 10¹⁹ g yr⁻¹ or 2 × 10 ¹² g s⁻¹. The minimum replenishment rate, corresponding to the longer P-R timescale, would be 16 times smaller. This rate of dust production is high compared to the inferred accretion rates of dust onto single WDs, but could be reconciled if the radius of a typical dust particle is larger than 1 μm.

The mechanism for creating the dust in the asteroidal scenario might be significantly different than for single WDs. It is unclear how an asteroid could tidally disrupt at the circumbinary disk inner radius, which is much larger than the Roche disruption radii of either the WD or the M dwarf. If the dust is from tidally disrupted bodies, multiple asteroids could have similar periastra and the dust forms from multiple collisions of these objects. The accretion rate would be lower if the disk were optically thick, since P-R drag could not act on all of the dust present. Assuming an inclination of 82°, an inclination consistent with the results of P09, the lack of reddening toward the system constrains the vertical height of the disk to be <0.1 R_☉.

Alternatively, low levels of mass transfer that might be present in this system could account for the presence of a dust disk as well. The combined photometric and spectroscopic behavior of the SDSS J0303+0054 system suggests low levels of wind material being accreted by the WD. The u and g light curves of the system show a brightening of the system coincident with the WD's eclipse. Parsons et al. (2010) interpreted this excess flux as due to the side of the M dwarf facing away from the WD being brighter. However, this is contradicted by the ∼96% depth of eclipse seen in the u band, despite a 10%–20% increase in flux leading up to the eclipse. If the M companion were responsible for the excess, the eclipse would be shallower than that predicted by the ratio of the M dwarf optical flux to the total flux from the system. Based on our best-fit model, the eclipse depth should be ∼90% in u. This difference suggests our model overpredicts the flux from the M dwarf, or underpredicts the flux from the WD. We instead interpret the brightening, which covers about half the orbit, as due to the side of the WD facing the secondary being hotter due to broad wind accretion on that hemisphere. The presence of multi-component Hα emission lines (Tappert et al. 2011a) is also indicative of accretion from the secondary wind.

Similar to CVs with infrared excesses, this wind accretion may provide enough "spillage" to create a low-mass dust disk. Using the above P-R drag mass-loss rate as an estimate of dust destruction, we can infer the rough mass injection required to obtain a steady-state dusty disk. If we assume that all of the metals in the secondary's wind condense into dust, and that half of the wind must contribute to the "spillage," this suggests that the M dwarf must be losing mass at a rate of 5 × 10²¹ g yr⁻¹ with a minimum rate an order of magnitude lower for the longer P-R drag and assuming roughly solar metallicity. This could be uncertain by a factor of 10 if the disk is radially optically thick—this can slow the destruction of dust grains if some of the dust is shadowed. If we compare this to the mass-loss rate of other M dwarf winds that accrete onto WDs (Debes 2006), we find that the M dwarf in SDSS J0303+0054 must be losing mass at a maximum rate ∼100–1000 times that of typical M dwarfs and more than the solar wind mass-loss rate (Debes 2006; Tappert et al. 2011b; Pyrzas et al. 2012). If the M dwarf is tidally distorted, as required by the presence of ellipsoidal variations, then this could enhance the mass loss in its wind. The escape velocity at the secondary's surface is smaller than the escape speed of the binary system, which could naturally trap a majority of the wind particles in stable orbits around the binary.

Neither scenario completely explains the origins of this unusual system. If there is significant dust production from low levels of mass transfer or from remnant planetary systems, it is possible that a significant number of other post-common envelope binaries contain dust as well. With over 1000 candidate WD+M systems from the WIRED survey, a more careful analysis of that population may yield other dusty post-common envelope binaries, and hopefully lead to a resolution of the origin of the dust, and how it fits into the evolution of compact binaries.

This work is based on data obtained from: (1) the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory (JPL), California Institute of Technology (Caltech), funded by the National Aeronautics and Space Administration (NASA); (2) the Two Micron All Sky Survey (2MASS), a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center (IPAC)/Caltech, funded by NASA and the National Science Foundation; (3) the SIMBAD database, operated at CDS, Strasbourg, France; and (4) the NASA/IPAC Infrared Science Archive, which is operated by JPL, Caltech, under a contract with NASA. M.C. thanks NASA for supporting his participation in this work through UCLA Sub-Award 1000-S-MA756 with a UCLA FAU 26311 to MIRA.