Interpreting Brightness Asymmetries in Transition Disks: Vortex at Dead Zone or Planet-carved Gap Edges?

We investigated the formation and evolution of a large-scale vortex in two scenarios. In the first scenario, the vortex forms at the edges of a gap opened by a giant planet. Note that we only model the vortex formation at the outer gap edge.⁵ In the second scenario, the vortex is excited by a steep transition in the viscosity, e.g., at the edge of a dead zone. Two-dimensional grid-based, global hydrodynamic disk simulations have been used to study the properties of the vortices under these conditions. We used the GPU-supported version of the code FARGO (Masset 2000), which numerically solves the vertically integrated continuity and Navier–Stokes equations, using locally isothermal approximation in a cylindrical coordinate system centered on the star.

We used the α prescription for the viscosity (Shakura & Sunyaev 1973) with ${10}^{-5}\leqslant \alpha \leqslant {10}^{-3}$. Initially, the gas density has a power-law profile ${\rm{\Sigma }}(R)={{\rm{\Sigma }}}_{0}{R}^{-p}$, with $p=0.5,1$, and 1.5 and ${{\rm{\Sigma }}}_{0}=1.45\times {10}^{-5},5.3\times {10}^{-4}$, and $1.45\times {10}^{-4}$, respectively, for which cases the disk mass is $0.02\,{M}_{* }$ independent of p. As the Toomre Q parameter (Toomre 1964) initially is higher than unity (${Q}_{\min }\simeq 13$ at the vortex radial distance) for the investigated density profiles or even in the center of the fully developed large-scale vortex (${Q}_{\min }\simeq 5$), the disk is gravitationally stable; therefore, we neglect the disk self-gravity.

The pressure scale height of the disk is assumed to have a power-law dependence on the radius, $H(R)={{hR}}^{1+\gamma }$, where h is the aspect ratio and γ is the flaring index. For the flat-disk assumption, we use h = 0.05 and $\gamma =0$, and for the flaring-disk models, we use h = 0.05 and $\gamma =2/7$ (Chiang & Goldreich 1997). We define a fiducial disk setup for which p = 0.5, $\gamma =0$, and h = 0.05 are used.

Assuming locally isothermal approximation, the disk's temperature profile is $T(r)\sim {R}^{2(\gamma -0.5)}$. In this approximation, the equation of state of gas is $P={c}_{{\rm{s}}}(R){\rm{\Sigma }}(R)$, where ${c}_{{\rm{s}}}(R)={\rm{\Omega }}(R)H(R)$ is the sound speed ${\rm{\Omega }}(R)$ and H(R) are the local Keplerian angular velocity and the local pressure scale height, respectively. The unit length is taken to be 1 au, and the unit mass is the stellar mass. Assuming that the unit time is the inverse of the Kepler frequency, the orbital period becomes 2π and the gravitational constant is unity.

The spatial extension of the computational domain is $3\mbox{--}50\,\mathrm{au}$, consisting of N_R = 512 logarithmically distributed radial and ${N}_{\phi }=1024$ equidistant azimuthal grid cells. At the inner and outer boundaries, a wave-damping boundary condition is applied (de Val-Borro et al. 2006). For this boundary condition, the disk mass is not conserved. Indeed, the disk mass increases with time in the simulation, but the total variation of the disk mass during the entire simulation is less than a percent. Note that setting the boundary conditions to open has no significant effect on the results but increases the computational time due to waves developed near the inner boundary. In planet-bearing models, we use $\epsilon H(a)$ as the smoothing of the gravitational potential of the planet with an appropriate value of $\epsilon =0.6$ (Kley et al. 2012).

In our simulations, the star is always at the center of the numerical domain. Therefore, the indirect potential is taken into account (see its importance in Mittal & Chiang 2015; Zhu & Baruteau 2016); i.e., two additional indirect potentials are added to the total potential due to (1) the shift of the barycenter by the disk material (i.e., large density asymmetry caused by the vortex) and (2) the shift of the barycenter due to the massive planet.

Simulations cover 200 × 10³ yr, which corresponds to ∼6350 orbits of the nonmigrating planet. Since the initial parameters are set such that the vortex forms at ∼18 au, the simulation time corresponds to ∼2600 Keplerian orbits at the distance of the vortex center.

To study the effect of planet mass on the formation and evolution of vortices developed at the gap edge, we set the planet-to-star mass ratio to $q=1.19,2.38,4.77,$ and $9.54\times {10}^{-3}$, corresponding to a planet mass of $1.25,2.5,5,$ and $10\,{M}_{\mathrm{Jup}}$, assuming a solar-mass central star. We assume a constant semimajor axis of ${a}_{{\rm{p}}}=10\,\mathrm{au}$ for the planet. We also run simulations in three different viscosity regimes for a $q=4.77\times {10}^{-3}$ giant planet, where $\alpha ={10}^{-3},{10}^{-4}$, and 10⁻⁵ are assumed. In the following, we refer to these models as "GAP models."

A gap is opened by the planet in all simulations because the gap-opening criterion is satisfied as long as $\alpha \lesssim 0.01$ for the assumed disk aspect ratio and planet-to-star mass ratio (Crida et al. 2006). The RWI is excited in all models in about 100 orbits, initially with mode number m = 5 being the fastest growth mode, as predicted by analytic theory for relatively low viscosity ($\alpha \leqslant {10}^{-3}$; Lovelace et al. 1999).

