DUST EVOLUTION CAN PRODUCE SCATTERED LIGHT GAPS IN PROTOPLANETARY DISKS

We make the assumption that solid particles will only fragment efficiently in the inner region of a disk, where impact velocities are high. This should generally be the case, since those velocities scale with either ${c}_{{\rm{s}}}$ (in the case of Brownian motion or turbulent mixing) or ${V}_{{\rm{k}}}$ (for vertically settling and radial or azimuthal drift), both of which decrease with r. We define ${r}_{{\rm{f}}}$ as a critical radius in the problem: fragmentation is suppressed beyond it. As highlighted above, there are two possible reasons for that suppression. First, the maximum collision speeds between particles could be lower than required for fragmentation; i.e., Δv_max < ${v}_{{\rm{f}}}$. Second, and more important, the particles could achieve high enough drift speeds that they move inwards before they can reach the sizes (and thereby relative velocities) at which fragmentation is effective. This latter case is the "drift limit," where ${a}_{{\rm{d}}}$ < ${a}_{{\rm{f}}}$; it can be reformulated using Equations (1) and (2) as

$$\begin{eqnarray}&&\epsilon \;{\alpha }_{{\rm{t}}}\lt \displaystyle \frac{| \gamma | }{3}{\left(\displaystyle \frac{{v}_{{\rm{f}}}}{{V}_{{\rm{k}}}}\right)}^{2},\end{eqnarray} \tag{ 4 }$$

where $\epsilon$ = Σ_d/Σ_g is the dust-to-gas mass ratio. γ, $\epsilon$, and ${\alpha }_{{\rm{t}}}$ are not expected to be uniform throughout the disk, but by assuming they are, we can roughly estimate the maximum radius that fulfills this condition,

$$\begin{eqnarray}&&\displaystyle \frac{{r}_{{\rm{f}}}}{100\;\mathrm{AU}}\sim \displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\;\displaystyle \frac{{\alpha }_{{\rm{t}}}}{{10}^{-3}}\;\displaystyle \frac{\epsilon }{0.01}\;{\left(\displaystyle \frac{| \gamma | }{2.75}\right)}^{-1}\;{\left(\displaystyle \frac{{v}_{{\rm{f}}}}{10\;{{\rm{m\; s}}}^{-1}}\right)}^{-2}.\end{eqnarray} \tag{ 5 }$$

Radial drift can significantly reduce $\epsilon$ in the outer disk, thereby pushing ${r}_{{\rm{f}}}$ well below 100 AU. Lower limits on ${r}_{{\rm{f}}}$ are likely given by sintering (e.g., Sirono 2011) or by a reduced ${v}_{{\rm{f}}}$ inward of the water evaporation front where particles loose their ice mantles.

The particle size distribution inside ${r}_{{\rm{f}}}$ (region I) will quickly reach a coagulation-fragmentation equilibrium, for which analytical fits have been derived by Birnstiel et al. (2011). Here we approximate these as

$$\begin{eqnarray}{{\rm{\Sigma }}}_{{\rm{d}}}(r,a)\;=\left\{\begin{array}{ll}C(r)\cdot {a}^{1/2}\phantom{\rule{5em}{0ex}}\mathrm{for} & a\leqslant {a}_{\mathrm{BT}}\\ {{\rm{\Sigma }}}_{{\rm{d}}}(r,{a}_{\mathrm{BT}})\cdot {\left(\displaystyle \frac{a}{{a}_{\mathrm{BT}}}\right)}^{-5/8}\;\mathrm{for} & {a}_{\mathrm{BT}}\lt a\lt {a}_{{\rm{f}}}\end{array},\right.\end{eqnarray} \tag{ 6 }$$

where a_BT (itself a function of radius) is the transition from Brownian motion to turbulent relative velocities (cf. Equation (37) in Birnstiel et al. 2011), and C(r) is a normalization constant such that ${{\rm{\Sigma }}}_{{\rm{d}}}(r)=\int {{\rm{\Sigma }}}_{{\rm{d}}}(r,a)\;{da}.$

Even though fragmentation is inefficient outside ${r}_{{\rm{f}}}$, larger radii can still be supplied with small particles (from region I) due to outward radial diffusion. If we consider such particles to be passive tracers of the gas, we can apply the analytical treatment of Pavlyuchenkov & Dullemond (2007) to predict their concentration outside ${r}_{{\rm{f}}}$. In this scenario, the surface density of a given particle size can be approximated in steady state with a power-law distribution,

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}(r,a)}{{{\rm{\Sigma }}}_{{\rm{g}}}(r)}\propto {\left(\displaystyle \frac{r}{{r}_{{\rm{f}}}}\right)}^{-\displaystyle \frac{3}{2}\;\mathrm{Sc}}\;{\rm{for}}\quad {r}_{{\rm{f}}}\lt r\lt {r}_{0},\end{eqnarray} \tag{ 7 }$$

where $\mathrm{Sc}$ is the Schmidt number, the ratio of the particle diffusivity to the gas diffusivity. Here we assume $\mathrm{Sc}$ = 1 for all particle sizes. r₀ is the radius at which the particle size a equals the drift limit ${a}_{{\rm{d}}}$.

