DEBRIS DISTRIBUTION IN HD 95086—A YOUNG ANALOG OF HR 8799

Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004) photometry was first published by Chen et al. (2012) where fluxes were extracted using point spread function (PSF) fitting designed for point sources. We examined the 24 μm image, and found that the source is slightly extended (FWHM of 598 × 587, compared to a typical point source's FWHM of 549 ± 003 × 545 ± 003). Therefore, we adopted measurements with photometry through a large aperture (aperture radius of 149 and sky annulus of 299–423), resulting in 52.95 ± 0.55 mJy as the total 24 μm flux for the system (∼16% higher than the PSF extracted flux). At 70 μm, the source is point-like with a total flux of 655 ± 33 mJy. The quoted flux uncertainties include instrumental repeatability, which is 1% and 5% at 24 and 70 μm, respectively.

Herschel Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRES; Griffin et al. 2010) observations were first published by Moór et al. (2013), and showed that the disk was resolved at 70 and 100 μm with an estimated inclination angle of 25°. Because the Herschel observations were reduced with an older pipeline (HIPE v9.2), we re-reduced the data with an updated pipeline (HIPE v13.0) that incorporated many new calibration improvements detailed in Balog et al. (2014). Specifically, we included the experimental tool for pointing reconstruction based on the gyroscope information that considerably reduces the pointing jitter and increases the stability of the PSF (termed as gyro correction). Finally, we mosaicked the final maps with a default "pixfrac" parameter 1.0, equivalent to drizzle parameter of 0.8, in the instrumental orientation so the PSF structures were always in the same location. No rotation/interpolation is then needed in further PSF subtraction analyses. Detailed analysis of the PACS images is presented in Section 3.1. We used an aperture size of 22'' at all three bands to measure the source flux, at which radius the encircled flux reached the maximum and flattened afterward (after proper aperture correction). We did not apply color correction for the PACS fluxes. The final PACS fluxes are 688 ± 48, 688 ± 48, and 489 ± 35 mJy at 70, 100, and 160 μm, respectively, in good agreement with the Moór et al. (2013) published results.

Inspecting the SPIRE maps, the source is clearly detected at all three wavelengths as a point source, and the surrounding field has no bright source or large-scale extended structure. Therefore, we used the level 2 point-source product (Jy beam⁻¹) provided by the Herschel Science Center (HIPE ver. 11) to measure fluxes. We used the peak value of the source to estimate the point-source flux. After correcting the pixelization effect (0.951, 0.931, 0.902 at 250, 350, and 500 μm, respectively; SPIRE Handbook, version 2.5, 2014), the source submillimeter fluxes are 186 ± 16 mJy, 112 ± 16 mJy, and 67 ± 17 mJy at 250, 350, and 500 μm, respectively. At 500 μm, there are a couple of faint peaks within 40'' of the target, which might contaminate the submillimeter fluxes. No color correction was applied for the SPIRE fluxes either. Our SPIRE fluxes are lower than the values published by Moór et al. (2013) who extracted source fluxes using aperture photometry, but consistent within uncertainties. Finally, HD 95086 was detected at 870 μm with the submillimeter telescope APEX by Nilsson et al. (2010). Its nominal beam size is 192 (FWHM), but the effective resolution is 27'' after filtering and smoothing in the data reduction (Nilsson et al. 2010). The source at 870 μm appears to be elongated (Figure 2 in Nilsson et al. 2010) and surrounded by faint peaks. Therefore, it is likely that the measured 870 μm flux (41.3 ± 18.4 mJy) is overestimated and contaminated by background galaxies.

A Spitzer InfraRed Spectograph (IRS; Houck et al. 2004) spectrum of the system was first published by Chen et al. (2006) where the excess spectrum was characterized by blackbody emission of 80 ± 30 K. The same spectrum was reduced and analyzed by Moór et al. (2013), who found the spectral shape of the excess along with long-wavelength photometry is better fitted with two blackbody temperatures: 187 ± 26 K and 57 ± 1.5 K. Finally, Chen et al. (2014) characterized the excess to be consistent with two temperatures: 225$^{+10}_{-7}$ K and 57 ± 5 K. The warm temperature discrepancy might be due to how the two IRS modules were joined and scaled relative to the photosphere. For this reason, we also performed the IRS data reduction, similar to what was done in Moór et al. (2013) but with the SMART software (ver. 8.2.7; Higdon et al. 2004). The two low-resolution modules (short-low (SL): 5.2–14.5 μm; long-low (LL): 14.0–38.0 μm) were extracted and combined independently using SMART's optimal 2 nod extraction (Lebouteiller et al. 2010). The SL spectrum agrees very well (<1%) with the expected Kurucz model flux for wavelengths shorter than 6.5 μm, while the LL spectrum was scaled down (by 0.963) to match the MIPS 24 μm large aperture photometry. The two modules then joined smoothly without further scaling.

