Effect of Dust Radial Drift on Viscous Evolution of Gaseous Disk

Here we consider structures of gas and dust grains in steady state. First, we ignore the vertical structures and simply consider a 2D disk. We use the polar coordinate ($R,\phi$). We treat only the surface densities of the gas and the dust, ${{\rm{\Sigma }}}_{{\rm{g}}}$ and ${{\rm{\Sigma }}}_{{\rm{d}}}$, respectively, instead of densities of the gas and dust, ${\rho }_{{\rm{g}}}$ and ${\rho }_{{\rm{d}}}$. Kretke et al. (2009) derived the formulae of velocities of the gas and the dust grains in the 2D disk. Very recently Dipierro & Laibe (2017) also derived similar formulae. From Equations (14) and (15) of Kretke et al. (2009), the radial and azimuthal velocities in the case of single-sized dust grains are given by

$$\begin{eqnarray}{v}_{{\rm{R}}}^{2{\rm{D}}} & = & -\displaystyle \frac{2{{St}}_{2{\rm{D}}}^{\prime }}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{\eta }_{2{\rm{D}}}{V}_{K}\\ & & +\displaystyle \frac{1}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{V}_{\mathrm{vis}}^{2{\rm{D}}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}{v}_{\phi }^{2{\rm{D}}} & = & {V}_{K}+\displaystyle \frac{1}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{\eta }_{2{\rm{D}}}{V}_{K}\\ & & -\displaystyle \frac{{{St}}_{2{\rm{D}}}^{\prime }}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{V}_{\mathrm{vis}}^{2{\rm{D}}},\end{eqnarray} \tag{ 2 }$$

and those for the disk gas are given by

$$\begin{eqnarray}{V}_{{\rm{R}}}^{2{\rm{D}}} & = & \displaystyle \frac{2{{St}}_{2{\rm{D}}}^{\prime }}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{\eta }_{2{\rm{D}}}{V}_{K}\\ & & +\left(1-\displaystyle \frac{1}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}\right){V}_{\mathrm{vis}}^{2{\rm{D}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}{V}_{\phi }^{2{\rm{D}}} & = & {V}_{K}+\left(\displaystyle \frac{1}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}-1\right){\eta }_{2{\rm{D}}}{V}_{K}\\ & & +\displaystyle \frac{{{St}}_{2{\rm{D}}}^{\prime }}{{{{St}}_{2{\rm{D}}}^{\prime }}^{2}+1}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}}{V}_{\mathrm{vis}}^{2{\rm{D}}},\end{eqnarray} \tag{ 4 }$$

where ${{St}}_{2{\rm{D}}}^{\prime }={{\rm{\Sigma }}}_{{\rm{g}}}/({{\rm{\Sigma }}}_{{\rm{g}}}+{{\rm{\Sigma }}}_{{\rm{d}}}){{St}}_{2{\rm{D}}}$. The Stokes number of the dust grains in the 2D disk, ${{St}}_{2{\rm{D}}}$, is given by

$$\begin{eqnarray}&&{{St}}_{2{\rm{D}}}={t}_{\mathrm{stop}}^{2{\rm{D}}}{{\rm{\Omega }}}_{{\rm{K}}},\end{eqnarray} \tag{ 5 }$$

where ${t}_{\mathrm{stop}}^{2{\rm{D}}}$ is the stopping time of dust grains in the 2D disk and ${{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{{{GM}}_{* }/{R}^{3}}$ is the Keplerian angular velocity at the midplane, where G and ${M}_{* }$ are the gravity constant and the mass of the central star, respectively. The Keplerian rotation velocity at the midplane is ${V}_{K}=R{{\rm{\Omega }}}_{{\rm{K}}}$. In the Epstein regime, the stopping time is written by (Takeuchi & Artymowicz 2001)

$$\begin{eqnarray}&&{t}_{\mathrm{stop}}(R,z)=\displaystyle \frac{{\rho }_{p}{s}_{{\rm{d}}}}{\sqrt{8/\pi }{\rho }_{{\rm{g}}}{c}_{s}},\end{eqnarray} \tag{ 6 }$$

where ${\rho }_{p}$, ${s}_{{\rm{d}}}$, and ${c}_{s}$ are the size and the internal density of dust grains and the sound speed, respectively. In the 2D disk, using the surface density, we can write the stopping time as

$$\begin{eqnarray}&&{t}_{\mathrm{stop}}^{2{\rm{D}}}(R)=\displaystyle \frac{\pi {\rho }_{p}{s}_{{\rm{d}}}}{2{{\rm{\Sigma }}}_{{\rm{g}}}{{\rm{\Omega }}}_{{\rm{K}}}}.\end{eqnarray} \tag{ 7 }$$

The pressure gradient force is parameterized by

$$\begin{eqnarray}&&{\eta }_{2{\rm{D}}}=-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{h}_{{\rm{g}}}}{R}\right)}^{2}\left(\displaystyle \frac{d\mathrm{ln}{{\rm{\Sigma }}}_{{\rm{g}}}}{d\mathrm{ln}R}+\displaystyle \frac{d\mathrm{ln}{c}_{s}^{2}}{d\mathrm{ln}R}\right),\end{eqnarray} \tag{ 8 }$$

and ${V}_{\mathrm{vis}}^{2{\rm{D}}}$ is the viscous velocity of gas without dust grains, which is given by (e.g., Lynden-Bell & Pringle 1974)

$$\begin{eqnarray}&&{V}_{\mathrm{vis}}^{2{\rm{D}}}=-\displaystyle \frac{3\nu }{R}\displaystyle \frac{d\mathrm{ln}(\nu {{\rm{\Sigma }}}_{{\rm{g}}}{R}^{1/2})}{d\mathrm{ln}R},\end{eqnarray} \tag{ 9 }$$

where ν is the kinematic viscosity of radial diffusion, given by $\nu =\alpha {c}_{s}{h}_{{\rm{g}}}$, where the α-prescription (Shakura & Sunyaev 1973) is used.

In the inviscid case ($\nu =0$ and thus ${V}_{\mathrm{vis}}^{2{\rm{D}}}=0$), only the first terms on the right-hand side (rhs) of Equations (1)–(4) remain, which are the same as Equations (2.11)–(2.14) of Nakagawa et al. (1986). These terms originated from gas–dust friction. On the other hand, if the dust surface density is very small (${{\rm{\Sigma }}}_{{\rm{d}}}\to 0$), the gas velocities in radial and azimuthal directions correspond to ${V}_{{\rm{R}}}^{2{\rm{D}}}={V}_{\mathrm{vis}}^{2{\rm{D}}}$ and ${V}_{\phi }^{2{\rm{D}}}=(1-{\eta }_{2{\rm{D}}}){V}_{K}$, respectively.

We have employed 2D approximation, meaning that the vertical structures of the gas and dust are similar to each other. The gas and dust disks are assumed to have the same scale height. However, the scale height of the dust disk can be much smaller than that of the gas disk, if the dust grains are settled in the midplane. In the next subsection, we discuss the gas and dust velocities, taking into account the vertical structure.

