DISCOVERY OF RESOLVED DEBRIS DISK AROUND HD 131835

Planet formation and evolutionary history is imprinted on the nature and distribution of circumstellar debris. Circumstellar debris disks are composed of dust and planetesimals that are in orbit around main-sequence stars, analogous to the asteroid belt and Kuiper Belt in our solar system. The materials in the disks are thought to be the second generation where the dust is constantly being replenished (Backman & Paresce 1993) from collisions between planetesimals and sublimation of comets (Harper et al. 1984). Most of the debris disks are discovered through detecting their infrared (IR) excess and characterized by their spectral energy distributions (SEDs). However, studying their SEDs alone provides only limited information. While the disk temperature and the total IR flux can be well determined from the SEDs, the spatial information of the systems is inaccessible. In addition, there is a possibility that the excess actually originates from an unrelated foreground or background source contaminating the beam.

In addition, with the SED information only, the location of the grains will be degenerate with the grain properties. For example, small grains farther away from the star might be heated to the same temperature as large grains that are closer toward the center. Thus, a single-temperature debris disk does not necessarily indicate that all the grains are located in a thin ring. A similar argument applies to the opposite scenario where small grains might be heated to a higher temperature compared to large grains at the same location. Therefore, for a multitemperature debris disk, the SED information alone is typically insufficient to distinguish between the cases where grains are located at the same spatial location (i.e., a ring) and those where grains are located at multiple spatially distinct locations (i.e., belts or extended disks). Kennedy & Wyatt (2014) argued that this degeneracy is small for A-type stars, owing to the truncation on the small end of the grain size distribution by radiation pressure. They proposed that most two-temperature disks around A-type stars probably arise from multiple spatial components. To confirm that the disk indeed has multiple spatial components, resolved images are required.

If the disks can be spatially resolved, we will be able to characterize the spatial distribution of the grains and the geometry of the disks. We can use distribution of grains to probe the underlying physics responsible for dust distribution and constrain the locations and nature of the grain parent bodies. In addition, the morphology of the planetesimal belts might be closely related to the presence of unseen planets dynamically sculpting the disk. For example, as seen in mid-IR thermal emission, HR 4796A has an inclined ring with one lobe being ∼5% brighter than the other (Telesco et al. 2000). Wyatt et al. (1999) argued that the secular perturbations from a planet of mass >0.1 M_J located close to the inner edge of the disk could introduce this brightness asymmetry. The planet would impose the forced eccentricity on the ring, causing one side of the ring to be closer to the star and thus showing the "pericenter glow." Such incidence demonstrates that the level of information provided from the direct imaging far surpasses that from the SED alone.

Here we have attempted to spatially resolve the disk around HD 131835, a system characterized solely by its SED prior to this work. Rizzuto et al. (2011) assign a 91% membership probability to the Upper Centaurus Lupus (UCL) moving group (a subgroup of the Sco–Cen association) based on its Galactic location and velocity. HD 131835 is a young star (∼15 Myr, based on the age of the UCL estimated by Mamajek et al. 2002; ∼16 Myr, based on an analysis of F-type pre-main-sequence members of the group by Pecaut et al. 2012) with a spectral type of A2IV (Houk 1982) at a distance of d = 122.7$_{-12.8}^{+16.2}$ pc (van Leeuwen 2007). The IR excess emission of HD 131835 was first reported by Moór et al. (2006) through searching the IRAS and Infrared Space Observatory (ISO) databases for the stars in the vicinity of the Sun. The MIPS data for this source were first published in Chen et al. (2012), which shows the MIPS 24 and 70 μm data for all of the de Zeeuw et al. (1999) Sco–Cen A-type stars and helps to put the excess emission for this star into context (e.g., it is one of only four UCL/Lower Centaurus Crux A-type stars with ${{L}_{{\rm IR}}}/{{L}_{*}}$ commensurate to beta Pic, $\gt {{10}^{-3}}$). Chen et al. (2014) found that the SED for HD 131835 can be described using a two-temperature model.

In contradiction to the traditional disk evolution scheme, CO gas in debris disks around A-type stars may be fairly common. Typically, the primordial gas dissipates in a few-million-year timescale during the protoplanetary phase. By the time it transforms into a debris system, the CO gas level will be undetectable. However, there are several exceptional cases. For example, 49 Cet and HD 21997 are 40 (Zuckerman & Song 2012) and 30 (Torres et al. 2008) Myr old A-type stars hosting gas-rich disks (Zuckerman et al. 1995; Moór et al. 2011). The gas is believed to be the second generation, possibly due to violent comet collisions (Zuckerman & Song 2012). In the case of HD 21997, there is an alternative explanation based on a recent Atacama Large Millimeter/submillimeter Array (ALMA) observation that indicates that the gas may be of primordial origin in this system (Kóspál et al. 2013). Another possible production mechanism involves constant resurfacing of the parent bodies and sublimation or photodesorption of the CO ice. For the β Pic system, its CO distribution is particularly interesting because it is not the same as the dust, is even more highly asymmetric, and may imply the presence of a second undetected planet (Dent et al. 2014). Not only has CO been detected and characterized for β Pic, 49 Cet, and HD 21997, but it has also been detected around five A-type stars in Upper Sco (J. Carpenter 2014, private communication).

