A TREND BETWEEN COLD DEBRIS DISK TEMPERATURE AND STELLAR TYPE: IMPLICATIONS FOR THE FORMATION AND EVOLUTION OF WIDE-ORBIT PLANETS

We searched the Spitzer observers log for main-sequence stars that were observed with the IRS Long Low (LL) module (both orders) in staring mode and with MIPS at 24 μm and 70 μm, and we accumulated a sample of 546 targets. The stellar properties of our target list are summarized in Table 1. After reducing and analyzing the data, we refined the sample as described in Section 2.5. This yielded 225 stars with significant excess (of which 174 had cold components), which are detailed in Table 2.

It is important to note that our sample comprised stars from a variety of Spitzer observational programs, each having targets selected in a different manner (e.g., nearby stars, stars with previously detected IR excess, stars with RV-detected planets, etc.). Hence, our sample was not selected to be statistically representative in any one sense; rather, it was designed to include as many relevant targets as possible.

We calculated the stellar temperature, T_⋆, from the star's V − K_s color using tabulated values for main-sequence stars (or interpolations between these values) from Cox (2000), and we estimated a 200 K 1σ uncertainty on these values. T_⋆ is listed for targets with significant excess in Table 2. Over the range of spectral types and metallicities of our sample of stars, V − K_s is an accurate indicator of stellar effective temperature, with a peak-to-peak scatter of less than ±1.5% and weak metallicity dependence (Masana et al. 2006). We use T_⋆ to parameterize stellar type, and this degree of accuracy is sufficient for our purpose. Some of the uncertainty in T_⋆ reflects the effect of interstellar reddening on V − K_s, although the great majority of stars in this sample are within the Local Bubble and therefore should not be strongly reddened (of the 174 targets for which we detect cold components, only 29 are beyond 100 pc, 19 are beyond 120 pc, and 4 are beyond 150 pc).

Ages for these stars were estimated from a combination of chromospheric activity measurements, X-ray emission, placement on the Hertzsprung–Russell (HR) diagram, surface gravity, membership in clusters and associations, and gyrochronology collected from the literature. Sierchio et al. (2013) describe the intercomparison of these methods and how they are applied, and the quoted ages are on the scale calibrated by Mamajek & Hillenbrand (2008). The ages are listed in Table 1, along with references for the measurements used to derive the ages and a quality flag for the age accuracy (ranging from 0 if no age could be determined to 3 if there were three or more age measurements with good agreement).

Although we built our target list around the Spitzer IRS spectra, we supported these spectra in our analysis with a suite of photometric data. The properties of the photometric systems we used are summarized in Table 3 and the photometric data for the targets with significant excess (see Section 2.5) are given in Table 2.

MIPS photometry at 24 μm provided an important calibration reference for the IRS data (Section 2.3), while 70 μm photometry provided a crucial constraint on the temperature of cold debris disks (Section 2.5). We used our in-house debris disk pipeline to reduce and extract photometry for the MIPS data as part of the effort to preserve the legacy of Spitzer measurements on debris disk studies (Su et al. 2010; K. Su et al. 2013, in preparation). Basic reduction (up to the post-basic calibrated data (post-BCD) mosaics) and calibrations of the MIPS data follow the descriptions by Engelbracht et al. (2007) and Gordon et al. (2007). In addition, the extraction of the 24 μm photometry was briefly described in Urban et al. (2012), where the source position at 24 μm is determined by a combination of point-spread function (PSF) fitting and two-dimensional (2D) Gaussian fitting. Both PSF fitting and aperture photometry³ were performed, and we preferentially used the aperture photometry at 24 μm because this was used to calibrate the photosphere model at 24 μm, as described in Section 2.4 and the Appendix (the PSF and aperture photometry agree to within a few percent). We then used the 24 μm source position to perform PSF fitting for data at 70 μm by minimizing the residual signal at the source position (for details, see K. Su et al. 2013, in preparation). For faint sources located in areas with structured background, the resultant 70 μm photometry can be negative, which reflects non-detection. The 1σ photometry uncertainty (listed in Table 2) includes the pixel-to-pixel variation near the source of interest, and detector repeatability (1% and 5% of the source flux at 24 and 70 μm, respectively). MIPS photometry data for all of our targets can be found in Gáspár et al. (2013), Sierchio et al. (2013), and K. Su et al. (2013, in preparation).