As the vortices are subject to a coalescing process, a large-scale vortex develops in a couple of additional planetary orbits, which is due to the domination of lower mode numbers (large vortices) after the saturation of the initially linear growth of the m = 5 mode (Godon & Livio 1999). However, the vortex decays due to the viscous evolution of the disk (see, e.g., Crida et al. 2006). For $\alpha ={10}^{-3}$ and 10⁻⁴, the lifetime of the vortex, ${\tau }_{{\rm{v}}}$, is shorter than 11 × 10³ and $110\times {10}^{3}\,{\rm{y}}{\rm{r}}$, corresponding to about 350 and 3500 planetary (or 144 and 1440 vortex) orbits, respectively.

We run several models to explore the effect of the disk parameters p, γ, or planetary mass on ${\tau }_{{\rm{v}}}$ in the $\alpha ={10}^{-4}$ model. The value of ${\tau }_{{\rm{v}}}$ decreases for steeper surface density distributions as it is found to be ∼90 × 10³ and $\sim 60\times {10}^{3}\,\mathrm{yr}$ (corresponding to ∼1100 and ∼785 vortex orbits) for p = 1.0 and 1.5, respectively. It is known that the indirect potential (emerged due to the barycenter shift) has the side effect of increasing the vortex strength and therefore its lifetime (Mittal & Chiang 2015; Zhu & Baruteau 2016; Regály & Vorobyov 2017). Since the mass confined inside the vortex decreases with increasing p, if the disk mass is kept constant, the indirect potential proportional to the mass of the asymmetry also decreases.

For the flaring-disk model, where $\gamma =2/7$ and p = 0.5, ${\tau }_{{\rm{v}}}\simeq 6.5\times {10}^{3}\,\mathrm{yr}$ (corresponding to about 200 planetary and ∼85 vortex orbits). This can be explained by the viscosity ${\nu }_{\mathrm{flared}}\simeq 5{\nu }_{\mathrm{flat}}$ at $R\simeq 18\,\mathrm{au}$, which results in a weaker vortex and faster viscous evolution of the gap edge.

Assuming the p = 0.5 flat-disk model, ${\tau }_{{\rm{v}}}$ is found to be correlated with the planetary mass, i.e., ${\tau }_{{\rm{v}}}\simeq 30,55,110$, and $150\times {10}^{3}\,\mathrm{yr}$ (corresponding to $\sim 393,720,1440$, and 1964 vortex orbits) for a $1.25,2.5,5$, and $10\,{M}_{\mathrm{Jup}}$ planet, respectively. This can be explained by the fact that the density gradient, and therefore the pressure gradient, is sharper for a higher planet mass, which results in a stronger vortex and longer vortex decay time.

Contrary to the previous cases, the vortex lifetime is longer than the simulations for $\alpha ={10}^{-5}$ models with p = 0.5 and 1, while the vortex is dissipated by the end of the simulation for p = 1.5. The left panel of Figure 1 shows the density distribution of the disk in the GAP model at the end of the simulation for $\alpha ={10}^{-5}$ assuming p = 0.5. We emphasize that the vortex can survive for a sufficiently long time to be observed in the submillimeter (let us say ≫100 × 10³ yr by the time the transitional disk phase presumably develops) only for a nearly inviscid assumption ($\alpha \leqslant {10}^{-5}$) in conjunction with the results of previous investigations (Godon & Livio 1999; de Val-Borro et al. 2007; Ataiee et al. 2013; Fu et al. 2014a).

Zoom In Zoom Out Reset image size

To model the formation of vortices in a disk with a dead zone, we reduced the value of α smoothly within a certain radius such that $\alpha (R)=\alpha {\delta }_{\alpha }(R)$. The viscosity reduction factor is given by

$$\begin{eqnarray}&&{\delta }_{\alpha }(R)=1-\displaystyle \frac{1}{2}(1-{\alpha }_{\mathrm{mod}})\left[1-\tanh \left(\displaystyle \frac{R-{R}_{\mathrm{dze}}}{{\rm{\Delta }}{R}_{\mathrm{dze}}}\right)\right],\end{eqnarray} \tag{ 1 }$$

where ${\alpha }_{\mathrm{mod}}=0.01$ is the depth of turbulent viscosity reduction. To quantify the dead-zone edge radius, ${R}_{\mathrm{dze}}=24\,\mathrm{au}$ was used, where we adopt the results of Matsumura & Pudritz (2005), who found that R_dze lies between 12 and 36 au, depending on the density of the disk. For the dead-zone edge width, we assumed ${\rm{\Delta }}{R}_{\mathrm{dze}}=1$ and $1.5{H}_{\mathrm{dze}}$, where ${H}_{\mathrm{dze}}={R}_{\mathrm{dze}}h$ (assuming a flat-disk model) is the disk scale height at the outer edge of the dead zone. Although these parameters correspond to 1.8 and 2.7 au assuming h = 0.05, according to Equation (1), the gross widths of the regions where the viscosity transitions occur are 3.6 and 5.4 au. Note that we did not study the effect of the inner edge of the dead zone. In the following, we refer to these models as "DZE models."