Besides having the IR excess from the dusty debris, HD 131835 hosts a detectable amount of carbon monoxide gas (Moór et al. 2013; A. Moór et al., in preparation). For HD 131835, Kastner et al. (2010) first reported the nondetection of CO emission with the 30 m Institut de Radio Astronomie Millimetrique telescope. However, Zuckerman & Song (2012) argued that if the comet collision model is correct, then the H₂/CO ratio is unconstrained and thus the upper mass limit of the CO gas for HD 131835 is $4.06\times {{10}^{-3}}$ M_E. Later, Moór et al. (2013) announced the discovery of detected submillimeter CO emission with the Atacama Pathfinder EXperiment (APEX) radio telescope through a survey. Once the CO gas is well characterized, it can provide us valuable information on the disk environment and the dust–gas associations of such systems with the coexistence of gas and debris at an age ≳10 Myr.

Here we report the discovery of the spatially resolved debris disk around HD 131835 at mid-IR wavelengths, tracing the thermal emission from the grains. Currently, these images are the only resolved images of the system. In this paper, we present the analysis, modeling, and characterization of the debris disk around HD 131835. Our mid-IR images are based on two epochs of observation as described in Section 2. The data processing details, including image reduction, photometry, and point-spread function (PSF) subtraction, are in Section 3. In Section 4 we analyze the stellar properties and show that a three-component model is required to describe the system. In Section 5 we characterize the grain temperatures in the system and investigate the possible grain compositions. We also briefly discuss how the current generation of imaging instruments and telescopes can improve our understanding of debris disks in the near future. Finally, we summarize our modeling results in Section 6.

We used the T-ReCS instrument on the Gemini South telescope to obtain mid-IR imaging of HD 131835 in two programs, GS-2008A-Q-40 and GS-2009A-Q-19. As part of these programs, the star HD 136422 was also observed as a photometric and PSF calibration standard. The observations of the reference star were interleaved between the observations of the target. All the images are taken in the Si-5 (λ_c = 11.66 μm, Δλ = 1.13 μm) and Qa (λ_c = 18.30 μm, Δλ = 1.51 μm) filters. In these bands, the T-ReCS pixel size is $0\buildrel{\prime\prime}\over{.} 08976\;\pm 0\buildrel{\prime\prime}\over{.} 00021$ on the sky.⁵ Owing to the high sky background and instrumental thermal background, we used the standard mid-IR ABBA chop-nod observing technique with a 15'' throw between chop-nod positions.

Since the disk had not been resolved prior to the 2008 observations, the first-epoch data were obtained with an arbitrary chop position angle. This turned out to be very close to the apparent position angle of the disk. To confirm that the feature was truly from the emission of the resolved disk and not from artifacts due to imperfect chopping and nodding, the second epoch of observation was obtained in 2009 with the chop position angle roughly perpendicular to the previous one. Observation conditions were generally very good, with diffraction rings visible in the calibration images. The sky was very transparent, with low water vapor ($\lt 1$ mm) measured toward zenith for most of the nights. The processed images of HD 136422 gave diffraction-limited resolution of $\sim 0\buildrel{\prime\prime}\over{.} 39$ and $\sim 0\buildrel{\prime\prime}\over{.} 54$ at 11.7 and 18.3 μm, respectively. HD 131835 (and HD 136422) were observed in two nights in 2008 and four nights in 2009, with the total on-source time of 6168 (2056) s at 11.7 μm and 5734 (1303) s at 18.3 μm. Details of the observations are summarized in Table 1.

3.1. Image Reduction

We reduced all the images following the standard mid-IR data reduction procedures. We linearly combined the chop and nod frames within each image to remove the instrumental and sky background. We then subtracted the average row background and corrected for the column offset due to small bias variations among different channels. The cores of the PSF standard star HD 136422 are fairly circular (as shown in Figures 1(c) and (d)) with the FWHM of $\sim 0\buildrel{\prime\prime}\over{.} 39$ and $\sim 0\buildrel{\prime\prime}\over{.} 54$ at 11.7 and 18.3 μm, respectively. The target images in both wavelengths show the extended and elongated emission compared to the PSFs. Before the PSF subtraction process, these mid-IR images already show prominent disk emission extending beyond $1^{\prime\prime}$ along the semimajor axis in both bands.

Zoom In Zoom Out Reset image size

3.2. Photometry

We performed aperture photometry on HD 131835 using HD 136422 as the photometric reference. The flux of HD 136422 was taken from Cohen et al. (1999). We applied aperture correction to our photometry measurement on HD 131835. Based on the observations of HD 136422, the PSF structure was fairly stable over the course of each night; the enclosed flux in an aperture radius of 1'' usually varied only from 1% to 3%. Thus, we adopted one relationship of aperture radius versus enclosed flux fraction for each night. Owing to the nonuniform residual background, the enclosed flux fluctuated at large aperture radii. We characterized these fluctuations as the uncertainties of the radius-enclosed flux relationships.

Next, we applied atmospheric extinction correction and color correction. The atmospheric extinction varied from 0.09 to 0.19 mag per airmass among different nights in both bands. These values are consistent with the extinction observed on Mauna Kea (Krisciunas et al. 1987), where the conditions are similar to Cerro Pachón. For color correction, since the 11.7 and 18.3 μm fell within the wavelength coverage of the IRS spectrum, we used its shape to compute the correction factors. Detailed information about the IRS data is provided in Section 4. The final flux of HD 131835 is 49.3 ± 7.6 mJy for 11.7 μm and 84 ± 45 mJy for 18.3 μm.