We now consider the gas and dust motions in a 3D disk. We use the cylindrical coordinate (R, ϕ, z). The gas and dust velocities are ${\boldsymbol{V}}=({V}_{{\rm{R}}},{V}_{\phi },{V}_{z})$ and ${\boldsymbol{v}}=({v}_{{\rm{R}}},{v}_{\phi },{v}_{z})$, respectively. The equations of motion for the dust grains are given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{{\rm{R}}}}{\partial t}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){v}_{{\rm{R}}}-\displaystyle \frac{{v}_{\phi }^{2}}{R}=-\displaystyle \frac{{{GM}}_{* }}{{R}^{2}}-\displaystyle \frac{{v}_{{\rm{R}}}-{V}_{{\rm{R}}}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{\phi }}{\partial t}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){v}_{\phi }+\displaystyle \frac{{v}_{{\rm{R}}}{v}_{\phi }}{R}=-\displaystyle \frac{{v}_{\phi }-{V}_{\phi }}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {v}_{z}}{\partial t}+({\boldsymbol{v}}\cdot {\rm{\nabla }}){v}_{z}=-\displaystyle \frac{{{GM}}_{* }}{{R}^{3}}z-\displaystyle \frac{{v}_{z}-{V}_{z}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 12 }$$

where we assume axisymmetric gravitational potential ${\rm{\Psi }}=-{{GM}}_{* }/\sqrt{{R}^{2}+{z}^{2}}$, adopting $R\sqrt{1+{(z/R)}^{2}}\simeq R$. The equations of the motion of the gas can be written by

$$\begin{eqnarray}\displaystyle \frac{\partial {V}_{{\rm{R}}}}{\partial t}+({\boldsymbol{V}}\cdot {\rm{\nabla }}){V}_{{\rm{R}}}-\displaystyle \frac{{V}_{\phi }^{2}}{R} & = & -\displaystyle \frac{{c}_{s}^{2}}{{\rho }_{{\rm{g}}}}\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial R}-\displaystyle \frac{{{GM}}_{* }}{{R}^{2}}\\ & & +\displaystyle \frac{{f}_{R}}{{\rho }_{{\rm{g}}}}-\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\displaystyle \frac{{V}_{{\rm{R}}}-{v}_{{\rm{R}}}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {V}_{\phi }}{\partial t}+({\boldsymbol{V}}\cdot {\rm{\nabla }}){V}_{\phi }+\displaystyle \frac{{V}_{{\rm{R}}}{V}_{\phi }}{R}=-\displaystyle \frac{{f}_{\phi }}{{\rho }_{{\rm{g}}}}-\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\displaystyle \frac{{V}_{\phi }-{v}_{\phi }}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}\displaystyle \frac{\partial {V}_{z}}{\partial t}+({\boldsymbol{V}}\cdot {\rm{\nabla }}){V}_{z} & = & -\displaystyle \frac{{c}_{s}^{2}}{{\rho }_{{\rm{g}}}}\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial z}-\displaystyle \frac{{{GM}}_{* }}{{R}^{3}}z\\ & & +\displaystyle \frac{{f}_{z}}{{\rho }_{{\rm{g}}}}-\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\displaystyle \frac{{V}_{z}-{v}_{z}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 15 }$$

where f_R, ${f}_{\phi }$, and f_z are the viscous forces in radial, azimuthal, and vertical directions, respectively, which are originated by the gas turbulence. Assuming the axisymmetric structure, we can express the viscous forces as

$$\begin{eqnarray}{f}_{R} & = & \displaystyle \frac{2}{R}\displaystyle \frac{\partial }{\partial R}\left[{\nu }_{{rr}}{\rho }_{{\rm{g}}}R\left(\displaystyle \frac{\partial {V}_{{\rm{R}}}}{\partial R}-\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{V}}\right)\right]\\ & & +\displaystyle \frac{\partial }{\partial z}\left[{\nu }_{{rz}}{\rho }_{{\rm{g}}}\left(\displaystyle \frac{\partial {V}_{z}}{\partial z}+\displaystyle \frac{\partial {V}_{{\rm{R}}}}{\partial R}\right)\right]+\displaystyle \frac{{\nu }_{\phi \phi }}{R}\left[2\displaystyle \frac{{V}_{{\rm{R}}}}{R}-\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{V}}\right],\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{f}_{\phi }=\displaystyle \frac{1}{{R}^{2}}\left[\displaystyle \frac{\partial }{\partial R}\left({\nu }_{r\phi }{\rho }_{{\rm{g}}}{R}^{3}\displaystyle \frac{\partial {{\rm{\Omega }}}_{{\rm{g}}}}{\partial R}\right)+\displaystyle \frac{\partial }{\partial z}\left({\nu }_{\phi z}{\rho }_{{\rm{g}}}{R}^{3}\displaystyle \frac{\partial {{\rm{\Omega }}}_{{\rm{g}}}}{\partial z}\right)\right],\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}{f}_{z} & = & \displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left[{\nu }_{{rz}}{\rho }_{{\rm{g}}}R\left(\displaystyle \frac{\partial {V}_{z}}{\partial R}+\displaystyle \frac{\partial {V}_{{\rm{R}}}}{\partial z}\right)\right]\\ & & +2\displaystyle \frac{\partial }{\partial z}\left[{\nu }_{{zz}}{\rho }_{{\rm{g}}}\left(\displaystyle \frac{\partial {V}_{z}}{\partial z}-\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{V}}\right)\right],\end{eqnarray} \tag{ 18 }$$

where ${{\rm{\Omega }}}_{{\rm{g}}}={V}_{\phi }/R$ and ${\nu }_{{ij}}$ indicates the kinetic viscosity associated with the term ${v}_{i}{v}_{j}$ of the Reynolds stress. The efficiency of the turbulence may be different in direction. Hence, we distinguish each component of ${\nu }_{{ij}}$ in Equations (16)–(18). In particular, ${\nu }_{r\phi }$ and ${\nu }_{\phi z}$ would be important, which are related to transport of the angular momentum due to radial and vertical shear motions, respectively. For simplicity, however, we just adopt ${\nu }_{{ij}}=\nu$ in the following. The equations of the motions (10)–(15) are the same as those adopted in Nakagawa et al. (1986), excepting the terms associated with the gas viscosity. However, a parameter range (e.g., dust-to-gas mass ratio, Stokes number of the dust grains) in which the basic equations are valid may not be obvious. We discuss the validity of these equations in Section 4.5.

The continuity equation of the gas is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{g}}}}{\partial t}+{\rm{\nabla }}({\rho }_{{\rm{g}}}{\boldsymbol{V}})=0.\end{eqnarray} \tag{ 19 }$$

For the dust grains, the continuity equation is expressed as

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rho }_{{\rm{d}}}}{\partial t}-{\rm{\nabla }}({\rho }_{{\rm{d}}}{\boldsymbol{v}}+{\boldsymbol{j}})=0,\end{eqnarray} \tag{ 20 }$$

where ${\boldsymbol{j}}$ is the mass flux due to the turbulence of the dust grains (e.g., Cuzzi et al. 1993; Dubrulle et al. 1995; Youdin & Lithwick 2007). From the analogy of the molecular diffusion, assuming the axisymmetric structure, we may obtain ${\boldsymbol{j}}=({j}_{R},{j}_{\phi },{j}_{z})$ as (e.g., Takeuchi & Lin 2002, 2005)

