THE 0.5–2.22 μm SCATTERED LIGHT SPECTRUM OF THE DISK AROUND TW Hya: DETECTION OF A PARTIALLY FILLED DISK GAP AT 80 AU^*

The surface brightness profiles of TW Hya reveal structure on top of a smooth decrease in surface brightness with radial distance; this structure deviates from what one would expect from a simple flared disk. We investigate the radial surface brightness profiles in our medium band NICMOS data and compare the behavior in these wavelengths to that seen in the shorter wavelength data.

Figure 2 shows the F171M through F222M azimuthally averaged surface brightness profiles. In all of them, we see the characteristic shift in behavior between ∼80 and 130 AU that is seen in the visible, 1.1, and 1.6 μm, namely, a change in the slope of the surface brightness. This feature is now seen in all the wavelengths of light in which TW Hya is observed, indicating that this is a feature caused by some change in morphology rather than a compositional feature localized in wavelength. We overplot simple flared disk model surface brightness profiles consistent with those used to describe the TW Hya disk in the submillimeter (Andrews et al. 2012), which clearly do not match the behavior of the surface brightness profiles for TW Hya.

Zoom In Zoom Out Reset image size

To highlight the structure in the disk, we scaled each pixel of the disk images by R², where R is the distance from the central star. Figure 3 shows the STIS, F160W, F171M, and F222M scaled images, which clearly shows a depression in the disk structure at about 80 AU, coincident with the slope changes in the surface brightness profiles. The structure is most pronounced in the higher spatial resolution STIS images, but is still visible at longer wavelengths. Additionally, the STIS image shows a possible arc structure exterior to the gap at a position angle (P.A.) of ∼270°.

Zoom In Zoom Out Reset image size

The depression, or deficit, of surface brightness relative to the surrounding disk material can be caused by several factors, including a drastic change in dust opacity at ∼80 AU, a sudden change in turbulence of the disk, shadowing due to structure just interior to 80 AU, or the presence of non-axisymmetric structures like spiral arms (Weinberger et al. 1999; Clampin et al. 2003; Fukagawa et al. 2004). It could also be caused by a protoplanet accreting material and opening a gap structure within the disk (Bate et al. 2003; Bryden et al. 1999). Current planet formation models have difficulty forming large planets far from the central star. However, examples of such systems, like HR 8799 (Marois et al. 2008), suggest that large bodies can form around stars at distances of several tens of AU.

In Section 4.3, we investigate how a physical gap in material might change the observed surface profiles and use models to determine limits on the depth of the gap as constrained from our surface brightness profiles as a function of wavelength. From these limits, we make a preliminary estimate to the limiting mass of a protoplanet capable of opening such a gap in the TW Hya disk.

Roberge et al. (2005) measured the surface brightness as a function of azimuth for the STIS image of TW Hya's disk and found it possessed a significant asymmetry between 65 and 140 AU consistent with the asymmetry caused by a disk inclined to the line of sight with forward scattering dust grains. The magnitude and position of the asymmetry corresponded to a measured P.A. of the brightness maximum of 2336 ± 4° and an inferred inclination of 15$^{+8.7\circ }_{-6.4}$ under the assumption of a Henyey–Greenstein scattering phase function for the grains with an asymmetry parameter g of 0.5. They found no significant asymmetry for larger radii. Recent CO line data place the P.A. of the disk at 150°–155° (Andrews et al. 2012; Rosenfeld et al. 2012), consistent with a P.A. for the side of the disk closest to Earth being at either 60° or 240°. Rosenfeld et al. (2012) determined that they could reproduce this apparent brightness asymmetry with a warped inner disk at 4 AU with an inclination of 8°, rather than 6°–7° as inferred by CO line data for the outer disk. This particular configuration also helps to explain larger than expected CO line velocities in the inner disk regions.

In this paper, we re-analyze the STIS image as well as the other NICMOS NIR images to search for similar asymmetries. For each passband, we took the final subtracted images and constructed photometric apertures in the shape of concentric annuli as a function of increasing radius with equal spacings of 0357 (20 AU), to improve signal-to-noise ratio (S/N). Each annulus was then further equally subdivided into azimuthal bins with width 20°. For each image a total of six annuli were used, from 088 (47 AU) to 25 (134 AU). For the STIS and F110W images, there was enough signal to noise to extend out to 44 (236 AU). We therefore had a measure of the surface brightness distribution between 50 and 130 AU for all passbands, and extending to 230 AU for STIS and F110W.

In the NICMOS images, certain areas were completely masked by either the coronagraphic spot or diffraction spikes. For the STIS image, masked areas included diffraction spikes and the coronagraphic wedge. For apertures where masked pixels constituted more than 2/3 of the aperture, no surface brightness was calculated. Each annulus was then scaled to the median of the brightest annulus to provide equal weight to the surface brightness in any one azimuthal bin at each radius. The final azimuthal surface brightness distribution was then constructed by taking a median of the first six annuli for all images and another distribution was constructed for the STIS and F110W image of the disk from the final four annuli. For the uncertainty at each point, we took the quadrature sum of calculated photometric uncertainty and the standard deviation of the azimuthal points at that angular position.