We have considered uncertainties from the process of measuring the target's flux, the standard star's flux, the aperture correction, and the extinction correction. However, by comparing the photometry measurement from each image (33 images at 11.7 μm and 14 images at 18.3 μm), we noticed the systematic errors dominated over the random errors. Thus, being conservative, we quoted the sample standard deviations of photometry measurements from the ensemble of images as the uncertainties of our mean fluxes. Our measured fluxes are consistent with the values from the IRS spectrum. Table 2 lists the flux measurements of HD 131835 from this work and from other published literature, and the SED is shown in Figure 2.

Zoom In Zoom Out Reset image size

Table 2.  Summary of the Measured Fluxes of HD 131835

			Color Corrected
Wavelength (μm)	Source	Flux (mJy)	Flux (mJy)	References
0.43	Hipparcos	2571 ± 48	⋯	(1)
0.55	Hipparcos	2728 ± 26	⋯	(1)
1.24	2MASS	1448 ± 29	⋯	(2)
1.66	2MASS	965 ± 35	⋯	(2)
2.16	2MASS	652 ± 11	⋯	(2)
3.4	WISE	237.2 ± 5.8a,b	⋯	(3)
4.6	WISE	166.8 ± 4.7a	165.3 ± 4.7	(3) (4)
11.7	Gemini South/T-ReCS	49.7 ± 7.7	49.3 ± 7.6	(5)
12	WISE	49.1 ± 2.2a	49.4 ± 2.2	(3) (4)
18.3	Gemini South/T-ReCS	83 ± 44	84 ± 45	(5)
22	WISE	160.5 ± 9.4a,c	⋯	(3)
24	Spitzer/MIPS	153.1 ± 3.1	161.7 ± 3.3	(6) (7) (8)
25	IRAS	186 ± 34	224 ± 40	(9)
60	IRAS	684 ± 62	681 ± 61	(9)
70	Spitzer/MIPS	659.2 ± 44.7d	710.0 ± 48.1	(7) (8) (10)
90	AKARI /FIS	560 ± 39	583 ± 41	(11) (12)
870	APEX/LABOCA	8.5 ± 4.4	⋯	(13)
Spectrum
5.2–37.9	Spitzer/IRS	⋯	⋯	(6)

References. (1) Høg et al. 2000; (2) Cutri et al. 2003; (3) Cutri et al. 2012; (4) Wright et al. 2010; (5) this paper; (6) Chen et al. 2014; (7) MIPS Instrument Handbook (Version 3.0, 2011 March) by MIPS Instrument and MIPS Instrument Support Teams; (8) Rieke et al. 2004; (9) Moór et al. 2006; (10) Chen et al. 2012; (11) Yamamura et al. 2010; (12) Liu et al. 2014; (13) Nilsson et al. 2010.

^aThe quoted WISE data are based on aperture photometry. The error has included the uncertainty from the rms scatter in the standard calibration stars. ^bThis is a lower limit due to saturation. ^cThis is an upper limit due to source confusion. ^dThe error is calculated by considering the source photon counting uncertainty, the detector repeatability uncertainty, and the absolute calibration uncertainty.

Download table as:  ASCIITypeset image

3.3. PSF Subtraction

The IR excess of HD 131835 has been measured in multiple wavelengths. Table 2 shows the summary of the measurements taken from the published literature. HD 131835 was detected in all four Wide-field Infrared Survey Explorer (WISE) bands (Cutri et al. 2012). However, the WISE w1 measurement represents a lower limit due to saturation, and the w4 measurement is an upper limit due to source confusion, as flagged in the WISE catalog by Cutri et al. (2012). The IRAS photometry measurements are taken from Moór et al. (2006). The Spitzer/MIPS 24 and 70 μm measurements are taken from Chen et al. (2012, 2014). The object 1456545-354138 in the AKARI catalog (Yamamura et al. 2010) was identified as HD 131835 by Liu et al. (2014). Among the four AKARI/FIS bands, only the source measured in 90 μm is confirmed and reliable (Yamamura et al. 2010). Finally, the longest wavelength point on the SED was measured in 870 μm by Nilsson et al. (2010).

We apply color correction to good photometry measurements between 3 and 100 μm. Moór et al. (2006) had computed the color-corrected fluxes for the IRAS points. For other points that fall within the wavelength coverage of the IRS spectrum, we use the shape of the IRS spectrum to determine the correction factors. The other photometry points are corrected based on the best-estimated local SED shape. Color-corrected fluxes and references used to compute the correction factors are listed in Table 2. We do not consider correcting the 870 μm point since it is on the Rayleigh–Jeans tail, so the correction factor should be close to unity. If a point is color corrected, the corrected flux is used in the following modeling and analysis. Figure 2 shows the SED (with the color-corrected fluxes). The SED shows a single-hump IR excess with the flux peaking around 70 μm. The system has ${{L}_{{\rm IR}}}/{{L}_{*}}\ \sim 2\times {{10}^{-3}}$.