Due to the sudden change in viscosity (i.e., in the accretion rate) at the outer dead-zone edge, density and pressure enhancements appear that are found to be unstable to RWI (Varniere & Tagger 2006; Terquem 2008). The fastest-growing mode is m = 5; thus, five anticyclonic vortices initially form that later coalesce (Lovelace et al. 1999; Li et al. 2000). As a result, a large-scale vortex develops inside the dead zone. Note that the large-scale vortex forms at about the same distance, $R\simeq 18\,\mathrm{au}$, as in the GAP models assuming ${R}_{\mathrm{dze}}=24\,\mathrm{au}$ (see Figure 1). The DZE vortices are also subject to slow decay due to the disk viscous evolution, similar to the GAP models. However, the lifetimes of the large-scale vortices are found to be longer than the simulation time in all models; i.e., DZE vortices are long-lived structures for these cases (Regály et al. 2012; Lin 2014).

To compare the morphologies of the full-fledged vortices in the GAP and DZE models, we determine the vortex strengths characterized by the Rossby number,

$$\begin{eqnarray}&&{R}_{{\rm{o}}}(R,\phi )=\displaystyle \frac{{\rm{\nabla }}\times ({\boldsymbol{u}}(R,\phi )-R{\rm{\Omega }}(R))}{2{\rm{\Omega }}({R}_{{\rm{v}}})},\end{eqnarray} \tag{ 2 }$$

i.e., the z component of the vorticity measured in the local frame of the vortex divided by the global vorticity of the Keplerian disk. The vortex strength, R_o, defined by the mean value of ${R}_{{\rm{o}}}(r,\phi )$ inside the vortex (set by the contour level of density at $1/e$ times the maximum density measured at the vortex center), is found to be −0.14 and −0.13 for p = 0.5 and 1 in the GAP models and −0.068, −0.076, and −0.041 for p = 0.5, 1, and 1.5 in the DZE models at the end of the simulations.

For a steady-state vortex with uniform vorticity,

$$\begin{eqnarray}&&{R}_{{\rm{o}}}^{\mathrm{Kida}}=-\displaystyle \frac{3}{4}\displaystyle \frac{{\chi }^{2}+1}{\chi (\chi -1)}+\displaystyle \frac{3}{4},\end{eqnarray} \tag{ 3 }$$

where χ is the vortex aspect ratio defined by the ratio of the azimuthal and radial axes of the vortex (Kida 1981; Chavanis 2000). For a vorticity field whose strength linearly depends on the distance from the vortex center, Goodman et al. (1987) proposed

$$\begin{eqnarray}&&{R}_{{\rm{o}}}^{\mathrm{GNG}}=-\displaystyle \frac{\sqrt{3}}{2}\displaystyle \frac{{\chi }^{2}+1}{\chi \sqrt{{\chi }^{2}-1}}+\displaystyle \frac{3}{4},\end{eqnarray} \tag{ 4 }$$

while Surville & Barge (2015) recently proposed a more elaborate Gaussian model, where

$$\begin{eqnarray}&&{R}_{{\rm{o}}}^{\mathrm{Gauss}}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\chi }^{2}+1}{{\chi }^{2}-1}\left(\displaystyle \frac{3}{2}-\sqrt{3}\right)+\displaystyle \frac{3}{2}.\end{eqnarray} \tag{ 5 }$$

We emphasize that Equations (4) and (5) give ${R}_{{\rm{o}}}=(1/4)(3\,-2\sqrt{3})\lt -0.116$ when $\chi \to \infty$. To compare vortex aspect ratios for the GAP and DZE vortices, we use the Kida approximation, as R_o is found to be $\geqslant -0.1$ in all DZE models.

The evolution of χ calculated according to Equations (2) and (3) is shown in Figure 2 for the DZE and GAP models. In the $\alpha ={10}^{-5}$ GAP models, χ is found to be $6\lesssim \chi \lesssim 8$ and shows only weak evolution in time. However, the vortex lifetime is shorter than the simulation in the $\alpha ={10}^{-4}$ GAP models, and $\chi \gtrsim 10$. Note that χ shows no clear trends with the planet mass, as we found that the strongest vortex, $\chi \simeq 6$, formed in models where a $q=2.5\,{M}_{\mathrm{Jup}}$ planet is assumed (dark blue curve in Figure 2). In contrast, χ evolves strongly in time for the DZE models: initially, the vortex strengthens, i.e., χ decreases up to ∼500 orbits; then it starts to weaken, i.e., χ increases, ${\chi }_{\max }\simeq 14$ and 10 for ${\alpha }_{\mathrm{dz}}={10}^{-4}$ and 10⁻⁵, respectively. Note that χ does not saturate by the end of the simulation for the DZE models. However, as we showed earlier, a quasi-stationary vortex mode can be reached after ∼5000 orbits of the vortex (see Figure 5 in Regály et al. 2012). In summary, $\chi \simeq 8$ in the GAP models where the vortex is maintained until the end of the simulations assuming $\alpha \leqslant {10}^{-5}$. For the DZE models, the vortex weakens by the end of the simulation and thus has a larger aspect ratio: $\chi \simeq 14$ and $\chi \simeq 10$ for ${\alpha }_{\mathrm{dz}}={10}^{-4}$ and ${\alpha }_{\mathrm{dz}}={10}^{-5}$, respectively.

Zoom In Zoom Out Reset image size