Figure 4 shows the STIS and F110W results. Between 50 and 130 AU, the asymmetry is similar in the two passbands. Beyond 130 AU, the azimuthal profile is of a lower signal to noise but roughly matches what is seen at shorter distances, contrary to what was reported in Roberge et al. (2005). The distance of 130 AU corresponds to a break in the radial surface brightness of the disk. For a comparison, we also overplot our best-fitting scattered light models of the disk, described in more detail in Section 4.5, assuming an inclination of 7°, with Henyey–Greenstein g parameters of 0 and 0.5. Even with no forward scattering, flared disks can present a brightness asymmetry when circular apertures are used due to a "foreshortening" effect. Even though the side of the disk closest to Earth is slightly brighter due to the line-of-sight angle to the observer from the disk surface, this occurs at a smaller stellocentric angular separation than that expected for a non-flared disk. In this scenario, the side of the TW Hya disk closest to the observer would be to the NE. For STIS and F110W, a model with moderate forward scattering can also reproduce the azimuthal brightness asymmetry, if the semi-minor axis of the TW Hya disk facing the observer is toward the SW.

Zoom In Zoom Out Reset image size

Figure 5 shows the F160W, F171M, F180M, and F222M azimuthal brightness profiles with the same models overplotted. We neglect the F204M data because they have lower signal to noise. For each of these four images, no significant azimuthal asymmetry is revealed. Whatever is the cause of the asymmetry at shorter wavelengths, it is not detected with significance in the longer wavelength observations. We discuss the potential causes of this in Sections 4.5 and 6.

Zoom In Zoom Out Reset image size

The scattered light spectrum of the disk around TW Hya is a combination of the input stellar spectrum and the intrinsic reflectivity of the disk. We removed the stellar spectrum of TW Hya by dividing the measured surface brightness by the flux density of the star in each filter. We divided the STIS spectrum of the disk by the point source spectrum of the central star.

TW Hya did not have direct photometric measurements in the broadband HST filters (F110W, F160W, and STIS CCD) and has a complicated SED in both the optical and NIR (see Section 3.4). To find TW Hya's broadband flux densities, we bootstrapped from its flux relative to the PSF reference star. Although the PSF does not have exactly the same spectrum as TW Hya, it is anchored to TW Hya by literature photometric measurements of both stars at V, J, and H bands and our measurement at F171M. To interpolate to STIS, F110W, and F160W, we fit stellar atmosphere models of a range of effective temperatures and gravities to the broadband photometry of the PSF, scaled to TW Hya, and propagated the systematic model uncertainty into the TW Hya photometry. We also compared our flux densities in the NIR to the spectrum of Vacca & Sandell (2011), and found general agreement within the uncertainties.

Figure 6 shows the total reflectance spectrum of the disk from 0.5 to 2.22 μm averaged over radial distances of 50–215 AU, the extent to which we detected the disk in the medium band filters. Also shown is the visible light spectrum of the disk over the same distances, normalized to the STIS photometric point, from Roberge et al. (2005). The overall spectrum is relatively neutral between 0.5 and 1.6 μm, becoming slightly blue at longer wavelengths. No absorption bands are seen within this broadband spectrum.

Zoom In Zoom Out Reset image size

Most previous determinations of TW Hya's spectral type, and hence its age and mass, have been based on optical spectral diagnostics (Webb et al. 1999; Alencar & Batalha 2002; Yang et al. 2005). The consensus value of the optical spectral type is K7 or an effective temperature of ∼4000 K and a mass of 0.5–0.8 M_☉. Recent NIR spectroscopy of TW Hya, however, implied a much later spectral type of M2.5 (Vacca & Sandell 2011), which corresponds to a lower mass and younger age.

We investigated this discrepancy by finding our own spectral type for TW Hya using the broad wavelength coverage of our G750L spectrum (5240–10266 Å) taken 2002 July 17. We compared our spectrum to those of K7–M2 stars in the STIS Next Generation Spectral Library (Heap & Lindler 2009, 2007; Gregg et al. 2004, 2006), whose spectra also incorporate STIS G750L data. We calculate a reduced $\chi _\nu ^2$ value from a comparison of the spectra at each wavelength with each of the six comparison stars listed in Table 2 as models, while excluding data in the H and Ca emission lines. The reduced χ_ν value was determined against the 961 wavelengths from the spectrum, for a total of 960 (958) degrees of freedom for a one- (two-)component model, respectively. For each model star, we determined a mean T_eff and uncertainty from the literature, computed over the range of published T_eff values. The closest match ($\chi _\nu ^2$ = 105) was to GJ 825, an M0V dwarf with an effective temperature of 3730 ± 160 K (Figure 7). The fit is rather poor at both wavelength extremes—TW Hya is dimmer in the blue and brighter in the red.