The Spitzer/IRS data were obtained and analyzed by Chen et al. (2014). It is generally presumed that the MIPS fluxes will have a lower absolute calibration uncertainty. Therefore, Chen et al. (2014) calibrated the IRS spectrum by tabulating the synthetic MIPS flux and comparing it to the measured MIPS flux. They presumed that the source did not vary with time and then scaled the IRS spectrum to match the MIPS 24 μm flux. More details of the spectral extraction and calibration can be found in Chen et al. (2014). The calibrated IRS spectra do not show any obvious solid-state emission or absorption features, suggesting that HD 131835 is not a system rich in small silicate grains and polycyclic aromatic hydrocarbons (PAHs). For the following analysis, we treat each point of the spectra as a single measurement that is directly comparable to photometry measurements. When including images in the fit, each pixel is considered as a data point and is treated equally as a point on the SED. For each image, only the central 68 by 68 pixels ($\sim 6^{\prime\prime}$ by $6^{\prime\prime}$) are considered in the fit. There are a total of 9595 data points, of which 10 points come from broadband photometry, 337 points come from the IRS spectrum, and 9248 points come from the images.

4.1. Stellar Properties

The measured flux of the system includes the contribution from the star and from the disk. The star has T_eff = 8770 K, logg = 4.0, and solar metallicity with A(V) = 0.187 mag (Chen et al. 2012). We use the IDL Astrolib routine ccm_unred.pro, which is based on Cardelli et al. (1989), to apply reddening to the Kurucz (1993) stellar atmosphere model. Since the disk emission is prominent in longer wavelengths, only the measurements with wavelength shorter than 5.8 μm are used to fit for the stellar flux to characterize the stellar contribution. We compare the Hipparcos and Two Micron All Sky Survey (2MASS) photometry to the synthetic fluxes of the model during the fitting process. The best fit is shown in Figure 2. The stellar fluxes at 11.7 and 18.3 μm of this best-fit model were used to scale the PSFs during the PSF subtraction process described in Section 3.3.

We parameterize a scaling factor as $\xi \equiv \ {{({{R}_{*}}/d)}^{2}}$ so that we can write

$$\begin{eqnarray}&&{{L}_{*}}=4\pi {{d}^{2}}\xi \sigma {{T}^{4}}=4\pi {{d}^{2}}\int _{0}^{\infty }{{I}_{\nu }}d\nu .\end{eqnarray} \tag{ 1 }$$

By integrating the unreddened stellar component of the best-fit SED model and solving the equation, we obtain $\xi =6.68\times {{10}^{-20}}$. When estimating the luminosity, the largest uncertainty comes from the distance measurement. Taking the distance to be $122.7\pm _{12.8}^{16.2}$ pc (van Leeuwen 2007), we find ${{R}_{*}}=1.41\pm _{0.15}^{0.19}{{R}_{}}$ and ${{L}_{*}}=10.5\pm _{2.2}^{2.8}{{L}_{}}$.

4.2. Model

Our goal is to find a disk model that can reproduce the observed SED and the images simultaneously. The SED sets strong constraints on the grain temperature distribution, whereas the images inform the spatial distribution of the grains. We aim to recover the extended mid-IR emission, with the disk flux peaking close to the center of the system, while simultaneously reproducing the broad IR excess hump on the SED. With the space and temperature distributions of the grains so constrained, we hope to further characterize the grain properties such as its size and composition. In our model, we assume the disk to be optically thin since its ${{L}_{{\rm IR}}}/{{L}_{*}}\ \sim 2\times {{10}^{-3}}$. To search for a model that accounts for the extended emission, we start with a simple two-dimensional continuous disk composed of only a single population of grains and then consider more complicated distributions.

4.2.1. A Continuous Power-law Disk Model

In contrast to perfect blackbodies, dust grains do not couple efficiently to radiation when the wavelength is much larger than the grain size. We consider a simple parameterization of the emissivity with a modified blackbody function, whereby the emissivity varies with the frequency as a power law with a positive index β. The grain temperature T_g at a distance r away from the star is found by balancing the energy intake and the energy output:

$$\begin{eqnarray}&&\frac{{{d}^{2}}}{{{r}^{2}}}\int {{F}_{\nu ,*}}\ {{\nu }^{\beta }}d\nu =4\pi \int {{B}_{\nu }}[{{T}_{g}}(r)]\ {{\nu }^{\beta }}d\nu .\end{eqnarray} \tag{ 2 }$$

${{B}_{\nu }}(T)$ here is the blackbody function. The grains are distributed with the surface density of ${{n}_{0}}\ {{(r/1\;{\rm AU})}^{{\Gamma }}}$, with n₀ being the two-dimensional number density of the grains at 1 AU. The density law applies to the region between the inner disk radius r_i and the outer radius r_o; it is zero elsewhere. The total disk flux ${{F}_{\nu ,{\rm disk}}}$ is

$$\begin{eqnarray}&&{{F}_{\nu ,{\rm disk}}}=C\int _{{{r}_{i}}}^{{{r}_{o}}}{{B}_{\nu }}[{{T}_{g}}(r)]\ {{\left( \frac{\nu }{1\ {\rm Hz}} \right)}^{\beta }}{{\left( \frac{r}{1\ {\rm AU}} \right)}^{{\Gamma }+1}}dr.\end{eqnarray} \tag{ 3 }$$

C is the scaling factor such that $C=2{{\pi }^{2}}{{n}_{0}}a_{0}^{2}(1\ {\rm AU}){{d}^{-2}}$. The characteristic grain radius is denoted as a₀.