In order to compare the spatial extensions of the vortices in the GAP and DZE models, the radial (left) and azimuthal (right) density profiles taken across the vortex center are shown in Figure 3. The radial and azimuthal extensions of the vortices, ${{\rm{\Delta }}}_{\mathrm{rad}}$ and ${{\rm{\Delta }}}_{\mathrm{azim}}$, are determined by measuring the spatial extension of the vortex at the $1/e$th of the density maximum. For the GAP model, ${{\rm{\Delta }}}_{\mathrm{rad}}\simeq 2.7\,\mathrm{au}$ and ${{\rm{\Delta }}}_{\mathrm{azim}}\simeq 25\,\mathrm{au}$ (corresponding to $\sim 244^\circ$ angular width), while for the DZE model, ${{\rm{\Delta }}}_{\mathrm{rad}}\simeq 4.2\,\mathrm{au}$ and ${{\rm{\Delta }}}_{\mathrm{azim}}\simeq 80\,\mathrm{au}$ (corresponding to $\sim 95^\circ$ angular width).

Zoom In Zoom Out Reset image size

Although the vortex excitation model (GAP or DZE) could be distinguished by measuring the vortex aspect ratios, the precise determination of χ would require a very high resolution (i.e., $\sim 0\buildrel{\prime\prime}\over{.} 01$ beam size for our model assuming a 100 pc source distance) to resolve the vortices radially. Moreover, in some cases, the vortex strengths are similar in the DZE and GAP models: $\chi \simeq 10$ for the $\alpha ={10}^{-5}$ GAP and ${\alpha }_{\mathrm{dz}}={10}^{-5}$ DZE models, respectively. Therefore, vortex aspect ratio is not an ideal parameter to use for distinguishing the two formation scenarios. However, vortices are about an order of magnitude wider in the azimuthal than in the radial direction, easing the determination of the vortex azimuthal width. The azimuthal extensions of the DZE vortices are found to exceed $\sim 180^\circ$, while they are only $\sim 90^\circ$ for the GAP vortices, in agreement with the literature (Lyra et al. 2014; Miranda et al. 2016, 2017; Zhu & Baruteau 2016). Note that Lyra & Mac Low (2012) and Flock et al. (2015) presented azimuthally concentrated DZE vortices presumably caused by MHD effects, but see Zhu & Stone (2014). Note also that the vortex azimuthal extension in the GAP models that includes disk self-gravity is found to be in the range of 90°–180° (Lin & Papaloizou 2011). However, the stretching requires a more massive disk than what we assumed. Therefore, we suggest that the vortex azimuthal extension is a key parameter that can be used to infer the vortex excitation mechanism as long as disk self-gravity is negligible.

Here we have to note that vertical shear instability (Nelson et al. 2013) is effective in transporting angular momentum in the disk causing $\alpha \simeq {10}^{-4}$ (Stoll & Kley 2014). Thus, the magnitude of viscosity, $\alpha ={10}^{-5}$, required for the long-term vortex in the GAP models is unreasonably small. However, ${\alpha }_{\mathrm{dz}}={10}^{-4}$ is a reasonable value for the disk viscosity in the dead zone. For this reason, we compare the observational properties of vortices formed in the $\alpha ={10}^{-5}$ GAP and ${\alpha }_{\mathrm{dz}}={10}^{-4}$ models.

In the GAP model, it is assumed that the dust is completely removed because the planet clears all the material inside its gap. However, a giant planet opens only a narrow gap whose width is ${\rm{\Delta }}{R}_{\mathrm{GAP}}\simeq 5-6{R}_{H}$, where ${R}_{H}=a{(q/3)}^{1/3}$ is the radius of the planetary Hill sphere (Crida et al. 2006). Therefore, a system of multiple planets is required to completely clear the inner disk (Dodson-Robinson & Salyk 2011; Zhu et al. 2011). In the hydrodynamic simulations of Isella et al. (2013), three equal-mass giant planets ($q=1.9\times {10}^{-3}$) orbiting at $a=21,34$ and 55 au were assumed to explain the observed brightness asymmetry of LkHα 330. They found complete gas clearing in their simulations due to overlapping of planetary gaps. We repeated their calculations with $\alpha =2\times {10}^{-3}$ and 2 × 10⁻⁴ (we use the same α assumptions for the viscosity as Isella et al. 2013) and confirmed the complete gas clearing. Note, however, that the vortex lifetime exceeds 10⁵ yr only for the nearly inviscid case, where α is on the order of 10⁻⁵.

We emphasize that, according to our simulations, the observed vortex structures in the gas, i.e., the azimuthal extension and contrast, for multiple- and single-planet setups are similar. This can be explained by the fact that RWI excitation at the outer edge of a common gap is driven by the outermost giant planet only, while the inner planets do not affect the vortex formation there. Although our hydrodynamic models with a single giant planet (presented in Section 2.2) cannot explain the inner disk clearing, assuming complete dust clearing inside the planetary orbit fairly models the formation and evolution of the large-scale vortex for multiple-planet systems.

In the DZE model, however, the inner disk is not cleared in the gas; therefore, to test the dust removal assumption, we investigated the dust drift in a viscously evolving 1D disk model. The radial drift of the solid particles in the viscously evolving dead zone is modeled in a representative particle approach, in which dust coagulation and fragmentation are not included. We followed the dust drift for 3 × 10⁵ yr with 10⁵ dust particles. In this simple 1D model, we assume that the disk is axisymmetric, thus the evolution of the surface density of the gas (${{\rm{\Sigma }}}_{g}(R,t)$) can be given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{g}(R,t)}{\partial t}=\displaystyle \frac{3}{R}\displaystyle \frac{\partial }{\partial R}\left[{R}^{1/2}\displaystyle \frac{\partial }{\partial R}(\nu (R){{\rm{\Sigma }}}_{g}(R,t){R}^{1/2})\right],\end{eqnarray} \tag{ 6 }$$