Zoom In Zoom Out Reset image size

Vacca & Sandell (2011) suggested that TW Hya is in fact a later type star with a hot spot. An alternative explanation is that TW Hya is a hotter star with substantial cool spot, such as the one imaged near the magnetic poles and seen in radial velocity measurements (Huélamo et al. 2008; Donati et al. 2011; Boisse et al. 2012). In order to test these hypotheses, we further fit our spectrum with combinations of two comparison stars of different spectral types. A significantly better ($\chi _\nu ^2$ = 35), though not perfect, fit is obtained with a combination of a K7 (45% of the flux) and M2 (55% of the flux) star. Our K7 comparison was HR 8086 (T_eff = 3990 ± 115 K), and our M2 comparison was HD 95735 (T_eff = 3600 ± 130 K). This combination provides a good fit to both the optical and NIR spectra of TW Hya. One limitation to our method is the lack of multiple comparisons for a finer temperature constraint, as well as disagreement on the T_eff of our comparison stars in the literature.

We find the effective radius of TW Hya by combining the known distances to the comparison stars and to TW Hya with interferometrically measured (van Belle & von Braun 2009) or model-inferred (Takeda et al. 2007) stellar radii of the comparison stars. We find a radius of 0.84 ± 0.2 R_☉ for the M0 comparison and 1.08 ± 0.15 R_☉ for the K7+M2 comparison, roughly in agreement with the findings of Vacca & Sandell (2011). We assume that the uncertainties in the radius are driven by the uncertainty in the parallax.

This complex situation involving multiple temperature components makes an inference of a mass, T_eff, and luminosity for TW Hya problematic. However, a comparison of our inferred temperature and radius with evolutionary models provides some insight into whether we have a low-mass star with accretion hot spots and a higher mass star with large cool spots. We overplot the inferred T_eff and radii from our spectral templates in Figure 8 with respect to isochrones for stellar masses of 0.4–0.8 M_☉ from Baraffe et al. (1998; inspection of other isochrones gives similar results within the uncertainties). If the coolest component represents the central star, then the inferred radius and T_eff are consistent with a 0.55 ± 0.15 M_☉ object at an age of ∼8 Myr, in line with other age and mass estimates of TW Hya. If the hot component dominates, however, the age would be closer to 20 Myr. Similarly, the intermediate case of a slightly higher T_eff object comparable to an M0 spectral type implies an older star (t ∼ 20–30 Myr).

Zoom In Zoom Out Reset image size

We conclude that the best interpretation of stellar spectrum is of a cooler star (T_eff ∼ 3600) with a mass of 0.55 ± 0.15 M_☉ and an age of 8 ± 4 Myr. There exists, therefore, a significantly warmed atmosphere due to accretion over portions of the star's surface. An estimation of the accretion luminosity can come from assuming the K7 component as the SED of the warmed photosphere due to accretion (additional emission may arise at shorter wavelengths due to the hotter shock itself). Under the assumption of blackbody emission, if we subtract off the underlying luminosity of the M2V component, the luminosity coming from accretion is estimated to be 0.03 L_☉, corresponding to an areal covering fraction of 10%. Assuming $L=GM\dot{M}/R$, this implies accretion rates for TW Hya of 2 × 10⁻⁹ M_☉ yr⁻¹, in line with other estimates of TW Hya's accretion rate (Batalha et al. 2002), which range from 4 × 10⁻¹⁰ to 5 × 10⁻⁹ M_☉ yr⁻¹.

Submillimeter measurements of the Keplerian rotation of gas in the TW Hya disk in principle provide an independent measure of the stellar mass. Given the near face-on inclination for the disk, even very precise measurements of Keplerian rotation in TW Hya provide a similar constraint to mass as to what we have obtained through our STIS spectroscopy: while CO measurements of TW Hya's disk constrain the inclination to 6°–7° and accuracies to within a degree, that still allows for a wide range of stellar mass, since the quantity M^1/2sin i is measured. If the inferred inclination is 6° and the mass of TW Hya is best fit by 0.6 M_☉ (Hughes et al. 2011), an uncertainty of 1° would correspond to masses that ranged from 0.85 M_☉ to 0.44 M_☉ (Qi et al. 2004; Hughes et al. 2011) Most of the prior SED and submillimeter spectral cube modeling have assumed M_* = 0.6 M_☉ (Calvet et al. 2002; Qi et al. 2008), while other investigators have assumed masses as high as M_* = 0.8 M_☉ (Donati et al. 2011; Andrews et al. 2012). In Section 4, we select our best-fitting mass of 0.55 M_☉ and determine that mass does not significantly impact most parameters of interest within the disk.

The disk models are generated as described in Jang-Condell & Turner (2012). The details of the radiative transfer modeling are also described in Jang-Condell & Sasselov (2003, 2004) and Jang-Condell (2008). Although we assume that the overall disk structure is axisymmetric, the calculation of the radiative transfer is done in 3D to include the curvature of the disk both in vertical height above the midplane (z) and in the azimuthal direction (ϕ).

