ACCRETING PLANETS AS DUST DAMS IN "TRANSITION" DISKS

Where, when, and how planets form is currently an unsolved problem in astrophysics. Central to this issue is how planets interact and sculpt the environment in which they are forming. The large number of observed exoplanets indicates that protoplanetary disks which churn out planets are the norm rather than a rare occurrence. However, connecting the properties of protoplanetary disks and exoplanets with a coherent picture of planet formation and evolution remains elusive.

Protoplanetary disks are made up of gas and dust particles; while the dust is only a minor component by mass (∼1%) it dominates the opacity of the disk material. Therefore, most observational indicators of protoplanetary disk evolution probe the evolution of the dust component rather than the gas which drives the dynamics. Protoplanetary disks appear to live for ∼1–10 Myr (Haisch et al. 2001; Hernández et al. 2007) until they are destroyed, most likely by photoevaporation (Clarke et al. 2001; Alexander et al. 2006, 2013; Owen et al. 2011, 2012). During their lifetime, disks are primarily optically thick out to large radii (≫10 AU), with accretion rates that decline with time in the range 10⁻⁶–10⁻¹⁰ M_☉ yr⁻¹ (e.g., Hartmann et al. 1998), and masses in the range 10⁻³–10⁻¹ M_☉ (e.g., Andrews & Williams 2005, 2007).

However, a small fraction of protoplanetary disks show a lack of opacity at NIR wavelengths, but return to levels comparable to those of primordial disks at MIR wavelengths (Strom et al. 1989; Skrutskie et al. 1990). The rarity of this population of protoplanetary disks, coupled with the fact that their inner disks appear to be cleared of dust, has led many authors to suggest that these disks may be caught in the act of evolving from a young star with a primordial disk to a disk-less star, and they have been aptly named "transition" disks (e.g., Kenyon & Hartmann 1995; Ercolano et al. 2011; Koepferl et al. 2013; Espaillat et al. 2014). Many mechanisms have been proposed in order to explain the properties of "transition" disks: tidal truncation by a companion (e.g., Calvet et al. 2005; Rice et al. 2006; Dodson-Robinson & Salyk 2011), photoevaporation (e.g., Clarke et al. 2001; Owen et al. 2011), grain-growth (e.g., Dullemond & Dominik 2005), photophoresis (Krauss et al. 2007), magneto-rotational instability (MRI-) driven winds (Suzuki & Inutsuka 2009), and magnetic winds (Armitage et al. 2013). Specifically, individual cases of these models can reproduce the observed characteristics of a given "transition" disk. However, no model or combination of models proposed can explain the observed population of "transition" disks. While photoevaporation will ultimately destroy the disk, clearing it from the inside out, (naturally creating a "transition" disk phase), it cannot explain the full population of observed "transition" disks (Alexander & Armitage 2009; Owen et al. 2011, 2012).

In fact, there is good evidence that observed "transition" disks do not represent a homogeneous population, but may contain several populations with different origins and lifetimes (e.g., Merín et al. 2010; Owen & Clarke 2012; Espaillat et al. 2014). Owen & Clarke (2012) showed that there are two distinct populations of "transition" disks, with different properties that can be separated by their millimeter flux. Owen & Clarke (2012) demonstrated there is a population of "transition" disks at very low millimeter fluxes (often the lowest millimeter fluxes of all protoplanetary disks); these disks have small (<10 AU) holes and low accretion rates (<10⁻⁹ M_☉ yr⁻¹). There is a second population of "transition" disks with high millimeter fluxes (often with the highest millimeter fluxes of all protoplanetary disks); these disks have large hole sizes (>10 AU) and high accretion rates (10⁻⁹–10⁻⁸ M_☉ yr⁻¹). The characteristics of the "transition" disks with low-millimeter fluxes are consistent with the concept of a disk in transition from primordial to cleared—if one assumes that millimeter flux is a proxy disk mass, which declines with time—as these disks appear to be entering a "transition" disk phase at the end of the their lifetimes. Furthermore, comparing this low-millimeter flux population alone to the X-ray photoevaporation scenario indicates good agreement between these "transition" disks and the model predictions (Owen et al. 2011, 2012; Owen & Clarke 2012).

However, it is the population of "transition" disks at high millimeter fluxes that still remains puzzling to understand. For photoevaporation to trigger disk clearing and create a transition disk, the accretion rate must drop below a threshold value (typically ∼10⁻⁹ M_☉ yr⁻¹) and the hole always develops around ∼1 AU (M_*/M_☉), meaning a "transition" disk with a hole at ∼20 AU and an accretion rate of 10⁻⁸ M_☉ yr⁻¹ is difficult to fit with a photoevaporative origin. Recently, it has been suggested that a significant dead zone (e.g., Morishima 2012; Bae et al. 2013) or an embedded planet (Rosotti et al. 2013) may trigger photoevaporation at higher accretion rates and larger hole sizes. Although, it still remains difficult to reconcile the required dust depletion in the inner disk with these ideas.

The original physical interpretation of a "transition" disk, inferred from the spectral energy distributions (SEDs)of an inner disk that is heavily dust depleted but still contains a significant gas reservoir (to power the observed accretion) while sharply switching to an outer disk with a significant mass in gas and dust, still remains today (e.g., Calvet et al. 2002, 2005). Many "transition" disks with large hole sizes have now been imaged in the sub-millimeter (e.g., Brown et al. 2009; Hughes et al. 2009; Andrews et al. 2011; van der Marel et al. 2013) with measured holes sizes that agree with those inferred from the SED, indicating the interpretation of a dust-depleted cavity is correct (at least for large dust particles with radii s ∼ 1 mm). Detailed SED modeling (e.g., Espaillat et al. 2007, 2008, 2010) and MIR imaging (Geers et al. 2007) have further confirmed that the inner disk¹ is also depleted of small dust particles, indicating the observed cavities must be dust-poor for particles with sizes s ≲ 1 mm. Furthermore, the radial transition from optically thin inner cavity to optically thick outer disk still remains unresolved, indicating the transition is sharp (<10 AU) and is inconsistent with the predictions of grain growth alone (Birnstiel et al. 2012). New ALMA observations of several "transition" disks have indicated that, as expected (from the measured accretion rates), the dust-poor cavities in "transition" disks contain a significant gas disk (van der Marel et al. 2013; Pérez et al. 2014) and in some cases the outer dust disk shows significant azimuthal asymmetries in the sub-millimeter images (van der Marel et al. 2013; Isella et al. 2013).