where $\nu (R)=\alpha {\delta }_{\alpha }(R){c}_{s}^{2}(R)$ is the kinematic viscosity of the gas in our dead-zone edge model, where ${\delta }_{\alpha }(R)$ is defined by Equation (1).

The coupling of the dust to the gas can be described by the Stokes number (St), which in the Epstein regime (where the interaction between particles and a single gas molecule becomes important, as the particle size is comparable to the mean free path of molecules) is

$$\begin{eqnarray}&&{St}=\displaystyle \frac{a{\rho }_{s}}{{{\rm{\Sigma }}}_{g}(R,t)}\displaystyle \frac{\pi }{2},\end{eqnarray} \tag{ 7 }$$

where a is the size (assumed to be mm and μm in two simulations, respectively) and the intrinsic density of the solid particle is set to ${\rho }_{s}=1.6\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. The radial velocity of the solid particles (${u}_{\mathrm{dust},{\rm{r}}}$) is given by

$$\begin{eqnarray}&&{u}_{\mathrm{dust},{\rm{r}}}=\displaystyle \frac{{u}_{\mathrm{gas},{\rm{r}}}}{1+{{St}}^{2}}+\displaystyle \frac{2}{{St}+{{St}}^{-1}}{u}_{\mathrm{drift}},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}&&{u}_{\mathrm{gas},{\rm{r}}}=-\displaystyle \frac{3}{{{\rm{\Sigma }}}_{g}(R,t){R}^{1/2}}\displaystyle \frac{\partial }{\partial R}({{\rm{\Sigma }}}_{g}(R,t)\nu (R){R}^{1/2})\end{eqnarray} \tag{ 9 }$$

is the radial velocity of the gas and

$$\begin{eqnarray}&&{u}_{\mathrm{drift}}=\displaystyle \frac{{c}_{s}^{2}(R)}{2{\rm{\Omega }}(R)R}\displaystyle \frac{d\mathrm{ln}P}{d\mathrm{ln}R}\end{eqnarray} \tag{ 10 }$$

is the maximum drift velocity of a particle, where ${\rm{\Omega }}(R)={R}^{-3/2}$ (see, e.g., Weidenschilling 1977; Nakagawa et al. 1986; Youdin & Lithwick 2007; Birnstiel et al. 2012b).

In our dust drift model, the disk extends from 0.1 to 50 au with 1500 radial cells. Similar to our fiducial 2D DZE hydrodynamic model, the outer dead-zone edge is set to ${R}_{\mathrm{dze}}=24\,\mathrm{au}$, ${\rm{\Delta }}{R}_{\mathrm{dze}}=1.5{H}_{\mathrm{dze}}$, and the initial surface mass density of the gas is ${{\rm{\Sigma }}}_{\mathrm{gas}}(R)={{\rm{\Sigma }}}_{0}{R}^{-0.5}$, where ${{\rm{\Sigma }}}_{{\rm{g}}}=1.45316\times {10}^{-5}\,{M}_{\odot }\,{\mathrm{au}}^{-2}$ is used.

Figure 4 shows the evolution of the semimajor axis of mm-sized (left panel) and μm-sized (right panel) dust particles in the fiducial DZE model. A negative pressure gradient results in sub-Keplerian gas velocity and causes rapid inward drift of dust particles, while the gas is super-Keplerian at positive pressure gradients, which causes outward drift of dust particles (see Figure 5). As a result, most of the mm-sized dust particles initially orbiting inside the dead zone drift toward the central star, while the dust initially orbiting outside the dead-zone edge is trapped at the pressure maximum. As a result, mm-sized dust particles are accumulated at the local pressure maximum. However, μm-sized dust particles are well coupled to the gas; therefore, the dead zone remains well populated with μm-sized dust.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The radial extension of the feeding zones for the pressure maxima evolves in time, which can affect their dust accumulation efficiency. Due to the viscous evolution of the disk, the pressure maximum shifts inward (see, e.g., Regály et al. 2012). We found that the mm-sized dust particles are completely absent in the dead zone by $\sim 250\times {10}^{3}\,\mathrm{yr}$, while the dead zone is still populated by μm-sized dust. Note that if we take into account the dust coagulation, the average dust size increases in time. Since the dust drift velocity increases with particle size, the dust clearing would be accelerated.

Note that Birnstiel et al. (2012a) found that the presence of the dead zone fails to explain the dust clearing of the inner disk observed in transition disks using a more elaborate model that also takes into account dust coagulation and fragmentation. However, they did not take into account the effect of the pressure maximum formed at the dead-zone outer edge. Since particles strongly coupled to the gas (less than several meters in size) cannot cross the pressure bump (Figure 5), the dead zone cannot be replenished with mm-sized dust from the outer disk.

In short, it is a plausible assumption that the inner disk is completely cleared of mm-sized dust within several 10⁵ yr for the DZE models assuming our fiducial disk parameters, while multiple planets are required to completely clear the inner disk of mm-sized dust for the GAP models.

Synthetic images are calculated by the 3D radiative transfer code RADMC-3D.⁷ Our imaging simulations contain three steps. First, the dust temperature is calculated in a thermal Monte Carlo simulation, then images at 880 μm are calculated using the ray-tracing module of RADMC-3D, and finally, ALMA observations are simulated using the Common Astronomy Software Applications (CASA) package⁸ v4.2.2.

