Effect of dust radial drift on viscous evolution of gaseous disk

TL;DR

This paper investigates how dust feedback influences the viscous evolution of protoplanetary disks, revealing conditions under which gas flows outward and impacts planetesimal formation environments.

Contribution

It introduces a 2D hydrodynamic model incorporating dust feedback and vertical settling, demonstrating its significant effects on gas flow and dust accumulation in disks.

Findings

Dust feedback can reverse gas flow from inward to outward.

Gas density within 10 AU can be significantly depleted.

Inner disk regions may have dust-to-gas ratios exceeding unity.

Abstract

The total amount of dust (or "metallicity") and the dust distribution in protoplanetary disks are crucial for planet formation. Dust grains radially drift due to gas--dust friction, and the gas is affected by the feedback from dust grains. We investigate the effects of the feedback from dust grains on the viscous evolution of the gas, taking into account the vertical dust settling. The feedback from the grains pushes the gas outward. When the grains are small and the dust-to-gas mass ratio is much smaller than unity, the radial drift velocity is reduced by the feedback effect but the gas still drifts inward. When the grains are sufficiently large or piled-up, the feedback is so effective that forces the gas flows outward. Although the dust feedback is affected by dust settling, we found that the 2D approximation reasonably reproduces the vertical averaged flux of gas and dust. We also performed the 2D two-fluid hydrodynamic simulations to examine the effect of the feedback from the grains on the evolution of the gas disk. We show that when the feedback is effective, the gas flows outward and the gas density at the region within $\sim 10\ \mbox{AU}$ is significantly depleted. As a result, the dust-to-gas mass ratio at the inner radii may significantly excess unity, providing the environment where planetesimals are easily formed via, e.g., streaming instability. We also show that a simplified 1D model well reproduces the results of the 2D two-fluid simulations, which would be useful for future studies.

Equations

$$v_{\mathrm{R}}^{2D}=-\frac{2S{t}_{2D}^{\prime }}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{g}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{\eta }_{2D}{V}_K$$

$$+\frac{1}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{g}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{V}_{\mathrm{vis}}^{2D},$$

$$v_{\varphi }^{2D}=V_K+\frac{1}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{g}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{\eta }_{2D}{V}_K$$

$$-\frac{S{t}_{2D}^{\prime }}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{g}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{V}_{\mathrm{vis}}^{2D},$$

$$V_{\mathrm{R}}^{2D}=\frac{2S{t}_{2D}^{\prime }}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{\eta }_{2D}{V}_K$$

$$+\left(1-\frac{1}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}\right)V_{\mathrm{vis}}^{2D},$$

$$V_{\varphi }^{2D}=V_K+\left(\frac{1}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}-1\right){\eta }_{2D}{V}_K$$

$$+\frac{S{t}_{2D}^{\prime }}{S{t}_{2D}^{\prime }{}^2+1}\frac{{\Sigma }_{\mathrm{g}}}{{\Sigma }_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}}{V}_{\mathrm{vis}}^{2D},$$

$$S{t}_{2D}=t_{\mathrm{stop}}^{2D}{\Omega }_{\mathrm{K}},$$

$$t_{\mathrm{stop}}(R,z)$$

$$t_{\mathrm{stop}}^{2D}(R)$$

$${\eta }_{2D}=-\frac{1}{2}{\left(\frac{h_{\mathrm{g}}}{R}\right)}^2\left(\frac{d\ln {\Sigma }_{\mathrm{g}}}{d\ln R}+\frac{d\ln c_s^2}{d\ln R}\right),$$

$$V_{\mathrm{vis}}^{2D}=-\frac{3\nu }{R}\frac{d\ln \left(\nu {\Sigma }_{\mathrm{g}}{R}^{1/2}\right)}{d\ln R},$$

$$\frac{\partial v_{\mathrm{R}}}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)v_{\mathrm{R}}-\frac{v_{\varphi }^2}{R}=-\frac{G{M}_{\ast }}{R^2}-\frac{v_{\mathrm{R}}-V_{\mathrm{R}}}{t_{\mathrm{stop}}},$$

$$\frac{\partial v_{\varphi }}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)v_{\varphi }+\frac{v_{\mathrm{R}}{v}_{\varphi }}{R}=-\frac{v_{\varphi }-V_{\varphi }}{t_{\mathrm{stop}}},$$

$$\frac{\partial v_z}{\partial t}+\left(\mathbf{v}\cdot \nabla \right)v_z=-\frac{G{M}_{\ast }}{R^3}z-\frac{v_z-V_z}{t_{\mathrm{stop}}},$$

$$\frac{\partial V_{\mathrm{R}}}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)V_{\mathrm{R}}-\frac{V_{\varphi }^2}{R}=-\frac{c_s^2}{{\rho }_{\mathrm{g}}}\frac{\partial {\rho }_{\mathrm{g}}}{\partial R}-\frac{G{M}_{\ast }}{R^2}$$

$$+\frac{f_R}{{\rho }_{\mathrm{g}}}-\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{V_{\mathrm{R}}-v_{\mathrm{R}}}{t_{\mathrm{stop}}},$$

$$\frac{\partial V_{\varphi }}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)V_{\varphi }+\frac{V_{\mathrm{R}}{V}_{\varphi }}{R}=-\frac{f_{\varphi }}{{\rho }_{\mathrm{g}}}-\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{V_{\varphi }-v_{\varphi }}{t_{\mathrm{stop}}},$$