The structure of the unperturbed planetless disk is generated in a two-step procedure. In step 1, we calculate a locally plane-parallel two-dimensional model for the entire disk. We use the same formalism developed by Calvet et al. (1991) and D'Alessio et al. (1998, 1999), with some simplifying assumptions. We calculate the disk from 0.25 to 256 AU, with radial bins spaced by factors of $\sqrt{2}$. At each radius, we assume that the disk is locally plane parallel. This quickly generates an estimate for the radial and vertical temperature of the disk and the surface density profile that we can then refine in step 2.

Then, in step 2, we remove the assumption of a local vertical plane-parallel surface, but keep the vertically integrated surface density profile fixed. We take a radially limited slice of this disk and refine its structure, this time taking the full three-dimensional (3D) curvature of the disk into account. We calculate radiative transfer in 3D as a numerical integration over the surface of the disk.

In both steps, we iteratively calculate the vertical density and temperature structure of the disk including radiative transfer and under the assumption of hydrostatic equilibrium. The main heating sources are stellar irradiation and viscous dissipation. The opacities used are the same as those used for calculating the disk brightness, although to speed up the calculations, the opacities are averaged over the stellar and disk thermal emission. To represent the optical depth of the disk to stellar light, we average the wavelength-dependent opacity over the Planck spectrum evaluated at T_eff and use the Rosseland mean opacity to represent the optical depth of the disk to its own radiation.

We base our disk structure modeling on our inferred stellar parameters from Section 3.4. We assume a stellar mass of 0.55 M_☉, radius of 1.08 R_☉, and total luminosity of 0.208 L_☉. The spectrum of the star is composed of two blackbodies, one at 3990 K and the other at 3600 K, contributing 55% and 45%, respectively, to the total luminosity. The opacities of the dust to stellar irradiation are calculated using this spectrum.

Given the uncertainties in mass and T_eff described in Section 3.4, we consider their impact on our models. Decreasing the stellar mass while keeping the total luminosity of the star fixed has the effect of making the disk puffier, due to the lower gravity. The disk is then able to intercept more light, making it appear brighter. Disk models differing only by the stellar parameters have the same basic shape and spectral behavior but change overall brightness by 5%–10%. Changing the stellar parameters also mainly impacts the inferred radial width of the gap, but not the depth.

The disk parameters are those derived by Calvet et al. (2002) from the SED of TW Hya. For the accretion rate, we assume a value of 1 × 10⁻⁹ M_☉ yr⁻¹ (see Section 3.4). We assume that the dust is well mixed with the gas and constant throughout the disk, with a dust-to-gas ratio of 0.0138 and a grain-size distribution of n(a)∝a^−3.5. Calvet et al. (2002) and Wilner et al. (2005) found that large maximum grain sizes (a_max) give good fits to the SED of TW Hya at long wavelengths (λ ≳ 1 cm), but there is a degeneracy between maximum grain size and viscosity parameter (α) in their models, where the viscosity is given by $\nu = \alpha c_s^2/\Omega _{\phi }$, where c_s is the thermal sound speed and Ω_ϕ is the Keplerian orbital angular velocity. Good fits to the SED of TW Hya are obtained with a_max = 1 mm, 1 cm, or 10 cm. This is because large grains are effectively invisible to optical and infrared wavelengths. The α parameter must be adjusted as the maximum grain size changes, because the overall mass of the disk increases with increasing grain size, assuming a fixed gas-to-dust ratio; however, the stellar mass accretion rate must be kept fixed.

Given that a_max, α, and $\dot{M}$ all have the primary effect of scaling the disk mass up or down and that the degeneracies inherently give limited constraints on these parameters, we assume for simplicity in all our disk models that $\dot{M} = 10^{-9}\, M_{\odot }\,\mbox{yr}^{-1}$, a_max = 1 cm, and α = 5 × 10⁻⁴. The effect of increasing the accretion rate would be to increase the total surface density of the disk, which could be offset by adjusting the viscosity parameter α downward. Our assumed parameters give a total disk mass within the simulation boundaries between 27 and 211 AU of 0.074 M_☉. This assumed disk mass is comparable to estimates of the disk mass from Herschel HD observations (Bergin et al. 2013).

The wavelength-dependent opacities are calculated using a Mie scattering code, miex, that includes distributions of particle sizes and can account for large dust particles (Wolf & Voshchinnikov 2004). The disk models rely on mean opacities calculated from these wavelength-dependent opacities. The Rosseland mean opacity evaluated at 100 K is used to represent the overall disk opacity to its own emission, and the Planck-averaged opacity evaluated at 100 K is used for calculating the thermal emission from the disk. The Planck-averaged opacity evaluated over the stellar spectrum is used to represent the opacity of the disk to stellar irradiation. We model the star as a two-temperature blackbody, with 55% of the luminosity from a 3990 K component, and 45% of the luminosity from a 3600 K component. For a model with a given minimum grain size, we calculate the entire disk model from the initial conditions in order to be completely self-consistent.