The photospheric contribution was subtracted using the best-fit Kurucz model (T_eff = 7500 K and log g = 4.0), and 2% of the photospheric emission was included in computing the excess flux uncertainty. Figure 1 shows the excess SED composed of broad-band photometry and the IRS spectrum. Our best-fit (weighted by the uncertainty) two blackbody temperatures are 55 ± 5 K (cold component) and 175 ± 25 K (warm component). The large error in the warm temperature is driven by the large uncertainty in the 7–10 μm excess flux after photospheric subtraction. The 175 K blackbody emission slightly (within 1σ) under predicts the excess flux shortward of 10 μm, suggesting either the system has another faint, hotter component or the emission from the warm component has a weak silicate feature from ∼μm size grains. We have used the fitting technique developed by Ballering et al. (2014), including various normalizations of the photosphere, to evaluate this emission. We find no evidence for silicate emission features, thus differing from those excesses identified by Ballering et al. (2014); instead, within the errors, a ∼300 K blackbody is an adequate fit. Adding this component has negligible effect on the fitted temperature of the 175 K component (reducing its temperature by about 10 K). Given its large uncertainty, we do not include this tentative, faint, and hotter component in our following discussion.

Zoom In Zoom Out Reset image size

The MIPS-SED data were obtained in 2006 Feb 20 (AOR key 16171264, PID 84, PI: Jura) with the 10 s × 10 cycles and 1' chop setting, resulting in a low-resolution (R = 15–25) spectrum from 55 to 95 μm with 600 s of integration on source. The data were reduced and calibrated as described by Lu et al. (2008). An extraction aperture of 50'' in the spatial direction was used, and the slit loss was corrected assuming a point source. We trimmed off the spectrum past 95 μm because of known contamination by the blue end of the second-order spectrum produced in this observing mode of the instrument that cannot be calibrated (Lu et al. 2008). The photosphere-subtracted MIPS-SED spectrum is shown in Figure 2, where a feature peaking near ∼70 μm can be seen. There is also one low discrepant point near 80 μm (by 2–3σ) that is problematic. We investigated the discrepancy and were not able to find an obvious reason outside the possible transient systematic errors in the array around 80 μm. However, the feature at ∼70 μm appears to be real since it consists of two high points that are ≳3σ above the overall continuum, and several adjacent points that are consistent with the instrumental resolution, unlike the single low point at 80 μm. We have also checked other MIPS-SED spectra for sources that have similar 70 μm fluxes as HD 95086, and none were found to have such an emission feature, suggesting that it is not due to instrumental artifacts. The shape of the MIPS-SED spectrum is consistent with the broad-band photometry measured by PACS and MIPS at 70 μm. Since the bandwidth of the MIPS 70 μm filter is wider than that of the PACS 70 μm filter, the slightly lower value of the MIPS 70 μm measurement is consistent with the presence of the expected feature. A Lorentzian profile fit to the MIPS-SED spectrum suggests that the feature peaks at 69.5 ± 0.5 μm with a linewidth (FWHM) of 4.2 μm, consistent with an emission-line feature that is not spectrally resolved in the MIPS-SED mode. Given the cold dust temperature (55 ± 5 K), the measured peak position is consistent with the laboratory measurement of crystalline olivine (Mg₂SiO₄, forsterite) at 50 K (Suto et al. 2006). The 69 μm band strength (integrated luminosity) is 8.3 × 10⁻¹³ erg cm⁻² s⁻¹, which is in the range of detected band strengths around protoplanetary disks (see Table 2 in Maaskant et al. 2014). This feature has also been detected in the β Pic debris disk using the Herschel/PACS spectrometer (de Vries et al. 2012). We discuss the possible implications of the presence of the 69 μm feature in Section 4.2.

Zoom In Zoom Out Reset image size

Fully understanding the instrumental PSFs is essential in characterizing resolved images; in particular, it has been shown that the high-pass filtering and drizzling methods for coadding images used in HIPE have a great impact on the resultant PSF and noise characteristics (Popesso et al. 2012). We, therefore, carried out some basic PACS PSF characterization by using PACS data for γ Dra, a routine calibrator, that were observed throughout Herschel's operation. We used the data obtained in the mini-scan map mode with two scan angles of 70° and 110° only (the same as the HD 95086 data). These data were processed with the same pipeline and procedure as our target star (same filtering width, masking radius, pixfrac, instrumental orientation, and gyro correction). We further inspected the quality of the data by measuring the FWHM of each observation and rejected data sets that have deviant values. Finally, we obtained a high signal-to-noise, co-added γ Dra image by median combining sets of 26, 5, and 26 observations at 70, 100, and 160 μm, respectively. These combined, high-quality γ Dra data were used as our point source standards when assessing the resolvability of an observation.

All the FWHM measurements were measured by fitting a two-dimensional Gaussian function on sub-regions centered at the target with fields of view (f.o.v.) of 211, 294, and 441 at 70, 100, and 160 μm, respectively. Note that different values of f.o.v. would result in slightly different FWHM values; therefore, it is important to be consistent with the f.o.v. used in comparisons. The average FWHMs for the γ Dra data are 575 ± 008 × 556 ± 004 with position angle (P.A.) of 32° ± 29° at 70 μm, 685 ± 004 × 670 ± 004 with P.A. of 33° ± 26° at 100 μm, and 1197 ± 016 × 1069 ± 015 with P.A. of 91° ± 5° at 160 μm. The beam (defined as the FWHM) is more elongated in the 160 μm channel and with a preferred P.A. compared with the beams at 70 and 100 μm. The variations in the measured FWHMs are all less than 2%, demonstrating the effectiveness of the new gyro correction.

