SPATIALLY RESOLVING THE HK Tau B EDGE-ON DISK FROM 1.2 TO 4.7 μm: A UNIQUE SCATTERED LIGHT DISK

HK Tau was observed in the thermal infrared on 2004 November 2 (UT) using NIRC2 behind the Laser Guide Star Adaptive Optics (LGS-AO) system (Wizinowich et al. 2006) on the W. M. Keck 10 m telescope. The observations were made while HK Tau was at an airmass between 1.16 and 1.09. The laser system was run at 6 W output power, which produced a guide star with an equivalent V-band magnitude of 9.9. The primary, HK Tau A, an R = 14.9 mag object 2''4 away from the science target, was used as the tip/tilt correction point source. The deformable mirror was run at 301 Hz.

NIRC2, a 1024× 1024 InSb array, was used with the smallest plate scale available (000993 ± 000005 pixel⁻¹; Ghez et al. 2004), providing a 102 square field of view. HK Tau was imaged using the L' (λ_o = 3.776 μm, Δλ = 0.700 μm) and Ms (λ_o = 4.670 μm, Δλ = 0.241 μm) filters. The largehex pupil mask was used throughout the observations. The HK Tau binary system was imaged in the L' filter with a 0.2 s integration co-added 100 times in three of the four quadrants of the array (avoiding one quadrant which suffers from an additional noise term that is difficult to calibrate). This observing cycle was repeated 27 times, providing a total integration time of 540 s. The same observing pattern was used for the Ms observations, using an integration time of 0.181 s and repeating the pattern 36 times for a total integration time of 651.6 s. The photometric standard HD 18881 (L' = 7.16 mag, M' = 7.17 mag; Leggett et al. 2003) was also observed during the run. The AO system generated a good correction at both wavelengths, with a median Strehl ratio of 0.71 and 0.75 for the L' and Ms data, respectively.

