CONSTRAINTS ON THE RADIAL VARIATION OF GRAIN GROWTH IN THE AS 209 CIRCUMSTELLAR DISK

The growth of submicron-sized dust grains into millimeter- and centimeter-sized particles is a fundamental component of the planet formation process (Beckwith et al. 2000; Natta et al. 2007). Grain growth directly affects the optical properties of dust particles (Henning & Stognienko 1996; D'Alessio et al. 2001), with composition and temperature having a lesser impact (Henning & Mutschke 1997; Semenov et al. 2003; Boudet et al. 2005). At long wavelengths, the slope α of the spectral energy distribution (SED, νS_ν∝ν^{α + 1}) can be related to the spectral index of the dust opacity β (κ_ν∝ν^β), in particular β = α − 2 for optically thin warm dust (Beckwith & Sargent 1991; Miyake & Nakagawa 1993). The millimeter dust opacity slope for small interstellar medium (ISM) dust grains is β_ISM ∼ 1.7 (Li & Draine 2001). A smaller value of β is measured if larger millimeter-sized grains are present (D'Alessio et al. 2001; Draine 2006).

Observations of dust in disks from submillimeter to centimeter wavelengths (Testi et al. 2001, 2003; Calvet et al. 2002; Natta & Testi 2004; Wilner et al. 2005; Rodmann et al. 2006; Lommen et al. 2009; Ricci et al. 2010a, 2010b) have provided strong evidence for grain growth in disks. These results, obtained mainly from spatially unresolved observations, assumed a constant dust opacity (i.e., constant β) throughout the disk. However, theoretical models of dust transport, fragmentation, and size evolution in protoplanetary disks predict that the average size of grains varies with the distance to the central star (hereafter referred to as radius; Dullemond & Dominik 2005; Birnstiel et al. 2010). Hence, radial variations in the dust opacity slope are expected.

The increased sensitivity and angular resolution of current interferometers have enabled the first studies of radial variations of grain growth within circumstellar disks. Observations of disks in Taurus at 1.3 and 3 mm with CARMA¹⁵ (Isella et al. 2010) and PdBI¹⁶ (Guilloteau et al. 2011) have constrained β as a function of radius. However, these results are only sensitive to radial variations larger than Δβ ∼ 0.6–0.7 at 3σ in the best cases, limited by the small wavelength separation between the two bands being observed. Increased wavelength coverage has been obtained by adding centimeter-wave observations (Banzatti et al. 2011), but the low signal-to-noise ratio, before the Karl G. Jansky Very Large Array (VLA) upgrade, made such studies challenging. To significantly improve these constraints, and effectively compare them to theoretical predictions of grain growth, increased wavelength coverage and enhanced sensitivity are required.

This Letter presents the first results of a program that combines spatially resolved observations from submillimeter to centimeter wavelengths, to constrain radial variations of grain growth within circumstellar disks. The pre-main-sequence star AS 209, also known as V1121 Oph, is a classical K5 T Tauri star with high accretion rate (1.3 × 10⁻⁷ M_☉ yr⁻¹; Johns-Krull et al. 2000). AS 209 is isolated from the main ρ-Ophiuchus cloud ($\alpha ={\rm 16^h49^m15\mbox{$.\!\!^{\mathrm{s}}$}3},\delta =-{\rm 14^d22^m08\mbox{$.\!\!^{\mathrm{s}}$}7}$) at a distance of 131  ±  50 pc (van Leeuwen 2007). Resolved 0.88 mm continuum observations reveal a smooth radial density structure down to 03 ∼ 40 AU (Andrews et al. 2009, 2010). We obtained additional spatially resolved observations at λ > 0.88 mm, to study grain growth inside the AS 209 disk.

Submillimeter Array (SMA) observations at 0.88 mm were obtained between 2006 May and 2007 June, in C, E, and V configurations (∼03 resolution), with total continuum bandwidth of 4 GHz. Its integrated flux density is 577  ±  60 mJy at 0.88 mm. Andrews et al. (2009) presented a detailed discussion of the observations and data reduction.