The final PACS images are shown in the upper panel of Figure 3. The measured FWHMs for the HD 95086 PACS data are 787 × 747 with P.A. of 112° at 70 μm, 891 × 843 with P.A. of 116° at 100 μm, and 134 × 122 with P.A. of 140° at 160 μm, suggesting that the source is resolved at ∼1.4 beams at both 70/100 μm, but only at ∼1.1 beams at 160 μm. This explains why the P.A. measured at 160 μm is quite different from the ones at 70/100 μm. Based on the resolvability (≲1.4 beam), we simulated model images at all three bands (for details see below) and at different inclinations, compared with the derived FWHMs and aspect ratio, and concluded that the source is inclined at 25° ± 5° at P.A. of 115° ± 10°. These values agree with the numbers derived by Moór et al. (2013).

Zoom In Zoom Out Reset image size

The source is resolved at ≲1.4 beams, consistent with the observation that the images at all three wavelengths are centrally peaked. To estimate the maximum flux of a point source in the data, we performed PSF subtraction using the γ Dra combined PSF by matching the peak value at the center. The PSF-subtracted images are shown in Figure 3, bottom row. The scaled point source fluxes are 397 mJy, 422 mJy, and 357 mJy at 70, 100, and 160 μm, respectively, which are ∼262, ∼574, and ∼1266 times the expected photospheric values at these wavelengths. The scaled point source fluxes are also much brighter than the expected values of the warm component (<10 mJy, see Figure 1). It is clear that a significant amount of dust emission (presumably from a planetesimal belt (P.B.)) is not resolved (up to ∼58%, ∼61%, and ∼73% at 70, 100, and 160 μm, respectively); i.e., ≳30%–40% of the total excess fluxes is in the resolved structure.

We adopted a similar strategy as in Su et al. (2009) and Matthews et al. (2014) to investigate simple structural parameters that would be consistent with the resolved images. The model images were constructed for a given model description, projected to the inclination angle of 25°, then convolved with the appropriate PSF. We then computed the azimuthally averaged radial profiles for both the observed and model data for comparison. Since PACS 70 μm and 100 μm channels have the better resolutions, we derived the best-fit parameters only for the 70/100 μm data, and verified them with the 160 μm data. The model surface brightness (SB) was normalized to the observed profile for data points within 7'' (ensuring good fits in the high signal-to-noise region), then we computed the χ² values for data points within 20'' at 70 μm, and 15'' at 100 μm, respectively. Limited by the marginally resolved images, we only explored a simple power-law distribution (SB(r) ∼ r^−α) with inner (r_in) and outer (r_out) radial cutoffs. With α = 3, we found r_in ∼ 13 and r_out ∼ 9''. With α = 4, we found r_in ∼ 18 and r_out ∼9''. Alternatively, we could also find reasonable fits for the case of α = 1, which is the power index for the P.B. derived for the resolved images of HR 8799 (Matthews et al. 2014). In summary, the power-law index (α) is unconstrained, and the resolved images are best described as a point source plus some form of extended structure. As we discuss in Section 3.2, the resolved structure in the HD 95086 system most likely arises from the disk halo component composed of only small grains that can be warm enough to emit efficiently at far-infrared wavelengths to account for the extended structure seen in the PACS images. We note that the disk halo component is found to have α = 3–4 for HR 8799 (Matthews et al. 2014). With this assumption, the derived r_in (∼120–165 AU) represents the boundary where the halo emission becomes dominant (i.e., the outer radius of the P.B.), while the r_out (∼800 AU) only represents the boundary where the faint disk halo can be well detected given the sensitivity. The best-fit SB distributions are shown in Figure 4 along with the observed profiles. Unlike the HR 8799 SB profiles where two power laws were required to fit the data, one single power-law model provides satisfactory fits to the HD 95086 data.

Zoom In Zoom Out Reset image size

To guide the SED modeling, we first computed the thermal equilibrium dust temperature distribution around HD 95086 (Figure 5) with assumed grain properties in an optically thin environment. We adopted two grain compositions: astronomical silicates and ice-coated silicates. We used the optical constants from Laor & Draine (1993) and Mie theory code to compute the grain absorption and scattering efficiency for compact, spherical astronomical silicates. For the absorption and scattering efficiency of icy grains, we used both multi-layer Mie calculations as well as standard effective medium theory rules (Garnett and Bruggeman equation; e.g., Voshchinnikov et al. 2005). For the inner warm belt, we adopt compact astronomical silicates as the composition for the SED fitting since the majority of the ice should be sublimated at the warm belt temperature. For the outer cold disk, we use the icy grains composed of compact silicate cores with icy mantles that occupy 30% of the grain volume. The thickness of the icy mantles is chosen so that the resultant model emission produces a good overall fit in the 60–70 μm region. Figure 5 shows that in terms of these two grain compositions (non-icy versus icy) the resultant temperatures for grains larger than ∼10 μm are not very different, but they could have ≳10 K difference for μm-size grains closer in. The observed two characteristic temperatures (∼175 K and ∼55 K) suggest two distinct dust locations (∼7 AU and ∼60 AU) for blackbody-like grains.

Zoom In Zoom Out Reset image size