The scattered light image of a disk is modeled as in Jang-Condell (2009) and Jang-Condell & Turner (2012). The scattering surface of the disk is defined to be where the optical depth from the star at a given frequency ν is τ_ν = 2/3. The brightness at the scattering surface of the disk, including multiple scattering, is

$$\begin{eqnarray} I^{\rm{{scatt}}}_{\nu } &=& \frac{\omega _{\nu }\mu R_*^2B_{\nu }(T_*)}{4r^2(\mu +\cos \eta)}\nonumber\\ &&\times \left\lbrace 1 + \frac{\omega _{\nu }}{1-g^2\mu ^2} \left[\frac{(2+3\mu)(\mu +\cos \eta)}{(1+2g/3)(1+g\cos \eta)} - 3\mu ^2 \right] \right\rbrace,\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where ω_ν is the albedo, $g=\sqrt{3(1-\omega _{\nu })}$, μ is the cosine of the angle of incidence to the scattering surface, B_ν is the stellar brightness, r is the total distance to the star, and η is the angle between the line of sight to the observer and the normal to the scattering surface. The observed brightness of a star at distance from the observer d is

$$\begin{equation} F_{\nu,{\rm {obs}}} = \pi B_{\nu } \left(\frac{R_*}{d}\right)^2 \end{equation} \tag{ 2 }$$

so we can express the surface brightness in scattered light in units of the apparent brightness of the star per square arcsecond:

$$\begin{equation} I^{{\rm {scatt}}}_{\nu } = \frac{\omega _\nu \mu (1+f_{\rm mult})}{4\pi (\mu +\cos {\eta })} \left(\frac{d}{\mbox{pc}}\right)^2 \left(\frac{r^2+z_s^2}{\mbox{AU}^2}\right)^{-1} \left(\frac{F_{\nu,{\rm {obs}}}}{\mbox{asec}^2}\right), \end{equation} \tag{ 3 }$$

where

$$\begin{equation} f_{\rm mult} = \frac{\omega _{\nu }}{1-g^2\mu ^2} \left[\frac{(2+3\mu)(\mu +\cos \eta)}{(1+2g/3)(1+g\cos \eta)} - 3\mu ^2 \right] \end{equation} \tag{ 4 }$$

represents the fractional increase in brightness caused by multiple scattering.

Finally, we need to choose a composition and size distribution (especially minimum grain size) of the dust to calculate the final scattered light that is emitted from the disk. We use the wavelength-dependent opacities calculated using Mie theory, as described above, to determine the scattered brightness. Many compositions and possible grain sizes are available, and it is not immediately clear whether a single composition is preferred for optically thick disks, or whether each individual disk has its own properties. The neutral color of TW Hya in the optical compared to the quite red colors of HD 100546 suggests that a range of compositions and size distributions exist amongst disks.

The gap in the disk is modeled as an ad hoc axisymmetric density perturbation parameterized by a Gaussian with adjustable width w and depth d centered at 80 AU. To model the truncation of the disk, we introduce an exponential cutoff at a knee k, as appears in the self-similar solution for an accretion disk as derived in Hartmann et al. (1998). If Σ₀(r) is the unperturbed disk surface density, then the new surface density profile is

$$\begin{equation} \Sigma (r) = \Sigma _0(r) \lbrace 1-d\exp [-(r- 80\mbox{ AU})^2/(2w^2)] \rbrace \exp [ -(r/k) ]. \end{equation} \tag{ 5 }$$

Since the primary heating source for the disk is stellar irradiation, the effects of shadowing and illumination must be accounted for in determining the vertical structure of the disk. The 3D density and temperature structure of the perturbed disk are calculated iteratively according to Jang-Condell (2008) and Jang-Condell & Turner (2012). That is, the illumination at the surface of the disk is determined by ray tracing and the disk temperatures are calculated accordingly. Once the disk temperatures are determined, the vertical density profile of the disk is iteratively recalculated assuming hydrostatic equilibrium as above, keeping the vertically integrated surface density profile constant.

How predictive can models of optically thick disks be for composition and size distribution? While we discuss below (Section 4.5) the best fit to our measurements of TW Hya's disk, here we consider the uniqueness of composition within such fits. We consider a fixed disk geometry by taking a single disk model from the parameter space we explored for TW Hya—a disk with a gap located at 80 AU with width w = 30 AU and depth d = 0.3, a truncation knee at k = 100 AU, and a maximum grain size of 1 cm. We then made six independent models of the disk structure that vary composition between pure water ice and astronomical silicates (Warren 1984; Laor & Draine 1993), and vary minimum grain sizes between 5 × 10⁻³, 1, and 10 μm. The dust is well mixed with the gas so that the dust density is proportional to the gas density. We calculated albedo and opacity to find the surface brightness of scattered light as a function of wavelength.