Outside of ${r}_{{\rm{f}}}$, particles grow toward a limiting size, ${a}_{{\rm{d}}}$, defined by equality between the drift and growth timescales. Most of the dust mass then resides close to ${a}_{{\rm{d}}}$ (see Figure 1). However, radial diffusion effectively broadens the particle size distribution: even if particle growth were to produce a hypothetical monodisperse size distribution like δ(a − ${a}_{{\rm{d}}}$) at a given r, radial diffusion will move particles from larger r inwards to that radius. Those particles will be smaller due to the radial dependence of ${a}_{{\rm{d}}}$, and therefore the size distribution will necessarily broaden in the presence of turbulent diffusion. This inward mixing is balanced by sweep-up collisions with drifting particles of size ${a}_{{\rm{d}}}$. To approximate these effects, we will parametrize the radial dependence of the surface density for a given grain size as a power law with index p_d. This index can be estimated from balancing drift, diffusion, and sweep-up collisions (T. Birnstiel et al. 2015, in preparation); it generally depends on various physical parameters of the disk, but in the following, we take it to be p_d = 2.5 (the results of interest here are rather insensitive to p_d).

The same principles also apply to the outward diffusion of the drifting particles. In this case, radial drift is opposing that outward diffusion. Similar to the case in region II, we can find a steady state solution for this scenario, but since radial drift velocities are much higher than the diffusion velocities, this outward motion of drifting particles is highly ineffective. Solving for the stationary solution analogous to Pavlyuchenkov & Dullemond (2007; see also Jacquet et al. 2012) then defines the second part of our approximation, where

$$\begin{eqnarray}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}(r,a)}{{{\rm{\Sigma }}}_{{\rm{d}}}({r}_{0},a)}=\left\{\begin{array}{ll}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{{p}_{{\rm{d}}}} & \mathrm{for}\;r\leqslant {r}_{0}\\ \displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}(r)}{{{\rm{\Sigma }}}_{{\rm{g}}}({r}_{0})}\cdot \mathrm{exp}\left({\displaystyle \int }_{{r}_{0}}^{r}\displaystyle \frac{\mathrm{St}\;\gamma }{{\alpha }_{{\rm{t}}}}{{dr}}^{\prime }\right) & \mathrm{for}\;r\gt {r}_{0}\end{array}.\right.\end{eqnarray} \tag{ 8 }$$

Here $\mathrm{St}$ = a ${\rho }_{{\rm{s}}}$ π/2 ${{\rm{\Sigma }}}_{{\rm{g}}}$ denotes the Stokes number at the disk mid-plane.

Given a set of disk parameters, Equations (4)–(8) can be used to reconstruct the particle size distribution as a function of radius. The results for an illustrative example⁵ are shown in Figure 2. Panel (a) shows the full simulation after 1 Myr (using the Birnstiel et al. 2010 code) in grayscale. The framework explained above was used to reconstruct the size distribution, given ${{\rm{\Sigma }}}_{{\rm{g}}}$ and ${{\rm{\Sigma }}}_{{\rm{d}}}$ from the simulation. The results are overlaid as contour lines at the same levels; colors correspond to the respective regions in Figure 1. The key features are reproduced by this approximation: copious amounts of small grains reside inside ${r}_{{\rm{f}}}$ and extend to larger radii, corresponding to regions I and II, respectively. The densities for small grains in region II decrease with r (reflecting the limits of turbulent diffusion), while in region III they increase with r (due to longer growth timescales and smaller ${a}_{{\rm{d}}}$). This gives rise to a small-grain density minimum at intermediate radii (region IV). The radial segregation of particle sizes (${a}_{{\rm{d}}}$ decreases with r) means that the location of this minimum decreases with a.

Zoom In Zoom Out Reset image size

Panel (b) aims to connect this density minimum feature for small particles to a more observationally oriented diagnostic. It shows the radial distribution of (vertically integrated) optical depths for the simulated and reconstructed (approximated) models at wavelengths of 1.65 μm and 1.3 mm, where spatially resolved data are relatively common for protoplanetary disks. The reconstruction generally approximates the more detailed simulation remarkably well in terms of this key diagnostic.

To evaluate the observable consequences of this rather general feature in current models of dust evolution, we used the simulation outputs to calculate images and radial intensity profiles using the RADMC-3D radiative transfer code (Dullemond 2012). To do that, we assumed the gas density is in vertical hydrostatic equilibrium, with the dust density of each particle size distributed according to a mixing-settling equilibrium (e.g., Schräpler & Henning 2004). The vertical density structure was iterated based on the output thermal structure from RADMC-3D following the approach of Kama et al. (2009), allowing a maximum dust density deviation of 10% in each iteration to avoid numerical instabilities. We also employed the partial diffusion approximation and modified random walk procedures advocated by Min et al. (2009). We adopted the dust opacities prescription from Ricci et al. (2010), and assumed isotropic scattering for simplicity.