V, J, H, and K_s photometry were used to model the stellar photosphere SED (Section 2.4) and to estimate T_⋆ (Section 2.1). We obtained Hipparcos V-band (ESA 1997) and Two Micron All Sky Survey (2MASS) J-, H-, and K-band photometry (Cutri et al. 2003) from the VizieR online database. Many of the targets in the sample are nearby stars and are severely saturated in the 2MASS data. To overcome this, we used heritage aperture photometry, which we transformed to match the 2MASS system. When both 2MASS and transformed heritage photometry of high quality were available, we averaged them. The references for the heritage photometry are in Table 2.

IRAC photometry (3.6 μm) was also used to model the stellar photosphere SED (Section 2.4). These data were taken in Spitzer cycle 7 (PID 70076, PI: Su). All the IRAC data were taken in subarray mode with four dithered positions to avoid saturation. We used the BCD products provided by the Spitzer Science Center (pipeline version S18.18), and performed the necessary steps (pixel solid angle correction and pixel phase correction; K. Su et al. 2013, in preparation) to extract the photometry. Aperture photometry was used for each individual data frame (64 frames per dithered position); and the final quoted flux is the median value of all measurements per star.

Wide-field Infrared Survey Explorer (WISE) data (Cutri et al. 2012) were obtained from VizieR, but are not included in Table 2. Because WISE photometry is less accurate than Spitzer photometry, we did not use it quantitatively. Instead, we used it as a qualitative confirmation of our Spitzer data.

Our IRS reduction started with the Level 1 BCD products, downloaded from the Spitzer Heritage Archive. In addition to IRS LL data, we also reduced Short Low (SL) module data, when available, and included them in our analysis. The Astronomical Observation Requests (AORs) used for each target are listed in Table 1. The basic reduction was performed using the Spectroscopic Modeling Analysis and Reduction Tool (SMART) software package (Higdon et al. 2004), scripted into a series of automated routines. For each spectral order of each IRS AOR, three files (2D spectra, uncertainty, and mask) were combined into a single three-plane file. Bad and rogue pixels were removed by the routine IRSCLEAN with the clean parameter set to 4096 (pixels with this value or higher were included in the clean). Next, when available, multiple Data Collection Events for the same nod position were combined, and then the background was removed from each 2D spectrum by subtraction of the opposite nod. The 2D spectra were then converted into one-dimensional (1D) spectra using SMART's optimal 2 nod extraction (Lebouteiller et al. 2010). The 1D spectra from both nods were combined. The result was a wavelength, flux, and uncertainty vector of each spectral order for each AOR. The "bonus" third order data were not used. The remainder of the data reduction and analysis was performed with the MATLAB software package.

IRS data at wavelengths near the order edges often were of poor quality, so data longward of 38 μm and shortward of 20.75 μm were discarded from the LL first order, data shortward of 14.3 μm were discarded from the LL second order, data longward of 14.7 μm and shortward of 7.55 μm were discarded from the SL first order, and data shortward of 5.25 μm were discarded from the SL second order. Systematic offsets in flux existed between IRS orders, which we fixed by applying a multiplicative correction factor to the LL first order spectrum and to both SL order spectra (if they existed), in order to align them with the LL second order flux. Initially, we determined the order correction factors using an automated routine, but due to the variety of shapes of the spectra and the presence of outlying data points, we found that fine-tuning these factors by eye was more reliable. The data from all available orders and modules were then combined into a single spectrum.

Next, we cut outlying data points from the spectrum in an iterative process. The spectrum was fit with a polynomial and the standard deviation of residual flux values around the fit was calculated. Points lying more than three standard deviations from the fit were discarded. This process was iterated six times. Because the scatter of the data generally increased toward longer wavelengths, we applied this process separately to the short and long ends of our spectrum; if SL data were available, the two sections were divided at 14 μm and each section fit with a fourth degree polynomial, whereas if no SL data were available, the sections were divided at 25 μm and each section was fit with a second degree polynomial. This procedure, on average, cut 10 data points from each spectrum, with the first iteration typically cutting five points, and the sixth iteration only cutting zero, one, or two points (nearly 90% of spectra had no points cut in the sixth iteration). Before cutting, the spectra had 181, 296, or 373 points, depending on the available orders. After the outliers were cut, the spectrum was smoothed by binning to a wavelength resolution of 0.7 μm.