CARMA observations at 2.8 mm were obtained between 2009 December and 2010 March, in C, B, and A configurations spanning baselines of 30–1800 m (∼05 resolution), with total continuum bandwidth of 8 GHz. Its integrated flux density is 40 ± 6 mJy at 2.8 mm. Atmospheric conditions were good, with zenith opacities τ_230 GHz ⩽ 0.15, except for B configuration observations where τ_230 GHz ∼ 0.4. Data were calibrated using the Multichannel Image Reconstruction, Image Analysis and Display software package (Sault et al. 1995).

VLA observations at 8.0 and 9.8 mm (Ka band) were obtained between 2011 January–May (project AC982), in CnB and BnA configurations, spanning baselines of 350 m to 36.4 km (∼02 resolution). Two 1 GHz basebands centered at 30.5 and 37.5 GHz provided 2 GHz of continuum bandwidth. At 8.0 and 9.8 mm, the integrated flux density of AS 209 is 1.2 ± 0.1 and 0.7 ± 0.1 mJy, respectively. Data were calibrated using the Common Astronomy Software Applications (CASA).

For all data sets, quasar observations interspersed between science target observations were used to calibrate the complex antenna gains. Absolute flux calibration was obtained by observing Neptune for CARMA, Uranus and Titan for SMA, and 3C286 for VLA. For each telescope, calibrated visibilities coming from different configurations were corrected for proper motion (van Leeuwen 2007) and compared over the spatial frequencies where these configurations overlap; these were found to agree within the absolute flux scale uncertainty (∼10%–15%). Fourier inversion of the visibilities and cleaning of the dirty image to produce dust emission map were performed using CASA.

Figure 1 presents the normalized real/imaginary part of the correlated emission from AS 209 as a function of spatial frequency (uv-distance), for each wavelength. The real part of the visibility will be constant with increasing spatial frequency for an unresolved point source. In contrast, for AS 209 we measure a declining profile with increasing uv-distance, indicative of spatially resolved emission at all wavelengths. Furthermore, short-wavelength emission (0.88, 2.8 mm) diminishes faster than long-wavelength emission (8.0, 9.8 mm); therefore, the structure at long wavelengths is more compact than at short wavelengths.

Zoom In Zoom Out Reset image size

Before interpreting these observational results in terms of grain properties variations, we explore other explanations for the differing amplitude profiles: (1) seeing affecting high-frequency observations would have to be >12 at 0.88 mm (10 × larger than observed seeing; Andrews et al. 2009), (2) optically thick emission at short wavelengths: modeling of these observations (Section 4) indicates τ < 1 over the scales studied. Consequently, we conclude that the long-wavelength emission comes from a smaller region than the short-wavelength emission.

Before proceeding, we account for continuum emission not arising from dust, such as free–free emission from photo-evaporative stellar winds. Such emission originates near the star (within 10 AU for TW Hya; Pascucci et al. 2011), appearing as a point source at the resolution of our observations. Indeed, in addition to the resolved emission detected with VLA, we detect a point-like component (at uv-distances >900 kλ) with a flux density of 0.08  ±  0.02 mJy (8 mm) and 0.10  ±  0.02 mJy (9.8 mm), corresponding to ∼7% and ∼14% of the total emission at these wavelengths. However, VLA observations at 5.2 cm from the Disks@EVLA Collaboration¹⁷ (C. J. Chandler et al. in preparation) place a 3σ upper limit of 72 μJy, suggesting that the point-like emission at 8 and 9.8 mm may be dust emission. Nonetheless, we conservatively include this point-like emission as a free–free component in models described next.

Observations were analyzed using the disk emission model described in Isella et al. (2009). The disk structure follows the two-layer disk approximation (Chiang & Goldreich 1997), defined by the boundary between the disk surface layer (optically thin to stellar radiation) and the disk interior (optically thick to stellar emission). The disk is passively heated by AS 209 (L_⋆ = 1.5 L_☉, T_⋆ = 4250 K, M_⋆ = 0.9 M_☉; Andrews et al. 2009). Dust in the disk surface absorbs the stellar radiation and re-radiates this energy at longer wavelengths. Because this hot dust layer will be optically thin to its own thermal radiation, half of the surface layer emission will be radiated outward while the rest will be radiated toward the disk, heating up its interior.

For the density structure, we employ the self-similar solution for a viscous accretion disk (Lynden-Bell & Pringle 1974)