The data were reduced in the following manner. Individual dark current subtracted images, with the binary system located in different quadrants were subtracted from each other to provide a zeroth-order sky subtraction. These subtracted images were then flat fielded and any bad pixels interpolated over. NIRC2 suffers from a field distortion which has been measured. Using code made available by Keck, the cleaned images were distortion-corrected. The residual sky level in the frames is measured by taking the median value of the quadrants which do not have a star in the field. This value (∼200 DN in the L' images, corresponding to 0.05% of the initial sky value) is then subtracted from the cleaned images. The images were shifted using a cubic polynomial algorithm so that the primary is moved to the center of the field, and averaged together. Because the code to undistort the NIRC2 images does not accurately conserve flux, the reduced images described above are used only for astrometry and measuring the disk structure. Images reduced in the same manner, but without the distortion correction step, are used to provide photometry. The Kp, L', Ms three-color image is shown in Figure 1, using an LGS-AO observation of HK Tau at K' band (2.12 μm) taken during engineering time. The primary has an FWHM of 0081 and 0096 at 3.8 and 4.7 μm, respectively.

Zoom In Zoom Out Reset image size

Aperture photometry of the observations of HD 18881 were used to calculate the zero points of the filters. The observations, taken at airmasses of 1.25–1.55, were airmass-corrected using the median atmospheric extinction coefficients from UKIRT (0.093 and 0.22 mag airmass⁻¹ at L' and M, respectively). The encircled energy of the stars (flux as a function of radius) was measured and used to determine the aperture correction. The measured zero points are listed in Table 1 and are used to calibrate the photometry for HK Tau.

On 2002 November 18, we observed HK Tau in the near-infrared using the Nasmyth Adaptive Optics System (NAOS) and the CONICA instrument (Rousset et al. 2003; Lenzen et al. 2003) installed on the Yepun 8.2 m Unit Telescope at ESO's Very Large Telescope. The observations were part of the NAOS Guaranteed Time Observing program (program 70.C-0565). We used the 00133 pixel scale for the J and H imaging and the 00270 plate scale for the Ks imaging. The observations were conducted in good seeing conditions (06–07). HK Tau A was used as the adaptive optics guide star with the visible wave front sensor. The FWHM of the central core of the point-spread function (PSF) was measured to be roughly 0075 at all wavelengths.

We obtained deep images with the broadband J, H, and Ks filters, immediately followed by short, unsaturated exposures with the associated narrowband filters (centered at 1.24, 1.64, and 2.12 μm). Total integration times of 180 s, 300 s and 180 s were recorded at J, H, and Ks, respectively, split into 2–5 independent images. Total integration times with the short narrowband filters ranged from 14 s to 40 s. When only two images were recorded with a filter, they were subtracted from each other for background subtraction, cleaned for bad pixels and cosmic rays, and shift-and-added. When more images were available, a sky was estimated by medianing all images, which was then subtracted from each image prior to cosmetic cleaning and shift-and-adding.

In the deep broadband images, HK Tau A was heavily saturated, and we used the short narrowband images to replace the central part of the PSF after appropriate scaling. The short and long images were first aligned using the location of the brightest point in the HK Tau B nebula. The scaling factor between the two images was estimated as the median of the ratio of the images in the region where the non-linearity limit was not reached in the long image and the S/N was higher than 5 in the short image. The final images are shown in Figure 2, as a three-color JHKs image.

Zoom In Zoom Out Reset image size

HK Tau was observed with the Gemini North Near InfraRed Imager and spectrograph (NIRI; Hodapp et al. 2003) on 2006 January 7 in the K band and 2006 January 18 in the L band. NIRI has a 1024× 1024 ALADDIN array that we used in f/6 mode, providing a plate scale of 01162 pixel⁻¹.

For the K-band spectroscopic observations, the 2 pixel wide "center" slit was used with the K grism, producing wavelength coverage from 1.9 to 2.49 μm at R ∼ 1300. The observations were made in position angle (P.A.) mode with the P.A. of the array set to 169 deg; the long axis of the slit was therefore aligned along the binary axis to allow spectra of both stellar components to be obtained. Four frames, each with an integration time of 75 s, were taken in an ABBA nod pattern (nod separation of 6''), providing a total integration of 300 s on source. On 2006 January 18, the L-band spectra of HK Tau B were obtained using the 2 pixel wide "blue" slit, which produced spectra from 2.85 to 4 μm with a spectral resolution R ∼ 1100. The P.A. of the array was again set to 169 deg. Using the same ABBA nod pattern as before, with a single frame integration time of 1 minute, HK Tau was observed for 80 minutes. A telluric standard (HIP 26734) was observed immediately after the HK Tau observations with a total integration time of 40 s and 360 s, and an airmass difference of 0.11 and 0.07, at K and L bands, respectively.

The data were reduced using the IRAF GNIRS and NIRI reduction packages. Each frame was flat-fielded and subsequent frames in the nod-pattern subtracted from each other to provide sky-subtracted spectra. The positive and negative spectra for all the frames were then aligned and co-added. Spectra of Argon lamps were taken during the K-band observations and these were reduced in a similar manner and used to provide a wavelength and spatial calibration for the K-band observations. For the L-band data, the telluric lines in the atmospheric standard were used for wavelength calibration. Using the wavelength and spatial distortion solutions, the data were rectified and resampled to a uniform wavelength scale, and the spectra were extracted within an 8 pixel aperture, corresponding to 1.2*FWHM and 1.0*FWHM of the K and L spectral profiles, respectively. The telluric standard spectra were divided by a T = 9000 K blackbody continuum, producing an atmospheric absorption profile. In the K-band telluric standard, the 2.166 μm Br γ line is interpolated over. The HK Tau spectra then have the atmospheric profile divided out.

We did not attempt to use the telluric standard as a spectrophotometric standard, as slit losses and atmospheric conditions were uncertain. Instead, we renormalized the spectra of HK Tau A so that the integrated photometry in the broadband K and L filters matched the photometry reported by Kenyon & Hartmann (1995). The spectra of HK Tau B were corrected for the same scaling factor, i.e., the ratio of the components in both bands was preserved. While HK Tau A is variable in the near-infrared at the 0.2–0.3 mag level, we note that with our choice of normalization the spectrum of HK Tau A fits nicely to the overall spectral energy distribution (SED) of that component. The resulting final K- and L-band spectra are shown in Figure 3, overlaid on the near-infrared SED of the two sources. Only the 2.03–2.42 μm and 2.95–3.9 μm ranges of the spectra are shown, to avoid the sections of the data that are strongly affected by telluric absorption. The spectra were re-binned into 0.03 μm wide bins to improve the signal to noise. The rms of the measurement within each bin is used as an estimate of the uncertainty.

Zoom In Zoom Out Reset image size

3.1.1. Photometry and Astrometry for the HK Tau Binary System

Fluxes for both HK Tau A and B, as well as the relative separation and P.A. of the binary system, are calculated for each bandpass in the thermal infrared. In the L' and Ms bands, fluxes are calculated within 50 pixel and 10 pixel radii (05 and 01 radius apertures). The smaller aperture is necessary in the more background-dominated Ms-band data. The fluxes are airmass- and aperture-corrected, and the magnitudes are listed in Table 1. The uncertainties listed are the rms flux values from the 3(4) individual images that were combined to make the median at L' (Ms). These uncertainties therefore do not include any systematic effects from, for example, the aperture correction. The AO fluxes listed in Table 1 are generally consistent with those found for HK Tau A and B when the flux ratios are combined with the Rebull et al. (2010) IRAC 3.6 and 4.5 μm fluxes. The VLT data were not observed with obtaining absolute fluxes in mind; we provide the flux ratio of the systems components for the VLT J, H, and Ks filters in Table 1. The observed binary separation and P.A., calculated from centroids on the stellar positions, are recorded in Table 1.

The SED for HK Tau is shown in Figure 4, combining the photometric results found here with those from the literature. The spatially resolved photometry for HK Tau A (squares) and HK Tau B (circles) are displayed, in comparison to the fluxes for the unresolved system, HK Tau AB (crosses). The photometry has been taken from Kenyon & Hartmann (1995) at 0.44–4.7 μm, Duchêne et al. (1999) at 0.44 μm, Woitas et al. (2001), S98, and Moneti & Zinnecker (1991) at 1.2–2.2 μm, White & Ghez (2001) at 2.2–3.5 μm, Spitzer IRAC and MIPS data from Rebull et al. (2010) at 3.6–70 μm, Keck LWS data from McCabe et al. (2003), IRAS data at 12–100 μm fluxes at 350 μm and 1.3 mm from Andrews & Williams (2005) and 1.4 and 2.7 mm fluxes from D03. The Spitzer IRS spectrum of HK Tau A+B from 5 to 35 μm (Furlan et al. 2006) is also included in this figure. We used the reduced IRAC images from the Taurus Spitzer Survey (Güdel et al. 2007; Rebull et al. 2010) to see if we could recover resolved IRAC photometry for the two components through PSF fitting and subtraction. Even constraining the P.A. and separation of the components, the binary was too tight for a stable solution to be found. Even without the IRAC photometry, however, it is clear that the slope of the SED for HK Tau B follows that of a normal photosphere, although one with significantly reduced fluxes for its spectral type and distance and slightly bluer than that seen for HK Tau A. This is typical for edge-on disks seen in scattered light (Stapelfeldt et al. 1997). The SED appears to be dominated by scattered light out to 11.7 μm, a fact born out by the spatially resolved image of the disk at 11.7 μm (McCabe et al. 2003).

Zoom In Zoom Out Reset image size

3.1.2. 3 μm Spectra of HK Tau B: Non-detection of PAH Emission and a Variable Water-ice Feature

Given the low mass, low UV, nature of the system, we did not broach the possibility of the 11.7 μm resolved disk being partly due to polycyclic aromatic hydrocarbon (PAH) emission in the McCabe et al. (2003) paper. A combined HK Tau A+B 2–3 μm spectrum by Whittet et al. (1988) showed no feature, although the primary dominates this spectrum. The NIRI spectra, resolving the spectra from both primary and secondary (Figure 3), show that neither HK Tau A or B show any evidence of the 3.3 μm PAH emission feature. An independent L-band spectra of HK Tau B from Terada et al. (2007) also does not detect the 3 μm PAH feature toward the disk. This is not surprising; in a Spitzer sample of 38 T Tauri stars, some with spectral types as low as M6, the lowest mass object detected with PAH emission has a spectral type of G8 (Geers et al. 2006).

Figure 3 shows that both components have a Br γ emission feature at 2.17 μm. Both components also show Na (2.20 μm) and Ca (2.26 μm) in absorption and have CO band heads in absorption. However, the only feature present in the L-band spectra of HK Tau B is the broad 3 μm water-ice absorption band, confirming the detection by Terada et al. (2007). Note that this feature is not detected in HK Tau A, ruling out that we are observing water-ice absorption along the line of sight to the binary system. We have converted the flux calibrated spectra into an optical depth plot in Figure 5; we do not believe that the feature is an absorption profile, but rather a reflectance one (see Section 5.3) and we only use τ here as a means to quantify the feature depth and compare it to previous work (Terada et al. 2007). The spectra has been rebinned into 0.03 μm bins, and the source continuum is measured by a blackbody fit to the 2.15–2.25 and 3.8–3.9 μm windows, in order to provide a direct comparison with Terada et al.'s analysis. We find a maximum optical depth of τ ∼ 0.6, with an rms of 0.1 dex due to uncertainties in defining the continuum level. This is considerably different from the τ = 1.32 measured by Terada et al., and we argue that this feature is variable, just like that for HV Tau C (Terada et al. 2007).

Zoom In Zoom Out Reset image size

3.1.3. The Scattered Light Disk around HK Tau B

3.1.4. An Asymmetric Disk

As originally noted by S98, the HK Tau B disk is not symmetric. Analysis of the horizontal profile of the brightest nebula shows that there appears to be two components to the scattered light profile, a bright "central" peak and a fainter, more extended, component, which appear to have centers that are laterally offset from each other. To quantify this, the horizontal profile of the bright nebula is fitted simultaneously with two Gaussian components at each wavelength, each of which have independent centers and widths. For ease of discussion we call these two components the core (narrower, brighter peak) and halo (fainter, wider peak). The results of these two-component fits to the disk profile are seen in Figure 6. In Table 2 we report the FWHM of the core and halo components, the ratio of the peak strength of the two components and the displacement, in arcsecs, between the center of the core component and the center of the halo component, with a positive shift indicating that the broad halo component is displaced in the southwest direction along the disk midplane from the core component. Both core and halo components are seen to decrease in size with increasing wavelength, becoming "peakier," although the ratio between the peak values of the two components appears steady. The asymmetry of the disk shows up in the displacement between the central position of the two components. The largest measured asymmetry is in the visible and it becomes less asymmetric (smaller displacements) at longer wavelengths. At 0.6 μm, the disk has a low level wing centered on the southwest side of the core. By 3.8 μm the more extended wing is centered slightly to the northeast. A similar behavior is observed in the fainter nebula as well, in that it appears to be off-center relative to the bright nebula from 0.6 to 2.12 μm.

Zoom In Zoom Out Reset image size

Table 2. Scattered Light Profile Fitting

Wavelength	FWHM		Core/Halo	Halo–Core
(μm)	Core	Halo	Flux Ratio	Displacementa
	('')	('')		('')
0.606	0.26	0.60	2.05	0.08
0.814	0.24	0.55	1.45	0.06
1.24	0.19	0.51	1.77	0.04
1.64	0.15	0.43	1.77	0.02
2.12	0.15	0.41	2.22	0.00
3.78	0.11	0.30	1.31	−0.01

Note. ^aA positive displacement indicates that the center of the halo component is located in the +x direction from the center of the core component, where x is the axis along the midplane of the disk.

Download table as:  ASCIITypeset image

The observed scattered light disk clearly changes morphology as a function of wavelength, both in terms of disk shape and flux. Figure 7 shows a 256 AU square region centered on HK Tau B, imaged at 0.6 and 0.8 μm (S98) with WFPC2 on the Hubble Space Telescope (HST), 1.2, 1.6, 2.12, 3.8, and 4.7 μm (this work). Each image has been rotated so that the disk axis is horizontal. Not only do the two scattered light nebulae get closer together, the flux ratio between the nebulae also increases. We quantify this in Table 3, where the FWHM of the long axis of the brightest scattered light nebula, measured through a two-dimensional Gaussian fit, is recorded, along with the separation and peak-to-peak flux ratio of the two nebulae. The nebula separation is measured by taking a cut across the mid-disk, perpendicular to the long axis, and fitting it simultaneously with two Gaussians. Given that we are fitting the vertical profile through the disk nebulae, which are unresolved in this data set, we keep the width of the Gaussians fixed. The best-fitting peak positions are subtracted to calculate the nebula separation. In all instances where the counternebula is clearly detected, the peak-to-peak flux ratio is the flux at the centroid of the bright nebula divided by the peak flux at the same location in the counternebula. For the L' data where the counternebula is confused with the diffraction ring, we take the flux ratio estimated from the primary nebula model subtraction.

Zoom In Zoom Out Reset image size

At all wavelengths presented here, the measured FWHM is significantly larger than that expected for a point source. We find that the FWHM of the scattered light nebula declines by a factor of two over the 0.6–3.8 μm range of observed wavelengths. The separation between the two nebulae, which is directly proportional to the disk-mass–opacity product, also declines from 024 to 017. The flux ratio between the nebulae has the strongest dependence with wavelength; the flux ratio steadily increases throughout the optical and near-infrared, with a factor of three increase between 0.6 and 2.2 μm, and increasing by another factor of 2.5 between 2.2 and 3.8 μm.

We can compare the observed wavelength dependence in HK Tau B with that seen in other scattered light edge-on disks. Both HV Tau C (Duchêne et al. 2010) and HH 30 (Watson & Stapelfeldt 2004) have inclinations greater than 83° and have been spatially resolved at more than three separate wavelengths in the optical and near-infrared. All of the original data was obtained from the authors, and the peak-to-peak flux ratio, nebula distance, and FWHM calculated in exactly the same manner. The observed chromaticity of the dark lane in HK Tau B is typical (e.g., Duchêne et al. 2010; Watson & Stapelfeldt 2004), and is a response to the wavelength dependence of the dust opacity (Watson et al. 2007); on the other hand, the chromaticity of the nebula flux ratio is quite unusual. Figure 8 plots the measured peak-to-peak nebula flux ratio and FWHM for all three disks as a function of wavelength. HV Tau C and HH 30 show a nebula flux ratio that declines by a factor of two from 0.8 to 2.2 μm, in comparison to the increase by a factor of ∼3 seen in HK Tau B. To date, this is the only edge-on disk system that has an increasing flux ratio with wavelength. Even without the information at λ>2 μm, it is clear that this disk has an unusual scattering behavior. The FWHM of the HV Tau C disk is near constant with wavelength (Duchêne et al. 2010), unlike the significant reduction in FWHM that is seen in HK Tau B with wavelength. This is true even when the nebula is broken down into separate "core" and "halo" components; both components have an FWHM that becomes smaller as the wavelength of observation increases. Both the FWHM and the nebula flux ratio are affected by the dust phase function, which is frequently quantified by the dust asymmetry parameter, g = 〈cos  θ〉 = ∫_4πS(ω)cos θdω, where θ is the scattering angle and S is the phase function. The observed wavelength dependence in these two observables suggests that the dust in this disk is quite forward throwing, and remains so throughout the mid-infrared.

Zoom In Zoom Out Reset image size

In this section we describe a grid of radiative transfer models of the disk in which both the disk geometry and dust properties are varied in order to investigate whether we can recreate the observed wavelength dependence of scattered light properties. Given that the SED of the HK Tau B component is lacking any information between 11.7 μm and the resolved millimeter data, we do not attempt to include any flux constraints in our modeling efforts. In addition, the significant discrepancy between the previous best-fit disk model to the scattered light data (S98) and that needed to replicate the spatially resolved millimeter data (D03) highlights that a concerted effort to reconcile the scattered light images with the millimeter data is beyond the scope of this work. We therefore restrict our efforts to recreating the observed disk morphology from 0.8 to 3.7 μm as a first step toward a global modeling of this disk.

Within this range of wavelengths, the observed disk morphology could contain (1) stellar photons, (2) stellar photons that have scattered off the disk, (3) thermal emission from the disk, either in LTE, or non-LTE emission from small grains, and (4) thermal emission photons that are scattered off the disk. The first option, that stellar photons from the star may start to penetrate the disk and contribute to flux in the brightest nebula, could potentially explain both the decrease in FWHM and increase in nebula flux ratio; however, the SED of HK Tau B (Figure 4) rules this out—the flux of the system falls throughout these wavebands and is not consistent with the starlight penetrating the disk. One potential way stellar photons could affect the scattered light distribution is if the star were particularly active, and hot starspots were contributing to the light eventually being scattered off the disk. Clear evidence for starspots altering the observed scattered light distribution in a circumstellar disk is provided by HH 30 (e.g., Stapelfeldt et al. 1999; Watson & Stapelfeldt 2007). Unlike HH 30, HK Tau B is a fairly inactive T Tauri star (EW Hα = 12.5 Å; e.g., McCabe et al. 2006) and we also do not have any evidence for time variation in the observed scattered light disk. The presence of a large starspot, especially if located in the northern hemisphere of the star and on the stellar meridian, could add a "beaming" effect that would make the FWHM of the nebula narrower and make the dust look more forward throwing. Critically, however, the wavelength dependence of the such a beaming effect acts in the opposite sense of the behavior we observe. For example, a 10,000 K starspot on a 3800 K stellar surface will make the disk scattered light distribution look narrower at I band than it will at K band, whereas we observe the disk to become narrower and "peakier" at longer wavelengths. We therefore do not include the presence of starspots in our model. We explore the remaining options 2–4 in Section 4.2.1.

While radiative transfer modeling of edge-on disk systems has been shown to be limited in its ability to find a unique fit to the disk geometry (e.g., Watson & Stapelfeldt 2004; Watson et al. 2007), the disk inclination, mass–opacity product, and scale height can be well constrained. HK Tau B has been modeled before (e.g., S98) and rather than use an iterative convergence fitting routine to find a single model that best fit the available data, we run models over a uniformly sampled grid in parameter space and take a Bayesian approach to find the range of dust properties that can recreate the observed wavelength dependence of the scattered light morphology. Monochromatic models are run at the central wavelengths corresponding to the I, H, Ks, and L' filters. We quantify the scattered light morphology for both the models and the data in the same manner, using three of the disk observables measured in Table 3: the FWHM of the brightest nebula, the distance between the top and bottom nebulae, and the peak-to-peak nebula flux ratio. Simultaneous multi-wavelength chi-square estimates of the models are done and Bayesian probability distributions are used to constrain the grain sizes present in the disk.

Following previous models of the disk (e.g., S98; McCabe et al. 2003), we simulate an azimuthally symmetric, flared circumstellar disk with

$$\begin{eqnarray*} &&\rho (r,z) = \rho _{0}(r)\,\exp \left(- \frac{z^{2}}{2H(r)^{2}}\right). \end{eqnarray*}$$

The scale height of the disk varies with radius as a power law, $H(r) = H_{0}(\frac{r}{r_{0}})^{\beta }$ and the surface density of the disk follows $\Sigma (r) = \Sigma _{0}(\frac{r}{r_{0}})^{\alpha }$. Given the observed asymmetry of the disk, the assumption of an azimuthally symmetric system is not correct. However, comparing the data and models using a chi-square based on the disk observables, rather than through a pixel-by-pixel determination, will minimize this bias.

We use the MCFOST radiative transfer code (Pinte et al. 2006, 2009), which has successfully modeled a number of other disks (Pinte et al. 2008; Duchêne et al. 2010; Glauser et al. 2008; Bouy et al. 2008; Benisty et al. 2010). The model grid cells are logarithmically spaced in the radial direction and linearly in the vertical direction. The disk is cut-off vertically at 10 scale heights, has an inner edge set to 1 AU and outer radius of 105 AU, a value previously determined from modeling of the WFPC2 scattered light images (S98). The disk is populated with dust in a power-law grain size distribution, such that dN(a) ∝ a⁻ⁿda, with the composition of the dust fixed to the porous dust composition of Mathis & Whiffen (1989).

For each wavelength, a monochromatic Monte Carlo model was run to determine the scattered specific intensity for each cell in the disk model. Photon packets were emitted from the central point source in an isotropic fashion; when a packet intercepts a cell, MCFOST chooses if the packet is absorbed or scattered depending on the average absorption and scattering opacities which are calculated under the assumption of the Mie theory. Each packet is traced through multiple scattering events until it is either absorbed⁷ or escapes the disk. The final scattered intensity in each direction is saved for each cell in the disk model. Ray tracing, with rays originating from the observer and integrating along the line of sight to each part of the disk, is then done to generate the final scattered light image at any specified inclination. This combination of Monte Carlo and ray-tracing methods has the distinct advantage of producing a high, uniform S/N model image. The only contributing uncertainty is the noise from the Monte Carlo process; making sure that this part of the computation is run with enough photon packets so that all scattered paths are followed and that each cell has a well defined specific intensity is important, especially for disks that are flat (low β) or have steep density profiles.

4.2.1. Investigating the Effects of Including Thermal Emission from the Disk

4.2.2. The Scattered Light Model Grid

The best-fitting model to the observables has a reduced $\chi ^{2}_{IHKsL^{\prime }}$ = 7.8, and is shown in comparison to the data in Figure 9. The first column of Figure 9 displays the data, the second column the simultaneous IHKsL' fit. A wide range of $\chi ^{2}_{IHKsL^{\prime }}$ are produced in the grid, with a maximum $\chi ^{2}_{IHKsL^{\prime }}$ of 2.1×10¹¹; clearly we are probing a wide enough range of parameter space to produce models that do not remotely fit the data. In addition to the $\chi ^{2}_{IHKsL^{\prime }}$ fit, we also show the χ²_I, χ²_H, χ²_K, and χ²_L single bandpass best fits. While these single-wavelength fits do a better job at matching the data at that individual wavelength, none of the single-wavelength models can recreate the other wavelength data as satisfactorily as the multi-wavelength fit.

Zoom In Zoom Out Reset image size

To estimate the most likely disk and dust parameters we employ a Bayesian approach, wherein each model is assigned a relative probability of exp(χ² /2). We assert that the range of parameters investigated here covers the full range of possible values such that the sum of the probabilities for each parameter is unity. The probability distribution for each model parameter is shown in Figure 10, and is plotted for each individual wavelength as well as the combined IHKsL' model (black thick dashed line). The values for each parameter found in the best-fitting IHKsL' model are also plotted on the probability distributions (black squares).

Zoom In Zoom Out Reset image size

All of the disk parameters are well constrained; all parameters but β have a probability distribution peaking above 0.6. There is only a 0.2% probability that the disk has a dust mass larger than 5 × 10⁻⁶M_☉, with the probability distribution strongly peaking (P = 0.76) at a disk mass of 1.55 × 10⁻⁶M_☉. Assuming a gas:dust mass ratio of 100, this suggests a total disk mass of <5 × 10⁻⁴M_☉ (with a 99.8% probability). Large grains are present in the HK Tau B disk; there is only a 0.03% probability that the maximum grain size is 3 μm in size or smaller. However, it appears that the population of large grains is a small fraction of the overall dust population, with a 99.99% probability that the size distribution has a power law of − 4.2 or smaller. Both H₀ and i are strongly constrained by the models, with a scale height of 3.8 AU at a radius of 50 AU and an inclination of 85 deg. The probability distributions for α, β, and dust mass, all show a double peaked distribution; it is possible that more than one family of solutions are capable of finding a good fit.

The best matching IHKsL' model, for our assumed dust composition and power-law size distribution, indicates the presence of large grains. However, the probability distributions in Section 4.3 favor a very steep power-law index (99% probability of n ⩾ 4.2), one that is not commonly inferred for dust size distributions and suggests that we may be observing a more complex size distribution than we have modeled. In this case, the indirect measure of a_max is not useful, and we instead focus on the fundamental dust properties and compare those to theoretical dust models to investigate the dust evolutionary state.

Figure 11 shows the dust properties for the best-fit IHKsL' multi-wavelength models of HK Tau B; the left panel shows the dust asymmetry parameter, g (blue solid line), while the right panel shows the albedo (green solid line) and dust opacity (blue solid line). The best-fitting model does a good job at reproducing both the FWHM of the primary nebula and the flux ratio between the two nebulae at all four wavelengths, both of which are primarily determined by the dust scattering phase function so we focus on the dust asymmetry parameter.

Zoom In Zoom Out Reset image size

In the g(λ) panel of Figure 11, we also plot the single-wavelength best fit for the four individual wavelengths fitted in this paper (blue circles), as well as g(0.8 μm) found by S98 and g(11.7 μm) found by McCabe et al. (2003). We find a slightly different g(0.8 μm) than S98, primarily because we use a slightly different disk parameterization and a coarser sampling in parameters. Overall, the models find that the dust in the HK Tau B disk is highly forward-throwing (large g values) throughout the optical and near-infrared before gradually decreasing through the mid-infrared. The multi-wavelength fit clearly has to make some compromises to fit all four wavelengths and does not provide as good a fit as any single-wavelength fit. It does, however, predict a g∼ 0.5 at 11.7 μm, in the range of that found by McCabe et al. (2003), providing independent support for such a dust model.

How does the modeled g(λ) for HK Tau B compare to that of standard ISM dust models? Figure 11 shows the predicted g(λ) for both the Mathis & Whiffen (1989) ISM dust model (black dotted line), and the Weingartner & Draine (2001, hereafter WD01) R_V = 3.1 dust model (green dotted line). Both predict a steeply declining g(λ) throughout the optical and near-infrared, reaching isotropic scattering (g = 0) somewhere in the 3.5–7 μm range, depending on the model. The HK Tau B dust is clearly not ISM-like in nature. While most of the extended Taurus dark cloud has dust with an R_V ∼ 3, the denser regions of the cloud have R_V ⩾ 4.5 (Whittet et al. 2001). WD01's R_V = 5.5 model (green dashed line) includes a small amount of grain growth over their R_V = 3.1 model. The addition of slightly larger grains increases the g value substantially at λ>1.5  μm, and predicts a flat g of ∼ 0.3 from 1.6 to 4 μm; a trend that is in the right direction but still produces g values considerably smaller than those observed for HK Tau B.

Our understanding of ISM dust properties, primarily deduced from the observed extinction law, has evolved considerably as the wavelength coverage over which we have measured the extinction law increases. The original ISM dust models, such as the Mathis et al. (1977) model, or the Mathis & Whiffen (1989) model used here, were developed through modeling the UV–1 μm extinction law. Considerable work had been done to extend the ISM extinction law and associated dust models out to longer wavelengths, but the uncertainties, especially at λ>3.4 μm have historically been very large (Mathis & Whiffen 1989). One of the key results from the Spitzer Space Telescope has been the extension of the ISM extinction law out to 8–24 μm and the realization that it flattens considerably beyond 3 μm (e.g., Indebetouw et al. 2005; Flaherty et al. 2007; Roman-Zuniga et al. 2007). K. Pontoppidan (2008, private communication) has constructed a post-Spitzer ISM dust model, which we also plot for comparison in Figure 11 (dot-dashed black line). The post-Spitzer model predicts a much flatter g(λ) relation than either of the ISM models, but the predicted g values are significantly lower than those observed at all wavelengths presented here. To summarize, our current knowledge of the dust properties for the ISM and also for dust in dense regions of dark clouds do not produce large enough g values to reproduce the observed scattered light distribution in the HK Tau B disk, and we argue that the observed g(λ) points to grain growth within the disk itself.

How does the multi-wavelength dust model for HK Tau B compare to dust models from other T Tauri circumstellar disks? Wood et al. (2002) present a dust model for the HH 30 circumstellar disk found by fitting the SED of the system from 0.5 μm to 3 mm, taking into account the disk geometry constraints provided by HST scattered light images of the disk (Cotera et al. 2001). The best fit g(λ) for the Wood et al. model, which includes dust grains up to 50 μm in size, is represented by the black dashed line in the left panel of Figure 11. The g(λ) relation is clearly flatter than that predicted from ISM dust models, and has g values that are higher than those predicted from either the post-Spitzer dust model or WD01's R_V = 5.5 model. The gradient in g(λ) for HH 30 is very similar to that found for HK Tau B, however, the HK Tau B dust is overall much more forward throwing than that in HH 30, suggesting a larger amount of grain growth than in HH 30.

Further indications that the HK Tau B dust is evolving lie in the disk mass discrepancy. The 1.4 mm spatially resolved map of HK Tau B requires a dust mass of 2 × 10⁻⁴M_☉ (D03), two orders of magnitude larger than models need to recreate the observed optical scattered light (S98 and this paper). This led D03 to suggest that dust settling toward the disk midplane has occurred in this disk and that the scattered light study, tracing the τ = 1 outer surface of the disk, and the optically thin millimeter study, were actually tracing separate dust populations. Our multi-wavelength study of HK Tau B, extending the spatially resolved scattered light observations out to a factor of ∼ 5 larger wavelength, suggest only a 0.2% probability of a disk dust mass on order of 5 × 10⁻⁶ to1 × 10⁻⁵M_☉. The discrepancy between the disk mass expected from the scattered light images and that from the disks thermal emission remains, reinforcing the case for a settled disk. In addition, the steep power-law size distribution, n = 4.2, is difficult to justify from a physical standpoint. The widely assumed n = 3.5 slope is theoretically grounded in the Dohnanyi (1969) model of a steady-state collisional cascade, for instance. Instead, it is possible that such a steep power law results from our attempt to fit an overly simplistic model to a more complex reality, such as a non-simple power-law size distribution, or the presence of two or more distinct dust size distributions, or vertical dust settling, as Pinte et al. (2008) suggested in the case of the IM Lupi. Note that observational evidence for dust settling has been argued for a number of other T Tauri disks (e.g., Duchêne et al. 2004; Fedele et al. 2008). Further modeling will be necessary to explore these options for HK Tau B and investigate whether the mass discrepancy can be resolved with a single disk model.

The profile of the brightest scattered light nebula of HK Tau B shows clear asymmetry at most wavelengths, quantified in Table 2 as a shift in the center position between two fitted Gaussians. The data have been taken on three separate dates, spanning seven years, with the R and I-band data taken in 1997, the J, H, andKs data taken in 2002, and the L' and Ms data taken in 2004. The R- and I-band data show a shift in the asymmetry even though they are taken on the same day. The J, H, and Ks data show this too. Therefore the asymmetry appears to be primarily a function of wavelength, although the lack of time coverage means that we cannot rule out time variability. The asymmetry decreases in magnitude toward longer wavelengths; if the disk contains a local overdensity a few AU out in the disk that causes a shadow (at a large enough radii to not cause any observed time variability over the seven years we have been observing this disk), then the wavelength dependence of both the opacity and phase function could reduce the size of the observed asymmetry with increasing wavelength.

HK Tau B is not the first asymmetric disk to be identified. HH 30, also an edge-on disk that has been spatially resolved in scattered light, shows clear brightness variations over an 11 year time period (Watson & Stapelfeldt 2007; Cotera et al. 2001; Stapelfeldt et al. 1999; Wood & Whitney 1998; Burrows et al. 1996). With evidence for both static and variable asymmetric components, Watson & Stapelfeldt (2007) suggest that there is an inclined magnetic dipole which may be affecting the inner edge of the disk and the accretion flow. In addition to HH 30, LkHa 263C (C. McCabe et al. 2012, in preparation) shows a strong lateral brightness asymmetry, as does HV Tau C (Stapelfeldt et al. 2003; Duchêne et al. 2010).

While we do not currently have any evidence for time variability in the observed disk asymmetry, there does appear to be time variability in the strength of the water-ice feature seen at 3 μm. Terada et al. (2007) obtained 3 μm spectra of HK Tau B in 2002; 4 years prior to our observation. They find a water-ice feature with a peak optical depth of 1.32, significantly different from the 0.6 ± 0.1 that we detect. Terada et al. (2007) also presented two independent observations of the same feature for the disk around HV Tau C; variability in the depth of the feature is also seen for this edge-on disk.

Terada et al. (2007) take the detection of the water-ice feature at 3 μm as an absorption feature from the outer edges of the protoplanetary disk, external to the observed scattered light region (see their Figure 3), and use the optical depth to calculate a column depth. However, the optical constants for water ice produce dips in the albedo throughout the 2–10 μm region, with the deepest albedo reduction at 3 μm (e.g., Warren 1984). Reduced albedos due to water features have been detected in disks before; L'-band coronagraphic observations of HD 142527 in scattered light show a dip that Honda et al. (2009) ascribe to a reduction in albedo. HD 142527 is in a near face-on orientation so the possibility of absorption from within the disk is excluded. Here we argue that detecting scattered light water albedo features is feasible for edge-on disks as well. By definition, scattering in an optically thick medium such as protoplanetary disks predominantly occurs where τ = 1, independent of grain opacity. Thus, spectral features in the grain opacity may not be visible in a scattered light spectrum. The possibility that the 3 μm ice feature in HK Tau B, and other disks, is albedo related and not absorption has to be strongly considered.

Plotted in Figure 5 is the optical depth associated with the albedo for pure water-ice grains (optical constants from Warren 1984). The three lines correspond to single-size grains of 3 μm (dotted line), 10 μm (solid line), and 30 μm (dashed line) in size. No fits between the model and the data are done, but 10 μm sized water-ice grains produces an optical depth feature that adequately matches the observed one. We therefore suggest that rather than being an absorption feature, the dip observed in the 3 μm spectra of HK Tau B is actually caused by a dip in the dust grain albedo from the scattered light from the disk.

If the observed feature is caused solely by an albedo effect, the variability in the depth suggests that the water-ice grains dominating the scattering have an effective radius that varies considerably on the timescale of a few years.

The search for grain growth in disks has come a long way since Beckwith et al. (1990). Work to measure the millimeter spectral index continues (e.g., Rodmann et al. 2006; Ricci et al. 2010) in addition to other complementary methods, such as SED fitting (D'Alessio et al. 1999, 2001; Chiang et al. 2001), scattered light imaging (e.g., S98; McCabe et al. 2002; Duchêne et al. 2004, Watson & Stapelfeldt 2004), and fitting of the 10 μm silicate feature (e.g., Kessler-Silacci et al. 2006; Sargent et al. 2009). Each method has its own limitations, but evidence for grains larger than those found in the ISM being present in T Tauri disks has grown stronger. At this point, we are beginning to see not just evidence for grain growth, but also for diversity in the dust properties (e.g., Sicilia-Aguilar et al. 2007; Sargent et al. 2009).

Scattered light imaging, probing the τ = 1 surface layers of the disks rather than the disk midplane, also seems to show a diversity of dust properties. While the disk around CB 26 needs grains smaller than 2.5 μm in size, with an unusually graphite-rich dust composition (Sauter et al. 2009), objects like HH 30 (Watson & Stapelfeldt 2004) and IRAS 04158+2805 (Glauser et al. 2008) show moderate amounts of grain growth, whereas HK Tau B provides evidence for grains larger than 5 μm. However, it should be stressed that the sensitivity of the scattered light disk morphology to the size of the scatterer, the wavelength of observation, and the disk geometry, can make comparisons between different disks complicated. The differing disk geometries can mean that scattered light observations at the same wavelength can probe substantially different regions in different disks (i.e., extended envelopes versus surface layers, versus layers closer to the midplane). Scattered light observations of the GG Tau circumbinary disk, showing no evidence for grain growth at λ < 2 μm finds strong evidence for large grains at 3.8 μm (McCabe et al. 2002, Duchêne et al. 2004). It should also be noted that single-wavelength models cannot constrain the dust properties as well as simultaneous fits to multiple wavelengths (Glauser et al. 2008; Pinte et al. 2008; Duchêne et al. 2010).

Given these rather extensive caveats, there are three T Tauri systems that have been observed in scattered light over a similar range of wavelengths, and which have a similar spectral type: HV Tau C (Duchêne et al. 2010), GG Tau (Duchêne et al. 2004), and HK Tau B, M0, K7/M0.5 and M2 type stars, respectively. All three disks have undergone multi-wavelength scattered light modeling and the best-fitting g values for each bandpass are shown in Figure 11 (GG Tau values from Duchêne et al. 2004; HV Tau C values from Duchêne et al. 2010). How do these disks compare?

All three disks show significantly different g values in the near- and mid-infrared. Both HK Tau B and HV Tau C have dust that appears to be significantly evolved from the ISM or dark cloud material. The HV Tau C dust appears to be similar to that predicted for the HH 30 disk (Wood et al. 2002), whereas HK Tau B has larger g values (more evolved dust?) than either of those disks. While the g values observed in the GG Tau circumbinary disk are similar to those predicted for dark cloud material with an R_V of 5.5, at least at wavelengths shorter than 2 μm, at 3.8 μm the best-fit g is significantly higher than that predicted for any of the ISM models, and Duchêne et al. (2003) argue for grain growth in this disk. Disk geometry does play a role here; while HV Tau C and HK Tau B are both nearly edge-on disk systems, GG Tau has an inclination of 43° and as Duchêne et al. (2004) point out, the observations at 3.8 μm appear to be coming from a layer of the disk closer to the disk midplane than those at shorter wavelengths.

To summarize, the three T Tauri disks highlighted here, all of a similar age and around similar type stars, have significantly different dust properties, showing a range of evolutionary states. The dust in the surface layer of the GG Tau disk is fairly similar to dust found in dense dark clouds, whereas the elevated g values found for HV Tau C and HK Tau B compared to the models argue for dust evolution from that state.