The absolute calibration of IRS data is only accurate to ∼10%, whereas the MIPS calibration is accurate to ∼2% (Engelbracht et al. 2007). To reduce systematic errors in our IRS data, we multiplicatively scaled our spectrum to be consistent with the measured MIPS 24 μm flux density for each target. The IRS spectra were calibrated as point sources, so if a source was slightly extended, it would have an incorrect slit-loss correction. The stellar photospheres (point-like) are dominant in the IRS spectra for the majority of sources in the sample; therefore, this has no significant impact on the result. If the MIPS photometry for a source was contaminated by background emission, this contamination would be passed to the IRS spectrum of the source through this scaling (the IRS data may or may not have picked up the contamination, depending on the slit orientation). Nevertheless, an additional multiplicative factor was later applied to the IRS data (see Section 2.4), which corrected any lingering errors from this process.

Some targets were observed with more than one IRS AOR. We combined the spectra from these AORs, interleaving their data points, then re-smoothing to a spectral resolution of 0.7 μm.

Although the spectrum of a main-sequence star peaks in the visible or near-IR, the stellar photosphere can still be the dominant source of flux density in the mid-IR. Thus, to extract the thermal emission of the debris disk (the IR excess), an accurate model spectrum of the stellar photosphere must be generated and subtracted from the data.

It is common to use photosphere spectra from detailed numerical models of stellar structures, such as KURUCZ/ATLAS9 (Castelli & Kurucz 2004). However, these models are not well tested in the mid-IR and far-IR. Sinclair et al. (2010) compared various families of models (KURUCZ/ATLAS9, MARCS, and NEXTGEN/PHOENIX) and found that the choice of model family can affect whether or not a star is determined to have an IR excess.

We opted for a simpler model, a Rayleigh–Jeans relation given by F_ν ⋆(λ) = RJ/λ², where RJ is an amplitude scale factor. This model is appropriate for wavelengths beyond 5 μm because stars in this study, mostly A through K dwarfs, have virtually no spectral features in the mid-IR, behaving as blackbodies in the Rayleigh–Jeans regime. For example, in the range of most interest for this study, the departure of a reference A0 star spectrum (Rieke et al. 2008) from a Rayleigh–Jeans fit between 15 and 30 μm is no more than ±0.3%, and the ATLAS9 (Castelli 2013) solar model follows Rayleigh–Jeans behavior to within ±0.24% from 20 to 40 μm. However, the theoretical spectra show minor systematic differences from the observations (Rieke et al. 2008), so use of them to improve the knowledge of the SEDs would be risky.

Modeling the photosphere thus came down to finding the appropriate scale factor (RJ) for each star in our sample. We used a set of color relations (derived from a sample of stars known to have no IR excess) to predict the MIPS 24 μm magnitude of the photosphere, [24], from measured V, J, H, K_s, and IRAC magnitudes—bands where IR excesses due to very close-in dust are very rare. The details of these relations are described in the Appendix. This procedure estimates the 24 μm photospheric outputs to about 2% rms, based solely on data at wavelengths short of 4 μm. We then converted [24] to flux density units and used it to find the scaling factor for the photosphere model, which is given in Table 2 for targets with significant excess (see Section 2.5).

We finally applied a small multiplicative factor to the IRS data if it clearly looked slightly offset from the photosphere model. Then, we subtracted the photosphere model from the IRS spectra and MIPS 70 μm point to obtain the IR excess, F_ν excess(λ) = F_ν(λ) − F_ν ⋆(λ). The photosphere model was assumed to have an uncertainty of 2%, thus the uncertainty in the excess flux density was $\sigma _{\rm excess}\left(\lambda \right) = \sqrt{\sigma \left(\lambda \right)^2 + (0.02F_{\nu \,\star }(\lambda))^2}$.