$$\frac{\partial V_z}{\partial t}+\left(\mathbf{V}\cdot \nabla \right)V_z=-\frac{c_s^2}{{\rho }_{\mathrm{g}}}\frac{\partial {\rho }_{\mathrm{g}}}{\partial z}-\frac{G{M}_{\ast }}{R^3}z+\frac{f_z}{{\rho }_{\mathrm{g}}}$$

$$-\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{V_z-v_z}{t_{\mathrm{stop}}},$$

$$f_R=\frac{2}{R}\frac{\partial }{\partial R}\left[{\nu }_{rr}{\rho }_{\mathrm{g}}R\left(\frac{\partial V_{\mathrm{R}}}{\partial R}-\frac{1}{3}\nabla \cdot \mathbf{V}\right)\right]$$

$$+\frac{\partial }{\partial z}\left[{\nu }_{rz}{\rho }_{\mathrm{g}}\left(\frac{\partial V_z}{\partial z}+\frac{\partial V_{\mathrm{R}}}{\partial R}\right)\right]+\frac{{\nu }_{\varphi \varphi }}{R}\left[2\frac{V_{\mathrm{R}}}{R}-\frac{1}{3}\nabla \cdot \mathbf{V}\right],$$

$$f_{\varphi }$$

$$f_z$$

$$+2\frac{\partial }{\partial z}\left[{\nu }_{zz}{\rho }_{\mathrm{g}}\left(\frac{\partial V_z}{\partial z}-\frac{1}{3}\nabla \cdot \mathbf{V}\right)\right],$$

$$\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\nabla \left({\rho }_{\mathrm{g}}\mathbf{V}\right)=0,$$

$$\frac{\partial {\rho }_{\mathrm{d}}}{\partial t}-\nabla \left({\rho }_{\mathrm{d}}\mathbf{v}+\mathbf{j}\right)=0,$$

$$j_R$$

$$j_{\varphi }$$

$$j_z$$

$${\rho }_{\mathrm{g}}(R,z)$$

$$h_{\mathrm{g}}(R)$$

$${\Omega }_{\mathrm{g}}(r,z)={\Omega }_{\mathrm{K}}\left[1+\frac{1}{2}{\left(\frac{h_{\mathrm{g}}}{R}\right)}^2\left(p+q+\frac{q}{2}\frac{z^2}{h_{\mathrm{g}}^2}\right)\right].$$

$$\eta (R,z)$$

$${\Sigma }_{\mathrm{g}}(R)$$

$$\frac{1}{R}\frac{\partial }{\partial R}\left[{\rho }_{\mathrm{d}}{v}_{\mathrm{R}}-{\rho }_{\mathrm{g}}\frac{{\nu }_{r\varphi }}{Sc}\frac{\partial }{\partial R}\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\right)\right]$$

$$+\frac{\partial }{\partial z}\left[{\rho }_{\mathrm{d}}{v}_z-{\rho }_{\mathrm{g}}\frac{{\nu }_{\varphi z}}{Sc}\frac{\partial }{\partial z}\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\right)\right]=0.$$

$${\rho }_{\mathrm{d}}{v}_{\mathrm{R}}-{\rho }_{\mathrm{g}}\frac{{\nu }_{r\varphi }}{Sc}\frac{\partial }{\partial R}\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\right)=F_{m,R},$$

Full text

1. Introduction

Terrestrial planets are formed by the accumulation of dust grains in protoplanetary disks (e.g., Youdin & Goodman, 2005; Tanaka et al., 2005; Okuzumi et al., 2012). In core-accretion scenario (e.g., Mizuno, 1980; Kanagawa & Fujimoto, 2013), the accumulation of dust grains is also critical for giant planet formation. How dust grains evolve in protoplanetary disk is directly connected with planet formation. Moreover, the evolution of dust grains may explain ring structures in protoplanetary disk (e.g., Zhang et al., 2015; Okuzumi et al., 2016), which are observed in several disks (e.g., ALMA Partnership et al., 2015; Momose et al., 2015; Akiyama et al., 2015; Nomura et al., 2016; Tsukagoshi et al., 2016). The evolution of gas and dust grains is one of the most important topics both from theoretical and observational point of view.

Many authors have investigated evolution of dust grains in protoplanetary disk using semi-analytical models (e.g., Youdin & Shu, 2002; Dullemond & Dominik, 2005; Tanaka et al., 2005; Birnstiel et al., 2012; Okuzumi et al., 2012) and using two-fluid hydrodynamic simulations (e.g., Zhu et al., 2012; Pinilla et al., 2012; Dipierro et al., 2015; Picogna & Kley, 2015; Jin et al., 2016; Rosotti et al., 2016). Dust grains radially drift due to gas–dust friction, and disk gas feels feedback from drifting dust grains. Most of previous studies considered gas–dust friction only for dust grains and ignored the dust feedback on the disk gas, because grains are not large and the amount of them is negligible as compared with gas. If the dust grains are highly accumulated or the grains grows up, however, the feedback from dust grains may not be negligible (e.g., Fu et al., 2014; Gonzalez et al., 2015; Taki et al., 2016; Gonzalez et al., 2017; Dipierro & Laibe, 2017).