Figure 9 shows a comparison between the six predicted STIS surface brightness profiles for the different models. Differing grain sizes and compositions indeed affect the resulting surface brightness profiles of a given disk structure. If the composition of the disk is incorrect, the structure of the disk can be incorrectly interpreted.

Zoom In Zoom Out Reset image size

In general, pure water ice disks are dimmer than those that possess pure silicates. The reason for this difference is primarily due to the higher opacity of a pure silicate disk, demonstrated in Figure 10. In this figure, we have plotted the height of the τ_ν = 2/3 surface (H) divided by the radius (R) for varying wavelengths. From Equation (3), the brightness of the disk is proportional to the angle of incidence μ at the surface,

$$\begin{equation} \mu \approx \frac{d}{dR}\left(\frac{H}{R}\right) - \frac{H}{R} = R \frac{d}{dR}\left(\frac{H}{R}\right). \end{equation} \tag{ 6 }$$

Note that H represents the wavelength-dependent penetration depth of stellar photons, not a thermal scale height. If H∝R^β, then μ = (β − 1)H/R. Therefore, the brightness scales roughly with H/R. The H surface occurs higher up in the disk for silicates, resulting in a higher surface brightness. The different structure of the surface brightness profiles in Figure 9 is also striking, and is clarified by looking at d(H/R)/dR for the different models (see Figure 11).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Water ice has strong absorption features at 1.5 and 2.0 μm  which are probed by our F160W/F171M and F204M observations. Given the ∼10% accuracy of the disk photometry measurements, we can in principle detect absorption features that have a depth of ∼15%–30%. Figure 12 shows the resulting reflectance spectra for each of our models. Silicates show no absorption lines in the visible to NIR, but are mostly neutral or red over this wavelength range for minimum grain sizes >1 μm. Water ice can show noticeable features, and disks with pure water ice can show detectable features at 20%–40% depth.

Zoom In Zoom Out Reset image size

These tests demonstrate that composition can play an important role in the observed scattered light spectrum and surface brightness profiles of a disk. For our final model of TW Hya, an exhaustive search of compositions, grain sizes, and structures is computationally prohibitive and beyond the scope of this paper.

To fit the observations of TW Hya's disk, we used dust opacities and scattering efficiencies for a composition with the same relative abundances as used to model the SED (D'Alessio et al. 2001): 29.6% organics, 40.4% water ice, 24.5% astronomical silicates, and 5.5% troilite. We test whether such a composition also fits the scattered light observations. We ran a suite of models, varying the following disk parameters as follows:

$$\begin{eqnarray} a_{{\rm min}} & \in & \lbrace 0.005, 0.5, 5\rbrace \mu \rom{m}, \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} w &\in & \lbrace 5, 10, 20, 30\rbrace \mbox{ AU}\ \ \ \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \, d &\in & \lbrace 0.0, 0.3, 0.5\rbrace\qquad \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \qquad \quad k &\in & \lbrace 60, 80, 100, 120, 150\rbrace \mbox{ AU}. \end{eqnarray} \tag{ 10 }$$

The surface brightness profiles are fit simultaneously spatially and across the seven wavelengths for which we have photometric data. To account for uncertainties in the overall normalization such as might arise from uncertainty in our chosen parameters, we allow the overall brightness to vary by a constant factor and calculate the reduced χ² value for each model with respect to the wavelength-dependent surface brightness profiles from 56 to 160 AU. The center of the gap is fixed at 80 AU.

The parameters for the best fits as measured by the minimum reduced $\chi ^2_{\nu }$ from the fits above are tabulated in Table 3, where we found that a disk gap of 30% depth and 20 AU width are preferred. The profiles are well fit by a truncation knee of 100 AU. We found that we could fit the observed brightness of the disk only with small a_min (0.005 μm), but if we were willing to scale the overall brightness, we could equally fit (similar χ²) the disk with a_min = 0.5 μm. However, larger grain sizes failed to reproduce the neutral-to-blue color observed over large portions of the disk, resulting instead in a redder color. Since the overall flux density of the disk is an observed quantity, we favor those models that naturally reproduced the observed surface brightness.

Table 3. Best-fit Disk Model Values

Model Parameter	Model Value
amin	0.005 μm
Gap Width	20 AU
Gap Depth	0.3
Truncation	100 AU
$\chi ^2_{{\rm nu}}$	3.29
Fixed Parameters
L⋆	0.208 L☉
Teff	3600 K
R⋆	1.08 R☉
$\dot{M}$	10−9 M☉ yr−1
α	5 × 10−4
amax	1 cm

Download table as:  ASCIITypeset image

A few of our assumptions could impact the expected brightness of the disk. For example, the luminosity of the star would differ by ∼20% if we instead used the luminosity determined by Vacca & Sandell (2011), rather than that by Webb et al. (1999) or our own determination. Since the surface brightness of TW Hya's disk due to scattered light depends on the square of the distance, uncertainties of up to 23% are possible, given the uncertainties in TW Hya's parallax. A different distance would also imply revised stellar parameters. A higher or lower mass of the central star will result in a puffier (brighter) or more vertically compact (dimmer) disk, respectively.