To interpret the excesses in a general way, we needed to assign them fiducial temperatures. The emission of disk grains is a function of their optical constants; however, as discussed in Section 1, debris disks radiate like blackbodies at one or two temperatures. The simplest possible description of the emission is in terms of these temperatures. In this section, we describe how we determined if each target had significant excess and whether the excess was best described by one or two blackbodies.

We fit F_ν excess(λ) with a one-component blackbody model:

$$\begin{equation} F_{\nu \,{\rm model1}}(\lambda)=c_{\rm one}F_{\nu \,{\rm BB}}(\lambda ,T_{\rm one}), \end{equation} \tag{ 1 }$$

where F_ν BB(λ, T) is the Planck function. The best-fit parameters, c_one and T_one, were found by minimizing the reduced χ²:

$$\begin{eqnarray} \chi _{\nu }^2 &=& \frac{1}{\nu } \left\lbrace \sum _i^N \frac{(F_{\nu \,{\rm excess}}\left(\lambda _i\right)-F_{\nu \,{\rm model}}\left(\lambda _i\right))^2}{\sigma _{\rm excess}\left(\lambda _i\right)^2}\right. \nonumber\\ &&\left. + 28\frac{(F_{\nu \,{\rm excess}}\left(70\, {\mu m}\right)-F_{\nu \,{\rm model}}\left(70\, {\mu m}\right))^2}{\sigma _{\rm excess}\left(70\, {\mu m}\right)^2} \right\rbrace , \end{eqnarray} \tag{ 2 }$$

where ν = N + 28 − n − 1, N is the number of data points in the IRS excess data, and n is the number of free parameters in the fit (n = 2). The MIPS 70 μm data point was weighted in the fit as 28 IRS data points, equivalent to the number of IRS data points (0.7 μm resolution) that would fit inside the equivalent width of the MIPS 70 μm spectral response function (19.65 μm). The fit was performed using MATLAB's lsqcurvefit algorithm. T_one was constrained to between 0 and 500 K.

While c_one represents the amplitude of the debris disk emission, a more useful measure of a debris disk's brightness is the fractional excess, f ≡ L_excess/L_⋆. We calculated f according to

$$\begin{equation} f=\left(\frac{F_{\nu \,{\rm excess\ max}}}{F_{\nu \,\star \,{\rm max}}}\right) \left(\frac{\lambda _{\star \,{\rm max}}}{\lambda _{\rm excess\ max}}\right), \end{equation} \tag{ 3 }$$

from Equation (2) of Wyatt (2008). F_{ν excess max} and F_{ν ⋆ max} are the peak flux density values of the disk emission and stellar photosphere emission, respectively, which occur at wavelengths λ_excess max and λ_⋆ max. While the peak flux densities for the disk components were easily calculated from our best-fit model, the Rayleigh–Jeans stellar photosphere model had no maximum. To overcome this, we created a blackbody function using the temperature of the star from Table 1, and scaled it to match the flux density of our Rayleigh–Jeans model at 24 μm. From this representation of the photosphere, we found the maximum flux density, and calculated f.

With our best fits in hand, we refined our sample to only those targets with statistically significant excess. First, we inspected all fits by eye and discarded targets with poor quality data that resulted in clearly unrealistic fits. Second, we discarded targets with excess that was too faint, f < 10⁻⁵. Third, we inspected the fits and identified all targets whose excess relied solely on the MIPS 70 μm data point (i.e., there was no excess in the IRS data). For these cases, we required this MIPS point to represent a significant excess, so we discarded targets where F_ν excess(70 μm)/σ_excess(70 μm) < 3 or where F_ν excess(70 μm)/F_ν ⋆(70 μm) < 1. This process resulted in 321 of the original 546 stars having no significant excess.

Of the remaining 225 targets with significant excess, we next determined whether the excess consisted of 1 or 2 components. To do this, we fit the excess SED of each target with a model consisting of the sum of two blackbodies:

$$\begin{equation} F_{\nu \,{\rm model2}}(\lambda)=c_{\rm cold}F_{\nu \,{\rm BB}}(\lambda ,T_{\rm cold}) + c_{\rm warm}F_{\nu \,{\rm BB}}(\lambda ,T_{\rm warm}). \end{equation} \tag{ 4 }$$