To explain the gas-rich, dust-poor, inner regions of "transition" disks requires that the small dust particles are depleted by around 10⁻⁴ from standard primordial values (e.g., Zhu et al. 2011). One of the most commonly invoked mechanisms is an embedded massive planetary companion, which creates a gap in the disk. Invoking a planetary companion is appealing since it naturally produces a leaky barrier to the gas (Calvet et al. 2005; Rice et al. 2006; Lubow & D'Angelo 2006; Dodson-Robinson & Salyk 2011; Gressel et al. 2013) allowing ongoing accretion, as well as driving azimuthal asymmetries in the outer disk such as the Rossby wave instability (e.g., Lovelace et al. 1999; Lin 2012; Lyra & Lin 2013). Furthermore, the pressure bump caused by a gap-opening planet is strong enough to trap dust particles with non-dimensionless stopping times (τ_s) of order unity (e.g., Rice et al. 2006; Pinilla et al. 2012; Zhu et al. 2012, 2013). Hole sizes of ∼20 AU and accretion rates of ∼10⁻⁸ M_☉ yr⁻¹ imply that particles with τ_s ∼ 1 are approximately 1 mm in size, indicating that millimeter-sized particles will easily be trapped by a planetary gap (Rice et al. 2006; Zhu et al. 2012) and can explain the observed millimeter images of "transition" disks (Pinilla et al. 2012). However, dust particles with smaller sizes (in particular, those that dominate the NIR opacity s < 1 μm) which have non-dimensional stopping times τ_s ∼ 10⁻³–10⁻⁴ (indicating they will be tightly coupled to the gas) will not be trapped in the pressure bump created by the planet, instead following the gas through the gap and into the inner disk. Detailed two-fluid simulations by Zhu et al. (2012) demonstrated this explicitly and concluded while the planetary scenario could explain the observed hole radii, accretion rates, and millimeter images, it would still have an optically thick inner disk (due to the small dust) and present with an SED consistent with that of a primordial disk. Thus, in order to rescue the planetary scenario one must include an additional effect (other than pressure trapping) in order to remove the small dust particles from the inner disk and thus explain all the observed features of a "transition" disk, including the SED. Furthermore, models that invoke multiple planets to carve out a large gap (e.g., Dodson-Robinson & Salyk 2011; Zhu et al. 2011) fail to explain all the characteristics. Zhu et al. (2011) showed that in the multiple planet scenario one could use planets to reduce the surface density in the inner disk, reproducing the observed SED and millimeter image, but failed to reproduce the observed accretion rates; modifying the model to match the observed accretion rates, Zhu et al. (2011) found that it could no longer reproduce the NIR dip in the SED.

A hint to the possible solution is that the observed accretion rates onto the star are high. This also implies a comparable accretion rate onto the embedded planet. In fact, this is backed up by the exciting discovery of a low-mass accreting object inside the gap of the "transition" disk HD 142527 (Close et al. 2014), which exhibits a significant accretion luminosity. As we will discuss, this additional accretion luminosity (from the forming planet) can, in certain cases, be the dominant force on small dust particles and provide the missing mechanism to remove the small dust particles from the inner disk. This article is organized as follows. In Section 2, we lay the theoretical basis for the role radiation pressure from an accreting planet may play on the dynamics of small dust particles. In Section 3, we develop a coupled one-dimensional (1D) secular gas and dust model and present the results of numerical calculations in Section 4. In Section 5, we use our results to compute synthetic observations, and discuss our results and model along with the caveats in Section 6, finally summarizing in Section 7. This article also includes two appendices: Appendix A discusses how to include planetary radiation pressure in a secular model and Appendix B covers the tests of the numerical method.

One of the most difficult things to explain about "transition" disks with large (>10 AU) cavities is their observed accretion rate, which is comparable to that of primordial disks (e.g., Espaillat et al. 2014). An appealing aspect of explaining "transition" disks with an embedded planet is that they naturally produce a leaky gap, which is known to trap millimeter-sized dust particles (Rice et al. 2006). The accretion rate through the gap and into the inner disk is comparable to the accretion rate into the gap, and by construction comparable to the accretion rate onto the planet (e.g., Lubow & D'Angelo 2006; Gressel et al. 2013). Following Lubow & D'Angelo (2006) we define the mass-loss rate onto the planet ($\dot{M}_p$) in terms of the accretion rate into the inner disk ($\dot{M}_{\rm inner}$), and therefore onto the star ($\dot{M}_*$) as

$$\begin{equation} \dot{M}_{p}=E\dot{M}_{\rm inner}. \end{equation} \tag{ 1 }$$

Conservation of mass across the gap implies that the accretion rate onto the planet and into the inner disk can be given in terms of the accretion rate into the gap from the outer disk ($\dot{M}_{\rm out}$) as

$$\begin{eqnarray} \dot{M}_p&=&\frac{E}{1+E}\dot{M}_{\rm out}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \dot{M}_{\rm inner}&=&\frac{1}{1+E}\dot{M}_{\rm out}. \end{eqnarray} \tag{ 3 }$$

Since simulations suggest E ∼ 3–10 (e.g., Lubow & D'Angelo 2006) and given the observed accretion rates onto the star in "transition" disks are $\dot{M}_*\sim 10^{-9}\hbox{--}10^{-8}$ M_☉ yr⁻¹, then a planetary origin for "transition" disks would also suggest $\dot{M}_p\gtrsim 10^{-9}\hbox{--}10^{-8}$ M_☉ yr⁻¹. In the following, we will make the assumption that "transition" disks with large cavities contain a massive planet (M_p ≳ M_J) and that E ≳ 1 as the simulations suggest.

Such a high accretion rate onto the planet will necessarily lead to a high accretion luminosity, which will be larger than the planet's intrinsic luminosity. For the scenario considered here this implies an accretion luminosity—estimated assuming the accreting material is free-falling onto the planet—of

$$\begin{eqnarray} L_p &=&\frac{GM_p\dot{M}_p}{R_p}=7\times 10^{-3}{\,L}_\odot \, E \left(\frac{M_p}{3\,{M_J}}\right)\nonumber \\ && \times \left(\frac{\dot{M}_*}{10^{-8}\,{\,M}_\odot \,{\rm yr}^{-1}}\right)\left(\frac{R_p}{10^{10}\, {\rm cm}}\right)^{-1}, \end{eqnarray} \tag{ 4 }$$

where M_p is the planet and R_p the planet radius. In deriving this expression we have made use of Equation (2) and implicitly assumed $\dot{M}_{\rm inner}=\dot{M}_*$ (which is true for a disk in steady state). Thus, if we compare the bolometric flux received at the gap edge (taken to be a Hill radius—R_H = a(M_p/3M_*)^1/3—from the planet, where a is the separation and M_* the stellar mass) compared to the star, assuming spherical dilution of the radiation we find

$$\begin{eqnarray} \frac{F^p}{F^*}&=& 0.7 E \left(\frac{M_p}{3\,{M_J}}\right)^{1/3}\left(\frac{\dot{M}_*}{10^{-8}\,{M}_\odot \,{\rm yr}^{-1}}\right)\left(\frac{R_p}{10^{10}\, {\rm cm}}\right)^{-1}\nonumber \\ && \times \left(\frac{L_*}{L_\odot }\right)^{-1}\left(\frac{M_*}{M_\odot }\right)^{2/3}, \end{eqnarray} \tag{ 5 }$$

where L_* is the star's bolometric luminosity. Therefore, at the gap edge we find that the accretion flux from the planet is comparable to, if not in excess of, the bolometric flux from the star. An obvious consequence of such a radiative flux is added feedback from radiation pressure. Since the dominant opacity source in protoplanetary disks is from dust particles, then this radiation pressure will act on the dust particles. The flux-mean opacity for an individual spherical dust particle is given by

$$\begin{equation} \kappa =\frac{3Q}{4\rho _ds}, \end{equation} \tag{ 6 }$$

where ρ_d is the dust particle density, s the dust particle radius, and Q is the radiative efficiency. For a perfect blackbody dust grain, Q = 1 when λ ≪ s and Q = (s/λ)² when λ ≫ s. We can estimate the magnitude of the radiation pressure due to a planet's accretion luminosity as a^rad = κF^p/c, where c is the speed of light and compare it to the other sources of acceleration on a dust particle. Therefore, to determine the appropriate value of Q we must estimate the radiation temperature of the accretion luminosity. In the viscous boundary layer model, magnetospheric accretion from a circumplanetary disk, or pure Bondi accretion scenario, we can roughly estimate the temperature of the radiation emerging from the accreting material at the planet's surface energetically as 3/2k_bT_acc = GM_pm_H/R_p (where k_b is the Boltzmann constant and m_H is the mass of a hydrogen atom), which yields T_acc ∼ 3 × 10⁵ K for the planets considered here. Wien's displacement law then gives an estimate of the photon wavelength of λ ∼ 0.01 μm. Ultimately, this radiation maybe reprocessed to longer wavelengths by extinction material; however, for the sake of simplicity, for the initial calculations presented here we assume the wavelength of the radiation is shorter than the size of the particles of interest. Thus, for all the dust particles, we set Q to unity with respect to the planetary accretion flux.

We now want to compare the magnitude of radiation pressure on a small (s ∼ 0.1 μm, ρ_d = 2 g cm⁻³) dust particle orbiting at the gap edge to the other forces that govern the dynamics of a dust particle. Evaluating the radiation pressure at the gap edge, we find a radiative acceleration of

$$\begin{eqnarray} a^{\rm rad}&=&3\times 10^{-3}{\, \rm cm\,\,s^{-1}}E\left(\frac{s}{0.1{\rm \, \mu m}}\right)^{-1}\left(\frac{M_p}{3{\,M_J}}\right)^{1/3} \left(\frac{M_*}{M_\odot }\right)^{2/3}\nonumber \\ && \times \left(\frac{\dot{M_*}}{10^{-8}{\,M_\odot \,{\rm yr}^{-1}}}\right) \left(\frac{R_p}{10^{10}{\,\rm cm}}\right)^{-1}\left(\frac{a}{20{\rm \, AU}}\right)^{-2} \end{eqnarray} \tag{ 7 }$$

compared with the radiation pressure from the star:

$$\begin{eqnarray} a^{\rm rad}_*&=&8\times 10^{-5}\,{\rm cm\,s}^{-2} \left(\frac{s}{0.1\,\mu {\rm m}}\right)\left(\frac{T_*}{4000{\rm \,K}}\right)^2\nonumber \\ &&\times \left(\frac{L_*}{{L_\odot }}\right)\left(\frac{a}{20{\, \rm AU}}\right)^{-2}, \end{eqnarray} \tag{ 8 }$$

where since the stellar blackbody peaks at λ > s we have used the perfect spherical blackbody approximation to suitably scale the radiative efficiency. The instantaneous radiative acceleration from the planet is large; thus, in order to assess when it might be dominant, we must compare it to the other main "acceleration" acting on the dust particle, namely, that due to drag from gas advection.

Adopting the Epstein regime, the acceleration on the dust particle due to dust advection is approximately

$$\begin{equation} a^{\rm adv}=\frac{\dot{M_*}c_s}{4\pi ^2 (H/R)a^2\rho _ds}, \end{equation} \tag{ 9 }$$

where c_s is the local sound speed, and H is the disk scale height at a radius R = a. Additionally, noting that both the acceleration due to dust advection and radiation pressure scale identically with dust density, accretion rate, and particle size, we can place a constraint on the value of E such that the radiation pressure is larger than the acceleration due to dust advection:

$$\begin{eqnarray} E &\gtrsim & \frac{1}{2}\left(\frac{H/R}{0.1}\right)^{-1}\left(\frac{a}{20\,{\rm AU}}\right)^{-1/4}\left(\frac{M_p}{3\,{M_J}}\right)^{1/3}\nonumber \\ && \times \left(\frac{R_p}{10^{10}\,{\rm cm}}\right)^{-1}\left(\frac{M_*}{M_\odot }\right)^{2/3}, \end{eqnarray} \tag{ 10 }$$

where since the disk is passively heated (e.g., Chiang & Goldreich 1997; D'Alessio et al. 2001), we adopt a T∝R^−1/2 temperature profile (e.g., Kenyon & Hartmann 1987). Thus, for the values of E typically found from simulations we see that the instantaneous radiation pressure from an accreting planet can dominate over the advection of small dust particles by a large factor. We strongly caution that this analysis is a very rough guide and ignores two important (and in practical considerations dominant) additional considerations. First, that the instantaneous radiation pressure only acts over part of the dust particles' orbit, while the dust drag due to advection acts over the entire orbit, necessarily weakening the role of radiation pressure. Second, the planetary gap will modify the surface density profile in the gap away from the simple steady-state ($\dot{M_*}=3\pi \nu \Sigma$) form we have used to estimate the advective drag force in Equation (9). This tends to weaken the strength of the advective drag force by reducing the coupling between the dust and gas through a reduction in the gas density close to the planetary gap (a very important consideration we will discuss in detail in Section 6). Ultimately, the competition between these two additional considerations determines whether the radiation pressure can stop the accretion of the small dust particles. This issue is the main aim of this work and we will answer this question through numerical simulations presented in Section 3. Finally, we note that for larger particle sizes, with non-dimensional stopping times τ_s ∼ 1, that dust drag due to the differential azimuthal velocity between the gas and dust becomes dominate over dust drag due to gas advection and radiation pressure. Thus, we expect that for dust particles with τ_s ∼ 1 that their dynamics will be governed by the pressure distribution of the gas disk.

Given the strength of the radiation pressure, this gives hope that for sufficiently high accretion rates, radiation pressure may be able to hold back the dust particles, and warrants further study in this work. Since the small dust particles are tightly coupled, they quickly transfer their excess momentum to the gas. Thus, one needs to check whether this will have a dynamical consequence on the gas itself. Adopting a dust particle number density distribution of the form n(s)ds∝s^−pds, the radiative acceleration on the gas will be given by

$$\begin{equation} a^{\rm rad}_{\rm gas}\approx \left\lbrace \begin{array}{@{}ll}Xa^{\rm rad}_{\rm dust}(s=0.1\,\mu {\rm m})\left(\frac{s_{\rm max}}{0.1{\rm \, \mu m}}\right)^{-1} & {p < 3}\\ Xa^{\rm rad}_{\rm dust}(s=0.1\,\mu {\rm m})\left(\frac{s_{\rm max}}{0.1{\rm \, \mu m}}\right)^{-1}\left(\frac{s_{\rm min}}{s_{\rm max}}\right) & {3 < p < 4}\\ Xa^{\rm rad}_{\rm dust}(s=0.1\,\mu {\rm m})\left(\frac{s_{\rm min}}{0.1{\rm \, \mu m}}\right)^{-1} & {p < 4} \end{array}\right.\!\!\!, \end{equation} \tag{ 11 }$$

where X is the dust-to-gas mass ratio and s_min and s_max are the minimum and maximum dust particle sizes, respectively. Therefore, we see for any dust distribution with p > 4 (namely the mass in the dust distribution is dominated at large particle sizes), radiation pressure will not inject significant momentum into the gas, and to first order the radiation pressure from the planet will only affect the dynamics of the dust distribution and not the gas provided the dust distribution has grown to sizes s_max ≫ 0.1 μm. Given protoplanetary disks are expected to have s_max ≳ 1 mm and p > 4 (e.g., Birnstiel et al. 2010; Pinilla et al. 2014) we can at this stage safely ignore the effect of the radiation pressure on the gas.

2.1. Consequences of Vertical Structure

Having seen that radiation pressure from an accreting planet may have a dynamical consequence on dust particles close to the gap edge, a concern is whether any circumplanetary disk may shield the outer regions of the protoplanetary disk from direct lines of sight to the planet. This can be considered by comparing the scale heights of any circumplanetary disk (H_cd) to the scale height of the protoplanetary disk (H). The protoplanetary disk will be passively heated by the central star at large radii (Chiang & Goldreich 1997; D'Alessio et al. 2001) so the scale height of the protoplanetary disk is given by (e.g., Kenyon & Hartmann 1987)

$$\begin{equation} \frac{H}{R}=0.04\left(\frac{R}{1{\rm \, AU}}\right)^{1/4}. \end{equation} \tag{ 12 }$$

In contrast, the circumplanetary disk is an active disk, heated by accretion, with a scale height of H_cd/R ≈ 0.3, which is relatively insensitive to the determining parameters in the range of interest (Martin & Lubow 2011). Furthermore, the circumplanetary disk is truncated near the orbit crossing radius of 0.4R_H (Martin & Lubow 2011). Thus, the height of the circumplanetary disk at its outer edge is

$$\begin{equation} H_{\rm cd}=0.2 {\rm \, AU}\left(\frac{a}{20{\rm \, AU}}\right)\left(\frac{M_p}{3{M}_{J}}\right)^{1/3}. \end{equation} \tag{ 13 }$$

Comparing this with the scale height of the protoplanetary disk at the orbit of the planet,

$$\begin{equation} H=1.7{\rm \, AU}\left(\frac{a}{20{\rm \, AU}}\right)^{5/4}, \end{equation} \tag{ 14 }$$

we see the circumplanetary disk is unable to shade the protoplanetary disk from direct lines of sight with the planet. Thus, the photons produced by accretion are able to directly impinge upon the edge of the protoplanetary disk.

In addition, it is well known that dust particles can sediment toward the mid-plane of the protoplanetary disk while conversely being lofted by turbulence, thus, it is important to check that the small grains will remain well mixed vertically in the disk. At a height of z ∼ H ≪ a the vertical components of gravity and radiation pressure are given by

$$\begin{eqnarray} F^{g}_z&=&-\frac{GM_*m_d}{a^3}z, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} F^{\rm rad}_z&=&\frac{m_d\kappa L_p}{4\pi c\left(R_H^2+z^2\right)^{3/2}}z. \end{eqnarray} \tag{ 16 }$$

Comparing gravity and radiation pressure to drag force vertically and assuming the dust particles are tightly coupled to the gas (hence adopting the terminal velocity approximation), we find a settling timescale of

$$\begin{equation} t_{\rm set}=\frac{\exp (-z^2/2H^2)}{\Omega \tau _s(1-\beta)}, \end{equation} \tag{ 17 }$$

where β represents the fractional reduction in vertical gravity due to radiation pressure given by

$$\begin{equation} \beta =\frac{\kappa L_p}{4\pi GcM_*[(M_p/3M_*)^{2/3}+(z/a)^2]^{3/2}}, \end{equation} \tag{ 18 }$$

Comparing this settling time to the turbulent lofting timescale t_loft ≈ z²/αc_sH, we can estimate the height z_dust of the dust layer in the disk by equating t_set and t_loft as

$$\begin{equation} \frac{\alpha }{2\tau _s(1-\beta)}=\left(\frac{z_{\rm dust}^2}{2H^2}\right)\exp \left(\frac{z_{\rm dust}^2}{2H^2}\right), \end{equation} \tag{ 19 }$$

which can be expressed in closed form as

$$\begin{equation} z_{\rm dust}=H \sqrt{2W_0\left(\frac{1}{2S(1-\beta)}\right)}, \end{equation} \tag{ 20 }$$

where S is the ratio of the viscosity to dimensionless stopping time (τ_s/α, e.g., Jacquet et al. 2012) and W₀ is the Lambert W function.² For the small particles we are interested in here S ≪ 1, z_dust > H. Thus, even without the added help of radiation pressure, the small dust is well mixed to many scale heights and is not settled into the mid-plane. Therefore, the photons produced by accretion are able to directly impinge upon the small dust particles at the edge of the disk. We note in passing that at large heights the small dust particles are no longer tightly coupled to the gas. Under certain circumstances, radiation pressure from the planet may be able to drive a "dust-particle wind" from several scale heights. Such a "wind" will have obvious implications for scattered light observations of "transition" disks such as those resulting from the SEEDs project (e.g., Dong et al. 2012).

In this section, we have set the ground work and seen that the radiation pressure from an accreting planet may have a dynamical consequence on the small dust particles at the outer edge of the planetary gap. In order to assess whether, for any sensible scenario, this can resolve the "transition" disk conundrum, we must turn to numerical calculations.

In order to model the "transition" disk problem, we must build a numerical model that allows us to assess whether the addition of planetary accretion luminosity can help explain the observed features. We choose to model the problem in a simple way, adopting a 1D radial model as has been used in previous studies of gas and dust disks with special regard to "transition" disks (e.g., Alexander & Armitage 2007, 2009; Zhu et al. 2011; Owen et al. 2011; Alexander & Pascucci 2012; Birnstiel et al. 2012; Pinilla et al. 2012).

3.1. Secular Models

The governing equations for the gas and dust disks with an embedded planet are given by (e.g., Clarke & Pringle 1988; Takeuchi & Lin 2002; Lodato & Clarke 2004; Alexander & Armitage 2007, 2009; Owen et al. 2011)

$$\begin{equation} \frac{\partial \Sigma _g}{\partial t}=\frac{3}{R}\frac{\partial }{\partial R}\left[R^{1/2}\frac{\partial }{\partial R}\left(R^{1/2}\nu \Sigma _g\right)-\frac{2\Lambda \Sigma _gR^{3/2}}{\sqrt{GM_*}}\right] \end{equation} \tag{ 21 }$$

and

$$\begin{equation} \frac{\partial \Sigma _d^i}{\partial t}=-\frac{1}{R}\frac{\partial }{\partial R}\left[R\Sigma _d^iv_d^i-\frac{\nu }{{\rm Pr}}{ R}\Sigma _g\frac{\partial }{\partial R}\left(\frac{\Sigma _d^i}{\Sigma _g}\right)\right], \end{equation} \tag{ 22 }$$

where the superscript i refers to a dust particle of size sⁱ, ν is the turbulent viscosity, and Pr is the Prandtl number describing the ratio of the turbulent viscosity to the dust diffusion due to the turbulence. Λ is the torque resulting from the planet, where we use the symmetric form used by Trilling et al. (1998), Armitage et al. (2002), Lodato & Clarke (2004), Alexander & Armitage (2009), and Alexander & Pascucci (2012); for a planet-to-star mass ratio of q, Λ is given by

$$\begin{equation} \Lambda = \left\lbrace \begin{array}{@{}ll}-\frac{q^2GM_*}{2R}\left(\frac{R}{\max (H,|R-a|)}\right)^4 & {{\rm if}\; R < a } \\ \frac{q^2GM_*}{2R}\left(\frac{a}{\max (H,|R-a|)}\right)^4 & {{\rm if}\; R > a}. \end{array}\right. \end{equation} \tag{ 23 }$$

Since protoplanetary disks are primarily passively heated (e.g., Chiang & Goldreich, 1997; D'alessio et al. 2001), a constant "α" viscosity requires ν∝R (e.g., Hartmann et al. 1998; Alexander et al. 2006; Owen et al. 2011). Therefore, we set ν = ν₀R where ν₀ = α(H/R)²ΩR² in all our simulations. In all calculations, we set α = 0.0065; H/R is normalized to 0.04 at 1 AU and we choose to set the Prandtl number to unity (e.g., Clarke & Pringle 1988; Alexander & Armitage 2007). The radiation pressure on the dust is included through its effect on the dust velocity via Equation (27).

In order to determine $v_d^i$, one needs to account for the additional impact on the dust due to the radiation pressure. The equations governing the evolution of a dust particle in a gas disk including radiation pressure³ are

$$\begin{equation} \frac{dv_R}{dt}=\frac{v_\phi ^2}{R}-\Omega ^2R-\frac{1}{t_s}\left(v_R-u_R\right)+a_R^{\rm rad} \end{equation} \tag{ 24 }$$

and

$$\begin{equation} \frac{dRv_\phi }{dt}=-\frac{R}{t_s} (v_\phi -u_\phi)+Ra_\phi ^{\rm rad}, \end{equation} \tag{ 25 }$$

where t_s is the stopping time given by t_s = τ_s/Ω. Since we are primarily interested in tightly coupled particles τ_s < 1, then the net effect of the planetary radiation pressure will be to impart an impulse on the dust particle every orbit, knocking it onto a slightly different orbit. Therefore, we adopt a 1D orbit-averaged (we are interested in the secular evolution) approach and can average the impulse imparted by the radiation pressure from the planet over an entire orbit, and we use Equations (24) and (25) to calculate the secular evolution of the dust particles. We do this by replacing a^rad by its orbit-averaged value 〈a^rad〉 in Equations (24) and (25); noting by symmetry $\langle a^{\rm rad}_\phi \rangle =0$ we can then proceed and solve for v_R. Following standard derivations (e.g., Takeuchi & Lin 2002; Armitage 2010) we can include the radiation pressure term. One takes the dust particles to move through a succession of circular orbits, so we may simplify the azimuthal equation to

$$\begin{equation} v_\phi -u_\phi =-\frac{1}{2}t_sv_R\Omega. \end{equation} \tag{ 26 }$$

Thus, solving for the radial dust velocity we find

$$\begin{equation} v_R=\frac{u_R\tau _s^{-1}-\eta \Omega R + \left\langle a^{\rm rad}_R \right\rangle /\Omega }{\tau _s+\tau _s^{-1}}, \end{equation} \tag{ 27 }$$

where η is a measure of the gas pressure gradient given by

$$\begin{equation} \eta =-\frac{{d}\log P}{{d}\log R}\left(\frac{c_s}{v_K}\right)^2. \end{equation} \tag{ 28 }$$

We discuss in Appendix A how to perform the orbit averaging of the planetary radiation pressure and we evaluate $\langle a^{\rm rad}_R\rangle$ and also present a closed form solution in the optically thin limit.

3.2. Numerical Method

Operationally, Equations (21) and (22) are integrated explicitly using a scheme that is second order in space and first order in time; the flux terms are reconstructed using a van Leer limiter. Furthermore, to make the timestep numerically manageable we follow Lodato & Clarke (2004) and Alexander & Armitage (2009) and smooth both the planetary torque and the radiation pressure inside the planet's Hill sphere and we do not attempt to model the disk properties in this region (essentially, in the regions where the flow is no longer 1D). In order to model the flow across the planetary gap, and onto the planet, we adopt a "leakage" prescription similar to those previously used in the literature (e.g., Alexander & Armitage 2009; Alexander 2012; Alexander & Pascucci 2012). Since the "leakage" is applied outside the planet's Hill sphere and the smoothing takes place inside, the nature of this smoothing does not affect the calculations (Alexander & Armitage 2009). The "leakage" is included as a sink and source term outside and inside the planet's Hill sphere that moves dust and gas from one side of the planet's orbit to the other. In order to determine the gas leakage rate, we measure the steady-state viscous rate ($\dot{M}_{\rm out}$) at 3a (Alexander & Armitage 2009). We then assume that the mass flux across the gap and into the inner disk is given by Equation (3) and the accretion rate onto the planet is given by Equation (2). For the dust we weight the gas leakage rate ($\dot{M}_{\rm out}$) by the dust concentration at the gap edge, and use the same value of E to decide how much dust is accreted by the planet. It is also necessary to apply a timestep limit to the "leakage" source and sink terms so that they do not account for a >1% change in the surface density. To determine the orbit-averaged radiation pressure on the dust we proceed as follows: at the beginning of each timestep we determine the accretion rate onto the planet and calculate the accretion luminosity using Equation (4). For each grid cell, we determine the optical depth on each cell face by integrating the radiative transfer equation as (assuming the dust is vertically well mixed)

$$\begin{equation} \tau ^i_{\rm face}=\sum _{j=k_p+1}^{i-1}\frac{\Sigma _g^j}{\sqrt{2\pi }H_j}\sigma ^{j-1}\Delta R^{j-1}, \end{equation} \tag{ 29 }$$

where k_p is the position of the planet, σ^j is the cell-centered cross-section of the dust distribution, and ΔR^j is the radial cell size. We calculate σ^j independently for each cell assuming each dust particle to have a geometric cross-section, the local dust size distribution, and the local dust-to-gas mass ratio. We then evaluate the orbit averaging integral (Equation (A2)) numerically in order to determine the orbit-averaged radiation pressure.

Furthermore, in order to include sufficient resolution around the planet we make use of static mesh refinement, where we increase the number of grid cells, by splitting each cell of the "mother" grid into 10 extra cells within 3 Hill radii of the planet. The "mother" grid is logarithmically spaced between 0.02 and 200 AU and composed of 400 cells, which results in approximately 40–50 cells per Hill radii in the vicinity of the planet. Since this high resolution requirement makes an evolutionary calculation computationally infeasible, we model a "quasi-steady-state" version of the problem. We do this by neglecting the planet's migration and evolving the model until a steady state is achieved in the gas distribution. This is done by using inflow boundary conditions at the outer boundary, where gas and dust are injected into the grid at a constant rate. At the inner boundary, we apply a zero torque boundary condition and set Σ = 0. In all cases, the gas initially has a Σ∝R⁻¹ profile (appropriately scaled for the incoming accretion at the outer boundary) and no dust is present. We then evolve the gas only for 1 Myr with the planet to obtain a steady state. After 1 Myr we introduce dust at the outer boundary, allowing the simulation to evolve for a further 1 Myr.

While a steady state is always achieved for the gas, this will not be the case for trapped dust particles, where a "quasi-steady-state" is achieved after 1 Myr. The concentration of dust particles in the inner disk slowly increases as the concentration of dust trapped at the outer edge grows, increasing the concentration gradient at the gap. This is a slow diffusive process and all our measurements are made at 2 Myr. Allowing the simulation to run for several 100 Myr, a steady state (where the rate dust particles entering the grid at the outer boundary equals the rate exiting at the inner boundary, excluding the planet's accretion) can be reached; however, such a situation results in an unphysically large dust concentration (∼100–1000) at the gap edges. The tests we used to ensure the numerical scheme is behaving accurately are discussed in Appendix B.

In all models we adopt a stellar mass of 1 M_☉, the dust distribution injected at the outer boundary has a dust-to-gas mass ratio of 0.01, with an MRN power-law distribution (p = 3.5; Mathis et al. 1977), a minimum particle size of s_min = 0.005 μm, and a maximum particle size of s_max = 1 mm. Such a choice is a reasonable starting point and is often the grain size chosen for radiative transfer modeling (e.g., D'Alessio et al. 2001; Koepferl et al. 2013). It is well known that the dust distribution can evolve away from the MRN distribution (Birnstiel et al. 2010, 2012). As such, if the particle distribution becomes strongly peaked toward small particle sizes it will weaken the effect of radiation pressure; however, if the particle distribution becomes strongly dominated by large particle sizes, it will strengthen the role of radiation pressure. The dust distribution is followed using 25 size bins, logarithmically spaced between s_min and s_max. We note that we do not include dust growth or fragmentation (a caveat we discuss in Section 6) and the dust distribution evolves only due to radial motion/trapping of individual dust species; therefore, our 25 size bins evolve independently. Finally we pick the median value of E from Lubow & D'Angelo (2006) of E = 6 which is comparable to the value used in similar 1D models of planet–disk interaction (Alexander & Armitage 2009; Alexander & Pascucci 2012).

We calculate several models that have parameters comparable to observed "transition" disks, but unfortunately the high numerical overhead⁴ restricts this parameter range to a small subset of the parameter space. The simulated parameters are shown in Table 1. Taking $\dot{M}_*=10^{-8}$ M_☉ yr⁻¹, we simulate the evolution of gas and dust in a protoplanetary disk with a planet. Placing an embedded planet on a fixed circular orbit at 20 AU with various masses, we follow the evolution of the gas and dust. For our "standard" model, we consider a planet with a mass of 4.0 M_J, a radius of 10¹⁰ cm, and do not include the radiation pressure feedback (simulation A).

Table 1. Simulation Parameters

Simulation	Mp	Rp	$\langle a_R^{\rm rad}\rangle$ (Y/N)	$\dot{M}$	a	E
	(MJ)	(cm)		(M☉ yr−1)	(AU)
A	4.0	1010	N	10−8	20	6
B	4.0	1010	Y	10−8	20	6
C	2.0	1010	Y	10−8	20	6
D	2.5	1010	Y	10−8	20	6
E	3.0	1010	Y	10−8	20	6
F	3.5	1010	Y	10−8	20	6
G	4.5	1010	Y	10−8	20	6
H	4.0	1010	Y	10−8	20	0.5

Note. $\langle a_R^{\rm rad}\rangle$ (Y/N) refers to whether radiation pressure from the planet was included in the simulation.

Download table as:  ASCIITypeset image

In Figure 1 we show the results of simulation A. This shows that the pressure gradient induced by the planet can trap large millimeter-sized dust particles, in agreement with previous studies (e.g., Rice et al. 2006; Zhu et al. 2012, 2013; Pinilla et al. 2012). However, the pressure gradient is unable to trap small sub-micron dust particles that freely flow across the gap and into the inner disk, essentially tracing the gas distribution. This simulation shows results similar to those presented by Zhu et al. (2012) using multi-dimensional two-fluid simulations. Essentially, the time to advect small dust particles into and across the gap is much shorter than the timescale on which they will feel the presence of the pressure bump and migrate toward it.

Zoom In Zoom Out Reset image size

In simulation B, we repeat simulation A (M_p = 4.0 M_J), but this time include the effect of radiation pressure. The resulting gas and dust surface densities are shown in Figure 2; here we see that, as before, the millimeter particles are trapped in the outer disk to a similar level to that found in simulation A (without radiation pressure). This means, as in simulation A, the dynamics of the millimeter-sized particles are dominated by the pressure bump outside the planet, and the inclusion of radiation pressure has little effect on their dynamics. However, Figure 2 shows that the small sub-micron-sized dust particles are significantly suppressed in the inner disk (we find a suppression factor of <10⁻⁴). Here, the radiation pressure traps the dust particles at the outer edge of the gap, and the resulting concentration in the inner disk is set by a balance of advection, diffusion, and radiation pressure at the gap edge.

Zoom In Zoom Out Reset image size

Furthermore, to assess whether the planet needs to accrete a significant fraction of the gas entering the gap, we repeat simulation B, but with E = 0.5 in simulation H. The resulting surface density profiles are shown in Figure 3. With this lower value of E, we find that, as in the case with no radiation pressure, the millimeter-sized dust particles are trapped by the planetary pressure bump and the sub-micron grains still make it across the gap. Comparing simulations A, B, and H we find that while the sub-micron grains are suppressed in the inner disk for E = 0.5 compared to the model with no radiation pressure (simulation A), the level of trapping for the sub-millimeter grains is not enough for the disk to give rise to an NIR dip in the SED and would still be classified as a primordial disk. Thus, we require the planet to be accreting the majority (E ≳ 1) of the gas entering the gap for the disk to appear as a "transition" disk.

Zoom In Zoom Out Reset image size

4.1. Variations in Planetary Mass

Finally, to conclude the initial study of the simulations presented in this work, we investigate the role of the planet's mass (simulations B–G). We consider a range of planets between 2.0 and 4.5 M_J all with radii of 10¹⁰ cm, on circular orbits at 20 AU. The resulting surface density profiles for the sub-micron dust are shown in Figure 4 for simulations B–G.

Zoom In Zoom Out Reset image size

These simulations clearly show for planet's with masses above ∼4 M_J, the dynamics of small sub-micron-sized dust particles are governed by radiation pressure from the accreting planet and the inner disk becomes optically thin at NIR wavelengths, ultimately giving rise to the signature of a transition disk in the SED for an accretion rate of $\dot{M}=10^{-8}$ M_☉ yr⁻¹. As in the comparisons between simulations A and B we find that in all cases the dynamics of the larger millimeter-sized particles are governed by the properties of the pressure bump and large planets trap more millimeter-sized grains; however, even planets with masses ∼2 M_J are able to trap enough millimeter grains to create a cavity in the millimeter images, as we will see in Section 6. While the radiation pressure at the gap edge does increase with planet mass (Equation (4)) this is not the major reason why the planet mass has a strong effect on the level of small dust particles entering the inner disk. An important fact that we could not include in our simple discussion in Section 2 is that bigger planets carve deeper gaps. This has two effects: (1) the optical depth at the gap edge where the dust is trapped drops, and (2) the reduction in surface density reduces the dust–gas coupling, making it easier for the radiation pressure to overcome the gas advection. These processes are discussed further in Section 6.

In order to compare the simulations with current observations, we calculate both the disk's SED and millimeter image. Since our simulations give us the dust distribution as a function of radius, we can use it to calculate the opacity to the incoming stellar and re-radiated thermal emissions. Unlike the calculation of the radiation pressure term we cannot simplify the opacity of the dust particles and thus we use tabulated values from Laor & Draine (1993), assuming a 50/50 mixture of graphite and silicate grains.

5.1. SED Calculation

In all cases, we assume the disk is observed face-on and vertically isothermal. Since the goal of this section is to calculate representative SEDs, rather than accurate models for SED fitting, we do not use a full numerical radiative transfer approach, but rather estimate the disk's radial temperature profile and, hence, brightness analytically. We calculate a disk temperature, assuming both an active disk (heated only by accretion; Shakura & Sunyaev 1973) and passive disk (heated only by stellar irradiation; Chiang & Goldreich 1997), and choose the maximum of the two temperatures to be the disk's temperature. We also apply a minimum disk temperature of 10 K and a dust sublimation temperature of 1500 K. The passive temperature profile is taken to be

$$\begin{eqnarray} T(R)&=&290\,{\rm K}\left(\frac{R}{1{\,\rm AU}}\right)^{-1/2}\exp (-\tau _R^*)\nonumber \\ && + 200\,{\rm K}\left(\frac{R}{1\,{\rm AU}}\right)^{-1/2} \left[1-\exp \left(-\tau _R^* \right)\right], \end{eqnarray} \tag{ 30 }$$

where $\tau _R^*$ is the radial optical depth to stellar radiation used to smoothly move between the radially optically thin and thick limits. Thus, the luminosity of a face-on thin disk is calculated as

$$\begin{equation} L_\nu =\pi \int _0^{\infty }{d}R\, 2\pi RB_{\nu }(T(R))[1-\exp (-\tau _\nu)], \end{equation} \tag{ 31 }$$

where τ_ν is the vertical optical depth at frequency ν as a function of radius, and B_ν is the Planck function.

5.2. Millimeter Images

Most of the millimeter images to date have been imaged by the Submillimeter Array at 880 μm (e.g., Andrews et al. 2011), therefore we calculate our images at 880 μm. We use the same radial temperature distribution used to calculate the SEDs and use our dust distribution to calculate the opacity at 880 μm, and the associated vertical optical depth. The source function for a face-on, thin disk is then approximately

$$\begin{equation} S_\nu =B_\nu [1-\exp (-\tau _{880})]. \end{equation} \tag{ 32 }$$

In order to calculate the image we assume the dust distribution to be axisymmetric at all particle sizes. Furthermore, since the gap features occur on smaller scales than typical observational resolutions, we further degrade the image by convolving it with an axisymmetric Gaussian beam with a standard deviation of 5 AU.

5.3. Resulting Observational Properties

We calculate the synthetic observations for simulations A–C. These simulations contain a 4.0 M_j (simulations A and B) or 2.0 M_J (simulation C) planet orbiting at 20 AU and the accretion rate onto the star is 10⁻⁸ M_☉ yr⁻¹. Simulation A contains no radiative feedback from the planet, while simulations B and C do contain radiative feedback. The SEDs for simulation A (solid line), simulation B (dashed), and simulation C (dot-dashed) are shown in Figure 5, compared to the stellar spectrum (dotted). Simulation A (without radiative feedback) clearly shows an SED that would be classified as a primordial disk based on its profile and NIR/MIR colors. However, simulation B (with radiative feedback) shows a lack of emission at NIR wavelengths and a strong MIR emission bump. Furthermore, the 10 μm silicate feature is visible indicating the presence of optically thin dust in the cavity. The SED of simulation B looks characteristically like those SEDs seen in "transition" disks and can only be created by removing the small sub-micron dust from the inner disk. Even though simulation C does include radiative feedback, the planet mass (2.0 M_J) is not massive enough to make trapping efficient enough that the inner disk remains optically thick at NIR wavelengths giving rise to a primordial SED that is very similar to that of simulation A.

Zoom In Zoom Out Reset image size

The 880 μm images are shown in Figure 6 for simulation A (center), B (right), and C (left). As expected from the millimeter-sized particle distributions, all disks show evidence for a large cavity in the millimeter images. We find that the peak of the millimeter emission occurs at roughly 1.5 × the planet's separation (∼30 AU) in good agreement with the model of Pinilla et al. (2012). Simulation C therefore represents a possible new observational class of "transition" disk: one that shows a primordial SED, but a large millimeter cavity. Two such disks (WSB 60, MWC 758) were serendipitously discovered by Andrews et al. (2011), and we will discuss the implications of this class further in Section 6.

Zoom In Zoom Out Reset image size

Therefore, in order to reproduce both the SEDs and millimeter images of "transition" disks, we find it is necessary to include radiative feedback from an accreting planet. However, should this radiative feedback fail due to too low a planet mass <3 M_J or too low a planetary accretion rate, one would expect to see observational diagnostics similar to simulation A (a millimeter cavity but a primordial SED).

Dust particles with τ_s ∼ 1, which have sizes of ∼1 mm at hole radii found in "transition" disks, are trapped by pressure gradients outside the planet's orbit, consistent with the images of dust at millimeter wavelengths (Brown et al. 2009; Hughes et al. 2009; Andrews et al. 2011; van der Marel et al. 2013). Small dust particles (≲ 1 μm), which dominate the opacity at IR wavelengths and ultimately whose absence give rise to the classification as a "transition" disk, have τ_s ≲ 10⁻³. Our calculations and those of other works, e.g., Zhu et al. (2012), have shown that without radiation pressure feedback, small particles follow the gas across the gap and give observation signatures of a primordial disk at NIR wavelengths rather than a transition disk (Ward 2009). This was conclusively demonstrated using two-fluid simulations (Zhu et al. 2012), which showed that an embedded planet alone could not prevent the sub-micron dust particles from entering the inner disk. Thus, with only dust trapping due to pressure gradients, small sub-micron particles are able to cross the planet's gap entrained in the gas and give rise to an observation of a primordial disk.

Including radiation pressure from accretion onto an embedded planet allows, for the first time, the observational demographics (SED, millimeter image, and accretion rate) of "transition" disks to be reproduced. In this model, radiation pressure holds small dust particles outside the planet's orbit, while the pressure gradient from the planet's gap traps the larger, millimeter-sized particles such that dust-depleted gas proceeds to flow into the inner disk. Ultimately, the amount of small sub-micron-sized dust particles allowed to flow across the gap is set by a balance between advection, turbulent diffusion, and radiation pressure, whereas for the millimeter-sized particles, it is set by the balance between turbulent diffusion and radial drift due to the pressure gradient.

It is at this point we must consider the difference between this scenario and that presented by Chiang & Murray-Clay (2007), who suggested radiation pressure from the star may trap dust in an MRI active layer at large radii (≳ 10 AU). Ultimately, Dominik & Dullemond (2011), showed that (irrespective of whether you could ever set up such a system) the dust wall was not stable, and the dust always overcame the stellar radiation pressure on a timescale shorter than the disk's lifetime. They found that dust advection behind the radiation pressure-supported layer (whose width was set by the attenuation scale of the stellar radiation) built up a large enough concentration gradient such that diffusion was able to overcome the radiation pressure and drive the radiation supported layer inward. At first glance, the ratio of the advective acceleration to the radiation pressure in our scenario (Equation (10)) is not much larger than that adopted by Dominik & Dullemond (2011). Thus, we must wonder where the difference lies, specifically, the reason why radiation pressure from the star is unable to hold back the dust in a standard disk, but radiation pressure from an accreting, gap-creating planet is.

The answer lies in the fact that the planet creates a deep gap reducing the surface density significantly within distances of order (R − a) ≲ 5R_H, due to the planetary torque. This reduction in surface density has several effects: (1) the reduction in surface density reduces the coupling between the gas and dust, and hence the advective drag force (∝Σ_g), meaning that, in reality, at the gap edge, Equation (10) is an (often severe) underestimate of the ratio of the radiation pressure to the advective drag force. Since this depends on the reduction in the surface density then it is obviously sensitive to the planet mass and this is one of the main causes of the planetary mass dependence found in the results rather than the increase in accretion luminosity. (2) The planet keeps the density in the gap low enough that the radiation from the planet remains optically thin until the radiation pressure and other forces become comparable, meaning the width of the radiation-pressure-supported layer is not set by sharp attenuation of the radiation field but rather by the properties of the gap (which is of a much large scale than the attenuation scale found in Dominik & Dullemond 2011). This allows the dust to extract more integrated momentum from the radiation field, while the increased width of the radiation-supported layer lowers the concentration gradient, weakening the diffusion that ultimately overcame the radiation field in Dominik & Dullemond's case. (3) The torque due to the planet acts to counter-balance the advective force of the gas flow, thus the "snow-plow" formed in Dominik & Dullemond's calculation is considerably wider, weaker, and forms on a much longer timescale, meaning turbulent diffusion is unable to overcome radiation pressure. Finally, we note that Dominik & Dullemond (2011) assumed a different (and somewhat unrealistic) dust distribution than that adopted in our calculation; they took all the dust to be of a single size of s = 0.1 μm. While the opacity of an individual dust grain is similar in both calculations, the attenuation of the radiation field in the Dominik & Dullemond (2011) calculation is much stronger than in our calculations. Our calculations use a more realistic dust distribution that has grown to approximately millimeter sizes, as observed in protoplanetary disks (e.g., Birnstiel et al. 2010; van der Marel et al. 2013; Pinilla et al. 2014), rather than one dominated by small (hence high-opacity) dust particles. Thus, in our calculations, we find that, unlike Dominik & Dullemond's (2011) problem, in the planetary hypothesis for "transition" disks, radiation pressure is able to overcome turbulent diffusion and trap small dust particles at the outer edge of the gap.

6.1. Caveats and Limitations

6.2. Possible Evolutionary Scenario for a "Transition" Disk

In this work, we have presented an update to the standard planetary hypothesis for the origin of "transition" disks by including feedback on the small dust particles from radiation pressure due to an accreting, embedded planet. By adopting the standard picture of planet–disk interaction for massive ≳ M_J planets, which carve deep gaps, while allowing ongoing accretion onto the star, we find that the observed accretion rates in "transition" disks (that cannot be explained by photoevaporation) imply a high accretion luminosity originating from the forming planet ≳ 10⁻³ L_☉. At the gap edge, radiation pressure from the planetary accretion luminosity can be the dominant force on small (≲ 1 μm) dust particles.

Using a simple 1D secular coupled gas and dust model for the disk, we find, in agreement with previous studies, that, without radiative feedback, the planet cannot explain all the observed features of "transition" disks. Without radiative feedback, massive planets can open deep enough gaps to trap the millimeter-sized dust particles while allowing the small dust particles to follow the gas into the inner disk, giving rise to a primordial SED but millimeter cavity. However, by including radiative feedback from the accreting planet, we find that above a planet mass threshold of ∼3–4 M_J radiation pressure is able to hold back the small s ≲ 1 μm dust particles, allowing dust-free gas to accrete across the gap and into the inner disk. We require that the planets accrete at least half of the material flowing into the gap in order to trap sufficient sub-micron dust particles. By computing synthetic SEDs and millimeter images, we find we are able to explain the observed NIR dip, millimeter cavity, and accretion rate within a single model.

The fact that this process possesses a planetary mass threshold suggests that "transition" disks with large holes and high accretion rates are not in fact disks rapidly transitioning from a primordial to diskless state, but rather, a rare and relatively long-lived state (0.5–1 Myr). Assuming a scenario where planets accrete and grow above the required mass threshold to appear as a "transition" disk, allows us to construct an evolutionary scenario where a planet-hosting disk would first show a millimeter cavity but primordial SED while the planet mass is low. Once it grows above the mass threshold it presents as a "transition" disk with a millimeter cavity and dip in the SED at NIR wavelengths. Finally, once the accretion rate in the disk (and hence onto the star and planet) drops or onto the planet the small dust then refills the inner cavity, again giving rise to a disk that has a millimeter cavity but primordial SED, before photoevaporation finally disperses the disk.

Finally, this model suggests that "transition" disks with large holes and high accretion rates should all have heavily accreting planetary objects. These planets would be detectable at close to face-on inclinations using AO imaging at Hα, as demonstrated by the detection of an accreting planetary object in the "transition" disk HD 142527 using the Magellan Adaptive Optics VisAO camera.

J.E.O. thanks the referee for comments that improved the manuscript. J.E.O. is grateful to Cathie Clarke, Ruobing Dong, Barbara Ercolano, Giovanni Rosotti, Ilaria Pascucci, Yanqin Wu, and Zhaohuan Zhu for helpful discussions. The calculations were performed on the Sunnyvale cluster at CITA which is funded by the Canada Foundation for Innovation.

The instantaneous acceleration on a dust particle due to radiation pressure is given by

$$\begin{equation} {a}^{\rm rad}=\frac{\kappa {F}^p}{c}. \end{equation} \tag{ A1 }$$

However, since we are interested in the long-term evolution of small dust particles, which are tightly coupled to the gas (τ_s ≪ 1), they will move slowly through a succession of Keplerian circular orbits (Takeuchi & Lin 2002). Therefore, we can consider this slow secular evolution by averaging the radiation pressure over an entire orbit. We describe the geometry of this setup in Figure 7, where we have transformed into a frame co-rotating with the planet. Considering an orbit-averaged radiation pressure

$$\begin{eqnarray} \left<a^{{\rm rad}}_R\right>&=&\frac{\frac{\kappa }{c}\int _0^TF^p_R{d}t}{T}, \end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray} \left<a^{{\rm rad}}_\phi \right> & = &\frac{\frac{\kappa }{c}\int _0^TF^p_\phi {d}t}{T}=0, \end{eqnarray} \tag{ A3 }$$

where T is the orbital period at a radius R in the co-rotating frame, i.e., $T=2\pi [\Omega (R)^{-1}-\Omega _p(a)^{-1}]=2\pi /\tilde{\Omega }$. Clearly the orbit azimuthal component is zero by symmetry. The radial component of the planet's flux in the radial direction is given by

$$\begin{equation} F^p_R=\frac{L_p\exp \left(-\frac{\tau _R}{\cos \delta }\right)}{4\pi d^2}\cos \delta, \end{equation} \tag{ A4 }$$

where τ_R is the radial (mid-plane) optical depth from the planet to a cylindrical radius, R (i.e., $\tau _R=\int _a^R\kappa \rho {d}R$) and we have used attenuation in slab geometry to include optical depth effects (valid provided d < a). One can find an analytic solution in the optically thin limit by proceeding as follows: defining r = R − a as the cylindrical radial distance between the disk material at radius, R, and the planet, and d is the distance between the planet and a region of the disk at a given R, ϕ position, then, using the rule of cosines, one can show

$$\begin{equation} d=\sqrt{a^2+(a+r)^2-2a(a+r)\cos \phi } \end{equation} \tag{ A5 }$$

and

$$\begin{equation} \cos \delta =\frac{a+r-a\cos \phi }{ \sqrt{a^2+(a+r)^2-2a(a+r)\cos \phi }}. \end{equation} \tag{ A6 }$$

Thus, the equation for the orbit-averaged acceleration due to radiation pressure is

$$\begin{equation} \left<a^{\rm rad}_R\right>=\frac{2\kappa L_p}{4\pi c a^2 T}\int _0^{t_{\rm crit}}\frac{1+\tilde{r} -\cos \phi }{[1+(1+\tilde{r})^2-2(1+\tilde{r})\cos \phi ]^{3/2}}{d}t, \end{equation} \tag{ A7 }$$

where $\tilde{r}=r/a$. Now, replacing dt with ${d}\phi /\tilde{\Omega }$, and using $T=2\pi /\tilde{\Omega }^{-1}$, the equation simply becomes

$$\begin{equation} \left<a^{\rm rad}_R\right>=\frac{2\kappa L_p}{8\pi ^2 c a^2}\int _0^{\phi _{\rm crit}}\frac{1+\tilde{r} -\cos \phi }{[1+(1+\tilde{r})^2-2(1+\tilde{r})\cos \phi ]^{3/2}}{d}\phi \end{equation} \tag{ A8 }$$

where ϕ_crit is given by

$$\begin{equation} \phi _{\rm crit}={\rm arcos}\left(1-\frac{R_H}{a}\right)+{\rm arcos}\left(\frac{1-R_H/a}{1+r/a}\right). \end{equation} \tag{ A9 }$$

The integral can be evaluated and formally has the following solution:

$$\begin{eqnarray} \left<a^{\rm rad}_R\right> &=& \frac{2\kappa L_p}{8\pi ^2 c a^2}(\tilde{r}^2 \sqrt{(2 + 2\tilde{r} + \tilde{r}^2 - 2(1 + \tilde{r}) \cos \phi _{\rm crit})/\tilde{r}^2} \nonumber\\ && \times \,E(\phi _{\rm crit}/2, -2\sqrt{1 + \tilde{r}})/\tilde{r}) \nonumber \\ && + \tilde{r} (2 + \tilde{r}) \sqrt{(2 + 2\tilde{r} + \tilde{r}^2 - 2(1 + \tilde{r}) \cos \phi _{\rm crit})/\tilde{r}^2} \nonumber\\ && \times F(\phi _{\rm crit}/2, -2\sqrt{1 + \tilde{r}}/\tilde{r}) \nonumber \\ && + 2(1 + \tilde{r}) \sin \phi _{\rm crit})/(\tilde{r} (1 + \tilde{r}) (2 + \tilde{r}) \nonumber\\ && \times \sqrt{2 + 2 \tilde{r} + \tilde{r}^2 - 2 (1 + \tilde{r}) \cos \phi _{\rm crit}}), \end{eqnarray} \tag{ A10 }$$

where E(θ, a) and F(θ, a) are the elliptic integral of the first and second kind. What is not immediately obvious from this expression is that the result is effectively independent of ϕ_crit and, thus, our choice of where to place the gap edge. This arises since most of the impulse comes when the dust particle is close to the planet in its orbit rather than at large angles. In the numerical calculations, we include the attenuation of the radiation due to absorption by the dust, and thus evaluate the orbit averaging integral numerically.

Zoom In Zoom Out Reset image size

We test whether the numerical scheme is behaving as expected by considering several test problems. First, to test that the "leakage" implementation across the gap does not provide artificial dust trapping, we consider a simple but highly relevant test problem. A very small dust particle (s = 10⁻⁵⁰ cm for the test), which will behave as a passive scalar, should exactly follow the gas distribution, as at such a small size it does not feel dust drag. Therefore, in the absence of radiation pressure, the dust-to-gas mass ratio should be radially fixed for the entire simulation range (excluding the region inside the planet's Hill sphere where the source and sink terms dominate) and the dust particles should be perfectly advected across the gap by the "leakage" scheme. The result of this test calculation is shown in Figure 8, where we find our scheme behaves accurately and is suitable for the required calculation. We also perform a resolution study and check that our scheme conserves mass to machine precision as expected.

Zoom In Zoom Out Reset image size

Furthermore, we can test out that we are accurately capturing the effect of radiation pressure on the small dust particles. We can construct a steady-state approximate analytic solution to the dust and gas problem. Assuming the dust particles are well coupled, τ_s ≪ 1, such that we can safely neglect dust drag, and the radial velocity of the dust particles is approximately given by

$$\begin{equation} v_R\approx u_R+\frac{\tau _s\kappa F^p_R}{c\Omega }. \end{equation} \tag{ B1 }$$

If we consider a steady disk, with constant influx of dust and gas at the outer boundary, then we may write

$$\begin{equation} v_R=-\left[\frac{\dot{M_*}}{2\pi R \Sigma _g}-\left(\frac{3\pi }{8c\Sigma _g\Omega }\right)F^p_R\right]. \end{equation} \tag{ B2 }$$

Given a steady disk, we can express the surface density of the gas as (e.g., Pringle 1981)

$$\begin{equation} \Sigma _g=\frac{\dot{M}_*}{3\pi \nu }\left(1-\sqrt{\frac{R_{\rm in}}{R}}\right), \end{equation} \tag{ B3 }$$

where R_in is the inner radius of the disk. Noting for ν∝R, if $F^p_R\propto R^{-2.5}$ then the dust and gas velocity have the same radial dependence (specifically, they are constant with a radius of ν∝R), then for R ≫ R_in, the dust concentration will be constant. Therefore, we can ignore the dust diffusion term and simply find the dust surface density using the following equation:

$$\begin{equation} \Sigma _d=\frac{X\dot{M}_*}{2\pi R \left[\frac{\dot{M_*}}{2\pi R \Sigma _g}-\left(\frac{3\pi }{8c\Sigma _g\Omega }\right)F^p_R\right]}, \end{equation} \tag{ B4 }$$

where X is the dust-to-gas mass ratio of the material injected into the disk at the outer boundary and Σ_g is given by Equation (B3). We simulate this steady problem with an accretion rate in the disk of 10⁻⁸ M_☉ yr⁻¹, an outer boundary of 100 AU, inner boundary of 0.04 AU, viscous α of 0.01, and H/R of 0.04 at 1 AU, where ν∝R throughout the entire grid. We use a dust particle of ρ_d = 2 g cm⁻³, a radius of 0.1 μm, and a dust-to-gas mass ratio of 0.01. We assume the flux scales as $F^p_R\propto R^{-2.5}$ and choose the flux (measured at 10 AU) to have values of 0, 5 × 10⁴, 7.5 × 10⁴, and 1 × 10⁵ erg s⁻² cm⁻². The comparison between the analytic solution (dashed lines) for the dust surface density and velocity and simulations (solid lines) are shown in Figure 9. We find excellent agreement between the code and analytic solutions, indicating the radiation pressure routine in our numerical method is behaving as expected.

Zoom In Zoom Out Reset image size