$$\begin{eqnarray}&&{j}_{R}={\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{r\phi }}{{Sc}}\displaystyle \frac{\partial }{\partial R}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{j}_{\phi }=0,\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{j}_{z}={\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{\phi z}}{{Sc}}\displaystyle \frac{\partial }{\partial z}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right).\end{eqnarray} \tag{ 23 }$$

Putting $\partial /\partial t=0$ in the equations of the motion and the continuity equations of the gas and the dust grains described above, we consider the structure of the gas and the dust grains. The radial and vertical gas velocities can be much smaller than ${V}_{K}$, and the deviation of ${V}_{\phi }$ from ${V}_{K}$ is also small, since the gravity of the central star dominates the gas motion. We can put f_r = 0 and f_z = 0 in Equations (13) and (15).

First, we consider the vertical structure of the gas and the dust grains in steady state. Neglecting small advection terms in Equations (12) and (15), we obtain ${v}_{z}-{V}_{z}=-{St}^{\prime} {{\rm{\Omega }}}_{{\rm{K}}}z$ and $d{\rho }_{{\rm{g}}}/{dz}=-{({\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}})({{\rm{\Omega }}}_{{\rm{K}}}/{c}_{s})}^{2}z$, as shown by Nakagawa et al. (1986). Assuming that the vertical gas structure is in hydrostatic equilibrium (${V}_{z}=0$), we obtain ${v}_{z}=-{St}^{\prime} {{\rm{\Omega }}}_{{\rm{K}}}z$. When ${\rho }_{{\rm{g}}}\gg {\rho }_{{\rm{d}}}$, moreover, the vertical structure of the gas density can be assumed by

$$\begin{eqnarray}&&{\rho }_{{\rm{g}}}(R,z)={\rho }_{{\rm{g}}}(R,0)\exp \left(-\displaystyle \frac{{z}^{2}}{2{h}_{{\rm{g}}}^{2}}\right),\end{eqnarray} \tag{ 24 }$$

where ${h}_{{\rm{g}}}$ is the scale height of the gas disk defined as

$$\begin{eqnarray}&&{h}_{{\rm{g}}}(R)={c}_{s}/{{\rm{\Omega }}}_{{\rm{K}}}.\end{eqnarray} \tag{ 25 }$$

Note that when ${\rho }_{{\rm{d}}}\sim {\rho }_{{\rm{g}}}$, the thickness of the gas disk may be smaller than ${h}_{{\rm{g}}}$ given by Equation (25) because the dust grains drag the gas as they sediment toward the midplane (Nakagawa et al. 1986). Since Equation (24) underestimates the gas density in the midplane, we may overestimate the effect of the dust feedback in this case. For simplicity, however, we use Equations (24) and (25) even if ${\rho }_{{\rm{d}}}\gt {\rho }_{{\rm{g}}}$. The validity of this assumption is discussed in Section 4.5. Assuming that the radial variations of physical quantities are given by a power law, as ${\rho }_{{\rm{g}}}(R,0)\propto {R}^{s}$, ${c}_{s}\propto {R}^{q/2}$, the angular velocity of gas (without the dust feedback) is described by (Takeuchi & Lin 2002)

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{g}}}(r,z)={{\rm{\Omega }}}_{{\rm{K}}}\left[1+\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{h}_{{\rm{g}}}}{R}\right)}^{2}\left(p+q+\displaystyle \frac{q}{2}\displaystyle \frac{{z}^{2}}{{h}_{{\rm{g}}}^{2}}\right)\right].\end{eqnarray} \tag{ 26 }$$

Hence, the deviation from ${{\rm{\Omega }}}_{{\rm{K}}}$ is expressed as ${{\rm{\Omega }}}_{{\rm{g}}}\,={{\rm{\Omega }}}_{{\rm{K}}}\sqrt{1-2\eta }$, where

$$\begin{eqnarray}&&\eta (R,z)=-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{h}_{{\rm{g}}}}{R}\right)}^{2}\left(p+q+\displaystyle \frac{q}{2}\displaystyle \frac{{z}^{2}}{{h}_{{\rm{g}}}^{2}}\right).\end{eqnarray} \tag{ 27 }$$

The surface density of the gas is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}(R)={\int }_{-\infty }^{\infty }{\rho }_{{\rm{g}}}{dz}=\sqrt{2\pi }{\rho }_{{\rm{g}}}(R,0){h}_{{\rm{g}}}(R).\end{eqnarray} \tag{ 28 }$$

We obtain the terminal vertical velocity of the dust grains as $-{St}^{\prime} {{\rm{\Omega }}}_{{\rm{K}}}z$ (Nakagawa et al. 1986). The dust grains are settled by the terminal velocity, while the grains diffuse owing to the turbulence. At the steady state, the dust settling is balanced with the turbulence diffusion (Cuzzi et al. 1993; Takeuchi & Lin 2002, 2005; Youdin & Lithwick 2007). In steady state, Equation (20) gives us

$$\begin{eqnarray}&&\displaystyle \frac{1}{R}\displaystyle \frac{\partial }{\partial R}\left[{\rho }_{{\rm{d}}}{v}_{{\rm{R}}}-{\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{r\phi }}{{Sc}}\displaystyle \frac{\partial }{\partial R}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)\right]\\ &&\quad +\displaystyle \frac{\partial }{\partial z}\left[{\rho }_{{\rm{d}}}{v}_{z}-{\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{\phi z}}{{Sc}}\displaystyle \frac{\partial }{\partial z}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)\right]=0.\end{eqnarray} \tag{ 29 }$$

Hence, we obtain

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}{v}_{{\rm{R}}}-{\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{r\phi }}{{Sc}}\displaystyle \frac{\partial }{\partial R}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)={F}_{m,R},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}{v}_{z}-{\rho }_{{\rm{g}}}\displaystyle \frac{{\nu }_{\phi z}}{{Sc}}\displaystyle \frac{\partial }{\partial z}\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)={F}_{m,z},\end{eqnarray} \tag{ 31 }$$

where ${F}_{m,R}$ and ${F}_{m,z}$ are mass fluxes of the dust grains in radial and vertical directions, respectively. Considering the situation that the dust settling is balanced with the diffusion, we put ${F}_{m,z}=0$. Using the terminal velocity of the dust grains, we obtain the vertical distribution of the dust density as (Takeuchi & Lin 2002)

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}(R,z)={\rho }_{{\rm{d}}}(R,0)\exp \left[-\displaystyle \frac{{z}^{2}}{{h}_{{\rm{g}}}^{2}}-\displaystyle \frac{{{St}}_{\mathrm{mid}}}{({\alpha }_{\phi z}/{Sc})}\exp \left(\displaystyle \frac{{z}^{2}}{{h}_{{\rm{g}}}^{2}}-1\right)\right],\end{eqnarray} \tag{ 32 }$$

where ${{St}}_{\mathrm{mid}}$ is the Stokes number of the dust grains at the midplane and ${\alpha }_{\phi z}={\nu }_{\phi z}/({h}_{{\rm{g}}}^{2}{{\rm{\Omega }}}_{{\rm{K}}})$. When the dust grains are relatively settled, we expand the expression of ${h}_{{\rm{d}}}$ with respect to $z/{h}_{{\rm{g}}}\ll 1$ and take the leading term. We obtain