There are seven free parameters total in this model, five of which are needed to calculate the total disk flux: β, Γ, C, r_i, and r_o. Two additional parameters, inclination i and position angle ϕ, are needed to generate model images. We took several steps to initialize the starting parameters for the fit. First, we set β to zero as if the emission came from perfect blackbodies. Then, we set r_i and r_o to a range of a few pixels to several 100 AU. Afterward, values for ${\Gamma }$ and C were obtained by fitting the model to the SED data only. Finally, i and ϕ were initialized to the best-fit inclination and position angle when fitting a ring to the images only. Once the parameters are initialized, we used a Levenberg–Marquardt least-squares fitter. The fitting process tends to drive r_o to unbounded values. The reason that r_o tends to drift to unbounded values is due to the model's underprediction of the flux beyond 35 μm. During the fitting process, after the model settles with the parameters that fit the majority of the data, the model will try to extend r_o to an arbitrarily large value so that it will contain slightly more cold dust. The contribution from this cold dust does make the SED fit better, but the amount is negligible. Thus, we fix r_o to 400 AU, beyond which the image fit does not get better. The best-fit parameters, along with the total chi-square (${{\chi }^{2}}$) and the reduced chi-square ($\chi _{\nu }^{2}$), are listed in Table 3. The contributions to the total chi-square from the broadband photometry, IRS spectrum, and images are 3%, 31%, and 66%, respectively.

Table 3.  One-, Two-, and Three-population Model Fits

Parameter	A Continuous Disk	A Continuous Disk + A Ring	A Continuous Disk + Two Rings	Unit
${{\chi }^{2}}$	9920	6255	5958
$\chi _{\nu }^{2}$	1.03	0.65	0.62
Continuous Component
β	0.76	$1.53_{-0.01}^{+0.02}$	1.64 ± 0.02
C	3.7 × 10−31	$(5\pm 1)\times {{10}^{-42}}$	$(9\pm 2)\times {{10}^{-44}}$	(AU−1)
Γ	1.0	0.48 ± 0.05	$0.53_{-0.07}^{+0.06}$
i	49	73.8 ± 1.0	$74.5_{-1.0}^{+0.9}$	(deg)
ϕ	74	$62.0_{-1.1}^{+1.0}$	$61.2_{-0.9}^{+1.0}$	(deg)
ri	0.24	37± 3	35± 3	(AU)
ro	400	400	$310_{-20}^{+30}$	(AU)
Ring Component I
β		$0.27_{-0.04}^{+0.05}$	0.59 ± 0.02
Cr		$(7.4\pm 1.5)\times {{10}^{-15}}$	$(5.9_{-0.7}^{+0.8})\times {{10}^{-16}}$
r		$61_{-6}^{+8}$	105 ± 5	(AU)
Ring Component II
β			0.32 ± 0.06
Cr			$(4\pm 1)\times {{10}^{-14}}$
r			220 ± 40	(AU)

Note. The chi-square (${{\chi }^{2}}$) and the reduced chi-square ($\chi _{\nu }^{2}$) are the minimum values estimated by least-squares fitting. The one-component parameters are from the least-squares fitting; uncertainties are not estimated since the model does not describe this system well at all. The two- and three-component model parameters and uncertainties are quoted from the marginal distributions of MCMC results. For the two-component model, the parameters (from top to bottom) corresponding to the quoted least squares are (β, C, Γ, i, ϕ, ri, ro, β, C_r, r) = (1.53, $5.0\times {{10}^{-42}}$, 0.47, 74.0, 61.8, 37.5, 400, 0.27, $7.2\times {{10}^{-15}}$, 62). For the three-component model, the parameters corresponding (from top to bottom) to the quoted least squares are (β, C, Γ, i, ϕ, ri, ro, β, C_r, r, β, C_r, r) = (1.64, $1.1\times {{10}^{-43}}$, 0.48, 74.3, 61.1, 36, 316, 0.59, $6.2\times {{10}^{-16}}$, 105, 0.33, $4.3\times {{10}^{-14}}$, 232).

Download table as:  ASCIITypeset image

This one-component model does not produce a good fit. The modeled flux is too concentrated in the central region of the disk, as the residual images show strongly negative values in the center owing to oversubtraction. As a result, the best-fit i and ϕ are not trustworthy. In addition, modeled SED underpredicts the flux shorter than 16 μm and beyond 33 μm. The drive to better fit the short-wavelength part of the SED is inconsistent with the spatial location of flux in the images. This suggests that the power-law spatial distribution is inconsistent with the resolved outer disk flux and the contribution to the warm SED flux. Since this model clearly does not fit the data, we did not investigate further for estimating the uncertainties of the best-fit parameters. This inconsistency could be potentially resolved by considering a multicomponent model. We therefore move on to a more sophisticated disk model by adding a second component.

4.2.2. A Continuous Power-law Disk + A Ring

This two-grain-population model is composed of a continuous disk and a thin ring. The first component is as described in the section above. The second component assumes a single population of grains located in a narrow ring at a single radius r away from the star. We adopt the same emissivity law in the form of ${{\nu }^{\beta }}$ for the ring component. However, this β parameter could have a different value than in the continuous component since the optical properties might be different for different grain populations. The flux from a ring model is

$$\begin{eqnarray}&&{{F}_{\nu ,{\rm ring}}}(r)=\ {{C}_{r}}\ {{B}_{\nu }}[{{T}_{g}}(r)]\ {{\left( \frac{\nu }{1{\rm Hz}} \right)}^{\beta }}.\end{eqnarray} \tag{ 4 }$$