$$\begin{eqnarray} \Sigma (R) =&&\Sigma _t\left(\frac{R_t}{R}\right)^{\gamma } \nonumber \\ &&\times \exp \left[ -\frac{1}{2(2-\gamma)}\left[\left(\frac{R}{R_t}\right)^{(2-\gamma)} -1\right]\right], \end{eqnarray} \tag{ 1 }$$

where Σ_t corresponds to the surface density at radius R_t. This Σ(R) prescription behaves as a power law for small radii, decreasing exponentially at large radii.

We compute the dust opacity κ_ν using Mie theory, assuming a population of compact spherical grains in a power-law size distribution, n(a)∝a^−q, for a_min < a < a_max with a_min = 0.01 μm. Optical constants for astrophysical grain components were obtained from Semenov et al. (2003), Zubko et al. (1996), and Warren (1984). Fractional abundances were adopted from Pollack et al. (1994; in volume: 8% silicates, 30% organics, 62% water ice). Recent revisions of the solar oxygen abundance (Asplund et al. 2009) would reduce the amount of water ice. However, the effect on the disk modeling is minimal, with a difference in the disk structure constraints well within their uncertainties as shown by Isella et al. (2010).

The best-fit model to a single wavelength observation is found through χ² minimization with three free parameters that describe the surface density (Σ_t, R_t, γ), since the disk geometry is fixed based on molecular line observations (i = 38°, P.A. = 86°; Andrews et al. 2009). The χ² probability distribution is sampled by varying the free parameters, generating a model state and comparing this state to the data in Fourier space, following the Markov chain Monte Carlo (MCMC) method described in Isella et al. (2009).

For this initial model fit, we select the maximum grain size (a_max) and particle-size distribution slope (q) that best fit the unresolved SED from submillimeter to centimeter wavelengths. For AS 209, these correspond to q = 3.5, a_max = 1.3 mm, resulting in β = 1.0 between 0.88 and 8.0 mm. Then, assuming a constant dust opacity throughout the disk, each wavelength is fitted separately with the aforementioned physical model. We obtain best-fit Σ(R) and T(R) profiles, which will be different for each wavelength if the assumption of a radially constant dust opacity is not valid (Isella et al. 2010).

Figure 2 shows the best-fit and 3σ constraints on T(R) and the disk optical depth τ_λ(R) = κ_λΣ(R), inferred from modeling each observation separately. As can be seen in this figure, a differing Σ(R) and T(R) was inferred for each wavelength when κ_λ is assumed to be constant with radius. We constrain R_t to be large for the short-wavelength emission (R_t = 61 ± 2 AU from 0.88 mm and R_t = 60 ± 6 AU from 2.8 mm), compared to R_t inferred at long wavelengths (R_t = 24 ± 3 AU from 8 mm, R_t = 26 ± 4 AU from 9.8 mm). We note that the derived temperature profiles differ by <4 K for R < 70 AU, and by ≲ 10 K for R > 70 AU. The best-fit model visibility profiles at each wavelength are shown as continuous lines in Figure 1. Maps of the observed emission, with corresponding best-fit model and residual emission maps are shown in Figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Since the disk is optically thin (right panel, Figure 2), the observed emission will depend directly upon the dust opacity and disk mass: S_λ∝κ_λ × Σ × B_λ(T). Hence, our observations constrain the product κ_λ × Σ × B_λ(T), where all these physical parameters (κ_λ, Σ, T) depend on the radius, but only the dust opacity depends on the observed wavelength. The modeling presented above, assuming a constant dust opacity with radius and fitting each wavelength separately, has resulted in a wavelength-dependent Σ_λ(R) and T_λ(R), which is unphysical but telling: the assumption of a constant κ_λ with radius is likely incorrect. Radial variations of the dust opacity are required to reconcile the differing visibility profiles (real part versus uv-distance, Figure 1) and the wavelength-dependent Σ_λ(R) and T_λ(R) (Figure 2).

5.1. Radial Variations of β

Expanding on the work presented in Isella et al. (2010), if the true but unknown disk physical quantities are $\overline{\kappa }_{\lambda }(R)$, $\overline{\Sigma }(R)$, and $\overline{T}(R)$, then at each wavelength

$$\begin{equation} \kappa _{\lambda } \Sigma _{\lambda }(R) B_{\lambda }(T_{\lambda }(R))=\overline{\kappa }_{\lambda }(R) \overline{\Sigma }(R) B_{\lambda }(\overline{T}(R)), \end{equation} \tag{ 2 }$$