Since fine debris is generated through collisional cascades by breaking up large parent bodies, we expect a full spectrum of particle sizes in a debris disk, generally in a power-law form, n(a) ∼ a^−q with n being the number density, a for the grain radius and q = 3.5–3.7 (Wyatt et al. 2011; Gáspár et al. 2012). Thus, the majority of particles in a collisional dominated disk are grains just slightly larger than the minimum grain size in the distribution (i.e., the average grain size is 〈a〉 ∼ (q − 1/q − 2) a_min where a_min is the minimum size of grains existing in a steady-state collisional disk). In a debris disk (no/very little gas), in addition to the gravitational force of the star, dust grains are subject to radiation pressure. Therefore, grains that are smaller than the blowout size (a_bl) are removed from the system, i.e., a_min ∼ a_bl. The blowout size is a function of stellar mass (M_*), luminosity (L_*), and grain density (ρ) given by a_bl = 1.15 (L_*/L_☉)(M_☉/M_*)(g cm⁻³/ρ) in units of μm. Given the stellar parameters for HD 95086 (L_* ∼ 7 L_☉ and M_* ∼ 1.6M_☉), the radiation blowout radius is 1.5 and 1.8 μm for densities of 3.3 g cm⁻³ (for compact silicates) and 2.6 g cm⁻³ (for icy silicates), respectively. Therefore, the average grain sizes in such a collisionally dominated disk are ∼3 μm. We also set the maximum grain radius (a_max) to be 1000 μm; grains larger than this size have negligible contribution to the model emission.

In our simple SED model, we adopt a uniform particle size distribution at all astro-centric distances and do not consider the effect of grain size segregation due to stellar radiation pressure. The grain-size segregation in debris disks is closely tied to the complex coupling between dynamical and collisional effects (e.g., Thebault et al. 2014); therefore, the impact of using a uniform size distribution is difficult to evaluate without multi-wavelength resolved images. For simplicity, we adopt q = 3.5 for the cold disk component because the cold disk, as we show later, is quite extended (∼63–190 AU), and the particle distribution for a unperturbed wide disk is most likely to have q = 3.5 distribution based on the computation of Thebault et al. (2014). For the inner warm component, on the other hand, we adopt a steeper power law, q = 3.65, based on the hypothesis that the warm component is a narrow belt closer to the star where collisional cascades are expected to reach collisional equilibrium (Gáspár et al. 2012; Thebault et al. 2014). Although the output far-infrared and submillimeter fluxes would be slightly different depending on the q values, the adopted q value has little impact on our model parameters for the warm component since we have no far-infrared/submillimeter constraint.

Two kinds of simple density distributions were used to construct the disk SED. One is a flat disk (constant surface density) with adjustable inner (R_in) and outer (R_in) radii, and the other one is a Gaussian profile ring with a peak radius (R_p) and a width (R_w as the FWHM). With the grain parameters (composition and size distribution) set, we searched for the best-fit two parameters in each density distribution. Because the warm and cold components are not spatially resolved and the (expected) peak emission of the warm component (see Figure 1) is in the region where emission from the cold component also contributes significantly, our SED fitting strategy is to obtain reasonable fits to the cold component first, then derive parameters for the warm component. For the cold component, we only used the data points longward of 35 μm and shortward of 500 μm to derive the best-fit (in χ² sense) size parameters. In the flat disk model, the best-fit disk extension is: R_in = 63 ± 6 AU and R_out = 189 ± 13 AU; while the Gaussian ring model gives R_p = 123 ± 3 AU with a width of R_w = 72 ± 11 AU. For the cold disk, the derived dust fractional luminosity is ∼1.5 × 10⁻³ with dust mass (up to 1000 μm) of 0.18 ± 0.01 M_⊕ for the flat disk, and 0.23 ± 0.01 M_⊕ for the Gaussian ring model. Our dust masses are lower than the value (0.5 ± 0.1 M_⊕) derived by Moór et al. (2013) who only used the 500 μm flux for mass estimation. As described in Section 2.1, faint background galaxies within 40'' of the source are likely to contaminate the submillimeter fluxes. A better estimate of the disk mass will have to come from high spatial resolution data like those from ALMA.

We used non-icy grains for the warm component with a_min = 2.2, a_max = 1000 μm, and q = 3.65. Because the warm component is less constrained, the uncertainty from the SED fitting is large. Assuming a flat disk density distribution, the warm component with radii ranging from ∼7 to ∼10 AU gave a consistent fit to the shape of the IRS spectrum. Alternatively, a Gaussian ring centered at 8 AU with a width of 2 AU resulted in a similar fit. The derived dust fractional luminosity is 1.5 × 10⁻⁴ with a total dust mass of ∼4–5 × 10⁻⁵ M_⊕. The model SED of the two (warm and cold) flat disks is shown in Figure 6(a).

Zoom In Zoom Out Reset image size