The grain temperature T_g is found by solving Equation (2), and C_r is a scaling factor for the ring model. There are five parameters in the ring component: β, C_r, i, ϕ, and r. Limited by the image resolution, we make the continuous disk and the ring sharing the same i and ϕ. Here, we again fix r_o to 400 AU, beyond which the image fit does not get better and the SED fit does not show any significant improvement. We first use least-squares fitting to find the best-fit model parameters. This model is a much better fit comparing to the previous single continuous disk model. The total ${{\chi }^{2}}$ decreases by more than 3000. The contributions to the total chi-square from the broadband photometry, IRS spectrum, and images are 3%, 14%, and 83%, respectively. The image fit improves significantly. The residual images do not suffer from the severe oversubtraction in the central regions anymore. Therefore, it is worth exploring the uncertainties of the fitting parameters for this model. However, the covariance matrix from the least-squares fitting does not always provide reasonable error estimations, especially for parameters that are degenerate. Degeneracies are prominent between β, ${\Gamma }$, C, and r_i for the continuous disk model and between β, C_r, and r for the ring model. For example, effects of temperature depression from increasing r can be compensated by having a larger value of β. To better quantify the uncertainties on the parameters, we use the ensemble MCMC method of Goodman & Weare (2010), as implemented in the emcee python package (Foreman-Mackey et al. 2013).

The emcee package is based on the Markov chain Monte Carlo (MCMC) method developed by Goodman & Weare (2010), where they utilize multiple "walkers" to propagate multiple chains simultaneously. We set the initial positions of the walkers by drawing random numbers based on Gaussian distributions with means equal to the best-fit parameters and the variances equal to 10% of the means. This MCMC algorithm adjusts the candidate step proposal distribution based on the walkers' current positions in the parameter space. We use 100 walkers and run for a few thousand steps after the burn-in period seems to be over. All of the marginalized probability density functions (PDFs) look fairly symmetric. We calculate the uncertainty on each parameter by measuring the 1σ interval on the marginalized PDFs, with the upper and lower bounds measured separately. The parameter values and their uncertainties are listed in Table 3.

The best fit shows that the continuous component contains hotter grains and is mainly contributing to the 11.7 μm emission while the ring component contains cooler grains and dominates the 18.3 μm emission. The hot continuous disk extends from a few tens of AU to hundreds of AU and contains grains as hot as 300 K. The cool ring locates around 60 AU, spatially overlapping with the continuous component. In comparison with the single-disk model, this two-grain-population model significantly improved the fit to the images and to the short-wavelength part of the SED. However, this model is still underpredicting the flux longer than 35 μm on the SED. The main reason for this underprediction is the outnumbered IRS spectral data points. The fit is driven by closely matching the spectrum in order to effectively lower the ${{\chi }^{2}}$. As a result, the model completely misses the photometry points at wavelengths longer than 35 μm. As a test, we considered excluding long-wavelength IRS points. We find that in order to produce a reasonably good fit with this two-component model, we have to cut out IRS points longer than ∼15 μm. Since we do not have a physical explanation to place such a cut, we include all IRS points in our analysis. To account for the long-wavelength emissions, we add another ring component in the following section.

4.2.3. A Continuous Power-law Disk + Two Rings

This three-grain-population model is composed of a continuous disk and two thin rings. The structures of these model components are described in the previous two sections. Although the inclination and the disk position angle are free parameters in the model, we constrain the three components to share the same values. As before, we use the least-squares fitting to find the best-fit model parameters and MCMC to quantify their uncertainties. The model parameters are listed in Table 3. The contributions to the total chi-square from the broad-band photometry, IRS spectrum, and images are 2%, 12%, and 86%, respectively. Figures 3 and 4 show the best-fit SED model and the image models, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This best-fit model is composed of a hot continuous power-law disk extended from 35 ± 3 AU to $310_{-20}^{+30}$ AU with temperature ranging from 330 to 150 K, a warm ring with temperature of 97 K located at 105 ± 5 AU, and the cold ring with temperature of 52 K located at 220 ± 40 AU. The middle panel of Figure 4 shows the PSF-convolved model disk in the two imaging wavelengths. Approximately 90% of the flux in the 11.7 μm model comes from the hot continuous component, while 10% comes from the warm ring. At 18.3 μm, about 2/3 of the flux comes from the warm ring and 1/3 of the flux comes from the hot disk. This model provides a reasonable fit to the images, despite that there are some small patches of oversubtraction and undersubtraction as seen from the image residuals (bottom panel in Figure 4). We characterize the variations in the resolution-element scales by dividing the sum of the smoothed residual flux values in a resolution-element-sized area by the sum of the corresponding smoothed noise. The variations in the resolution-element scales are within 1σ, suggesting that the model fits the images reasonably well. The cold ring component is too faint to affect either imaging channel. However, this cold ring is essential for contributing to the long-wavelength part of the SED and is a key reason why we can constrain r_o, unlike in the previous two models. We found that this three-component model is able to provide a sensible fit to the SED and the images simultaneously.

4.2.4. Other Disk Models

We also tried using a broken-power-law disk model. This model is composed of one continuous disk described by two power-law grain distributions: one between the inner radius and the intermediate radius and the other between the intermediate radius and the outer radius. The entire disk shares a single value of emissivity power index β, in which we assume that the grain composition is the same across the disk. The best-fit result improves slightly compared to the one from a single continuous disk model described by only one power-law grain distribution (Section 4.2.1). However, this best-fit broken-power-law model suffers from a similar inconsistency to the single-power-law model, where the model flux is too concentrated in the center of the disk in both bands. The reason lies under the assumption of a single grain population. Since the temperature of the entire disk is governed by a single emissivity power index, grains farther away will have lower temperatures. Thus, although a steeply rising density power law can make the model disk flux more extended, it will also drive the corresponding SED model too bright at longer wavelengths since the flux originated from low-temperature grains. In order to have warm dust grains be responsible for the extended emission instead, we must introduce a second emissivity power index to our model (as used in Section 4.2.2).