where the right-hand side of this equation corresponds to the true disk properties and the left-hand side encompasses the constraints found from our multi-wavelength observations. Taking the ratio of Equation (2) for λ₁ ≠ λ₂, and assuming that at long wavelengths κ_λ∝λ^−β, we obtain

$$\begin{eqnarray} &&\Delta \beta (R)=\log ^{-1}\left[\frac{\lambda _1}{\lambda _2}\right]\times \log \left[\frac{\Sigma _{\lambda _1}B_{\lambda _1}(T_{\lambda _1})/B_{\lambda _1}(\overline{T})}{\Sigma _{\lambda _2}B_{\lambda _2}(T_{\lambda _2})/B_{\lambda _2}(\overline{T})}\right], \end{eqnarray} \tag{ 3 }$$

where Σ_λ, T_λ, and $\overline{T}$ depend on the radius R. Equation (3) is a useful prescription to constrain radial variations of the opacity spectral slope (Δβ(R) = β(R) − β_C) in dual-wavelength observations. Here β_C corresponds to the spectral slope of the assumed radially constant dust opacity.

In the Rayleigh–Jeans regime (when hc/λ ≪ k_BT, k_B: Boltzmann constant), knowledge of $\overline{T}(R)$ to derive Δβ(R) is unnecessary since B_λ(T)∝T. However, in the cold outer disk and at short wavelengths, the Rayleigh–Jeans assumption is ill-founded (e.g., hν ∼ k_BT for T = 20 K at 0.88 mm). Consequently, we require an estimate of $\overline{T}(R)$. Since the temperature profiles inferred from each wavelength (T_λ(R)) are not very different, we average them to deduce $\overline{T}(R)$ (dashed line in Figure 2).

Equation (3) shows that in logarithmic space Δβ(R) is the slope of a line that goes through points {x = log (λ), $y=\log [\Sigma _{\lambda }B_{\lambda }(T_{\lambda })/B_{\lambda }(\overline{T})]$}, with y evaluated at radius R. Hence, to measure Δβ(R) we find the slope of a straight line—in a least-square sense—that fits our constraints on the product $\Sigma _{\lambda }\times B_{\lambda }(T_{\lambda })/B_{\lambda }(\overline{T})$ from our set of multi-wavelength observations. To set confidence intervals on Δβ(R), we make use of the Bayesian approach upon which our modeling is based. Since the result of the MCMC algorithm is a fully sampled posterior probability distribution function (PDF) for all the model parameters, we construct a PDF of the product $\Sigma _{\lambda }\times B_{\lambda }(T_{\lambda })/B_{\lambda }(\overline{T})$ at each radius R and for each wavelength. Random samples of these PDFs are taken at each wavelength, the slope of the line through points {x = log (λ), $y=\log [\Sigma _{\lambda }\times B_{\lambda }(T_{\lambda })/B_{\lambda }(\overline{T})]$} at a radius R, is computed. This slope corresponds to one random sampling of Δβ(R). Hence, the PDF for Δβ(R) is constructed by performing a large number of random samples. The peak of this PDF is the best-fit value of Δβ at radius R. Confidence intervals are derived from the region of the PDF that contains 68.3%, 95.5%, and 99.7% of all samples at equal probability (1σ, 2σ, and 3σ).

Figure 4 presents the constraints on the radial variation of β obtained from our multi-wavelength observations. The values of β allowed by our observations are significantly different than β_ISM ∼ 1.7, for R ≲ 70 AU. Furthermore, we find a gradient on β(R) inconsistent with a constant value at the 10σ level.

Zoom In Zoom Out Reset image size

5.2. Radial Variations of a_max

To derive Equation (3), the assumption of κ_λ∝λ^−β must be satisfied. We caution that for a_max ∼ 0.1–1 mm this approximation may break down (Draine 2006), as illustrated in Figure 4 (right). Therefore, rather than inferring a_max(R) from β(R), we constrain it directly by fitting a specific dust opacity $\overline{\kappa }_{\lambda }$ to the constraints on the product κ_λ × Σ_λ × B_λ(T_λ) at each radius. With a knowledge of $\overline{T}(R)$, and for a fixed set of dust properties (composition and grain-size distribution), we estimate the values of a_max and $\overline{\Sigma }$ that satisfy Equation (2), now written as