Although the two-component model gives an overall good fit to the observed SED, the model images using these parameters do not fit the observations. Figure 7(a) shows the difference image (observation–model) for the two-component disk model at 70 μm, where very negative values at the source position are evident. The extension of the cold disk derived from the observed SED (radius ≲ 200 AU) implies that the cold disk should not have been resolved with Herschel's resolution (FWHM of 56 at best). We can obtain reasonable SED fits to the cold disk by extending it up to ∼800 AU; however, the surface density distribution has to have a steeper drop-off (r^−1.5 compared to a flat disk of r⁰). This very wide planetesimal disk, perhaps analogous to the solar system's scattered disk of KBOs, could generate dust grains in situ from collisions. Using a uniform particle size distribution (a^−3.5) at all astro-centric distances up to ∼800 AU, the resultant model images from such a wide disk with a steep radial gradient of surface density (∼r^−1.5) are similar to the ones using much smaller outer disk radii, i.e., not resolved. In order to account for the detected far-infrared emission at large astro-centric distances up to ∼800 AU, an extended component composed of only small dust grains on unbound or barely bound orbits (i.e., a disk halo) is our best alternative. These small grains are most likely generated in the cold belt, which may include a very wide scattered disk structure. However, given that the source is only resolved at ≲1.4 beams, no meaningful constraint can be placed based on the images.

Zoom In Zoom Out Reset image size

Without some information on the brightness of the P.B., there is a wide range of parameters in this three-component SED model that can produce satisfactory fits. For example, if we assume the contribution of the P.B. in the SED is roughly equal to the point source fluxes that are scaled to match the peak fluxes of the PACS images for subtraction, we are able to construct a three-component model that fits the overall SED (Figure 6(b)) and the far-infrared images. In this three-component model, the warm component is assumed to be the same as the two-component model (it has negligible contribution at far-infrared), and the emission at far-infrared is partitioned into the P.B. (blue line in Figure 6(b)) and disk halo (red line in Figure 6(b)). The P.B. component has a dust fractional luminosity of ∼8 × 10⁻⁴ and a total dust mass of ∼0.2 M_⊕ while the disk halo component has one third (∼3 × 10⁻⁴) of the fractional luminosity in the P.B. component, but only ∼10% of the dust mass as a result of the disk halo being composed of ∼μm size grains. Figure 7(b) shows a detailed comparison between the observed and three-component model images at 70 μm. Overall, the residual suggests the three-component model agrees with the observation within ±2σ.

Due to the large distance of HD 95086, the warm excess is inferred from the SED only, unlike for other nearby debris systems (Vega and Fomalhaut (Su et al. 2013), and Eri (Backman et al. 2009)) where the warm and cold components are spatially separated. Recently, Kennedy & Wyatt (2014) presented an informative way to examine the properties of two-temperature disks, and tested the hypothesis that emission from a single narrow belt can account for the two observed temperatures by varying the particle size distribution. They conclude that the one-narrow-belt scenario is unlikely to produce two observed dust temperatures around early-type stars if the blowout size is preserved. Furthermore, the derived temperature and fractional luminosity ratios between the HD 95086 warm and cold components (R_T ∼ 3.2 and $R_{f_d}\sim$0.1–0.2) are located right at the center of the R_T versus $R_{f_d}$ plot (Figure 5 in Kennedy & Wyatt 2014), indicating that varying the particle size distribution in a narrow belt is unlikely to explain the observed SED.

Another possibility is that the warm excess emission could result from the small amount of particles being dragged in due to Poynting–Robertson (P-R) drag. The effect of P-R drag in collisional dominated debris disks has been studied analytically by Wyatt (2005). In the case of the HD 95086 disk (a P.B. at ∼100 AU with f_{d, cold} ∼ 10⁻³ producing excess emission at ∼10 AU), the maximum amount of material available due to P-R drag in terms of optical depth is <0.05f_{d, cold}. Furthermore, the analytical maximum value was found to be a factor of ∼7 higher than using a numerical model with a better treatment of collisional effects (van Lieshout et al. 2014). Thus, the observed value f_{d, warm} ∼1.5 × 10⁻⁴ for the warm component, is a factor of 20 higher than the maximum value that the cold belt can supply due to P-R drag, making this hypothesis unlikely. We, therefore, consider the warm component as an independent belt in the following discussion.

The 69 μm emission band, a solid state feature originating from magnesium-rich olivine crystals, was first observed astronomically by the Infrared Space Observatory around a Herbig Ae/Be star, HD 100546 (Malfait et al. 1998). The peak position and strength (FWHM) of the 69 μm feature strongly depends on the iron content of the crystals and their temperatures (Bowey et al. 2002; Koike et al. 2003; Suto et al. 2006). See Sturm et al. (2013) and Maaskant et al. (2014) for recent reviews on the iron-content and temperature dependence, and other secondary effects like grain sizes and shapes. The feature detected in the HD 95086 MIPS-SED data is not spectrally resolved, and we cannot put constraints on the fraction of iron composition as has been done for β Pic (de Vries et al. 2012); however, the iron fraction has to be less than a few percent to have the feature peaked at 69 μm at ∼50 K.

To further validate the detection of the forsterite 69 μm feature and derive a rough estimate of the crystalline abundance, we took the mass absorption coefficients of synthetic forsterite (Fo100 in Koike et al. 2003), and constructed simple model fits to the MIPS-SED data. The model has two components: (1) emission from the crystalline forsterite and (2) a blackbody emission representing the underlying dust continuum, with both at the temperature of 55 K derived for the cold component. The model spectrum was smoothed to match the MIPS-SED resolution, and the ratio between the two components was adjusted to fit the overall shape of the MIPS-SED data. Since we do not have information about the location where the 69 μm feature originates, the temperature of this component is not necessarily the same as that of the cold disk (underlying continuum). However, we can easily rule out the possibility that these crystalline grains originate from the warm belt with a dust temperature of ∼175 K since the blue part (50–60 μm) of the model spectrum would have a much bluer slope. Given that the continuum in both blue and red sides of the feature is relatively flat, we assume the crystalline grains are located in the cold disk with the same dust temperature as the underlying continuum.