The feedback from the dust grains affects gas structures. For instance, Kretke et al. (2009) showed that the inner edge of the dead zone is oscillated by the feedback. Taki et al. (2016) found that the pressure bump in gas is deformed by the feedback from the dust grains trapped in the gas bump. Recently, using two-fluid (gas and dust grains) SPH simulations, Gonzalez et al. (2017) demonstrated that the grains are trapped in the pressure bump formed by the feedback. In this paper, we examine how the feedback influences the global evolution of the gas disk induced by viscosity. In Section 2, we describe the structures of gas and dust grains in steady state, considering the vertical settling of the dust grains. In Section 3, we show the results of 2D two-fluid hydrodynamic simulations and a simplified 1D model, and discuss the effects of the dust feedback on the disk evolution and the planet formation in Section 4. Section 5 contains our summary.

2. Structure in steady state

2.1. Velocities of gas and dust in 2D 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, $\Sigma_{\rm g}$ and $\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 the similar formulae. From the 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

[TABLE]

and those for the disk gas are given by

[TABLE]

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

[TABLE]

where $t_{\rm stop}^{2D}$ is the stopping time of dust grains in the 2D disk, and $\Omega_{\rm K}=\sqrt{GM_{\ast}/R^{3}}$ is the Keplerian angular velocity at the mid-plane, where $G$ and $M_{\ast}$ is the gravity constant and the mass of the central star, respectively. The Keplerian rotation velocity at the mid-plane is $V_{K}=R\Omega_{\rm K}$. In the Epstein regime, the stopping time is written by (Takeuchi & Artymowicz, 2001),

[TABLE]

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

[TABLE]

The pressure gradient force is parameterized by

[TABLE]

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

[TABLE]

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

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

We have employed 2D approximation, meaning that the vertical structures of the gas and dust are similar with 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 mid-plane. In the next subsection, we discuss the gas and dust velocities, taking into account the vertical structure.

2.2. Gas and dust motions in a 3D disk

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

[TABLE]

where we assume axis-symmetric gravitational potential $\Psi=-GM_{\ast}/\sqrt{R^{2}+z^{2}}$, adopting $R\sqrt{1+\left(z/R\right)^{2}}\simeq R$. The equations of the motion of the gas can be written by

[TABLE]

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

[TABLE]

[TABLE]

where $\Omega_{\rm g}=V_{\phi}/R$, and $\nu_{ij}$ indicates the kinetic viscosity associated with the term of $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 with 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 these 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, etc.) 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

[TABLE]

For the dust grains, the continuity equation is expressed as

[TABLE]

where $\mathbf{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 $\mathbf{j}=(j_{R},j_{\phi},j_{z})$ as (e.g., Takeuchi & Lin, 2002, 2005)

[TABLE]

2.3. Gas and dust structures in steady state

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}\Omega_{\rm K}z$ and $d\rho_{\rm g}/dz=-(\rho_{\rm g}+\rho_{\rm d})(\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}\Omega_{\rm K}z$. When $\rho_{\rm g}\gg\rho_{\rm d}$, moreover, the vertical structure of the gas density can be assumed by

[TABLE]

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

[TABLE]

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 towards the mid-plane (Nakagawa et al., 1986). Since equation (24) underestimates the gas density in the mid-plane, 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}>\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 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)

[TABLE]

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

[TABLE]

The surface density of the gas is given by

[TABLE]

We obtain the terminal vertical velocity of the dust grains as $-St^{\prime}\Omega_{\rm K}z$ (Nakagawa et al., 1986). The dust grains are settled by the terminal velocity, while the grains diffuse due 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

[TABLE]

Hence, we obtain

[TABLE]

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),

[TABLE]