Another model we tried is composed of two continuous disks. We found that the two components can share the same values for ${\Gamma }$, r_i, r_o, i, and ϕ while still producing a good fit compared to having two independent sets of parameters. However, this model is very sensitive to the β values. A small adjustment in β will introduce a significant deviation in the resulting model. Thus, we must keep the two β parameters independent in order to maintain a good fit. This result again indicates that there is more than one population of grains in the system. We notice that one of the continuous components can be simplified into a ring without changing the best-fit result significantly ($\chi _{\nu }^{2}$ changes from 0.65 to 0.64). Therefore, we performed the detailed analysis of a continuous power-law disk plus a ring model (Section 4.2.2) instead of the two-continuous-component model.

Our models support the inference of multispatial components from a multitemperature SED. Recently, Kennedy & Wyatt (2014) argued that if the SED of a debris disk around an A-type star shows multiple temperature components, it is an indication that the system hosts multiple populations of grains in different locations. Our SED modeling of the disk around HD 131835 indicates that there are multiple temperatures. From modeling with the resolved images in mid-IR, we confirm that the system indeed has grains at multiple spatial locations: a hot continuous component extends from 35 ± 3 AU to $310_{-20}^{+30}$ AU, a warm ring located at 105 ± 5 AU, and a cold ring located at 220 ± 40 AU. Our model indicates that the two separated narrow rings are embedded in an extended disk component. Although not all the model components are completely spatially separated, we are confident that the dust is not concentrated in a single belt. Our modeling result agrees with the argument made by Kennedy & Wyatt (2014), adding to the small poll of confirmed cases with high-resolution observations.

Grains with effective temperatures hotter than blackbodies are responsible for the observed disk emission, since perfect blackbody grains at these spatial locations would not have the appropriate color temperatures. In our models, the grains are assumed to emit according to a modified blackbody function where the emissivity is proportional to ${{\nu }^{\beta }}$. In the case of the three-population model, the grains in the β values for the continuous component, the warm ring, and the cold ring are 1.64 ± 0.02, 0.59 ± 0.02, and 0.32 ± 0.06, respectively. β = 0 corresponds to blackbody grains. The positive β values indicate that the grains are small since they are inefficient in absorbing and emitting at long wavelengths, and β = 1 corresponds to small grains in the limit $2\pi a\ll \lambda$, where a is the grain size and λ is the observed wavelength. For β ≤ 1, real materials with complicated emissivities are required to explain the observations. Observational evidence points out that it is common for debris disks around A stars to have β values approximately within 2 (Booth et al. 2013), and our results fall within this region.

We searched for the potential grain compositions by first examining the equilibrium grain temperatures with our modified blackbody models and then compared them to the equilibrium temperatures of specific grain compositions. Using the model parameter r and β, we plotted the corresponding equilibrium temperatures at different stellar locations in Figure 5. Our goal is to use these temperature curves to identify the possible grain compositions. We consider some common compositions such as graphite, astro-silicate, and ice. We use the optical constants of graphite and silicate from Draine & Lee (1984) and Laor & Draine (1993). The optical property of ice is a function of its temperature owing to varying crystal structures. Warren (1984) provided the index of refraction of pure ice at −1°C, −5°C, −20°C, and −60°C. Although the optical properties are different at different crystal temperatures, the variation is quite small. Thus, to first order, we adopt the optical constants of ice assuming that the ice preserves a similar crystal structure to when it is at T = −60°C. We calculated their equilibrium temperatures using the Mie theory (Bohren & Huffman 1983) with the assumption that grains are spherical. The grain size contours are overplotted on the temperature–stellocentric distance map in Figure 5.

Zoom In Zoom Out Reset image size

Grains composed solely of graphite, astro-silicate, or ice in the equilibrium temperature cannot be used to explain all the observed disk properties. Among the three compositions, graphite grains have the highest temperature, and ice grains are the coolest. We have also considered using the other types of silicates such as olivine and pyroxene. However, their optical properties are similar to the astro-silicate, and thus the grain temperatures are on the same order of magnitude as the astro-silicate grains. The warm and cold rings could be made up by graphite and silicate grains closer to micron size. On the other hand, the hot continuous disk component would require the emitters to be really small. Since small grains become inefficient emitters at long wavelengths, these grains will sustain much higher temperatures than a blackbody. However, even though nanometer-sized graphite grains are significantly hotter than blackbodies, they still could not match the temperature of the hot continuous disk. In addition, grains this small are not likely to be present in the system as discussed below.