$$\begin{eqnarray}&&{\rho }_{{\rm{d}}}(R,z)={\rho }_{{\rm{d}}}(R,0)\exp \left(-\displaystyle \frac{{z}^{2}}{2{h}_{{\rm{d}}}^{2}}\right),\end{eqnarray} \tag{ 33 }$$

with the scale height of the dust disk given by

$$\begin{eqnarray}&&{h}_{{\rm{d}}}(R)={h}_{{\rm{g}}}(R)\sqrt{\displaystyle \frac{{\alpha }_{\phi z}/{Sc}}{{\alpha }_{\phi z}/{Sc}+{{St}}_{\mathrm{mid}}}}.\end{eqnarray} \tag{ 34 }$$

We set ${Sc}=1$, because we treat the dust grains with ${{St}}_{\mathrm{mid}}\lt 1$ (Youdin & Lithwick 2007). The surface density of the dust grains is given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}(R)={\int }_{-\infty }^{\infty }{\rho }_{{\rm{d}}}{dz}=\sqrt{2\pi }{\rho }_{{\rm{d}}}(R,0){h}_{{\rm{d}}}(R).\end{eqnarray} \tag{ 35 }$$

Using Equations (28), (34), and (35), we obtain the ratio of ${{\rm{\Sigma }}}_{{\rm{d}}}$ to ${{\rm{\Sigma }}}_{{\rm{g}}}$ as

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}=\displaystyle \frac{{\rho }_{{\rm{d}}}(R,0)}{{\rho }_{{\rm{g}}}(R,0)}\sqrt{\displaystyle \frac{\alpha }{\alpha +{{St}}_{\mathrm{mid}}}}.\end{eqnarray} \tag{ 36 }$$

Similar to above, we assume that ${v}_{{\rm{R}}},{v}_{z}$, ${V}_{{\rm{R}}},{V}_{z}$, and ${v}_{\phi }-{V}_{K}$, ${V}_{\phi }-{V}_{K}$ are much smaller than ${V}_{K}$. Leaving only the first-order terms with respect to these small values in Equations (10)–(14), we obtain the velocities of the dust grains as

$$\begin{eqnarray}&&{v}_{{\rm{R}}}(R,z)=-\displaystyle \frac{2{St}^{\prime} }{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}\eta {V}_{K}+\displaystyle \frac{1}{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}{V}_{\mathrm{vis}},\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}{v}_{\phi }(R,z) & = & {V}_{K}+\displaystyle \frac{1}{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}\eta {V}_{K}\\ & & -\displaystyle \frac{{St}^{\prime} }{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}{V}_{\mathrm{vis}},\end{eqnarray} \tag{ 38 }$$

and those of the gas as

$$\begin{eqnarray}{V}_{{\rm{R}}}(R,z) & = & \displaystyle \frac{2{St}^{\prime} }{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}\eta {V}_{K}\\ & & +\left(1-\displaystyle \frac{1}{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}\right){V}_{\mathrm{vis}},\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}{V}_{\phi }(R,z) & = & {V}_{K}+\left(\displaystyle \frac{1}{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}-1\right)\eta {V}_{K}\\ & & +\displaystyle \frac{{St}^{\prime} }{{St}{^{\prime} }^{2}+1}\displaystyle \frac{{\rho }_{{\rm{g}}}}{{\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}}{V}_{\mathrm{vis}},\end{eqnarray} \tag{ 40 }$$

where ${St}^{\prime} ={\rho }_{{\rm{g}}}/({\rho }_{{\rm{g}}}+{\rho }_{{\rm{d}}}){St}$, ${St}={t}_{\mathrm{stop}}{{\rm{\Omega }}}_{{\rm{K}}}$, and η is defined by Equation (27). The viscous velocity of gas (without dust feedback) is given by

$$\begin{eqnarray}&&{V}_{\mathrm{vis}}(R,z)=-\displaystyle \frac{1}{R}\left[3{\nu }_{r\phi }\displaystyle \frac{d\mathrm{ln}({\nu }_{r\phi }{\rho }_{{\rm{g}}}{R}^{1/2})}{d\mathrm{ln}R}-q{\nu }_{\phi z}\displaystyle \frac{d\mathrm{ln}({\rho }_{{\rm{g}}}z)}{d\mathrm{ln}z}\right].\end{eqnarray} \tag{ 41 }$$

If ${\nu }_{r\phi }={\nu }_{\phi z}=\nu$, we obtain ${V}_{\mathrm{vis}}$ as

$$\begin{eqnarray}&&{V}_{\mathrm{vis}}=-\displaystyle \frac{3\nu }{R}\left[p+\displaystyle \frac{2q}{3}+2{\left(\displaystyle \frac{z}{{h}_{{\rm{g}}}}\right)}^{2}\left(\displaystyle \frac{5q+9}{6}\right)\right].\end{eqnarray} \tag{ 42 }$$

We consider the net mass transfers of gas and dust grains in the disk. The net radial velocity of the gas is defined by (Takeuchi & Lin 2002)

$$\begin{eqnarray}&&\langle {V}_{{\rm{R}}}\rangle (R)=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{{\rm{g}}}}{\int }_{-\infty }^{\infty }{V}_{{\rm{R}}}{\rho }_{{\rm{g}}}{dz}.\end{eqnarray} \tag{ 43 }$$

Similarly, the net radial velocity of the dust grain is defined by

$$\begin{eqnarray}&&\langle {v}_{{\rm{R}}}\rangle (R)=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{{\rm{d}}}}{\int }_{-\infty }^{\infty }{v}_{{\rm{R}}}{\rho }_{{\rm{d}}}{dz}.\end{eqnarray} \tag{ 44 }$$

We adopt the following values as the fiducial case: ${M}_{* }=1{M}_{\odot }$, ${\rho }_{{\rm{g}}}(R=1\,\mathrm{au},z=0)=5.3\times {10}^{-10}\,{\rm{g}}\,{\mathrm{cm}}^{-3},$ ${h}_{{\rm{g}}}=2.8\times {10}^{-2}\,\mathrm{au}$ at $R=1\,\mathrm{au}$, $p=-2.25$, and $q=-0.5$. In this case, the gas surface density is obtained by ${{\rm{\Sigma }}}_{0}{(R/1\mathrm{au})}^{-1}$ with ${{\rm{\Sigma }}}_{0}\,=540\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, and the aspect ratio of the gas disk is proportional to ${R}^{1/4}$.