$$\begin{equation} \kappa _{\lambda } \Sigma _{\lambda }(R) \frac{B_{\lambda }(T_{\lambda }(R))}{B_{\lambda }(\overline{T}(R))} =\overline{\kappa }_{\lambda }(a_{{\rm max}}(R)) \overline{\Sigma }(R), \end{equation} \tag{ 4 }$$

where the right-hand side corresponds to our model (with parameters a_max and $\overline{\Sigma }$) and the left-hand side has been constrained by our multi-wavelength observations (i.e., we have a PDF for the product $\kappa _{\lambda }\times \Sigma _{\lambda }(R)\times {B_{\lambda }(T_{\lambda }(R))}/{B_{\lambda }(\overline{T}(R))}$).

At each radius R, we constructed a two dimensional grid of parameters $\lbrace a_{{\rm max}},\overline{\Sigma }\rbrace$. At each point in the grid, we compute the product $\overline{\kappa }_{\lambda }(a_{{\rm max}})\times \overline{\Sigma }$ at each wavelength. We then find the probability that such measurement will have occurred (given our observational constraints on the left-hand side of Equation (4)), and construct the likelihood function of the parameters $\lbrace a_{{\rm max}},\overline{\Sigma }\rbrace$. Best-fit values for a_max and $\overline{\Sigma }$ are found by maximizing the likelihood; confidence intervals are obtained from the marginalized likelihoods.

Figure 5 presents our constraints on a_max(R) and $\overline{\Sigma }(R)$ for two representative values of q. The same composition presented in Section 4 is assumed, however a different dust mixture will influence the derived a_max(R). For a composition that includes updated oxygen abundances (Asplund et al. 2009), a smaller a_max(R) is inferred (× 2 smaller), well within the uncertainties of this derivation given other unknown parameters (e.g., particle-size distribution slope). Across the disk, grains have grown at least up to ∼0.5 mm, with small grains present in the outer disk and large grains in the inner disk.

Zoom In Zoom Out Reset image size

We compare our observational constraints with theoretical models of grain growth, employing the approximations presented in Birnstiel et al. (2012) for the evolution of a_max with radius. Since the true $\overline{\Sigma }(R)$ profile constrained by our observations depends on the assumed value of q (see Figure 5), these theoretical prescriptions will depend on q as well. Our observational constraints on a_max are consistent with a radial drift limited grain population, where the head-wind felt by dust particles makes them spiral toward the star. A fragmentation-dominated population, where the turbulent relative motion of particles causes collisions that either grow or fragment these grains, seems incompatible for standard values of the turbulence parameter (α_t = 0.01; Shakura & Sunyaev 1973), fragmentation threshold velocity (u_t = 10 m s⁻¹; Blum & Wurm 2008), and 100:1 gas-to-dust ratio. These curves represent barriers that prevent further size increase, hence, the smaller of the two is considered the upper limit to growth. However, the parameters that go into deriving the fragmentation limit are very uncertain, while the physics and parameters in radial drift are better established. Either no fragmentation barrier exists (for a low-turbulence disk, α_t ≲ 0.001, the maximum collision speed never reaches the fragmentation threshold velocity in AS 209, making fragmentation impossible) or the gas-to-dust ratio is ∼10 × larger, allowing for the fragmentation-limited barrier to be consistent with our a_max(R) constraint. It remains to be demonstrated whether these models reproduce the observed β(R) profile and the millimeter fluxes. Such analysis requires numerical simulations starting from the disk parameters derived here and will be discussed in a future paper.

5.3. Summary

We acknowledge T. Birnstiel for useful discussions. A.I., J.M.C., L.M.P acknowledge support from NSF award AST-1109334. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Ongoing CARMA development and operations are supported by the National Science Foundation under a cooperative agreement, and by the CARMA partner universities. The SMA is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, funded by the Smithsonian Institution and Academia Sinica. Part of this research was carried out at the Jet Propulsion Laboratory, Caltech, under a contract with the National Aeronautics and Space Administration.

Facilities: CARMA - Combined Array for Research in Millimeter-Wave Astronomy, SMA - SubMillimeter Array, VLA - Very Large Array