The upper panel of Figure 8 shows the best-fit model where the total mass of the forsterite component is ∼6.2 × 10²⁵ g (0.01 M_⊕). Using the total dust mass we derived from the cold disk (0.2 M_⊕, Section 3.2), the abundance of the crystalline grains is ∼5%. This abundance is likely to be a lower limit because the laboratory measured absorption coefficients are based on finely ground powder (sizes of 0.1–0.5 μm), and the 69 μm feature is not sensitive to grain size (the band profiles are similar for grains less than ∼10 μm, Sturm et al. 2013).⁴ Since no prominent crystalline features are seen in the observed 20–40 μm range (while our best-fit model predicts strong features), the crystalline grains must have sizes ≳ a few micrometers. Therefore, the derived abundance should be taken as a lower limit.

Zoom In Zoom Out Reset image size

The presence of crystalline silicates at a low temperature and far from the star is surprising, given that relatively high temperatures are required to convert circumstellar amorphous silicates to the crystalline form (Fabian et al. 2000). The detection of the forsterite emission feature implies that there is a substantial transport of material from the inner zone of the system, or that material has been released from the disruption of a large body where high temperatures were reached during formation.

It is interesting to note that the 69 μm feature was not detected in the MIPS-SED data of β Pic (the bottom panel of Figure 8), while the crystalline abundance is estimated to be 3.6% ± 1.0% using the PACS spectroscopic data (de Vries et al. 2012). With the published PACS spectrum smoothed to match the MIPS-SED resolution, we found that the low equivalent width of the detected feature is completely washed out in the MIPS-SED mode (the inset of the bottom panel of Figure 8). This is also consistent with the derived mass values for the crystalline grains (2.8 × 10²³ g for the β Pic disk, while the mass we estimate for HD 95086 is ≳200 times more). The estimated crystalline abundances implied by the observed feature strengths in these two systems appear discrepant. The discrepancy lies in the total dust mass used to derive the abundance, which is model dependent. The SED models used for the β Pic system include grains up to 100 μm for all three-temperature components (de Vries et al. 2012), while our models for HD 95086 include grains up to 1000 μm, i.e., a factor of ∼3 in the derived dust mass in a typical a^−3.5 size power law ($M_d\sim \sqrt{a_{{\rm max}}}$). Indeed, the total non-crystalline, derived dust mass by de Vries et al. (2012) is a factor of ∼3 lower than the value (4.7 × 10²⁶ g) derived from recent ALMA observation for β Pic (Dent et al. 2014). Furthermore, the reported 3.6% crystalline abundance is only based on one of the temperature components (∼90 K with a total mass of 8 × 10²⁴ g). Using the ALMA derived dust mass, the crystalline abundance in the β Pic disk should be ≲1%. As we compare the disk properties between HD 95086 and HR 8799 in the following subsection, it is also noteworthy that no feature was detected in the MIPS-SED data around HR 8799 (Figure 8; Su et al. 2009), suggesting the crystalline abundance is <1%.

The eccentricity, e, of debris disks usually can be measured in two ways: (1) the offset between the ring center and star position, and (2) the degree of asymmetry seen in thermal images (i.e., pericenter glow). Both methods have been applied to the Fomalhaut debris disk (Kalas et al. 2005; Stapelfeldt et al. 2004) and result in a consistent value (e ∼ 0.1). Since the star contributes negligibly in the PACS wavelengths, the centroid we measured from the PACS images is the center of the disk. Based on the calibration observations of the PSF star, γ Dra, a typical centroiding uncertainty is 001. The position of the star is well measured by Two Micron All Sky Survey (2MASS) with a known proper motion from Hipparcos; therefore, the uncertainty in its position relative to the Herschel detection is dominated by the absolute accuracy of the World Coordinate Systems in Herschel/PACS maps. We used a bright galaxy (LEDA 3079352, 2MASS J10571665−6838171) present in the field of our PACS images to estimate the absolute pointing uncertainty. There appears to be an offset of ∼22 between the 2MASS and Herschel coordinate systems. The typical pointing uncertainty in Herschel observations (∼2'') makes it difficult to pin point the star position with high accuracy (i.e., the first method is inapplicable).