Since FARGO solves dimensionless equations, the simulations can be scaled in terms of spatial extent and also in mass in the locally isothermal approximation. Therefore, to make our synthetic images comparable to observations, we can scale the hydrodynamic simulations such that the disk extends from 9 to 150 au with the planetary orbital radius at 30 au. In this case, the vortex is located at about 30–45 au from the star, similar to the observations (see, e.g., Brown et al. 2009), and the simulation time corresponds to ∼1 Myr. For both the DZE and GAP models, we run two set of simulations, where the core of the vortex is optically thin and thick (by appropriately scaling the density).

In the RADMC-3D simulations, we use a 3D spherical mesh ($R,\phi ,\theta$). The radial grid extends between 9 and 150 au and contains N_R = 256 grid cells. Along the azimuthal and poloidal angular coordinates, we use ${N}_{\phi }=512$ and ${N}_{\theta }=180$ grid cells, respectively. While the azimuthal grid is equidistant, the grid cells in the poloidal coordinate are distributed such that we use 10, 160, and 10 grid points in the [0, $\pi /2-0.25$], $[\pi /2-0.25$, $\pi /2+0.25$], and $[\pi /2+0.25$, π] intervals, respectively. In each poloidal interval, we use a uniform, equidistant grid.

Hydrodynamic frames are interpolated to the RADMC-3D mesh and converted to volume density assuming a Gaussian vertical density distribution assuming a constant aspect ratio of h = 0.05, similar to the 2D hydrodynamic simulations. In the radiative transfer model, we use the optical constants of the astronomical silicate from Weingartner & Draine (2001) for the dust particles and calculate the absorption and scattering cross-sections using Mie theory. Dust grains in the model have a size distribution of $n(a)\propto {a}^{-3.5}$ between 0.1 and 1000 μm, corresponding to the steady-state mass distribution of the collisional fragmentation cascade (Dohnanyi 1969). The central radiation source in the model has parameters representative of that of a Herbig Ae star: ${R}_{\star }=2.5\,{R}_{\odot },{T}_{\star }=9500\,{\rm{K}},\ {M}_{\star }=2.0\,{M}_{\odot }$.

To determine the temperature of individual dust grains, we run a thermal Monte Carlo simulation using $1.2\times {10}^{8}$ photon packets. The continuum images at 880 μm are calculated with ray tracing consisting of 800 by 800 pixels, with a pixel size of 0.89 mas, assuming a source distance of 140 pc.

Simulated ALMA observations are calculated with the CASA simobserve task to generate synthetic visibilities and the clean task for imaging. The full 12 m array was used for the simulated observations in three different configurations, resulting in a spatial resolution of approximately 025, 018, and 01. We assume the full 7.5 GHz continuum bandwidth of the ALMA correlator and calculate the visibilities for a total integration time of 30 minutes. The simulations include thermal noise arising from a water vapor column of 0.913 mm representative of Band 7 observations.

In both the optically thin and thick cases, we calculated images at two inclination angles, $30^\circ$ and $60^\circ$, and five different azimuth angles for the vortex between 0° and 180° in 45° steps. To study the effect of the shape of the synthesized beam on the morphology of the images, we simulated observations assuming a source decl. of 25° and 70°, and the times of the observations are centered on 0, −1, and −2 hr angles.

The most striking difference between the vortices excited at the edge of a planet-carved gap and at the edge of a dead zone is their azimuthal extent (see Section 2.4). Our hydrodynamic simulations show that vortices formed at the edge of a planetary gap are significantly more compact in the azimuthal direction than their counterparts at the edge of a dead zone. In the following, we show that this property can be used in observations to infer the formation mechanism of the vortices.

To test the applicability of this technique, we determine the azimuthal extent of the vortices in the synthetic ALMA observations in the following way. As a first step, we deproject the images with a known inclination and position angle. In real observations of well-studied objects, both the inclination and the position angle are known to a reasonable (several degrees) accuracy from gas-line observations. Then, we converted the image from Cartesian to polar coordinates centered on the stellar position. Finally, we measured the azimuthal size of the region where the intensity is 66% of that of the peak intensity in the core of the vortex. The measured sizes of the vortices are presented in Figure 6 for different position angles of the vortex, hour angles the observations are centered on, and optically thick and thin vortex core assumptions.

Zoom In Zoom Out Reset image size

For low inclination angles ($30^\circ$; panels (a), (c), and (e) of Figure 6), the DZE and GAP models can be separated in most cases based on the measured azimuthal size of the vortex. While for the lowest resolution (∼025; see panel (a)), some confusion may occur between the two models, at the highest resolution (01; see panel (e)), the DZE and GAP models are clearly separated. For high inclination angles ($60^\circ$; panels (b), (d), and (f) of Figure 6), separating the two models becomes more challenging, especially for the lowest spatial resolution.

Based on our results presented in Figure 6, we conclude that if (1) the disk is optically thin, (2) the inclination is low (probably 45° or less), and (3) the resolution is sufficiently high (the beam size is at least twice as small as the projected distance from the star to the vortex core), the DZE and GAP models can be distinguished. As we show in the following, the morphology of the vortex in the two different models can be very similar for higher optical depths and inclination angles and lower spatial resolutions.

3.3.1. Effect of Optical Depth