Before presenting the vertical averaged velocities of the gas and the dust, we show typical vertical structures of the gas and the dust grains. Figure 1 shows the vertical distributions of the gas and dust densities (top) and the mass flux density of gas (bottom) in the cases of $\alpha ={10}^{-2}$ and $\alpha ={10}^{-3}$. In this figure, the Stokes number of the dust grain in the midplane is set to be 0.1, and the ratio of ${{\rm{\Sigma }}}_{{\rm{d}}}$ to ${{\rm{\Sigma }}}_{{\rm{g}}}$ is 0.01. As shown by Takeuchi & Lin (2002), the gas moves outward near the midplane even when the dust feedback is not considered. However, the dust feedback makes the outward velocity of gas faster. When $\alpha ={10}^{-2}$, the gas flows inward above the dust layer ($z\gtrsim {h}_{{\rm{g}}}$), while the gas in the dust layer flows outward ($z\lt {h}_{{\rm{g}}}$). Above the dust layer, the inward mass flux density is comparable to that assumed in the 2D disk given by ${{\rm{\Sigma }}}_{{\rm{g}}}{V}_{\mathrm{vis}}^{2{\rm{D}}}/(\sqrt{8\pi }{h}_{{\rm{g}}})$. When $\alpha ={10}^{-3}$, the thin dust layer is formed near the midplane. The gas moves outward near the midplane. The outward mass flux density near the midplane is very large. On the other hand, above the dust layer, the gas flows inward. The inward mass flux above the dust layer is the same level as that in the 2D disk, as in the case of $\alpha ={10}^{-2}$.

Zoom In Zoom Out Reset image size

Figure 2 shows the vertical averaged radial velocities of the gas and the dust grains given by Equations (43) and (44). For comparison, we also plot the gas and dust radial velocity in the 2D disk obtained by Equations (3) and (1). In the figure, we set ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}=0.01$. When $\alpha ={10}^{-2}$, the gas velocity is hardly affected by the grains if the Stokes number is small enough. As the Stokes number increases, however, the infall gas velocity decreases, and when ${{St}}_{\mathrm{mid}}=1$, the infall velocity becomes minimum. As the gas viscosity decreases, the infall velocity due to the viscosity ${V}_{\mathrm{vis}}$ decreases. As a result, the radial gas velocity becomes positive owing to the feedback from the grains, and the gas moves outward. For instance, when $\alpha ={10}^{-3}$, a relatively small grain with ${{St}}_{\mathrm{mid}}=0.1$ can make the gas flow outward, without any pileup of grains (${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\,=0.01$ in the figure). For a smaller viscosity (e.g., $\alpha ={10}^{-4}$), the thickness of the dust layer is very thin. Although the outward gas flux density in the dust layer is very large, the contribution from the inside of the thin dust layer is not significant. In fact, the outward velocity when $\alpha ={10}^{-4}$ is smaller than that when $\alpha ={10}^{-3}$, since the dust feedback is ineffective in the case of the very small viscosity. As seen in the bottom panel of Figure 2, the dust grains always flow inward and the infall velocity is hardly affected by the gas viscosity.

Zoom In Zoom Out Reset image size

When the size of the dust grain is sufficiently small (${{St}}_{\mathrm{mid}}\lt {10}^{-2}$), the gas and dust radial velocities in the 2D disk given by Equations (3) and (1) agree well with the vertically averaged velocities, because the dust layer is not significantly thin. As the size of grains increases, the velocity in the 2D disk deviates from the vertical averaged velocity. In this case, the dust grains are settled and ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}$ near the midplane is much larger than ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$. As a result, the effect of the dust feedback is enhanced owing to the dust settling for relatively large viscosity such as the cases with $\alpha ={10}^{-2}$ and $\alpha ={10}^{-3}$. On the other hand, if α is very small, the effect of the dust feedback is ineffective because the thickness of the dust layer is very thin. When $\alpha ={10}^{-4}$, for instance, the vertical averaged velocity is comparable to that of the 2D disk. Although the vertical dust settling affects the effect of the dust feedback as discussed above, Equations (1) and (3) reasonably reproduce the vertical averaged velocities given by Equations (44) and (43) within a factor of ∼2.

We show the dependence of $\langle {V}_{{\rm{R}}}\rangle$ on ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ and ${{St}}_{\mathrm{mid}}$ when $\alpha ={10}^{-3}$ in Figure 3. If the Stokes number of the dust grains is sufficiently small, $\langle {V}_{{\rm{R}}}\rangle$ is negative regardless of the value of ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$. As the Stokes number of the dust grain reaches unity, $\langle {V}_{{\rm{R}}}\rangle$ increases and the gas can move outward for small ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$. When ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is as large as $\gtrsim 0.5$, $\langle {V}_{{\rm{R}}}\rangle$ becomes smaller as ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ increases. In this case, ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}$ is much larger than unity at the midplane, and hence the gas velocity at the midplane becomes small (see Equation (39)).

Zoom In Zoom Out Reset image size

In the top panel of Figure 4, we illustrate the relation between the dust-to-gas surface density ratio and the Stokes number of the dust grains when $\langle {V}_{{\rm{R}}}\rangle =0$. If the dust-to-gas mass ratio is larger than this critical value, the gas velocity is positive. When $\alpha ={10}^{-4}$, the small dust grains with ${{St}}_{\mathrm{mid}}=0.01$ can make the gas move outward if ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\gtrsim 0.01$. For large dust grains with ${{St}}_{\mathrm{mid}}\simeq 1$, the gas can flow outward when the dust-to-gas surface density ratio is only $\sim {10}^{-4}$ and $\alpha ={10}^{-4}$. Even for the relatively large viscosity, the gas moves outward if relatively large dust grains are highly accumulated (e.g., when $\alpha ={10}^{-2}$, ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\,\sim 0.1$ for dust grains with ${{St}}_{\mathrm{mid}}\simeq 0.1$). If the size of the grains is small enough (e.g., ${{St}}_{\mathrm{mid}}\simeq 0.05$ in the case of $\alpha ={10}^{-2}$), on the other hand, the gas moves inward independent of the dust-to-gas surface density ratio.

Zoom In Zoom Out Reset image size

For comparison, we plot the dust-to-gas surface density when ${V}_{{\rm{R}}}^{2{\rm{D}}}=0$. As seen from the top panel of Figure 4, though the ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ when ${V}_{{\rm{R}}}^{2{\rm{D}}}=0$ is slightly larger than that when $\langle {V}_{{\rm{R}}}\rangle =0$, they reasonably agree with each other, if ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is smaller than unity. Setting ${V}_{{\rm{R}}}^{2{\rm{D}}}=0$ and assuming ${{\rm{\Sigma }}}_{{\rm{d}}}\ll {{\rm{\Sigma }}}_{{\rm{g}}}$, we obtain the condition from Equation (1) as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\right)}_{c}\simeq (1+{{St}}_{2{\rm{D}}}^{-2}){\left(1+\displaystyle \frac{2}{{{St}}_{2{\rm{D}}}}\displaystyle \frac{{\eta }_{2{\rm{D}}}{V}_{K}}{{V}_{\mathrm{vis}}^{2{\rm{D}}}}\right)}^{-1}.\end{eqnarray} \tag{ 45 }$$

As seen in the bottom panel of Figure 4, Equation (45) can reproduce the dust-to-gas surface density ratio when ${V}_{{\rm{R}}}^{2{\rm{D}}}=0$, if ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\lesssim 1$. Equation (45) also reasonably agrees with the dust-to-gas surface density when $\langle {V}_{{\rm{R}}}\rangle =0$, if ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is smaller than unity. Note that the condition of Equation (45) corresponds to the accretion rate of dust grains (${\dot{M}}_{{\rm{d}}}=2\pi R{{\rm{\Sigma }}}_{{\rm{d}}}{v}_{{\rm{R}}}^{2{\rm{D}}}$) being equal to the gas mass accretion rate owing to the viscous diffusion ($2\pi R{{\rm{\Sigma }}}_{{\rm{g}}}{V}_{\mathrm{vis}}^{2{\rm{D}}}$), in the case of ${{St}}_{2{\rm{D}}}\ll 1$.