We have also assumed that the dust isotropically scatters light. If the dust were instead forward scattering by a moderate degree, such as with a Henyey–Greenstein asymmetry parameter of 0.5, the approximate scattering angle of TW Hya would result in a dimming of the disk of ≈25%, based on our modeling. Due to flaring in the disk, our isotropic models predict a decrease in flux along the semi-minor axis of the disk pointing to the observer. Conversely, forward scattering dust, even in a flared disk, will produce a brightness peak along the semi-minor axis of the disk tilted toward the observer. The true orientation of the disk is degenerate between these two possibilities, so we orient our models such that the brightness maximum is aligned with the observed P.A. in the STIS and F110W data.

Closer to the star, the models with isotropic scattering do not completely reproduce the azimuthal brightness distribution, while the forward scattering model is more consistent. Our models suggest an inclination of ∼7° and forward scattering dust for the disk. Forward scattering grains in the outer disk are marginally favored as both the STIS and F110W data still show a trough and peak structure in the azimuthal surface brightness profile—while the S/N is lower, the expectation from our models is a larger brightness asymmetry at larger distances, which is not favored by the data. This is indicative of either a decrease in propensity to forward scatter or lack of forward scattering grains in the outer disk.

Figures 13 and 14 show the radial data compared to the best-fit model for all of the observations and radial spectral cuts of the disk. In general, the model fits the observed spectrum quite well, with the exception of the region around 80 AU, the position of the gap. This is also where the model fails to reproduce the observed surface brightness profile in some of the wavelengths observed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In comparison, models with no gap provide poor fits to the wavelength-dependent surface brightness profiles. In Figure 2, we show brightness profiles of disk models without a gap at 80 AU in comparison to the data at the F222M wave band. The disk models include truncation at k = 60, 100, and 150 AU and have been scaled by 1.05, 0.73, and 0.57, respectively, in order to fit the data across all wavelengths. The reduced $\chi ^2_{\nu }$ for the models are 20.8, 8.18, and 4.31, respectively. The model with truncation at 150 AU is the best fit, but still does not provide as good a fit as compared to the model with a 30% gap.

While our model successfully reproduces the observed properties of the TW Hya disk over our wavelength range, there remain some uncertainties that will require follow-up study. In particular, we further discuss the nature of the observed brightness asymmetry for the disk seen at short wavelengths, the more spatially extended nature of small dust grains in optical light compared to continuum submillimeter and millimeter observations, and alternative origins for the observed gap at 80 AU.

The brightness asymmetry first reported by Roberge et al. (2005) and further characterized in Section 3.2 is strange given the asymmetry is only apparent at wavelengths ⩽1.1 μm, whereas our models at any disk inclination would predict asymmetry at wavelengths >1 μm also. A warped inner disk, as suggested by Roberge et al. (2005) and Rosenfeld et al. (2012), also seems unlikely, as that would also be expected to produce asymmetry at long wavelengths. What could cause a wavelength dependence to the asymmetry?

For an inclined, optically thick disk with isotropically scattering grains, the near side of the disk appears slightly brighter because of the lower opacity along the line of sight. However, the near side also appears foreshortened. Our azimuthal profiles were taken assuming a circular disk, which would artificially decrease the surface brightness for P.A.s corresponding to the side of the disk pointing toward the observer. Therefore, our models do not uniquely predict the brightness asymmetry.

Forward scattering grains cancel the foreshortening effect by brightening the near side. If the disk grains are substantially more forward scattering at short wavelengths, the apparent brightness asymmetry could change with wavelength. However, forward scattering grains would also make the whole disk look dimmer than observed, since we image it nearly face-on.

TW Hya is known to have variable accretion (e.g., Dupree et al. 2012), and this could drive changes in the scale height of the inner disk. Variability in the inner disk structure is unlikely as the source of the changing asymmetry. Although the inner disk could change scale height due to variable accretion, the F110W, STIS images and STIS spectroscopy show the same asymmetry but were taken two, and then two more, years apart, respectively, while the F160W image does not show asymmetry and was taken only three months before the F110W image.

If portions of the outer disk were shadowed from direct starlight by a non-axisymmetric bump or warp in the inner disk, that might be able to produce asymmetry in the observed scattered light (Rosenfeld et al. 2012). The structure would become more optically thin at longer wavelengths. The problem is that the visible–NIR extinction would have to be nearly gray to match the observed photometry of the disk. That only happens from very large grains, whereas we determine that very small grains are plentiful in the outer disk.