We then turned to the second method to assess the degree of asymmetry we could have detected with Herschel observations. Since the resolved structure is mostly dominated by the disk halo (the P.B. is not resolved at Herschel's resolution), the presumption is that the disk halo shape follows closely with that of the P.B. This assumption is reasonable if the disk halo results from small grains generated in the parent-body belt that are ejected and/or placed in highly eccentric orbits due to radiation pressure. There appears to be some asymmetry seen in the PSF-subtracted images (bottom row of Figure 3); however, it is not significant (less than 3σ). The resolved component is consistent with the circular case with e = 0 (our baseline model shown in Figures 6(b) and 7(b)). We then constructed model images by placing the disk halo components at various eccentricities and compared them with our baseline model. We found that the difference is noticeable (more than 5σ) when e > 0.25. From this exercise, we conclude that the resolved structure (disk halo) around HD 95086 has e < 0.3.

Given that both HR 8799 and HD 95086 have bright debris disks and wide-orbit planets revealed via direct imaging, it is instructive to compare their properties in regard to different aspects of the planetary configuration. Table 1 lists the properties of the HR 8799 and HD 95086 systems. The debris structures revealed from excess SEDs and resolved images show that both systems have a warm (∼170 K) excess near the water-ice line and a cold excess at ∼50 K, surrounded by an extended disk halo. The HD 95086 system is ∼10 times dustier than the HR 8799 in terms of fractional dust luminosity, but the ratio of fractional luminosities of the warm and cold components is similar (∼0.1). The extended disk halo in HR 8799 can be traced up to ∼2000 AU with a sensitivity of 10⁻² mJy arcsec² (Matthews et al. 2014); interestingly, a similar size is seen in the HD 95086 system with the same detection limit, implying its disk halo is >5 times brighter given that it is ∼2.3 times more distant. Based on the resolved images, Matthews et al. (2014) estimated the HR 8799 disk halo accounts for 40%, 46%, and 52% of the total flux at 70, 100, and 160 μm, respectively. In the case of HD 95086, the marginally resolved images can only give a lower limit for the contribution from the disk halo, which is ≳42%, ≳39%, and ≳27% at 70, 100, and 160 μm, respectively. In other words, the debris structures in these two systems are very similar except that the debris in the HD 95086 system is ∼10 times more in dust luminosities and two to five times more in observed dust masses (derived from SED fitting) compared to those values in the HR 8799 system.

Table 1. Properties Comparison between HR 8799 and HD 95086

	HR 8799		HD 95086
Stellar Parameter
M* (M☉)	1.5		1.6
L* (L☉)	5.7		7.0
Teff (K)	7500		7500
R* (R☉)	1.4		1.6
Age (Myr)	∼30–90(1)		∼20(2)
Distance (pc)	39.4		90.4
SED parameter
	Warm	Cold	Warm	Cold
Td (K)	∼150	∼45	∼175	∼55
fd	1.8E −5	2.2E −4	1.5E −4	1.5E −3
Md (M⊕)	1E −6	0.1	5E −5	0.2–0.5
Disk orientation
i (deg)	25 ± 5		25 ± 5
P.A. (deg)	62 ± 10		115 ± 10
Planet	e, d, c, b		b
Mp (MJ)	9 ± 2, 9 ± 3, 9 ± 3, 7 ± 2(3)		5 ± 2(4)
ap (AU)	15.4, 25.4, 39.4, 69.1		61.5$^{+5.7}_{-4.9}$(5)

References. (1) Baines et al. 2012, see Section 4.4 for details; (2) Meshkat et al. 2013; (3) Goździewski & Migaszewski 2014; (4) Rameau et al. 2013b; (5) this work, assuming e_b = 0.

Download table as:  ASCIITypeset image

That both systems are 10–100 Myr in age is exciting because this is the epoch when interesting processes like dynamical settling, formation of rocky planets, and possibly even (late) formation of ice giants are expected to occur in planetary systems. The age estimate of HD 95086 (∼20 Myr) is quite robust given its association with the LCC moving group, but the age of HR 8799 has been widely debated. Most of the dating methods applicable for early-type stars favor an age younger than ∼150 Myr (Marois et al. 2008), while asteroseismic analysis based on the ground-based data favors an older age of ∼1 Gyr (Moya et al. 2010). Recent optical monitoring using the Microvariability and Oscillations in STars (MOST) space data found significant frequency and amplitude changes in HR 8799's pulsation, making precise asteroseismic analysis difficult (Sódor et al. 2014). The age of HR 8799 was revisited by Baines et al. (2012) with the directly measured stellar diameter from interferometric CHARA data and evolutionary models, who derived two age estimates: 33$^{+7}_{-13.2}$ Myr if the star is still contracting toward the zero-age main sequence (ZAMS), and 90$^{+381}_{-50}$ Myr if the star is evolving from ZAMS. Therefore, it is possible that HR 8799 could be as young as HD 95086 (20 Myr), but it is more likely older, by up to ≳10 times.

If both HR 8799 and HD 95086 are at the same age, the difference in the observed dust level must reflect different initial conditions. Based on the evolution models of debris disks developed by Kenyon & Bromley (2008), the relative disk luminosity (f_d) around a 2 M_☉ star at age of 20 Myr can be up to 10 times higher if one disk started with 10 times more mass than the other (Figure 17 in that paper). This suggests that HR 8799 was born with a less massive protoplanetary disk than that of HD 95086 if both have the same age. However, all four directly imaged planets in the HR 8799 system are ∼2 times more massive than the HD 95086 planet, implying the opposite trend in initial conditions. Alternatively, if both systems formed with similar initial conditions (mass and angular momentum), the difference in the dust level must be related to debris evolution. The models by Kenyon & Bromley (2008) predict the decline in disk luminosity is roughly a power law, f_d∝t⁻ⁿ with n ≈ 0.6–1.0. With possible ages of 40–200 Myr for HR 8799, the differences in f_d between HD 95086 and HR 8799 suggest n ∼1.0–3.3, a steeper decline than the model prediction.

Perhaps the biggest difference between these two systems is that HR 8799 hosts four massive (7–9 M_J) planets between the warm and cold belts, while there is only one ∼5 M_J planet currently known that is just interior to the HD 95086 cold belt. Because HD 95086 is more than two times more distant than HR 8799, the current direct imaging observations are only sensitive to ≳5 M_J as close as ∼30 AU to the star (Rameau et al. 2013a). Similar to the well-resolved two-belt systems (HR 8799, Vega, Fomalhaut, and Eri), the large gap between the warm and cold excesses might host multiple planets.

We can constrain possible planet configurations of the HD 95086 system based on two dynamical considerations: (1) multiple planets must be dynamically stable for a time span of the age of the HD 95086 system, ∼20 Myr, and (2) the planets' gravitational perturbations must clear and maintain the gap between the warm and cold belts. For the former, we make use of the numerical results of Faber & Quillen (2007) on the stability of multiple planet systems in coplanar configurations of roughly geometrically spaced orbits. For the latter we make use of numerical estimates of a planet's chaotic zone in which debris can be cleared and a gap can be maintained by each planet. The chaotic zone is owed to the overlap of first order mean motion resonances, as shown by Wisdom (1980), who estimated the chaotic zone width as Δa ≃ 1.3μ^2/7a_p, where μ = M_p/M_*, and M_p, a_p are the planet's mass and orbital semi-major axis. In this zone, initially circular test particle orbits exhibit strongly chaotic motion. Subsequent numerical studies have supported this estimate for planet-to-star mass ratios, μ, up to about 10⁻³ (Duncan et al. 1989). Morrison & Malhotra (2014) carried out numerical integrations for μ up to 10^−1.5, to determine the range in the planet's chaotic zone where particles are not merely chaotic, but actually cleared (by either collision or escape); they report cleared zone widths, both interior and exterior to a planet's circular orbit, as follows: Δa_{cl, int} = 1.2μ^0.28a_p and Δa_{cl, ext} = 1.7μ^0.31a_p, In this section, we use these "cleared zone" widths to estimate the effects of planets on the clearing of debris. However, we continue to refer to these as "chaotic zone" widths, to acknowledge the underlying dynamical mechanism of clearing. Below, we consider possible planetary systems in coplanar, non-resonant configurations.

4.5.1. Single-planet System

4.5.2. Two-planet System

Let us consider whether the gap could be maintained by two planets on low eccentricity (e ∼ 0.1) orbits. The known planet, HD 95086b, of mass ∼5 ± 2 M_J on a nearly circular orbit of radius ∼61.5 ± 5 AU would both shepherd the outer belt as well as clear debris inward of its orbit by about one chaotic zone width, i.e., down to about 42 AU at most. The second unseen planet must be located near the inner belt. Assuming that this planet is of mass ∼13 M_J (the conventional upper limit of planetary masses), and its pericenter is located one chaotic zone width (Δa_{cl, int}) away from the outer edge of the inner belt, we estimate that its semi-major axis is ∼13 AU. Such a putative planet could shepherd the inner belt and clear debris out to a maximum distance of about 18 AU. This is insufficient to account for the entire 8–80 AU observed gap. We can also consider a higher mass, unseen sub-stellar object near the inner belt. Currently, the best detection limit (5σ) for low-mass companions close to the star is ∼21 M_J at an astro-centric distance of 018 (∼16 AU) (Meshkat et al. 2014). Demanding that a 21 M_J object be able to shepherd the inner belt edge at ∼8 AU implies that its pericenter is located at about one chaotic zone width away, i.e., at about 14 AU. This means that it can clear debris out to about (14 + Δa_{cl, ext}) ≈ 20 AU. This is also insufficient to account for the entire 8–80 AU gap. We conclude that two planets on low eccentricity orbits are unable to account for the gap.

If we relax the low eccentricity condition, then it is possible to account for the gap provided one or both planets' orbits are of significant eccentricity. First, let us consider that HD 95086b is of low eccentricity but the inner unseen planet is of high eccentricity. Then, as we noted above, HD 95086b can clear debris inward of its orbit to about ∼42 AU astro-centric distance, so we require the inner unseen planet to account for debris clearing in the range 8–42 AU. Adopting the maximum mass for the inner planet, M₁ = 13 M_J, we can estimate its pericenter and apocenter distances: q₁ ≃ 8/(1–1.2(M₁/M_*)^0.28) ≃ 12 AU, and Q₁ ≃ 42/(1 + 1.7(M₁/M_*)^0.31) ≃ 30 AU. Therefore, the inner planet (of maximum planet mass 13 M_J) must have orbital eccentricity ∼0.5 to account for the observed gap. Solutions with somewhat lower planet mass and eccentricity are possible. For example, consider two planets of similar mass, 5 M_J, and similar orbital eccentricity, ∼0.3, with the inner planet's orbit of semi-major axis ∼16 AU and HD 95086b's orbit of semi-major axis ∼40 AU. Their pericenter and apocenter distances, augmented with one chaotic zone width each, could shepherd the warm and cold belt edges as well as clear the 8–80 AU gap; moreover, the planets' orbits are also marginally stable for the 20 Myr age of the system. This two-planet configuration is illustrated in Figure 9 as "system B."

Zoom In Zoom Out Reset image size

4.5.3. Three-planet System

4.5.4. Four-planet System