We have qualitatively considered having porous grains and using a broken emissivity power law in order to attempt to address the abnormally warm dust. However, neither consideration is in our favor. Porous dust grains have lower temperatures than compact spheres (Kirchschlager & Wolf 2013). Thus, introducing porosity will make the discrepancy worse for finding a physical grain composition for the observed dust temperatures. A similar qualitative result applies for using a broken emissivity power-law case. A broken emissivity power law assumes that the emissivity is 1 when $\nu \gt {{\nu }_{c}}$ and ${{({{\nu }_{c}}/\nu )}^{\beta }}$ when $\nu \lt {{\nu }_{c}}$, where ${{\nu }_{c}}$ is a free parameter indicating the critical frequency at which the emissivity function changes. Compared to the smooth emissivity law we used, the grains with the broken emissivity power law are less efficient in absorbing the stellar light. Therefore, grains with a broken emissivity power law will sustain lower temperatures. In addition to making the grains cooler, a broken emissivity power law introduces an additional model parameter, ${{\nu }_{c}}$, which makes the model more complex than a smooth emissivity power law.

Although the high disk temperature draws the connection to small grains, no features of small grains are detected in the system. Nanometer-sized grains are in the molecular region. For HD 131835, we do not observe any solid-state features in the 10 and 20 μm regions from the IRS spectrum, which indicates that the grains are probably larger and more likely to be in the micrometer-sized region (e.g., Draine & Lee 1984). Stochastically heating small grains such as PAHs could be one way to make grains have much higher temperatures than in the equilibrium states (Draine & Li 2001). However, the IRS spectrum does not show the spectral signatures of PAHs either.

Furthermore, small grains are unlikely to be present in the system when considering the radiation pressure from the star. The blowout grain size is calculated by balancing the radiation pressure with the gravity. Taking the mass of the star to be 1.9 M $_{}$ (Chen et al. 2012) and using the grain density of $\rho =3.5\;{\rm g}\;{\rm c}{{{\rm m}}^{-3}}$, we get a blowout grain radius for HD 131835 of

$$\begin{eqnarray}&&{{a}_{{\rm blowout}}}=0.91\ \mu {\rm m}\left( \frac{{{L}_{*}}}{10.5\ {{L}_{}}} \right)\left( \frac{1.9\ {{M}_{}}}{{{M}_{*}}} \right)\left( \frac{3.5\;{\rm g}\;{\rm c}{{{\rm m}}^{-3}}}{\rho } \right).\end{eqnarray} \tag{ 5 }$$

Grains smaller than this blowout size are unlikely to survive in the system. Although small grains can be trapped in resonance owing to gas drag while migrating, the effect is unlikely to be significant with the gas level present in this system. Since small grains are unlikely to be responsible for the hot continuous disk emission, the nature of these abnormally warm grains is not completely understood.

Identifying unique disk models around A stars can be challenging. For example, recently, with Herschel observations in 70, 100, and 160 μm, Booth et al. (2013) have shown that the disks around A stars have various morphologies, ranging from systems that can be fit with just a narrow ring to the ones that require wider or multiple rings. Although direct images can provide constraints on grain size distribution and dust properties, sometimes finding a proper model can be quite difficult. For example, a detailed modeling work on the β Leo debris disk with multiwavelength observations shows the degeneracy between one-, two-, and three-component models (Churcher et al. 2011). For HD 131835, the three-component model gives a reasonable fit, but more data are needed to set better constraints. We would like to better characterize the grain properties, understand the distribution of grains, and constrain the location of planetesimal belts in much greater detail.

Fortunately, with the new generations of high-resolution, high-sensitivity, and high-contrast instruments, detailed disk characterization is foreseeable. HD 131835 is located in the southern sky, making it a favorable target for GPI (Gemini Planet Imager), SPHERE (Spectro-Polarimetric High-contrast Exoplanet Research), and ALMA. GPI and SPHERE can potentially image the disk in scattered light in the near-IR, providing the currently uncharacterized scattering properties of the dust grains around HD 131835, in both total intensity and linearly polarized light. Designed specifically for high-contrast imaging, these instruments have the potential to probe dust-scattered light at small inner working angles (e.g., Perrin et al. 2014). With superior sensitivity and resolution, ALMA is capable of detailed mapping of the CO gas in this system. In addition, ALMA observations would trace the distribution of larger grains that are more tightly coupled to the locations of the large parent bodies. Future observations with these facilities hold great promise in further characterizing the dust distribution and dynamics.

HD 131835 shows strong IR excess, and here we present the discovery of the resolved debris disk in the mid-IR. The debris disk's properties can be constrained using all the available observations on HD 131835, including 15 photometry points, the IRS spectrum, and resolved images at 11.7 and 18.3 μm. From our modeling result, the disk clearly cannot be described by a single continuous population of modified blackbody grains. The images alone can be described by a two-grain-population model that is composed of a continuous power-law disk and a ring. The continuous component contains hotter grains and dominates the emission in the 11.7 μm image, whereas the ring component contains cooler grains and dominates the emission in the 18.3 μm image. However, in order to obtain a good fit to the SED simultaneously, an additional ring component is needed. This three-component model is composed of a hot continuous power-law disk, a warm ring, and a cold ring. In this model, the disk fluxes in the imaged wavelengths are contributed primarily from the hot continuous disk and the warm ring. Since the cold ring peaks at a longer wavelength, this third component does not show up in the mid-IR images, so its spatial location is relatively unconstrained. The excess emission in far-IR and submillimeter, however, can be well described by this third component. Starting with the simplest model, we found that a model with three components and therefore three grain populations can well describe the images and the SED simultaneously.

Work by L.-W.H. is supported by the National Science Foundation Graduate Research Fellowship number 2011116466 under Grant number DGE-1144087. P.G.K. and J.R.G. acknowledge support from NASA NNX11AD21G, NSF AST-0909188, and University of California LFRP-118057.