The vertical distribution of the gas flow is affected by the dust settling, which can enhance (and reduce) the effect of the dust feedback, as shown in previous subsections. However, as seen in Figure 2, when we consider the net flux integrated over the vertical direction, the approximation of the 2D disk may be reasonable. In the following, we focus on the net flux of gas and dust and discuss the evolution of the surface densities.

In steady state, the distribution of the gas surface density is given by ${\dot{M}}_{{\rm{g}}}=2\pi R{{\rm{\Sigma }}}_{{\rm{g}}}{V}_{{\rm{R}}}=\mathrm{constant}$ (Pringle 1981). The gas surface density is given by ${\dot{M}}_{{\rm{g}}}/(2\pi {{RV}}_{{\rm{R}}}^{2{\rm{D}}})$. For a distribution of grains, ${{\rm{\Sigma }}}_{{\rm{d}}}={\dot{M}}_{{\rm{d}}}/(2\pi {{Rv}}_{{\rm{R}}}^{2{\rm{D}}})$ (e.g., Testi et al. 2014). When ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\ll 1$ and ${{St}}_{2{\rm{D}}}\ll 1$, the radial velocity of dust grains is ${v}_{{\rm{R}}}^{2{\rm{D}}}=-2{{St}}_{2{\rm{D}}}({{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}){\eta }_{2{\rm{D}}}{V}_{K}$. For simplicity, we assume single-size dust grains. With ${{\rm{\Sigma }}}_{{\rm{g}}}\propto {R}^{-s}$, we obtain ${{\rm{\Sigma }}}_{{\rm{d}}}\propto {R}^{-(s+2f+1/2)}$, where f is the flaring index defined by ${h}_{{\rm{g}}}/R\propto {R}^{f}$. For ${V}_{{\rm{R}}}^{2{\rm{D}}}$, the first term of Equation (3) is approximately $2{{St}}_{2{\rm{D}}}({{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}){\eta }_{2{\rm{D}}}{V}_{K}$ and the second term is ${V}_{\mathrm{vis}}^{2{\rm{D}}}$. If the feedback is negligible, $s=2f+1/2$. The index of the power law of ${{\rm{\Sigma }}}_{{\rm{g}}}$ does not change if the feedback is effective, because the first term in ${V}_{{\rm{R}}}^{2{\rm{D}}}$ has the same dependence on the radius as ${V}_{\mathrm{vis}}^{2{\rm{D}}}$. Hence, in steady state, ${{\rm{\Sigma }}}_{{\rm{d}}}\propto {R}^{-4f-1}$ and ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\propto {R}^{-2f-1/2}$. In reality, the grain sizes are determined by the coagulation and fragmentation and are not necessarily constant over the disk. When turbulent fragmentation limits the size, however, larger dust grains dominate the total mass, and therefore the dust surface density is determined by the maximum size of the grains, which only weakly depends on the location (${s}_{{\rm{d}}}\propto 1/\sqrt{R}$) when $s=1,f=1/4$ (Birnstiel et al. 2012). This situation is not very different from the case with the single-size dust grains.

When the viscosity is sufficiently large but ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is not very large, the gas flows inward in the disk. As discussed above, the slope of ${{\rm{\Sigma }}}_{{\rm{g}}}$ is the same as that in the case that the feedback is negligible. Hence, the distribution of ${{\rm{\Sigma }}}_{{\rm{g}}}$ does not change so much, as long as the gas is supplied from outside. Note that in this case, the accretion rate of gas ${\dot{M}}_{{\rm{g}}}$ is reduced because the inflow velocity decreases owing to the feedback, but ${{\rm{\Sigma }}}_{{\rm{g}}}$ is not changed. On the other hand, when the viscosity is low or ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is sufficiently large as shown in Figure 4, the gas flows outward in the disk. In this case, the disk structure with ${\dot{M}}_{{\rm{g}}}\lt 0$ is no longer allowed. Since there is no gas supply from the inside of the disk in most cases, the disk may be depleted from the inside.

3.1.1. Two-dimensional Simulations

To examine how the feedback from the grains affects viscous evolutions of gaseous disks, we performed 2D ($R,\phi$) hydrodynamic simulations. We extended the publicly available FARGO code (Masset 2000) to include the dust grains. We simulated the evolutions of disk gas and dust grains by solving the two-fluid (gas and dust grains) equations of motion and continuity. Because we are simulating the 2D disk, we do not consider the settling of the dust grains in the simulations. However, as shown in Figure 2, the velocities of the gas and the dust grains in the 2D disk given by Equations (1)–(4) reasonably agree with the vertical averaged velocity in the 3D disk. Hence, the approximation of the 2D disk would be reasonably valid in this case. The computational domain ranges from 4 to 100 au from the central star. The resolution is 512 and 128 cells in radial and azimuthal directions, respectively.

The initial condition of gas surface density is set as ${{\rm{\Sigma }}}_{{\rm{g}}}={{\rm{\Sigma }}}_{0}{(R/1\mathrm{au})}^{-1}$ with ${{\rm{\Sigma }}}_{0}=570\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. We adopt a simple locally isothermal equation of state and assume a disk aspect ratio of ${h}_{{\rm{g}}}/R={H}_{0}{(R/1\mathrm{au})}^{1/4}$ with H₀ = 0.028. The initial angular velocity of gas is given by ${{\rm{\Omega }}}_{K}\sqrt{1-\eta }$, and the mass of the central star is assumed by $1{M}_{\odot }$. The radial velocity of gas is set by ${V}_{\mathrm{vis}}^{2{\rm{D}}}$. The initial surface density of dust grains is 0.01 of gas everywhere over the disk. The initial angular velocity of dust grains is ${{\rm{\Omega }}}_{K}$, and the initial radial velocity is zero. We adopt dust grains with 3 cm in size, corresponding to the Stokes number ${{St}}_{2{\rm{D}}}$ of 0.1 at $R=10\,\mathrm{au}$ for the initial state. For simplicity, we neglect coagulation and fragmentation of dust grains, and hence the grain size does not change during the simulations.

Since the gas velocity is very sensitive to the value of ${\eta }_{2{\rm{D}}}$ and ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ especially in the case of small viscosity, numerical instability occurs when small discontinuities of η and ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ exist at the innermost part of the computational domain. To avoid this instability, we introduced a "coupling damping region" at the innermost radii ($4\,\mathrm{au}\lt R\lt 6\,\mathrm{au}$). In this region the dust feedback on gas is gradually reduced by $\cos (\pi {x}^{2}/2)$, where $x={R}_{1}-R/({R}_{1}-{R}_{0})$ with ${R}_{0}=4\,\mathrm{au}$ and ${R}_{1}=6\,\mathrm{au}$, respectively. Hence, at the innermost annulus ($R=4\,\mathrm{au}$), the gas velocity is set to ${V}_{\mathrm{vis}}^{2{\rm{D}}}$. Moreover, in this region, we force all the physical quantities to be azimuthally symmetric by overwriting the quantities with their azimuthal average at every time step (see de Val-Borro et al. 2006). That is, in the coupling damping zone, ${{\rm{\Sigma }}}_{{\rm{g}}}$, ${V}_{{\rm{R}}}^{2{\rm{D}}}$, and ${V}_{\phi }^{2{\rm{D}}}$ are related to their azimuthally averaged values as

$$\begin{eqnarray}&&\displaystyle \frac{{dX}}{{dt}}=-\displaystyle \frac{X-{X}_{\mathrm{avg}}}{{\tau }_{\mathrm{damp}}}f(R),\end{eqnarray} \tag{ 46 }$$

where X represents ${{\rm{\Sigma }}}_{{\rm{g}}}$, ${V}_{{\rm{R}}}^{2{\rm{D}}}$, and ${V}_{\phi }^{2{\rm{D}}}$; ${X}_{\mathrm{avg}}$ denotes the azimuthally averaged values of them; and ${\tau }_{\mathrm{damp}}$ is the orbital period at the boundary. The function f is a parabola of the form $y={x}^{2}$, scaled to be zero at the boundary layer ($R=6\,\mathrm{au}$) and unity at the opposite edge ($R=4\,\mathrm{au}$). For dust grains, the velocity at the innermost boundary is given by Equations (37) and (38), respectively. The surface densities of gas and dust are set to be $\dot{M}=\mathrm{constant}$ (see Zhu et al. 2011). At the outer boundary, the velocities of gas and dust are given by Equations (37)–(40), and the surface densities are fixed at the initial values. Here the radial distribution of ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is expected to be smooth. The turbulent radial mass flux of j_R would be negligible in this case. Hence, we ignore the radial diffusion flux of Equation (21) for simplicity.

3.1.2. One-dimensional Simulations

Figure 5 illustrates evolutions of the radial velocity and surface densities of gas and dust grains, as well as the dust-to-gas surface density ratio in the viscous case with $\alpha ={10}^{-2}$. We compare the results with the 1D model in the figure. In this case, since the gas viscosity is large and the dust grains are not so concentrated, the gas flows to the central star in the entire region of the disk. As discussed in Section 2.5, the distribution of ${{\rm{\Sigma }}}_{{\rm{g}}}$ is hardly changed from the initial distribution, which is the steady state without the feedback. The ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is also distributed as $1/R$, as expected in Section 2.5. Since the velocities of the gas and the dust quickly converge to the values in steady state, the simplified 1D model is able to well reproduce the results of the 2D two-fluid hydrodynamic simulations.

Zoom In Zoom Out Reset image size

We now show results with different α and ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$. The agreement between the 1D model and 2D simulations is always good in the cases presented below. Figure 6 shows the evolution of the disk with $\alpha ={10}^{-2}$ and initial ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}=0.05$. Because the initial ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is larger than the case of Figure 5, the gas moves outward in the wide region of the disk. In this case, the gas surface density decreases in the inner region of the disk. Owing to the gas removal due to the feedback from the grains, the dust-to-gas mass ratio increases. The drift velocity of grains decreases as the dust-to-gas mass ratio increases, as pointed out by Dra̧żkowska et al. (2016) and Ida & Guillot (2016) (see also Equation (1)), which leads to the further increase of the dust-to-gas mass ratio. Owing to this positive feedback cycle, the inner region of the disk quickly becomes very dust rich. In this case, after $\sim {10}^{5}\,\mathrm{yr}$, the gas surface density is only $\sim 10 \%$ of the initial value at the inner region. The dust-to-gas mass ratio significantly exceeds unity, and the Stokes number of grains is ∼0.1. The dust-to-gas density ratio is also very high. Using Equation (36), we can estimate ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}$ at the midplane at $10\,\mathrm{au}$ as 50 (at $9.6\times {10}^{4}\,\mathrm{yr}$) and 150 (at $1.6\times {10}^{5}\,\mathrm{yr}$). Note that the decreases of ${{\rm{\Sigma }}}_{{\rm{d}}}$ near the inner boundary ($R=6\,\mathrm{au}$) originated from the difference of ${\eta }_{2{\rm{D}}}$ in the inside and outside of the coupling damping zone, since ${\eta }_{2{\rm{D}}}$ in the damping zone is fixed to the value in the steady state without the feedback.