The raytraced images and (azimuthally averaged) intensity profiles at wavelengths of 1.65 μm (H-band) and 1.3 mm for the simulation of Figure 2 are shown in Figure 3. Notably, the significant reduction in the optical depth profile produces a depression in the scattered light intensity, starting roughly at ${r}_{{\rm{f}}}$ (in this case ∼7.4 AU) and reaching a minimum at about 50 AU. The tendency of this depression is to move inward with time because radial drift is depleting ${{\rm{\Sigma }}}_{{\rm{d}}}$, thus reducing $\epsilon$, which in turn decreases ${r}_{{\rm{f}}}$ according to Equation (5). The depression noted in the scattered light intensity resembles the recently observed "gap" features seen in the disks around TW Hya (Debes et al. 2013; Akiyama et al. 2015) and HD 169142 (Quanz et al. 2013). Such a feature is often linked to a planetary mass companion, which can induce a similar pattern in the scattered light morphology (e.g., de Juan Ovelar et al. 2013). However, in this Letter the dip is purely caused by opacity effects: the total dust surface density is continuous.

Zoom In Zoom Out Reset image size

In principle, one could distinguish between these two potential origins for such a feature with resolved data at an optically thin wavelength like 1.3 mm. There, the "gap" seen in scattered light will be absent if it is caused solely by dust evolution. Still, the intensity profile at 1.3 mm may also exhibit some more complex behavior outside of ${r}_{{\rm{f}}}$, where the lack of fragmentation leads to a reduced surface density of mm-sized particles. A second drop in mm-wave intensity at even larger radii may be caused by a lack of large particles in the outer disk. The locations and shapes of such features depends strongly on the choice of disk parameters: special care must be taken in their interpretation (the example presented here is not necessarily representative behavior at 1.3 mm).

We have implicitly assumed that a stationary distribution is always reached, and thereby neglected some time-dependent effects. A region devoid of small particles only occurs if the population of larger particles has had enough time to remove small grains by sweep-up collisions. Relating the sweep-up timescale (Birnstiel et al. 2012a, Equation (A.6)) to the age of the disk ${t}_{\mathrm{disk}}$, we can derive a condition for the global behavior described above with respect to the turbulence parameter,

$$\begin{eqnarray}&&{\alpha }_{{\rm{t}}}\gt \displaystyle \frac{1}{{\left({{\rm{\Omega }}}_{{\rm{k}}}\;{t}_{\mathrm{disk}}\right)}^{2}}\displaystyle \frac{1}{\epsilon \left|\gamma \right|}{\left(\displaystyle \frac{H}{r}\right)}^{-2},\end{eqnarray} \tag{ 9 }$$

where H = ${c}_{{\rm{s}}}$/${{\rm{\Omega }}}_{{\rm{k}}}$ is the pressure scale height and ${{\rm{\Omega }}}_{{\rm{k}}}$ the Keplerian frequency. Another criterion for the existence of a "region IV" comes from the assumption that both small and large grains are always in a vertical mixing-settling equilibrium distribution. Without efficient vertical mixing, the inward drifting large grains, located near the mid-plane, would not be able to sweep up the small grain population in the disk atmosphere. Forcing the vertical mixing timescale to be less than the disk lifetime also constrains ${\alpha }_{{\rm{t}}}$, to

$$\begin{eqnarray}&&{\alpha }_{{\rm{t}}}\gt \displaystyle \frac{1}{{t}_{\mathrm{disk}}{{\rm{\Omega }}}_{{\rm{k}}}}.\end{eqnarray} \tag{ 10 }$$

Furthermore, Equation (4) can be recast as a limit on ${\alpha }_{{\rm{t}}}$,

$$\begin{eqnarray}&&{\alpha }_{{\rm{t}}}\lt \displaystyle \frac{\left|\gamma \right|}{3\epsilon }{\left(\displaystyle \frac{{v}_{{\rm{f}}}}{{V}_{{\rm{k}}}}\right)}^{2},\end{eqnarray} \tag{ 11 }$$

which ensures that turbulence is too weak to cause fragmentation at a given radius. Figure 4 shows these three conditions for ${\alpha }_{{\rm{t}}}$ as function of radius for the parameters of the example simulation (cf. Figure 2). In this case, the sweep-up timescale is not relevant for this simulation. Also, the lower bound on ${\alpha }_{{\rm{t}}}$ imposed by Equation (10) is irrelevant throughout most of the disk. The dashed line in Figure 4 corresponds to the adopted ${\alpha }_{{\rm{t}}}$ in the simulation. The intersection of ${\alpha }_{{\rm{t}}}$ and the upper limit is by definition at ${r}_{{\rm{f}}}$ and corresponds to the onset of the dip. In cases where ${\alpha }_{{\rm{t}}}$ < 10⁻⁴, the regions depleted in small grains can reach down to a few AU, while for ${\alpha }_{{\rm{t}}}$ > 10⁻² fragmentation is active throughout most of the disk. Figure 4 also shows that the largest possible size of the gap is only weakly dependent on ${\alpha }_{{\rm{t}}}$.

Zoom In Zoom Out Reset image size