The true azimuthal extent of the vortex can only be measured in submillimeter continuum images if the disk is optically thin, in which case the DZE and GAP models can easily be distinguished (Figure 7). Once the disk becomes optically thick, the surface brightness is determined by the temperature instead of the surface density. Since the highest surface density is in the core of the vortex, the disk optical thickness becomes highest there. The size of the region where the observed intensity drops below a certain fraction of the peak intensity (i.e., the azimuthal extension of the horseshoe-shaped asymmetry) also increases with the optical depth. Therefore, the azimuthal extension of the vortex is overestimated if the disk is optically thick at the vortex core. For the optically thick case, the azimuthal contrast across the vortex is also decreased by the saturation of the emission at the vortex core. The combination of these effects can make the appearance of a vortex at a planetary gap edge very similar to that of a vortex at the edge of a dead zone, but only for low resolution (Figure 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3.2. Effect of Beam Shape

In interferometric observations, the synthesized beam is not generally circular but rather elliptic due to the combination of the elevation of the target above the horizon and the antenna configuration. The ellipticity of the beam causes the structures in the image to look more elongated along the direction of the beam major axis. Such elliptical beams can cause characteristic changes to the apparent morphology of the observed vortices.

The azimuthal extent of the vortex in the GAP models depends on the position angle of the beam ellipse. The apparent azimuthal extent of the vortex is highest if the beam major axis is perpendicular to the vortex. For the optically thick case, the vortex in the GAP model can look azimuthally as much extended as a DZE vortex (panels (a) and (b) of Figure 9).

Zoom In Zoom Out Reset image size

At low inclination angles, the peak of the emission is located in the middle of the extended horseshoe-shaped structure, i.e., at the vortex center. At high inclination angles, however, the peak of the brightness can shift from the vortex center, depending on the relative position of the beam major axis and the vortex (panels (b) and (d) in Figure 8). We note that the horseshoe-shaped asymmetry is always azimuthally less extended in the GAP models than in the DZE models in the optically thin case (Figure 8).

In contrast to the GAP models, the vortex in the DZE models can split into two bright blobs if the beam is elongated. For optically thick DZE vortices, this is always the case. If, however, a DZE vortex is optically thin, the number of bright blobs and/or the azimuthal extent of the observed horseshoe depends strongly on the position angle of the beam ellipse with respect to that of the vortex. If the beam major axis is perpendicular to the vortex, an azimuthally extended horseshoe appears with two distinct bright peaks (panel (a) of Figure 10). In contrast, if the beam major axis is parallel to the vortex, a single-peaked asymmetry appears with a significantly smaller azimuthal extension (panel (b) of Figure 10).

Zoom In Zoom Out Reset image size

We note that this direction-dependent "smearing" effect also depends on the ratio of the structure size to the beam size. Beam-shape effects cannot change the morphology of structures in the image that are larger than a few beam sizes. Therefore, increasing the angular resolution in the observations can in most cases break the degeneracy, making it possible to distinguish DZE and GAP models based on the vortex azimuthal extent.

While there is now observational evidence for nonaxisymmetric surface brightness distribution in the submillimeter for more than a dozen transitional disks, the precise morphology of the asymmetry is rather uncertain in many cases due to insufficient spatial resolution and/or the low signal-to-noise ratio of the observation. However, in a few cases, if ALMA observations are available and the structure of the disk is known well enough, we can qualitatively compare the observed morphology of the disk—most importantly, the azimuthal extent of the observed vortices—to our models.

There are four well-known transitional disks with ALMA observations and evidence for a vortex in the disk: Oph IRS 48 (van der Marel et al. 2013), HD 142527 (Casassus et al. 2013), HD 135344B (Pérez et al. 2014), and SR 21 (Brown et al. 2009). Van der Marel et al. (2013) explained the extreme strong asymmetry in the disk of IRS 48 with a vortex excited by a massive ($\gt 10\,{M}_{\mathrm{Jup}}$) giant planet. Van der Marel et al. (2013) presented a numerical model of a gap-edge vortex showing remarkable similarities with the observations.

The azimuthal extent of the asymmetries in the other three sources, HD 142527, HD 135344B, and SR 21, seems to be significantly larger than that in the disk of Oph IRS 48 and thus may not necessarily be consistent with a planet-induced vortex. To fit the observed images, Pérez et al. (2014) assumed a superposition of a planet-induced vortex and a full ring for both HD 135344B and SR 21. While only a symmetric dust ring is expected to form at the outer edge of a planet-carved gap assuming canonical disk viscosity (${10}^{-3}\leqslant \alpha \leqslant {10}^{-2}$). Therefore, such azimuthally extended asymmetries can be better explained by a sudden jump in the viscosity (proposed by Varniere & Tagger 2006 and later compared to observations by Regály et al. 2012), rather than by putative giant planet(s).

Our models of a DZE vortex bear especially remarkable similarities with the images of SR 21. In Brown et al. (2008), the image of SR 21, obtained with the Sub-Millimeter Array (SMA), shows a horseshoe-shaped structure with more than $270^\circ$ in azimuthal extent and two bright blobs on opposite sides of the arc. In contrast, the ALMA image in Pérez et al. (2014) shows only a single maximum along the horseshoe-shaped asymmetry at the position between the two blobs seen in the SMA image. In Figure 10, we show that an azimuthally extended DZE vortex can reproduce the morphology of both the SMA and the ALMA images remarkably well if we take into account the differences in the shapes of the synthesized beams in the two observations. Note that since our goal was to show a qualitative comparison of our models with the observations, we did not use the exact same uv-coverage as in the real SMA and ALMA observations of SR 21. Instead, we changed the decl. of our target in the synthetic observations such that the resulting synthesized beam will have a similar shape to that in the real observations.

So far, horseshoe-shaped asymmetries in protoplanetary disks have been interpreted with vortices induced by planets. Here we suggest that in a few cases where the vortex is azimuthally extended, this may not be the case, and the vortex could be induced by a viscosity jump (Regály et al. 2012). Unfortunately, the images obtained so far do not have high enough spatial resolution that an azimuthally extended vortex could be distinguished from other models, e.g., a planet-induced vortex + ring proposed by Pérez et al. (2014). A significantly higher spatial resolution (by a factor of about 2–3 at least) compared to already existing ones is needed to unambiguously identify the presence of such extended vortices.

Here we have to mention some caveats of our model. Regarding the hydrodynamical simulations, the most obvious caveat is that we run simulations in 2D. However, Meheut et al. (2012c) found that the vortex formation and development in 3D are indeed very similar to that of 2D: the disk vertical stratification only slightly decreases the growth rate of vortices compared to that of 2D.

We neglect the effect of disk self-gravity, while Lin & Papaloizou (2011) found that in self-gravitating disks, the gap-edge vortex coagulation ends at mode m = 2 for relatively high disk masses ($\geqslant 0.031{M}_{* }$). As a result, two vortices may present simultaneously in the disk. Thus, in the GAP models, we should see a double-peaked submillimeter brightness distribution similar to that of the DZE models. Lovelace & Hohlfeld (2013) found that the disk self-gravity could be important regarding the vortex formation if $Q\lt R/h$. Recently, Zhu & Baruteau (2016) found that taking into account the disk's self-gravity results in a weakened vortex, especially for massive disks. Note, however, that they artificially imposed a density jump on the disk to form a large-scale vortex. Therefore, DZE or GAP vortex formation in self-gravitating disks assuming modest disk mass is worth investigating.

We neglect the planet growth in the GAP models. Recently, Hammer et al. (2017) found that the vortex azimuthal extension generated by a slowly growing giant planet (T_growth being on the order of several 1000 planetary orbits) can be larger by about a factor of two, similar to that of the DZE models. Note, however, that they prescribed planet growth by an artificial sinusoidal function (see their Equation (1)) neglecting the viscous dynamics inside the planet Hill sphere. Thus, it is worth modeling RWI excitation with a more elaborate accretion prescription, like what is described in Kley (1999).

We assume a locally isothermal disk, which is valid only if the cooling time of the gas is fast compared to the local orbital timescale. This assumption, however, can be plausible for a several-million-years-old transition disk. Nevertheless, it is worth investigating vortex formation and survival in adiabatic disks, as Lyra et al. (2009b) found stronger RWI excitation in adiabatic disks, which may have consequences for the vortex evolution.

We completely neglect the MHD effects in our simulations. However, Zhu & Stone (2014) showed that 2D simulations assuming a kinematic viscosity (approximated to the stresses of their 3D MHD models) show vortex formation that is very similar to that of their full MHD simulations. Note also that a long vortex lifetime in the 3D MHD simulations of Zhu & Stone (2014) was only observed if a nonideal effect, such as the ambipolar diffusion, was taken into account. The nonideal effect of ambipolar diffusion results in suppression of MRI turbulence, i.e., low disk viscosity, which helps to maintain the vortex.

In our dead-zone models, we assume quite sharp viscosity transitions to excite RWI. The simulations of Lyra et al. (2009b) and Regály et al. (2012) agree in that RWI excitation occurs only for sharp dead-zone edges in which the width of the viscosity transition is smaller than twice the local pressure scale height. However, Dzyurkevich et al. (2013) have shown that resistivity changes very smoothly (therefore, MRI generated viscosity parameter α too) in the outer dead-zone edge. As a result, the dead-zone edge in an α disk is about five times thicker than the required threshold value for RWI excitation. However, replacing the α prescription with resistive MHD, a smooth change in gas resistivity, is found not to imply a smooth transition in turbulent stress (Lyra et al. 2014). As a result, a large-scale vortex is stable even in smooth resistivity transition, and the abovementioned requirements for exciting RWI are only an inherent feature of α models.

Regarding the calculation of synthetic images, there are also some caveats in our model. The simple assumption that the dust density distribution follows that of the gas (i.e., the dust-to-gas ratio is constant throughout the disk) neglects the earlier results of Johansen et al. (2006), who observed strong enhancement of the dust-to-gas ratio in the presence of anticyclonic vortices. Birnstiel et al. (2013) also found that the dust-to-gas mass ratio could become higher than the global value if the viscosity parameter was larger than the particle's Stokes number. Moreover, we neglect the dust feedback to the vortex evolution too. Recently, Fu et al. (2014b) investigated the dust feedback on the vortex excited by a giant planet having a $q=5\times {10}^{-3}$ planet-to-star mass ratio and found that the dust feedback can destroy the vortex if a large amount of dust ($10\,{M}_{\oplus }$) is accumulated inside the vortex. Note that Regály et al. (2013) investigated the migration trapping of low-mass planets ($10\,{M}_{\oplus }$) in the presence of the dead-zone edge vortex and found similar vortex destruction as the planet enters the vortex. However, since the trapping of the low-mass planet is found to be only temporary, the vortex redevelops due to an existing sharp density jump formed at the dead-zone edge. Therefore, we think that the gap-edge vortex is vulnerable to and the dead-zone–edge vortex is stable against dust feedback. In any case, to give a more reliable model for the evolution of gap-edge and dead-zone–edge vortices, both the dust and gas dynamics should be incorporated in hydrodynamic models.