Zoom In Zoom Out Reset image size

The evolution with a lower viscosity ($\alpha ={10}^{-3}$) is shown in Figure 7. The initial ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is the same as in Figure 5. Since the viscosity is small, in this case the gas moves outward in the entire disk. The gas surface density also significantly decreases at the inner region of the disk. As in the case of Figure 6, the removal of gas due to the dust feedback makes the inner region dusty. After $\sim {10}^{5}\,\mathrm{yr}$, the dust-to-gas surface density ratio increases up to unity. The dust-to-gas density ratios at the midplane at $10\,\mathrm{au}$ are estimated as 1 (at $9.6\times {10}^{4}\,\mathrm{yr}$) and 10 (at $1.6\times {10}^{5}\,\mathrm{yr}$).

Zoom In Zoom Out Reset image size

The gas flows inward when the net radial velocity of the gas $\langle {V}_{{\rm{R}}}\rangle$ is positive. As shown in Figure 1, the mass flux density depends on the altitude. The dust grains are settled to the midplane, and the dust layer in which the dust grains are highly concentrated is formed near the midplane, while the dust density is quite small above the dust layer. When the net mass flux is positive, the gas in the upper part of the dust layer still flows inward, whereas the gas in the dust layer flows outward. In Figure 8, we illustrate the mass fluxes integrated in the inside and the upper part of the dust layer, in terms of the Stokes number of the dust grains, when ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}=0.01$ at $10\,\mathrm{au}$ (${\rho }_{{\rm{g}}}=3.0\times {10}^{-12}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, ${h}_{{\rm{g}}}/R=0.05$). For relatively large dust grains (e.g., ${{St}}_{\mathrm{mid}}\sim 0.1$), the net radial velocity of the gas is positive. However, above the dust layer ($z\gt 2{h}_{{\rm{d}}}$), the gas flows inward. When $\alpha ={10}^{-3}$ and ${{St}}_{\mathrm{mid}}=1$, for instance, the amount of the gas mass flux above the dust layer is $\sim -{10}^{-11}\ {M}_{\odot }\,\,{\mathrm{au}}^{-1}\,\,{\mathrm{yr}}^{-1}$, which is comparable to the mass flux without the dust feedback ($={{\rm{\Sigma }}}_{{\rm{g}}}{V}_{\mathrm{vis}}^{2{\rm{D}}}/2$). This indicates that even when the disk gas depletes from the inside of the disk, such as in Figures 6 and 7, the gas accretion onto the central star would not stop. If this disk is observed, we cannot find the depletion of the disk gas, from the viewpoint of the accretion rate onto the central star. We need to directly observe the disk gas to identify the physical condition of the protoplanetary disk.

Zoom In Zoom Out Reset image size