where $St_{\rm mid}$ is the Stokes number of the dust grains at the mid-plane and $\alpha_{\phi z}=\nu_{\phi z}/(h_{\rm g}^{2}\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

[TABLE]

with the scale height of the dust disk given by

[TABLE]

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

[TABLE]

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

[TABLE]

Similar to above, we assume $v_{\rm R},v_{z}$, $V_{\rm R},V_{z}$, and $v_{\rm\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), (11), (13), and (14), we obtain the velocities of the dust grains as

[TABLE]

and those of the gas as

[TABLE]

where $St^{\prime}=\rho_{\rm g}/(\rho_{\rm g}+\rho_{\rm d})St$, $St=t_{\rm stop}\Omega_{\rm K}$, and $\eta$ is defined by equation (27). The viscous velocity of gas (without dust feedback) is given by

[TABLE]

If $\nu_{r\phi}=\nu_{\phi z}=\nu$, we obtain $V_{\rm vis}$ as

[TABLE]

2.4. Radial net flows of gas and dust grains

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

[TABLE]

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

[TABLE]

We adopt the following values as the fiducial case: $M_{\ast}=1M_{\odot}$, $\rho_{\rm g}(R=1\ \rm{AU},z=0)=5.3\times 10^{-10}\rm{g/cm}^{3}$, $h_{\rm g}=2.8\times 10^{-2}\ \rm{AU}$ at $R=1\ \rm{AU}$, $p=-2.25$ and $q=-0.5$. In this case, the gas surface density is obtained by $\Sigma_{0}(R/1\ \rm{AU})^{-1}$ with $\Sigma_{0}=540\rm{g/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 mid-plane is set to be $0.1$, and the ratio of $\Sigma_{\rm d}$ to $\Sigma_{\rm g}$ is $0.01$. As shown by Takeuchi & Lin (2002), the gas moves outward near the mid-plane 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<h_{\rm g}$). Above the dust layer, the inward mass flux density is comparable with that assumed in the 2D disk given by $\Sigma_{\rm g}V_{\rm vis}^{2D}/(\sqrt{8\pi}h_{\rm g})$. When $\alpha=10^{-3}$, the thin dust layer is formed near the mid-plane. The gas moves outward near the mid-plane. The outward mass flux density near the mid-plane 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}$.

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 2D disk obtained by equation (3) and (1). In the figure, we set $\Sigma_{\rm d}/\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_{\rm mid}=1$, the infall velocity becomes minimum. As the gas viscosity decreases, the infall velocity due to the viscosity $V_{\rm vis}$ decreases. As a result, the radial gas velocity becomes positive due to the feedback from the grains and the gas moves outward. For instance, when $\alpha=10^{-3}$, a relatively small grain with $St_{\rm mid}=0.1$ can make the gas flow outward, without any pile-up of grains ($\Sigma_{\rm d}/\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 flows inward and the infall velocity is hardly affected by the gas viscosity.

When the size of the dust grain is sufficiently small ($St_{\rm mid}<10^{-2}$), the gas and dust radial velocities in the 2D disk given by equation (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 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 mid-plane is much larger than $\Sigma_{\rm d}/\Sigma_{\rm g}$. As a result, the effect of the dust feedback is enhanced due 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 $\alpha$ 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 with that of the 2D disk. Although the vertical dust settling affects the effect of the dust feedback as discussed above, equation (1) and (3) reasonably reproduce the vertical averaged velocities given by equation (44) and (43) within a factor of $\sim 2$.

We show the dependence of $\left<V_{\rm R}\right>$ on $\Sigma_{\rm d}/\Sigma_{\rm g}$ and $St_{\rm mid}$ when $\alpha=10^{-3}$ in Figure 3. If the Stokes number of the dust grains is sufficiently small, $\left<V_{\rm R}\right>$ is negative regardless of the value of $\Sigma_{\rm d}/\Sigma_{\rm g}$. As the Stokes number of the dust grain reaches unity, $\left<V_{\rm R}\right>$ increases and the gas can move outward for small $\Sigma_{\rm d}/\Sigma_{\rm g}$. When $\Sigma_{\rm d}/\Sigma_{\rm g}$ is large as$\gtrsim 0.5$, $\left<V_{\rm R}\right>$ becomes smaller as $\Sigma_{\rm d}/\Sigma_{\rm g}$ increases. In this case, $\rho_{\rm d}/\rho_{\rm g}$ is much larger than unity at the mid-plane, and hence the gas velocity at the mid-plane becomes small (see, equation 39).

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 $\left<V_{\rm R}\right>=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_{\rm mid}=0.01$ can make the gas move outward if $\Sigma_{\rm d}/\Sigma_{\rm g}\gtrsim 0.01$. For large dust grains with $St_{\rm 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}$, $\Sigma_{\rm d}/\Sigma_{\rm g}\sim 0.1$ for dust grains with $St_{\rm mid}\simeq 0.1$). If the size of the grains is small enough (e.g., $St_{\rm 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.

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

[TABLE]

As seen in the bottom of Figure 4, equation (45) can reproduce the dust-to-gas surface density ratio when $V_{\rm R}^{2D}=0$, if $\Sigma_{\rm d}/\Sigma_{\rm g}\lesssim 1$. Equation (45) also reasonably agrees with the dust-to-gas surface density when $\left<V_{\rm R}\right>=0$, if $\Sigma_{\rm d}/\Sigma_{\rm g}$ is smaller than unity. Note that the condition of equation (45) corresponds to that the accretion rate of dust grains ($\dot{M}_{\rm d}=2\pi R\Sigma_{\rm d}v_{\rm R}^{2D}$) is equal to the gas mass accretion rate due to the viscous diffusion ($2\pi R\Sigma_{\rm g}V_{\rm vis}^{2D}$), in the case of $St_{2D}\ll 1$.

2.5. Surface density distributions in gas and dust

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\Sigma_{\rm g}V_{\rm R}=\rm{constant}$ (Pringle, 1981). The gas surface density is given by $\dot{M}_{\rm g}/(2\pi RV_{\rm R}^{2D})$. For a distribution of grains, $\Sigma_{\rm d}=\dot{M}_{\rm d}/(2\pi Rv_{\rm R}^{2D})$ (e.g., Testi et al., 2014). When $\Sigma_{\rm d}/\Sigma_{\rm g}\ll 1$ and $St_{2D}\ll 1$, the radial velocity of dust grains is $v_{\rm R}^{2D}=-2St_{2D}(\Sigma_{\rm d}/\Sigma_{\rm g})\eta_{2D}V_{K}$. For simplicity, we assume single-size dust grains. With $\Sigma_{\rm g}\propto R^{-s}$, we obtain $\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}^{2D}$, the first term of equation (3) is approximately $2St_{2D}(\Sigma_{\rm d}/\Sigma_{\rm g})\eta_{2D}V_{K}$ and the second term is $V_{\rm vis}^{2D}$. If the feedback is negligible, $s=2f+1/2$. The index of the power law of $\Sigma_{\rm g}$ does not change if the feedback is effective, because the first term in $V_{\rm R}^{2D}$ has the same dependence on the radius as $V_{\rm vis}^{2D}$. Hence, in steady state, $\Sigma_{\rm d}\propto R^{-4f-1}$ and $\Sigma_{\rm d}/\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 $\Sigma_{\rm d}/\Sigma_{\rm g}$ is not very large, the gas flows inward in the disk. As discussed above, the slope of $\Sigma_{\rm g}$ is the same as that in the case that the feedback is negligible. Hence the distribution of $\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 due to the feedback, but $\Sigma_{\rm g}$ is not changed. On the other hand, when the viscosity is low or $\Sigma_{\rm d}/\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}<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. Evolution of gas disk

3.1. Model description

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 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 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 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\ \rm{AU}$ to $100\ \rm{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 $\Sigma_{\rm g}=\Sigma_{0}(R/1\ \rm{AU})^{-1}$ with $\Sigma_{0}=570\rm{g/cm}^{2}$. We adopt a simple locally isothermal equation of state and assume the disk aspect ratio as $h_{\rm g}/R=H_{0}(R/1\ \rm{AU})^{1/4}$ with $H_{0}=0.028$. The initial angular velocity of gas is given by $\Omega_{K}\sqrt{1-\eta}$, and the mass of the central star is assumed by $1M_{\odot}$. The radial velocity of gas is set by $V_{\rm vis}^{2D}$. The initial surface density of dust grains is 0.01 of gas everywhere over the disk. The initial angular velocity of dust grains is $\Omega_{K}$ and the initial radial velocity is zero. We adopt dust grains with 3 cm in size, corresponding to the Stokes number $St_{2D}$ of 0.1 at $R=10\ \rm{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 on the value of $\eta_{2D}$ and $\Sigma_{\rm d}/\Sigma_{\rm g}$ especially in the case of small viscosity, numerical instability occurs when small discontinuities of $\eta$ and $\Sigma_{\rm d}/\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\ \rm{AU}<R<6\ \rm{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\ \rm{AU}$ and $R_{1}=6\ \rm{AU}$, respectively. Hence, at the innermost annulus ($R=4\ \rm{AU}$), the gas velocity is set to $V_{\rm vis}^{2D}$. Moreover, in this region, we force all the physical quantities to be azimuthally symmetric by overwriting the quantities with their azimuthally average at every time step (c.f., de Val-Borro et al., 2006). That is, in the coupling damping zone, $\Sigma_{\rm g}$, $V_{\rm R}^{2D}$ and $V_{\phi}^{2D}$ are related towards their azimuthally averaged values as

[TABLE]

where $X$ represents $\Sigma_{\rm g}$, $V_{\rm R}^{2D}$ and $V_{\phi}^{2D}$, and $X_{\rm avg}$ denotes the azimuthally averaged values of them, and $\tau_{\rm 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\ \rm{AU}$), and unity at the opposite edge ($R=4\ \rm{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}=\rm{constant}$ (c.f., Zhu et al., 2011). At the outer boundary, the velocities of gas and dust are given by equation (37)-(40), and the surface densities are fixed at the initial values. Here the radial distribution of $\Sigma_{\rm d}/\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

For the purpose of confirming the validity and usefulness of our analytic formulas, we also calculated the evolution by the simplified 1D model with our analytic formulas and compared the results to these of the 2D fluid simulations. In the 1D model, we simultaneously solved the 1D continuity equation for the dust grains $\partial\Sigma_{\rm d}/\partial t=-(1/R)\partial(R\Sigma_{\rm d}v_{\rm R}^{2D})/\partial R$ and that for the disk gas $\partial\Sigma_{\rm g}/\partial t=-(1/R)\partial(R\Sigma_{\rm g}V_{\rm R}^{2D})/\partial R$. Here, we used equations (1) and (3) as the radial velocities of the dust grains and the disk gas ($v_{\rm R}^{2D}$ and $V_{\rm R}^{2D}$), respectively. The initial and boundary conditions are the same as the 2D simulations described above. For the coupling damping region, we just cut off the dust feedback to gas in $R<6{\rm AU}$.

3.2. Results

Figure 5 illustrates evolutions of the radial velocity and surface densities of gas and dust grains, and 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 is not so concentrated, the gas flows to the central star in the entire region of the disk. As discussed in subsection 2.5, the distribution of $\Sigma_{\rm g}$ is hardly changed from the initial distribution which is the steady state without the feedback. The $\Sigma_{\rm d}/\Sigma_{\rm g}$ is also distributed as $1/R$, as expected in subsection 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.

We now show results with different $\alpha$ and $\Sigma_{\rm d}/\Sigma_{\rm g}$. The agreement between the 1D model and 2D simulations are always good in the cases presented below. Figure 6 shows the evolution of the disk with $\alpha=10^{-2}$ and initial $\Sigma_{\rm d}/\Sigma_{\rm g}=0.05$. Because the initial $\Sigma_{\rm d}/\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); 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}\ \rm{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 over unity, and the Stokes number of grains is $\sim 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 mid-plane at $10\ \rm{AU}$ as $50$ (at $9.6\times 10^{4}\ \rm{yr}$) and $150$ (at $1.6\times 10^{5}\ \rm{yr}$). Note that the decreases of $\Sigma_{\rm d}$ near the inner boundary ($R=6\ \rm{AU}$) is originated from the difference of $\eta_{2D}$ in the inside and outside of the coupling damping zone, since $\eta_{2D}$ in the damping zone is fixed to the value in the steady state without the feedback.

The evolution with a lower viscosity ($\alpha=10^{-3}$), is shown in Figure 7. The initial $\Sigma_{\rm d}/\Sigma_{\rm g}$ is the same as 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}\ \rm{yr}$, the dust-to-gas surface density ratio increases up to unity. The dust-to-gas density ratios at the mid-plane at $10\ \rm{AU}$ are estimated as $1$ (at $9.6\times 10^{4}\ \rm{yr}$) and $10$ (at $1.6\times 10^{5}\ \rm{yr}$).

4. Discussion

4.1. Vertical structure of gas flow

The gas flows inward when the net radial velocity of the gas $\left<V_{\rm R}\right>$ is positive. As shown in Figure 1, the mass flux density depends on the altitude. The dust grains are settled to the mid-plane and the dust layer in which the dust grains are highly concentrated is formed near the mid-plane, while the dust density is quite small above the dust layer. When the net mass flux is positive, the gas in the upper 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 of the dust layer, in terms of the Stokes number of the dust grains, when $\Sigma_{\rm d}/\Sigma_{\rm g}=0.01$ at $10\ \rm{AU}$ ($\rho_{\rm g}=3.0\times 10^{-12}\ \rm{g/cm}^{2}$, $h_{\rm g}/R=0.05$). For a relatively large dust grains (e.g., $St_{\rm mid}\sim 0.1$), the net radial velocity of the gas is positive. However, above the dust layer ($z>2h_{\rm d}$), the gas flows inward. When $\alpha=10^{-3}$ and $St_{\rm mid}=1$, for instance, the amount of the gas mass flux above the dust layer is $\sim-10^{-11}\ M_{\odot}/\ \rm{AU}/\ \rm{yr}$, which is comparable with the mass flux without the dust feedback($=\Sigma_{\rm g}V_{\rm vis}^{2D}/2$). This indicates that even when the disk gas deplete from the inside of the disk such as 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.

4.2. Dust growth and fragmentation

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 which 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_{\rm frag}$ are given by

[TABLE]

where $v_{\rm 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\rm{m/s}$ – $10\rm{m/s}$ (Beitz et al., 2011). When $\alpha=10^{-3}$ and $c_{s}=10^{3}\rm{m/s}$, the Stokes number of the maximum dust grains with $v_{\rm frag}=10\rm{m/s}$ is about $\sim 0.1$. In this case, the dust feedback can influence on the entire evolution of the gas disk as shown in Figures 6 and 7, which leads formation of rocky planetesimals.

If $v_{\rm frag}\ll 10\rm{m/s}$,the dust feedback does not work, since the dust grains cannot grow to a sufficiently large size. If $v_{\rm frag}\gg 10\rm{m/s}$, 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$.

4.3. Implication of planet formation

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 Figure 6 and 7, when the dust feedback is effective, after $10^{5}\ \rm{yr}$, $\Sigma_{\rm d}/\Sigma_{\rm g}$ increases over unity within $10\ \rm{AU}$. Using equation (36), we can estimate $\rho_{\rm d}/\rho_{\rm g}$ at the mid-plane as $>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 (metalicity) and etc. The streaming instability can be caused in more early state of the disk evolution than $10^{5}\ \rm{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 decreases. 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 be still 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 decreases. 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.

4.4. Implication of disk observation

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 decreases 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 mid-plane which are observed by sub-millimeter may decreases. In this case, a hole structure can be found by the observations of dust continuum by sub-millimeter, 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\ \rm{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.

4.5. Validity of the model

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 motions (equations 10 – 12) (see, Appendix A). Although this formulation is widely used in the previous studies (e.g., Fu et al., 2014; Gonzalez et al., 2017; Dipierro & Laibe, 2017), this treatment of the dust fluid may not be appropriate in the cases where the dust density is comparable with 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>1/2$. At the mid-plane where the dust density is comparable with 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, b). We must take into account the momentum transfer due to the dust turbulent in this case, and 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 towards the mid-plane, so the gas density in the mid-plane is larger than that given by equation (24). As the dust-to-gas density ratio increases, the settling velocity of the dust grains is slow down ($v_{z}=-St^{\prime}\Omega_{\rm K}z\propto 1/(1+\rho_{\rm d}/\rho_{\rm g})$), which leads 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 $\Sigma_{\rm d}/\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 (mid-plane), 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 mid-plane 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 dust-to-gas mass ratio is less than $\sim 1$, and we expect that qualitative results are valid up to $\Sigma_{\rm d}/\Sigma_{\rm g}\sim 1$. If the dust grains are highly concentrated and dominate gas, it may be necessary to have more sophisticated formulation to describe the dynamics of gas and dust. If 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.

5. Summary

In this paper, we considered the effect of the dust feedback on the viscous evolution of the gas disk. Our results are summarized as follows:

We present the analytical expressions for the velocities of gas and dust grains in the 2D disk, considering viscous evolution of gas and gas–dust friction (equations 1 – 4). Considering the vertical dust settling, we also derived the formula of the radial velocities of the gas and dust grains (equation 43 with equation 39). For the net radial velocities of the gas and the dust, the analytical expression of the 2D disk reasonably agree with the formula considering the dust settling (see Figure 2). 2. 2.

We found that the feedback from the grains significantly affects the radial velocity of the gas (see, Figures 2 and 3). When the gas viscosity is sufficiently large and the dust grains are not concentrated, the gas infall slows down due to the dust feedback. As the viscosity decreases or the dust-to-gas mass ratio increases, the gas flows even outward due to the dust feedback. 3. 3.

We also demonstrated the 2D two-fluid hydrodynamic simulations, and showed how the feedback changes the evolutions of gas. As long as the viscosity is large and the initial dust-to-gas mass ratio is small, the gas flows inward. In this case, the gas disk evolves as in the case where the dust feedback is not effective, though the infall velocity of the gas decreases (see, Figure 5). When the viscosity is small or the initial $\Sigma_{\rm d}/\Sigma_{\rm g}$ is large, the feedback drastically changes the evolution of the disk. The gas flows to the outside of the disk, and the gas at the inner region is significantly depleted (see, Figures 6 and 7). The gas removal slows down the infall velocity of the dust grains, and the dust-to-gas mass ratio further increases. 4. 4.

We presented the idea of a simplified 1D model. We solve the continuity equations with the velocities given by equations (1) and (3). Since the velocities quickly converge to those in the steady state, the simplified 1D model is able to well reproduce the results of the 2D hydrodynamic simulations, as can be seen in Figures 5 – 7. 5. 5.

We also discuss the effect of the vertical structure on the disk evolution. If the disk gas is deposited, the gas accretion onto the star is possible. We need to observe the gas directly to detect the gas depletion. 6. 6.

In the inner region of the disk, the size of silicate dust grains would be strongly limited by the fragmentation. This situation is not very different from the cases assumed in our simulations. If the silicate grains remove the gas via the feedback, rocky planetesimals would be formed via the streaming instability.

In this paper, we pointed out that the dust feedback can significantly influence the evolution of gas disk. However, we must consider the size evolution of dust grains, because the feedback strongly depends on the grain size. In order to understand the evolution of gas disk and grains, we should take into account consistently the evolution of grain size and the evolution of gas, with the dust feedback.

We would like to thank the anonymous referee for his/her helpful comments. This work was supported by the Polish National Science Centre MAESTRO grant DEC- 2012/06/A/ST9/00276. Numerical computations were carried out on the Cray XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. S. O. is supported by Grants-in-Aid for Scientific Research (#15H02065, 16K17661, 16H04081) from MEXT of Japan. T. M. is supported by Grants-in-Aid for Scientific Research (#26800106, 15H02074, 17H01103) from MEXT of Japan.

Appendix A Validity of basic equations of the gas and dust grains

In order to describe flows including turbulence, the Reynolds-averaged Navier-Stokes equations of motion are commonly used. In this way, the density and all velocity components are devided into mean and fluctuating parts as $\rho=\overline{\rho}+\rho^{\prime}$ and $\mathbf{v}=\overline{\mathbf{v}}+\mathbf{v}^{\prime}$ (the former term is the mean part and the later term is the fluctuating part). The Reynolds-averaged Navier-Stokes equations of the gas and the dust grains in protoplanetary disks are derived by Cuzzi et al. (1993). Here we derive our basic equations (10)–(15) from the Reynolds-averaged equations (equations A.11–A.13 of Cuzzi et al. (1993)). As general forms of the Reynolds-averaged equations, Cuzzi et al. (1993) obtained

[TABLE]

where $P$ and $\Psi$ are a pressure and a gravity potential, respectively, and $\sigma_{ij}$ is the Reynolds stress tensor defined by

[TABLE]

and $F_{i}$ is the gas–dust friction force given by

[TABLE]

where the sign of equation (A3) is positive for the gas, while it is negative for the dust grains. As discussed later, if the fluctuating motion of the dust grains is dominated by the gas–dust fraction force, we may neglect the second term of equation (A3) (see also, Appendix B of Cuzzi et al. (1993)).

First let us consider the equations of the motion for the gas. If the mach number of the turbulence is not so very high, we can treat the gas fluid as weakly compressible fluid. In this case, since $\overline{\rho}_{g}\gg\rho^{\prime}_{g}$, we can neglect the time variation of $\overline{\rho^{\prime}_{g}v_{g,i}^{\prime}}$ (2nd term of LHS of equation A1) and terms related with the advection due to the turbulence (4th and 5th terms of RHS of equation A1). Using the Newtonian viscous stress tensor, instead of the Reynolds stress tensor, we obtain the equations of motions for the gas as equations (13) – (15).

Since we treat the dust fluid as pressureless fluid, $P_{d}=0$ in equation (A1) for the dust grains. If adopting the gradient diffusion hypothesis (Cuzzi et al., 1993; Takeuchi & Lin, 2002), the mass flux due to the turbulence can be written by

[TABLE]

where $D$ is the diffusion coefficient which may be given by

[TABLE]

where $\nu$ is the kinetic viscosity of the gas. In this case, the advection terms due to the turbulence (4th and 5th terms of RHS of equation A1) are proportional to $\nu\overline{\rho_{d}}/Sc$. As long as the distributions of $\overline{\rho_{d}}/\overline{\rho_{g}}$ and $\overline{v}$ are not steep, those terms are much smaller than $F_{i}$ given by equation (A3), because $\overline{\rho_{d}}\left(\overline{v_{g,i}}-\overline{v_{d,i}}\right)\sim\eta\overline{\rho_{d}}V_{K}$. We can neglect 4th and 5th terms of RHS of equation (A1).

The Reynolds stress of the dust grains may not be negligible. Here let us consider how the fluctuating motion of the dust grains responses to the fluctuating motion of the gas with the velocity of $v_{0}$. Assuming the decay time of the gas fluctuating motion is longer than the time scale that we now consider, we regard that $v_{0}$ is independent of time. If the fluctuation motion of the dust grain is dominated by the gas–dust friction, we may write a governing equation of the fluctuation motion of the dust grain as

[TABLE]

In this case, the fluctuating velocity of the dust grains is given by

[TABLE]

The difference between the fluctuating velocities of the gas and the dust grains decays as $\exp(-t/t_{\rm stop})$. For small grains, especially, the velocity of the dust fluctuating motion converges quickly to that of the gas fluctuating motion. In this case, the correlation $\overline{v_{d,i}^{\prime}v_{d,j}^{\prime}}$ would be approximately given by $\overline{v_{g,i}^{\prime}v_{g,j}^{\prime}}$ in equation (A2). Similarly, the correlation $\overline{\rho_{d}^{\prime}v_{d,i}^{\prime}}$ is also approximated as $\overline{\rho_{d}^{\prime}v_{g,i}^{\prime}}$ and they cancel out each other in the second term of equation (A3). When $\overline{\rho_{d}}\ll\overline{\rho_{g}}$, the components of the Reynolds stress tensor for the dust grains are much smaller than these for the gas becuase the dust density is small. However, when $\overline{\rho_{d}}\sim\overline{\rho_{g}}$, the components of the Reynolds stress tensor for the dust grains are comparable with these for the gas. In this case, we may need to consider the kinetic viscosity of the dust grains in the equation of the motion. However, since there is a large uncertainty about treatment of the dust fluctuating motion, we drop the term related with the Reynolds stress of the dust grains for simplicity. Ignoring a time-variation of the turbulence, and the Reynolds stress term and the advection terms due to the turbulence, we derive the equations of motions for the dust grains (10)–(12) from the Reynolds-averaged Navier-Stokes equation (A1).

Appendix B Vertical structures of gas and dust grains

The vertical structures of the gas and the dust grains are described by the Gaussian distributions as equations (24) and (33). However, when $\rho_{\rm d}\gtrsim\rho_{\rm g}$, the vertical distributions deviate from the Gaussian distributions.

Nakagawa et al. (1986) considered the vertical structure of the gas by adopting a simple dust vertical distribution as $\rho_{\rm d}(z)=const.=\rho_{\rm d}{}_{,0}$ for $|z|<h_{\rm d}$, otherwise $\rho_{\rm d}(z)=0$ (see Appendix A of that paper). The gas density at the mid-plane is given by $\rho_{\rm g}(0)\left(1+f\Sigma_{\rm d}/\Sigma_{\rm g}\right)$, where $\rho_{\rm g}(0)$ is the gas density at $z=0$ given by equation (24), and $f$ is a function of $h_{\rm d}/h_{\rm g}$ which is an order of unity (e.g., $f=0.5$ when $h_{\rm d}/h_{\rm g}=1$ and $f=1$ when $h_{\rm d}\ll h_{\rm g}$). If $\Sigma_{\rm d}/\Sigma_{\rm g}\ll 1$, hence, the vertical distribution of gas is given by the Gaussian distribution (equation 24). If $\Sigma_{\rm d}/\Sigma_{\rm g}\gg 1$, however, the gas density at the mid-plane is larger than that given by equation (24), and the vertical distribution also deviates from the Gaussian distribution.

Here we consider self-consistent vertical structures of the gas and the dust grains. The vertical gradient of the gas density is given by (see, section 2.3):

[TABLE]

From equation (31) with $F_{m,z}=0$, the vertical gradient of the dust density is given by

[TABLE]

Solving equations (B1) and (B2), we obtain the vertical structures of the gas and the dust grains. In Figure 9, we show the vertical structures of the gas and the dust grains if $\alpha=10^{-3}$ and $St_{\rm mid}=0.1$. For comparison, we plot the Gaussian distributions given by equations (24) and (33). When $\Sigma_{\rm d}/\Sigma_{\rm g}=1$, the gas distribution agrees with equation (24), whereas the gas density near the mid-plane is slightly enhanced. When $\Sigma_{\rm d}/\Sigma_{\rm g}=5$, the gas density near the mid-plane increases, while the gas density apart from the mid-plane slightly decreases. The thickness of the dust layer becomes thicker than that expected by equation (33), because the settling velocity is slow down by the factor of $\rho_{\rm g}+\rho_{\rm d}$ in $St^{\prime}$. The dust density at the mid-plane decreases, due to vertical swelling of the dust layer. The dust-to-gas density ratio is smaller than that expected by the Gaussian distributions of equations (24) and (33), while it is larger at the upper part of the dust layer.