Alternatively, the lack of asymmetry at longer wavelengths could be the result of a changing disk structure or grain population with depth, as longer wavelengths probe deeper in the disk. The forward scattering grains seen at wavelengths shortward of 1.1 μm could no longer be present deeper in the disk to be replaced by a different (perhaps larger) grain population. Our models predict where the changeover in grain properties occurs as defined by the scattering surface between F110W and F160W. At 80 AU, the scattering surface of the disk is at 19.9 AU above the midplane for 1.1 μm, while it is at 19.3 AU for 1.6 μm. If such a vertical transition exists, it would be within a very narrow layer of the disk. If any upper layer structure is the cause for the azimuthal asymmetries, they must reside only in layers >19.3 AU above the midplane. Such an abrupt change seems unlikely.

A more detailed examination of the effect of forward scattering grains, coupled with realistic models of forward scattering from aggregate grains, will be needed to make progress on interpreting the asymmetry.

We next discuss the overall size of the disk. As has been seen in CO and dust continuum data at submillimeter and millimeter wavelengths, the gas in the disk of TW Hya extends much further than submillimeter dust emission (Andrews et al. 2012). This dichotomy apparently does not extend to smaller grain sizes as the scattered light disk is nearly co-spatial with the extent of the CO gas.

We find that the best-fit model for the scattered light has a large truncation knee, at 100 AU, although we explored values down to 60 AU. Smaller values of k produce steep surface brightness profiles not warranted by our data (see Figure 2). In comparison, radio interferometric observations find much smaller truncation radii for similar models. Using the similarity solution (Lynden-Bell & Pringle 1974),

$$\begin{equation} \Sigma (R) = \frac{c_1}{R^{\gamma }} \exp \left[-\left(\frac{R}{k} \right)^{2-\gamma } \right]. \end{equation} \tag{ 13 }$$

Andrews et al. (2012) find k = 35–45 AU for values of γ between −1 and 1 fit SMA observations at 870 μm continuum and the CO J(3–2) line; Hughes et al. (2008) find k = 30 AU and γ = 0.7 for SMA observations at 1.3 mm and 870 μm (345 and 230 GHz) continuum, although Hughes et al. (2011) find k = 50 AU and γ = 1 is the best fit to the CO J(3–2); and Isella et al. (2009) find k = 17.5 AU and γ = −0.3 for CARMA observations at 1.3 mm continuum. Models by Gorti et al. (2011) find that a truncation knee of 35 AU and γ = 0.7 provides good fits to observations of gas emission lines in TW Hya.

The model fits to radio observations cited above predict that the disk density beyond ∼50 AU should be sharply cut off. As k and γ decrease, the more sharply the disk will be truncated. However, the scattered light observations clearly indicate the presence of disk material beyond 200 AU, so we do not expect to find a small truncation knee. The scattered light brightness profile of the TW Hya disk suggests that the disk truncation is at a larger radius than the radio observations imply. Since scattered light probes only the diffuse surface layer of the disk while radio observations probe the deeper interior structure of the disk, the difference in truncation radii may be attributable to layered disk structure, where the vertical profile of the disk varies with distance in a way not adequately modeled simply by vertically integrated surface density profiles. Conversely, there could be an overall lack of large grains beyond ∼60–80 AU, indicative of grain growth frustration beyond this radial distance in the disk.

Finally, we investigate other possibilities for the presence of the gap. It is interesting to note in Figure 14 that at 80 AU, the center of the gap, the spectral character of the disk is different from that seen in other portions of the disk. If this is real, then a different composition in this annulus may be causing the overall structure of the disk to change, or vice versa. Gaps in disks tend to receive larger amounts of UV photons due to isotropic scattering of Lyα photons from higher layers of the disk (Bethell & Bergin 2011). The excess UV photons can accelerate photodesorption of ices, which could in turn change the spectral character of the disk within the gap (Bergin et al. 2010; Hogerheijde et al. 2011).

Conversely, compositional changes in this gap might be caused by enhanced collisions within the disk or the presence of an accreting protoplanet causing a steep truncation to the disk. Both of these possibilities are supported in part by the recent discovery of a sharp cutoff in the millimeter emission from dust at 60 AU, which is not present in similar high-resolution imaging of CO gas (Andrews et al. 2012). Since we see scattered light from dust well beyond this cutoff, it is clearly due to a change in the bulk size distribution or millimeter emission properties of the dust, which in of itself may be a signature of a planetary companion or change in bulk dust composition.

The gap could instead be an unresolved spiral structure in the disk, as seen in NICMOS images of HD 141569, where spiral structure at lower spatial resolution mimicked a gap structure (Weinberger et al. 1999; Clampin et al. 2003). This could explain both the asymmetry as well as the decrease in surface brightness, and a hint of non-axisymmetric structure is seen in the higher spatial resolution STIS images of the TW Hya disk. Higher spatial resolution images might further solve this mystery.

The provocative nature of such a structure indicates that significant, large-scale physical processes at large separations occur for protoplanetary disks and can significantly impact their midplane structures. Whatever the specific origin of this feature within the TW Hya disk, the combined observations of optical/NIR scattering with longer wavelength photometry and spatially resolved imaging may lead to a new understanding of the conditions in which planetary systems form.