In reality, the grain size is strongly limited by the coagulation and fragmentation, though we did not consider this effect in this paper. However, for silicate grains, since the size is strongly limited by the fragmentation, the size of grains that locally dominates the grain density is not changed so much in the inner region (Birnstiel et al. 2012). In this case, the maximum size of the dust grains ${s}_{\mathrm{frag}}$ is given by

$$\begin{eqnarray}&&{s}_{\mathrm{frag}}\sim \displaystyle \frac{2}{3\pi }\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{g}}}}{{\rho }_{p}\alpha }{\left(\displaystyle \frac{{v}_{\mathrm{frag}}}{{c}_{s}}\right)}^{2},\end{eqnarray} \tag{ 47 }$$

where ${v}_{\mathrm{frag}}$ is the fragmentation threshold velocity, which depends on the composition (e.g., ice or silicate) of the dust grains. The fragmentation threshold velocity of the silicate dust grains may be obtained as 1–10 m s⁻¹ (Beitz et al. 2011). When $\alpha ={10}^{-3}$ and ${c}_{s}={10}^{3}\,{\rm{m}}\,{{\rm{s}}}^{-1}$, the Stokes number of the maximum dust grains with ${v}_{\mathrm{frag}}=10\,{\rm{m}}\,{{\rm{s}}}^{-1}$ is about ∼0.1. In this case, the dust feedback can influence the entire evolution of the gas disk, as shown in Figures 6 and 7, which leads to formation of rocky planetesimals.

If ${v}_{\mathrm{frag}}\ll 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$, the dust feedback does not work, since the dust grains cannot grow to a sufficiently large size. If ${v}_{\mathrm{frag}}\gg 10\,{\rm{m}}\,{{\rm{s}}}^{-1}$, the dust grains can grow up to planetesimals without fragmentation, as shown by Okuzumi et al. (2012), or the dust quickly falls to the star because they can quickly grow up to the size of ${St}\sim 1$.

When the dust grains are more concentrated as compared with the gas, the grains can grow quickly to planetesimals via the streaming instability (e.g., Youdin & Goodman 2005). Roughly speaking, when ${\rho }_{{\rm{d}}}\gtrsim {\rho }_{{\rm{g}}}$, the streaming instability may be significantly developed (e.g., Youdin & Johansen 2007; Dra̧żkowska & Dullemond 2014; Carrera et al. 2015). As shown in Figures 6 and 7, when the dust feedback is effective, after ${10}^{5}\,\mathrm{yr}$, ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ increases over unity within $10\,\mathrm{au}$. Using Equation (36), we can estimate ${\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}}$ at the midplane as $\gt 10$, because ${St}\simeq 0.1$. This ratio is sufficiently high for the streaming instability to occur, though the accurate condition of the streaming instability may depend on the viscosity, the total amount of the dust grains (metallicity), etc. The streaming instability can be caused in an earlier state of the disk evolution than ${10}^{5}\,\mathrm{yr}$. We should consider the planetesimal formation with the disk evolution.

Once the streaming instability occurs, since the dust grains grow up to the planetesimals, the dust-to-gas mass ratio will decrease. However, all dust grains are not converted to the planetesimals (e.g., Johansen et al. 2012; Simon et al. 2016). If the efficiency of the planetesimal formation is sufficiently large to remove the dust grains, the dust feedback will be ineffective. The gas density will be restored to that in steady state when the dust feedback is not considered. On the other hand, when the efficiency is moderate, the dust feedback can still be effective after the streaming instability turns on. In this case, the planetesimals are formed as the gas depletes, which indicates that the migration speed of the planetesimals would significantly decrease. If the planet is formed, the migration of the planet would also be slow, which is favorable for planet formation. This situation can continue until most of the dust grains fall to the central star.

Since the disk gas is depleted from the inner region, a hole structure of gas would be formed. For the dust grains, on the other hand, the density does not decrease in the gas-depleted region. However, if the dust grains are effectively translated to the planetesimals, as discussed in the previous subsection, a density of relatively large grains at the midplane that are observed by submillimeter may decrease. In this case, a hole structure can be found by the observations of dust continuum by submillimeter, as well as the observation of molecular lines. Recent ALMA observation done by Sheehan & Eisner (2017) has revealed that the young protoplanetary disk has the hole structure of the dust continuum within $\sim 10\,\mathrm{au}$. This hole structure may be explained by the gas depletion due to the dust feedback and the planetesimal formation triggered by the gas depletion.

We briefly comment on the validity of the treatment of the dust fluid. We treat the dust grains as pressureless fluid and ignore some effects of turbulence in the equations of motion (Equations (10)–(12); see Appendix A). Although this formulation is widely used in the previous studies (e.g., Fu et al. 2014; Dipierro & Laibe 2017; Gonzalez et al. 2017), this treatment of the dust fluid may not be appropriate in the cases where the dust density is comparable to the gas density, or the Stokes number of the dust grains is large. According to Hersant (2009), the pressure of the dust fluid is not negligible when ${St}\gt 1/2$. At the midplane, where the dust density is comparable to the gas density, the dust settling is prevented by the turbulence driven by the instability of the dust fluid (e.g., Sekiya 1998; Sekiya & Ishitsu 2000; Johansen et al. 2006; Chiang 2008; Lee et al. 2010a, 2010b). We must take into account the momentum transfer due to the dust turbulence in this case, and a more sophisticated formulation for dust and gas may need to be explored.

Moreover, the vertical structures of the gas and the dust grains may be different from what we have assumed in Equations (24) and (33) when the dust density is larger than the gas density. Gas experiences drag force from the dust grains that sediment toward the midplane, so the gas density in the midplane is larger than that given by Equation (24). As the dust-to-gas density ratio increases, the settling velocity of the dust grains slows down (${v}_{z}=-{St}^{\prime} {{\rm{\Omega }}}_{{\rm{K}}}z\propto 1/(1+{\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}})$), which leads to vertical swelling of the dust layer. As a result, the vertical distributions of the gas and the dust also deviate from the Gaussian distributions when ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}$ is larger than unity (see Appendix A of Nakagawa et al. 1986, or Appendix B). When adopting Equations (24) and (33), we may therefore have overestimated the dust feedback near a bottom of the dust layer (midplane), while underestimating the dust feedback in an upper region of the dust layer. If considering the vertical averaged quantities, the approximation of the Gaussian distribution may be reasonable. In any case, if the dust-to-gas density ratio in the midplane is of the order of unity and the Stokes number of the dust grains is of the order of unity, the drag force exerted on the gas is at most comparable to the (vertical component of) gravitational force from the central star, and therefore we expect that the qualitative outcomes presented in this paper are not very much affected. If the dust-to-gas density ratio further increases, it may be necessary to consider more sophisticated formulation of basic equations (e.g., onset of the streaming instability), which is beyond the scope of this paper.

In summary, our results may be safely used when the dust-to-gas mass ratio is less than ∼1, and we expect that qualitative results are valid up to ${{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\sim 1$. If the dust grains are highly concentrated and dominate gas, it may be necessary to have a more sophisticated formulation to describe the dynamics of gas and dust. Streaming instability or other hydrodynamic instabilities may come into play when dust grains dominate, in which case the dust-to-gas mass ratio may not increase too much after all.