Finding the optimal set of the four parameters c_cold, T_cold, c_warm, and T_warm was again done by minimizing the reduced chi-squared (Equation (2)), now with n = 4.

Our definition of "cold" was the coldest well-detected component of the excess that was below 130 K. Morales et al. (2011) found a continuous distribution of warm components above ∼130 K, centered around 190 K, whose temperatures were set by the water ice line. Thus, any component above 130 K was considered "warm" for the purposes of this study. A component with a temperature of 110 K, for example, would be considered warm if another colder temperature component was also detected, but would be considered cold if no colder component was detected. To implement this, we fit each target once using 100 K as the division between warm and cold, and again with the division at 130 K. We then selected the better of these two cases (based on reduced χ²) to represent the best two-component model. For many targets, these two cases yielded identical fits. In both cases, the minimum cold temperature allowed by the fit was 0 K, and the maximum allowed warm temperature was 500 K.

Deciding if each target warranted a two-component model was a two-step process. First, we required that the reduced χ² of the two-component fit be at least three times better than that of the one-component fit (i.e., all targets where $\chi _{\nu 1}^2$/$\chi _{\nu 2}^2 <3$ were deemed to be better fit by a single component). Second, for the remaining targets, we identified those for which the presence of the cold component relied entirely on the MIPS 70 μm point, meaning the IRS excess was entirely fit by the warm component. For these targets to have significant cold components, we required their MIPS 70 μm excess to fall more than 3σ above the excess predicted by the warm component. That is, if one of these targets had [F_ν excess(70 μm)  −  c_warmF_ν BB(70 μm, T_warm)]/σ_excess(70 μm) < 3, we concluded that it was best fit by one component. This step was analogous to the third step we performed when deciding if each target had excess or not. The targets that passed both of these criteria were deemed to have two components.

This process identified 100 single-component cold disks, 51 single-component warm disks, and 74 two-component disks. We found only one target (HIP45585) that required two warm components (both >130 K) to fit properly, and we discarded this target from our results as it had no cold components (it is counted in the 321 targets that we discarded). We found no targets requiring two components colder than 100 K. Our methods were insensitive to possible additional cold components below ∼50 K. We concluded that our fitting procedure allowed us to find all cold components above this limit. Equivalently, our method reliably fit the inner edges of the cold disks.

For the two-component disks, we calculated the fractional excesses for the warm and cold components, f_warm and f_cold, using Equation (3). One-component excesses were deemed warm or cold depending if they had temperatures above or below 130 K. The verdicts for all targets are given in Table 1 and the parameters of our best fits for targets with significantly detected components are listed in Table 2.

We estimated that the uncertainty in T_cold was 10 K. This estimate was conservative; warmer cold components had more IRS data points in excess of the photosphere, so their fits were more constrained, with a temperature uncertainty of 3–7 K.

We were unable to accurately constrain the temperatures of very cold components, which were detected only at 70 μm. With only one data point, blackbodies at a wide range of cold temperatures could be given the amplitude necessary to match the point, so a degeneracy existed between T and f in blackbody fits to the excess SED. We made this distinction at 45 K; above this temperature there was enough information in the IRS spectra to constrain the temperature. Cold temperatures above 45 K were considered true detections, while cold temperatures below 45 K were considered unconstrained and were replaced with upper limits at 45 K. By using the 70 μm data to determine if such targets had significant excesses, we ensured that these disks were truly very cold, rather than warm but very faint. Of the 174 cold components in our sample, 25 had temperature upper limits.

Some IRS spectra exhibit mineralogical emission features (e.g., HIP41081, HIP57971). Our model did not explicitly account for these features, but they peak sharply at ∼10 μm and do not resemble blackbody functions. These features did not influence our fits, except to raise the value of the minimum reduced χ².

The data and best fits for a small sample of our targets are shown in Figures 1–3 (for targets with two components, one cold component, and one warm component, respectively). The left panels in all of these figures show the measured data with the photosphere and blackbody models, while the right panels illustrate the excess data and models with the photosphere subtracted. Histograms of the significant warm and cold debris disk temperatures are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size