Revised description of dust diffusion and a new instability creating multiple rings in protoplanetary disks

TL;DR

This paper introduces a revised dust diffusion model conserving angular momentum, analyzes the secular gravitational instability, and discovers a new instability (TVGI) that can form multiple rings in protoplanetary disks, relevant to planetesimal formation.

Contribution

It formulates angular momentum conserving dust diffusion equations, identifies a new two-component viscous gravitational instability, and demonstrates its potential to form multiple rings in disks.

Findings

Secular gravitational instability is a monotonically growing mode.

The new TVGI efficiently accumulates dust and forms rings.

Both instabilities can produce rings about 10 au wide beyond 50 au in HL Tau.

Abstract

Various instabilities have been proposed as a promising mechanism to accumulate dust. Moreover, some of them are expected to lead to the multiple-ring structure formation and the planetesimal formation in protoplanetary disks. In a turbulent gaseous disk, the growth of the instabilities and the dust accumulation are quenched by turbulent diffusion of dust grains. The diffusion process has been often modeled by a diffusion term in the continuity equation for the dust density. The dust diffusion model, however, does not guarantee the angular momentum conservation in a disk. In this study, we first formulate equations that describe the dust diffusion and also conserve the total angular momentum of a disk. Second, we perform the linear perturbation analysis on the secular gravitational instability (GI) using the equations. The results show that the secular GI is a monotonically growing mode, contrary to the result of previous analyses that found it overstable. We find that the overstability is caused by the non-conservation of the angular momentum. Third, we find a new axisymmetric instability due to the combination of the dust-gas friction and the turbulent gas viscosity, which we refer to as two-component viscous gravitational instability (TVGI). The most unstable wavelength of TVGI is comparable to or smaller than the gas scale height. TVGI accumulates dust grains efficiently, which indicates that TVGI is a promising mechanism for the formation of multiple-ring-like structures and planetesimals. Finally, we examine the validity of the ring formation via the secular GI and TVGI in the HL Tau disk and find both instabilities can create multiple rings whose width is about 10 au at orbital radii larger than 50 au.

Equations

$$\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial t}+\frac{1}{r}\frac{\partial \left(r{\Sigma }_{\mathrm{d}}{v}_r\right)}{\partial r}=\frac{1}{r}\frac{\partial }{\partial r}\left(rD\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right),$$

$${\Sigma }_{\mathrm{d}}\left[\frac{\partial v_{\varphi }}{\partial t}+v_r\frac{\partial v_{\varphi }}{\partial r}\right]=-{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }{v}_r}{r}-{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}},$$

$$\frac{\partial \left({\Sigma }_{\mathrm{d}}{j}_{\mathrm{d}}\right)}{\partial t}$$

$$=-r{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}}+j_{\mathrm{d}}\frac{1}{r}\frac{\partial }{\partial r}\left(rD\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right).$$

$$\frac{\partial \Sigma }{\partial t}+\frac{1}{r}\frac{\partial \left(r\Sigma u_r\right)}{\partial r}=0,$$

$$\Sigma \left[\frac{\partial u_{\varphi }}{\partial t}+u_r\frac{\partial u_{\varphi }}{\partial r}\right]=-\Sigma \frac{u_{\varphi }{u}_r}{r}+{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}},$$

$$\frac{\partial \left(\Sigma j_{\mathrm{g}}\right)}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}_r\Sigma j_{\mathrm{g}}\right)=r{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}}.$$

$$\frac{\partial }{\partial t}\left(\Sigma j_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}{j}_{\mathrm{d}}\right)$$

$$=j_{\mathrm{d}}\frac{1}{r}\frac{\partial }{\partial r}\left(rD\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right).$$

$$\frac{\partial \left({\Sigma }_{\mathrm{d}}{j}_{\mathrm{d}}\right)}{\partial t}$$

$$=-r{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}}-D\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\frac{\partial j_{\mathrm{d}}}{\partial r}.$$

$$\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r\left(v_r-\frac{D}{{\Sigma }_{\mathrm{d}}}\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right){\Sigma }_{\mathrm{d}}\right]=0,$$

$${\Sigma }_{\mathrm{d}}\left[\frac{\partial j_{\mathrm{d}}}{\partial t}+\left(v_r-\frac{D}{{\Sigma }_{\mathrm{d}}}\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right)\frac{\partial j_{\mathrm{d}}}{\partial r}\right]=-r{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}}.$$

$$\frac{\partial \left({\Sigma }_{\mathrm{d}}{j}_{\mathrm{d}}\right)}{\partial t}$$

$$=-r{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}}.$$

$$\frac{\partial \Sigma }{\partial t}+\frac{1}{r}\frac{\partial \left(r\Sigma u_r\right)}{\partial r}=0,$$

$$\Sigma \left(\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}\right)$$

$$+\frac{\partial }{\partial x_j}\left[\Sigma \nu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_k}{\partial x_k}\right)\right]$$

$$+{\Sigma }_{\mathrm{d}}\frac{v_i-u_i}{t_{\mathrm{stop}}},$$

$${\nabla }^2\Phi =4\pi G\left(\Sigma +{\Sigma }_{\mathrm{d}}\right)\delta (z),$$

$$\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial t}+\frac{1}{r}\frac{\partial \left(r{\Sigma }_{\mathrm{d}}{v}_r\right)}{\partial r}=\frac{1}{r}\frac{\partial }{\partial r}\left(rD\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right),$$

$$\displaystyle\Sigma_{\mathrm{d}}\biggl{[}\frac{\partial v_{r}}{\partial t}+$$

$$=$$

$$-{\Sigma }_{\mathrm{d}}\frac{v_r-u_r}{t_{\mathrm{stop}}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{v}_{r}D\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right),$$

$$\displaystyle\Sigma_{\mathrm{d}}\biggl{[}\frac{\partial v_{\phi}}{\partial t}$$

$$=-{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }}{r}\left(v_r-\frac{D}{{\Sigma }_{\mathrm{d}}}\frac{\partial {\Sigma }_{\mathrm{d}}}{\partial r}\right)-{\Sigma }_{\mathrm{d}}\frac{v_{\varphi }-u_{\varphi }}{t_{\mathrm{stop}}},$$

$$\frac{\partial }{\partial t}\left(\Sigma j_{\mathrm{g}}+{\Sigma }_{\mathrm{d}}{j}_{\mathrm{d}}\right)$$

$$=\frac{1}{r}\frac{\partial }{\partial r}\left(r^3\Sigma \nu \frac{\partial \Omega }{\partial r}\right),$$

$$n\delta \Sigma +ik{\Sigma }_0\delta u_x=0,$$

$$n\delta u_x=2\Omega \delta u_y-\frac{c_{\mathrm{s}}^2}{{\Sigma }_0}ik\delta \Sigma -ik\delta \Phi -\frac{4}{3}\nu k^2\delta u_x+ϵ\frac{\delta v_x-\delta u_x}{t_{\mathrm{stop}}},$$

$$n\delta u_y=-\frac{\Omega }{2}\delta u_x-\nu k^2\delta u_y-ik\frac{3\nu \Omega }{2{\Sigma }_0}\delta \Sigma +ϵ\frac{\delta v_y-\delta u_y}{t_{\mathrm{stop}}},$$

$$n\delta {\Sigma }_{\mathrm{d}}+ik{\Sigma }_{\mathrm{d},0}\delta v_x=-D{k}^2\delta {\Sigma }_{\mathrm{d}},$$

$$n\delta v_x=2\Omega \delta v_y-\frac{c_{\mathrm{d}}^2}{{\Sigma }_{\mathrm{d},0}}ik\delta {\Sigma }_{\mathrm{d}}-ik\delta \Phi -\frac{\delta v_x-\delta u_x}{t_{\mathrm{stop}}},$$

$$n\delta v_y=-\frac{\Omega }{2}\left(\delta v_x-\frac{ikD}{{\Sigma }_{\mathrm{d},0}}\delta {\Sigma }_{\mathrm{d}}\right)-\frac{\delta v_y-\delta u_y}{t_{\mathrm{stop}}},$$

$$\delta \Phi =-\frac{2\pi G}{k}\left(\delta \Sigma +\delta {\Sigma }_{\mathrm{d}}\right),$$

Full text

diffusion — hydrodynamics — instabilities — turbulence — protoplanetary disks

††journal: ApJ

1 Introduction

Planets are thought to form through a process in which dust grains accumulate and undergo collisional growth in protoplanetary disks. Planetesimals are 10 km sized objects that form during the growth from dust grains to planets. Understanding the planetesimal formation is important in order to reveal the planet formation. Various instabilities are proposed as a possible mechanism to form planetesimals. The streaming instability is one example (e.g., Youdin & Goodman, 2005; Youdin & Johansen, 2007; Johansen & Youdin, 2007). The streaming instability is driven by friction between gas and drifting dust, resulting in the dust-clump formation. Planetesimals are expected to form if the resulting dust clumps self-gravitationally collapse (e.g., Johansen et al., 2007). Secular gravitational instability (GI) is another possible mechanism of the planetesimal formation (e.g., Ward, 2000; Youdin, 2011; Michikoshi et al., 2012; Takahashi & Inutsuka, 2014). The secular GI grows as a result of the decrease of the Coriolis force on the dust due to the dust-gas friction in a self-gravitationally stable disk. Since the growth of the secular GI also results in accumulation of dust grains, this instability has been proposed as the formation mechanism of planetesimals in outer disks and debris disks (Takahashi & Inutsuka, 2014; Tominaga et al., 2018). Moreover, the secular GI is one of the possible mechanisms to create multiple rings that recent observations with Atacama Large Millimeter/submillimeter Array (ALMA) have found in some protoplanetary disks (e.g., ALMA Partnership et al., 2015; Andrews et al., 2016; Tsukagoshi et al., 2016; Isella et al., 2016; Fedele et al., 2017, 2018). Based on the linear analysis, Takahashi & Inutsuka (2016) showed that the multiple rings in the HL Tau disk can form via the growth of the secular GI.

Protoplanetary disks are thought to be turbulent, resulting in the angular momentum transport in gas disks. Various hydrodynamic and magnetohydrodynamic instabilities have been proposed as mechanisms for driving the turbulence. These include the convective instabilitiy (e.g., Lin & Papaloizou, 1980; Klahr & Hubbard, 2014; Lyra, 2014), the vertical shear instability (e.g., Urpin & Brandenburg, 1998; Urpin, 2003; Nelson et al., 2013) and the magnetorotational instability (e.g., Balbus & Hawley, 1991, 1998). It is known that dust grains diffuse in such a turbulent gas disk, and the diffusion coefficient is evaluated based on numerical simulations and analyses using the Langevin equations (e.g., Carballido et al., 2005; Johansen & Klahr, 2005; Youdin & Lithwick, 2007). The dust diffusion is the most crucial process against the accumulation of dust grains. In previous work on the secular GI, the dust diffusion is modeled by introducing a diffusion term in the continuity equation for the dust, which is used for the linear analysis (e.g., Youdin, 2011; Shariff & Cuzzi, 2011; Takahashi & Inutsuka, 2014; Shadmehri, 2016; Latter & Rosca, 2017). Although this modeling has been widely used (e.g., Cuzzi et al., 1993; Goodman & Pindor, 2000; Dullemond & Penzlin, 2018), the modeling has a problem that the total angular momentum of the dusty gas disk is not conserved even when we consider the back-reaction from the dust to the gas and solve equations for both compnents. As described in Section 2, the non-conservation of the angular momentum affects mainly on the motion of the dust. Hence, we need to resolve this problem to precisely discuss the accumulation process of dust grains. Goodman & Pindor (2000) also pointed out this problem. They artificially modified the momentum equations for the dust and avoided the problem arising from the diffusion term. This modification, however, is not validated. In this work, we first formulate phenomenological equations that guarantee the conservation of the angular momentum. Our formulation is based on the Reynold-averaging, which was also used in Cuzzi et al. (1993) and Shariff & Cuzzi (2011). Next, we perform the linear analysis and investigate the stability of dusty-gas disks by using the formulated equations,.

This paper is organized as follows. In Section 2, we show the problem mentioned above and describe how to resolve this problem. Basic equations are presented in Section 3. We show results of the linear analysis for cases with and without the turbulent viscosity of gas. We find a new instability driven by the combination of the friction and the turbulent viscosity besides the secular GI. In Section 5, we discuss effect of the disk thickness and the multiple-ring formation via those instabilities.

2 Angular momentum transport due to turbulent diffusion

First of all, we show that equations used in the previous work (e.g., Takahashi & Inutsuka, 2014; Shadmehri, 2016; Latter & Rosca, 2017; Dullemond & Penzlin, 2018) do not conserve the total angular momentum of a dusty gas disk. We assume that a disk is infinitesimally thin and axisymmetric and solve the evolutionary equations of the vertically integrated disk. We adopt a cylindrical coordinate system $(r,\phi)$. The continuity equation and the azimuthal equation of motion for the dust are as follows:

[TABLE]

[TABLE]

where $\Sigma_{\mathrm{d}}$ is the surface density of the dust, $v_{r}$ and $v_{\phi}$ are the radial and azimuthal velocities of the dust, respectively. The diffusion coefficient is denoted by $D$, and the azimuthal velocity of the gas is $u_{\phi}$. The stopping time of a dust grain is $t_{\mathrm{stop}}$. The term on the right hand side of Equation (1) is often introduced to model the dust diffusion in a turbulent disk (e.g., Youdin, 2011). From these equations, we obtain an equation governing the evolution of the angular momentum of the dust $\Sigma_{\mathrm{d}}j_{\mathrm{d}}\equiv\Sigma_{\mathrm{d}}rv_{\phi}$:

[TABLE]

The first term on the right hand side stands for the angular momentum transport due to the friction between the dust and the gas. If we take into account the back-reaction from the dust to the gas and the gas motion, this friction term does not violate the conservation of the total angular momentum of a dusty gas disk, as shown in the following. The continuity equation and the azimuthal equation of motion for gas are the following:

[TABLE]

[TABLE]

where $\Sigma$ and $u_{r}$ are the surface density and the radial velocity of the gas, respectively. We obtain an equation describing the time evolution of the angular momentum of the gas $\Sigma j_{\mathrm{g}}\equiv\Sigma ru_{\phi}$ as follows:

[TABLE]

Equations (3) and (6) give the following equation:

[TABLE]

The term on the right hand side cannot be written in a form of the divergence of the angular momentum flux because it is proportional to the specific angular momentum. Thus, the volume integral of this term is not zero in general, meaning that the total angular momentum of the disk is not conserved. This is due to the treatment of the dust diffusion. The diffusion term introduced in the continuity equation for the dust (Equation (1)) unphysically changes the angular momentum of the dust disk and prevents the total angular momentum from being conserved. The gas motion is also affected by the unphysical change of the angular momentum of the dust because it is partly transported to the gas through the dust-gas friction. The effect on the gas motion is, however, smaller than that on the dust motion by a factor of the dust-to-gas mass ratio because the rate of the angular momentum transport due to the friction is proportional to the surface density of the dust (Equation (6)). This implies that the non-conservation of the angular momentum affects mainly on the motion of the dust.

In order to discuss how the non-conservation of the angular momentum of a dusty gas disk affects the motion of the dust, we rearrange Equation (3) as follows:

[TABLE]

We do not consider the effect of the second term on the left hand side since this term becomes the surface term and vanishes when it is integrated over all space. The second term on the right hand side changes the total angular momentum of the disk. If we consider that the disk rotates with the Keplerian velocity, that term represents negative (positive) torque exerted on the dust when the surface density gradient of the dust is positive (negative). In the dust accumulating region, for example, the dust in the inner part ($\partial\Sigma_{\mathrm{d}}/\partial r>0$) loses the angular momentum and goes inward, and vice versa. This unphysical torque prevents the accumulation of the dust. In other words, the previous studies underestimated the degree of dust accumulation. Hence, in order to discuss the accumulation process precisely, we need to re-formulate equations that conserve the angular momentum.

If the dust grains are small, and the stopping time is short compared to the Keplerian period, the diffusion is mainly driven by radial kicks from the turbulent gas (Youdin & Lithwick, 2007), meaning that the specific angular momentum of the dust grains does not change along the diffusion flow. Equation (1) is rearranged to yield

[TABLE]

which means that the advection velocity caused by the diffusion is written by $-D\Sigma_{\mathrm{d}}^{-1}\partial\Sigma_{\mathrm{d}}/\partial r$. Then, we consider a case that the time evolution of the specific angular momentum is governed by the following equation:

[TABLE]

(see Appendix A for the detailed derivation). From this equation, we obtain the following equation for the angular momentum of the dust:

[TABLE]

This equation shows that the angular momentum of the dust disk changes only through the friction. From Equations (6) and (11), we can show that the total angular momentum of the dusty-gas disk is conserved. In this way, we see that the equations conserve the total angular momentum if the advection velocity by the diffusion is considered in the equation for the specific angular momentum. Such an advection term is naturally derived by the Reynolds averaging of hydrodynamic equations for the dust. We also formulate a radial momentum equation consistent with Equations (1) and (10) based on the Reynolds averaging in Appendix A. We summarize a set of basic equations including the newly formulated dust equations in the next section.

3 Basic equations

In this work, we perform the linear stability analysis of an infinitesimally thin and axisymmetric disk with taking into account the conservation of the angular momentum in the disk as discussed in the previous section. We summarize the basic equations for gas and dust in this section. We use the following equations for the gas and the Poisson equation, which were also used in Takahashi & Inutsuka (2014, 2016):

[TABLE]

[TABLE]

[TABLE]

where $u_{i},v_{i}$ are the $i$-th component of the velocities of the gas and the dust, $c_{\mathrm{s}}$ is the sound speed, $\Phi$, $G$ and $M_{\ast}$ are the gravitational potential of the disk, the gravitational constant and the mass of the central star, respectively. We denote the coefficient of the turbulent viscosity by $\nu$ that is measured by the dimensionless parameter $\alpha\equiv\nu\Omega c_{\mathrm{s}}^{-2}$ (Shakura & Sunyaev, 1973), where $\Omega$ is the angular velocity of the gas disk.

Next, we summarize equations for the dust. We derived the following equations for the dust, which include Equation (1), based on the Reynolds averaging:

[TABLE]

[TABLE]

[TABLE]

where $c_{\mathrm{d}}$ is the velocity dispersion of the dust. The advection velocity $-D\Sigma_{\mathrm{d}}^{-1}\partial\Sigma_{\mathrm{d}}/\partial r$ that appears on the left hand side of Equation (14) and on the both sides of Equation (15) is not included in the previous studies. The last term on the right hand side of Equation (14) is also not included. See Appendix A for detailed derivation. From Equations (1), (4), (12) and (15), we can derive the evolutionary equation for the total angular momentum:

[TABLE]

The term on the right hand side shows the angular momentum transport due to the viscosity. Thus, our basic equations conserve the total angular momentum even if we take into account the dust diffusion.

We perform the linear analysis in local Cartesian coordinates $(x,y)=(r-R_{0},R_{0}(\phi-\Omega t))$ corotating with the angular velocity of the disk $\Omega$ at the radius $R_{0}$. As an unperturbed state, we assume that the surface densities of the gas and the dust are uniform, $u_{x,0}=v_{x,0}=0$, and $u_{y,0}=v_{y,0}=-(3/2)\Omega x$. Considering axisymmetric perturbations $\delta\Sigma,\delta\Sigma_{\mathrm{d}},\delta u_{x},\delta u_{y},\delta v_{x},\delta v_{y},\delta\Phi$ proportional to $\exp[nt+ikx]$, we linearize Equations (1), (4), (12)–(15). The linearized equations are given as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $n$ is the growth rate, $k$ is the wavenumber, and $\epsilon\equiv\Sigma_{\mathrm{d},0}/\Sigma_{0}$ is the dust-to-gas mass ratio. We do not consider the large-scale pressure gradient in the background state or the vertical motion of the gas and the dust in this work, which precludes the streaming instability. This setup is suitable to study growing modes different from the streaming instability.

4 Results

In this section, we present the results of the linear analyses with and without the turbulent viscosity and compare our results with those of the previous work. There are six modes in our analysis because we solve the equations for both gas and dust as described above. In the most simplified limit for the gas and the dust where we do not consider any of the dust-gas friction, the dust diffusion, or the turbulent gas viscosity, these are two density waves for each component and two static modes, and hence six modes in total (see Appendix B). A static mode is a steady solution of the perturbation equations. Figure 1 shows how the six modes change by adding three physical processes step by step and which mode becomes unstable. As shown in Figure 1, the two static modes become unstable, which is explained in more detail in the following subsections.

4.1 Without Turbulent Viscosity

First, we show results obtained when we neglect the viscosity term in Equation (12), which corresponds to Step (2) in Figure 1. We note that we relate $D$ and $c_{\mathrm{d}}$ to the dimensionless turbulent strength $\alpha$ as shown below (Equations (32) and (33)) although $\nu$ does not appear in this section. This is because it is easier to compare the results with those in Section 4.2.

We find one static mode ($n=0$) and one mode that can be unstable (the modes B’ and C in Figure 1). The static mode is the perturbed state where the dust has the same azimuthal velocity with the gas, and the radial force balance holds. The latter mode corresponds to the secular GI. By using the terminal velocity approximation ($t_{\mathrm{stop}}\ll n^{-1}$) 111The terminal velocity approximation sometimes refers to $v_{i}=u_{i}+t_{\mathrm{stop}}\Sigma^{-1}\partial\left(c_{\mathrm{s}}^{2}\Sigma\right)/\partial x_{i}$. This expression is valid only when the Coriolis force and the gas viscosity can be neglected. “The terminal velocity approximation” adopted here is more general one, which refers to neglecting the time derivative of the relative velocity by assuming it is small enough compared to the friction term. and assuming $t_{\mathrm{stop}}\ll\Omega^{-1}\ll n^{-1}$, we obtain the following approximate dispersion relation of the secular GI:

[TABLE]

where

[TABLE]

Figure 2 shows the dispersion relation of the secular GI obtained in this work and that of the previous work for $D=10^{-4}c_{\mathrm{s}}^{2}\Omega^{-1},c_{\mathrm{d}}=0,\epsilon=0.1,\tau_{\mathrm{s}}\equiv t_{\mathrm{stop}}\Omega=0.01$ and $Q\equiv c_{\mathrm{s}}\Omega/\pi G\Sigma_{0}=3$. As mentioned in Takahashi & Inutsuka (2014), the secular GI mode obtained in this work is also stabilized for the long wavelength perturbation by the Coriolis force exerted on the dust. In contrast, the short wavelength perturbation is stabilized by the turbulent diffusion. For wavelengths where the secular GI mode is unstable, the gas pressure gradient force dominates the Coriolis force acting on the gas as a stabilizing force, which forms an azimuthal zonal flow (see, Latter & Rosca, 2017), and the dust accumulates by the self-gravity of itself. The growth rate of the secular GI obtained in this analysis is several times larger than that obtained in the previous work. In addition, we find the secular GI is a monotonically growing mode, while the secular GI found in the previous work is overstable depending on parameters222The approximated dispersion relation, Equation (24), is linear for $n$, although that derived in the previous work (Equation (13) in Takahashi & Inutsuka, 2014) is a quadratic equation. This is why the overstable mode does not appear in our formulation.. This shows that the overstable mode obtained in the previous work is due to the unphysical torque acting on the dust. We derive an approximate condition for the growth of the secular GI from Equation (24). The condition is that at least two solutions for $n=0$ exist in the region $k>0$, which is equivalent to the condition that the equation $A_{0}/k^{2}=0$ has two distinct positive real solutions. We then obtain the following approximate condition for the instability:

[TABLE]

Equation (30) is equivalent to the condition derived in Latter & Rosca (2017) in the case of $D=0$ (see Equation (43) in their paper). If we assume $c_{\mathrm{d}}=0$ and $(1+\epsilon)D\ll\epsilon c_{\mathrm{s}}^{2}t_{\mathrm{stop}}$, Equation (30) becomes

[TABLE]

which is equivalent to the condition obtained in the previous work. For instance, when $(1+\epsilon)D\ll\epsilon c_{\mathrm{s}}^{2}t_{\mathrm{stop}}$ is satisfied, Equation (50) in Latter & Rosca (2017) is equivalent to the above condition (Equation (31)). In this way, we find that our formulation does not change the condition for the secular GI if $D$ is small enough to satisfy $(1+\epsilon)D\ll\epsilon c_{\mathrm{s}}^{2}t_{\mathrm{stop}}$.

We find that the mode exchange between the secular GI mode and the classical GI mode, which is an unstable density wave (see Appendix B), occurs at $k=k_{\mathrm{c},-},k_{\mathrm{c},+}$ where the growth rate of the classical GI mode becomes zero (Figure 3). The mode exchange is reconnection of the curves in the $k$ – $n$ plane of dispersion relations for two different modes. We note that the classical GI mode seen here is the dust GI mode mediated by the self-gravitationally stable gas disk. Hereafter, we simply refer to the classical GI mode as the dust GI mode in order to avoid confusion with the gas GI mode that does not appear since $Q$ is larger than unity in this paper. We can see the mode exchange even in the linear analysis only for the motion of the dust (see Figure 10 in Youdin, 2011). As shown in Figure 3, we designate the growing mode in the limited range of wavenumber where the dust GI remains unstable for $t_{\mathrm{stop}}\to\infty$ (no friction between gas and dust) as “the dust GI”. On the other hand, we designate the growing mode in the disconnected curves in the regions of wavenumber where the dust GI is stable for $t_{\mathrm{stop}}\to\infty$ as “the secular GI”. Figure 4 shows the maximum growth rate of the instabilities as a function of $\tau_{\mathrm{s}}$ and $\alpha$ for $\epsilon=0.1,Q=3$, where we use the following equations to calculate the diffusion coefficient and the velocity dispersion (Youdin & Lithwick, 2007):

[TABLE]

[TABLE]

We also show the upper limit of $\alpha$ obtained from the approximate condition for the instability (Equation (30)) in Figure 4. The exact upper limit of $\alpha$ is well represented by Equation (30). The secular GI is the fastest growing mode in the colored region above the long-dashed line, while the dust GI mode grows faster in the region below the long-dashed line. The most unstable wavenumber on the dashed line is $k_{\mathrm{c},-}$. The dust GI becomes the fastest growing mode for smaller $\alpha$ cases since the stabilizing effect of the diffusion and the velocity dispersion is smaller.

4.2 With Turbulent Viscosity: Two-Component Viscous Gravitational Instability

Next, we discuss the linear stability in the case with the turbulent viscosity acting on the gas (Step 3 in Figure 1). In this case, we find a new instability that is different from the secular GI. We designate this new instability as two-component viscous gravitational instability (TVGI). TVGI did not appear in Takahashi & Inutsuka (2014), in which the basic equations did not guarantee the conservation of the total angular momentum. The origin of TVGI is the static mode that appears in the case without the turbulent viscosity (the mode labeled (C) in Figure 1). The static mode is the steady solution where the gas and the dust have the same azimuthal velocity ($\delta u_{y}=\delta v_{y}$) as mentioned in Section 4.1 and Appendix B. In this steady state, the self-gravity balances mainly with the Coriolis force, that is, $2\Omega\delta u_{y}-ik\delta\Phi=2\Omega\delta v_{y}-ik\delta\Phi\simeq 0$. This radial force balance is not realized once we consider the gas viscosity. The viscosity decreases $\delta u_{y}$, which results in the relative velocity in the azimuthal direction between the dust and the gas. The relative azimuthal velocity causes the decrease of $\delta v_{y}$ through the friction term. The decrease of both $\delta u_{y}$ and $\delta v_{y}$ prevents the radial force balance since the Coriolis force in the radial direction decreases. Thus, both dust and gas accumulate by the self-gravity, which is the physical interpretation of TVGI. The decrease of the Coriolis force by both viscosity and friction is important for the growth of TVGI. If we neglect the friction, this unstable mode becomes a static mode that satisfies the radial force balance for both components and

[TABLE]

[TABLE]

[TABLE]

(see Equations (17) – (22)). This indicates that TVGI is different from the so-called viscous instability, which grows in the one fluid system and does not need the friction (e.g., Lynden-Bell & Pringle, 1974; Schmit & Tscharnuter, 1995; Gammie, 1996; Lin & Kratter, 2016).

By using the terminal velocity approximation and assuming $t_{\mathrm{stop}}\ll\Omega^{-1}\ll n^{-1}$, the dispersion relation is reduced to the following quadratic equation:

[TABLE]

[TABLE]

We here neglect the second and higher order terms of $\nu k^{2}$ by assuming that the turbulence is weak ($\alpha\ll 1$). Two modes obtained from Equation (37) are the secular GI and TVGI. In fact, if we assume $\nu=0$, the solutions of Equation (37) gives the static mode ($n=0$) and the secular GI mode obtained from Equation (24). Figure 5 shows the dispersion relations of TVGI and the secular GI for $\alpha=10^{-3},\epsilon=0.1,\tau_{\mathrm{s}}=0.3$ and $Q=5$. We calculate the diffusion coefficient and the velocity dispersion from Equations (32) and (33). In this case, the secular GI does not grow. Since $n$ is real for both TVGI and the secular GI, they are not oscillating modes. We find that TVGI can grow even in the case where the secular GI is stable. The growth rate of TVGI is very small at long wavelengths because the angular momentum transport by the viscosity becomes ineffective as $k$ decreases. We derive the condition for the growth of TVGI from the approximate dispersion relation (Equation (37)) in the case where the secular GI is stable. We assume that the disk is self-gravitationally stable and $\omega_{\mathrm{gd}}^{2}>0$, that is,

[TABLE]

We consider cases with $B_{2}>0$ below, which is satisfied when $\omega_{\mathrm{gd}}^{2}>0$ and $\epsilon\nu k^{2}/t_{\mathrm{stop}}\ll\Omega^{2}$. In this case, the condition is that for a certain wavenumber $k$ Equation (37) has one negative solution and one positive solution, which is equivalent to the condition for the existence of the wavenumber where $B_{0}<0$ is satisfied. From Equation (40), we obtain the following quadratic equation for $k>0$:

[TABLE]

For the existence of wavenumbers that satisfy $B_{0}<0$, the discriminant of the right hand side of Equation (42) must be positive. We then find the following condition for the instability:

[TABLE]

The left hand side of Euqation (43) does not depend on $\nu$, which is partly because we assume the turbulence is weak. Another reason is that the infinitesimally small viscosity is enough for TVGI to grow. To clarify the physical meaning of the above condition, we only consider the leading term in Equation (43) by assuming $\alpha\ll\tau_{\mathrm{s}}\ll 1$ and $c_{\mathrm{d}}^{2}/c_{\mathrm{s}}^{2}\sim D\Omega/c_{\mathrm{s}}^{2}\sim\alpha$, and obtain

[TABLE]

or

[TABLE]

The left hand side of Equation (45) represents a timescale in which the dust transverses the length $H=c_{\mathrm{s}}\Omega^{-1}$ with the terminal velocity. The right hand side represents a diffusive time scale of the surface density perturbation of the dust with the length scale $\sim H$. Thus, Equation (43) represents the condition that the dust grains accumulate with the terminal velocity by overcoming the turbulent diffusion. This physical picture is analogous to that of the one-component secular GI discussed in Youdin (2011). We can also estimate the most unstable wavelength when the higher order terms of $\nu k^{2}$ are negligible. The growth rate of TVGI is also small in this limit since it is determined by the efficiency of the angular momentum transport by the turbulent viscosity. Hence, in this case, the growth rate of TVGI is approximately given by $-B_{0}/B_{1}$ (see Equation (37)). By neglecting the higher order terms of $\nu k^{2}$, the growth rate is reduced as follows:

[TABLE]

We here neglect $c_{\mathrm{d}}^{2}k^{2}$ since this term has a smaller effect to stabilize TVGI than the diffusion. The most unstable wavenumber $k_{\mathrm{max}}$ is of the order of a wavenumber where $n/\nu k^{2}$ has the local maximum, which is

[TABLE]

The right hand side of the Equation (47) is about $4.5H^{-1}$ for $\alpha=10^{-3},\epsilon=0.1,\tau_{\mathrm{s}}=0.3$ and $Q=5$, which is consistent with the most unstable wavenumber seen in Figure 5. We note that when the higher order terms of $\nu k^{2}$ are not negligible, the condition for the instability depends on $\nu$ although Equations (43) and (45) do not.

Next, we consider the other case where the secular GI is also unstable. Figure 6 shows the dispersion relations for $\alpha=10^{-4},\epsilon=0.1,\tau_{\mathrm{s}}=0.03$ and $Q=3$. Although TVGI and the secular GI are unstable in this case (Equations (30) and (43)), only one growing mode appears. This is because the mode exchange occurs between TVGI and the secular GI (Figure 7). If we neglect the viscosity, eigenvalues and eigenfunctions of the static mode and the secular GI mode degenerate at wavenumbers where the growth rate of the secular GI is zero. If we include the finite viscosity, the dispersion relations of the destabilized static mode (TVGI) and the secular GI reconnect at the wavenumbers, which results in one growing mode. As shown in Figure 7, at wavenumbers where the secular GI is unstable for $\nu=0$, we designate the growing mode as the secular GI even when $\nu$ is not zero. On the other hand, in a band of wavenumber $k$ where the secular GI is stable for $\nu=0$, we call the mode TVGI. Figure 8 shows the maximum growth rate of the instabilities for $\epsilon=0.05$ and $Q=10$. In the colored region above the dotted line, TVGI is the fastest growing mode. The dotted line almost coincides with the maximum $\alpha$ determined by Equation (30). The secular GI is the fastest growing mode in the region enclosed by the dotted line and the long-dashed line. Figure 8 also shows that, compared to the secular GI, TVGI can grow in a region where the turbulence and the friction are strong. We expect that TVGI grows earlier than the secular GI since the stopping time becomes larger as the dust grows in protoplanetary disks (see Figure 8). As described above, the dust grains accumulate through the growth of TVGI. Therefore, TVGI should be a promising mechanism to form planetesimals.

The turbulent viscosity makes two of the modified density waves overstable (cf., Schmit & Tscharnuter, 1995). We do not discuss the properties of this viscous overstability in this work since the growth time is found to be longer than the typical disk lifetime (see also Section 5.2).

5 Discussion

5.1 Effect of Disk Thickness

We assume that the disk is razor thin in the above analysis, while a real disk has a finite thickness. The thickness of the disk reduces the self-gravity estimated for the infinitesimally thin disk case. We here show the effect on the dispersion relation and the maximum growth rate. The self-gravitational potential $\delta\Phi$ reduced by the disk thickness is approximately given as follows (Vandervoort, 1970; Shu, 1984):

[TABLE]

where $H_{\mathrm{d}}$ is the thickness of the dust disk that is represented by the following (Youdin & Lithwick, 2007):

[TABLE]

Figure 9 shows how the dispersion relations of the secular GI and TVGI change due to the thickness for $\epsilon=0.05$ and $Q=10$. The normalized stopping time and the strength of turbulence ($\tau_{\mathrm{s}},\alpha$) are set to be ($0.1,2.5\times 10^{-5}$) for the left panel of Figure 9, while the right panel is for ($\tau_{\mathrm{s}},\alpha$)$=$($1,2.5\times 10^{-4}$). The most unstable mode is the secular GI and TVGI for the former and latter cases, respectively. The growth rate decreases by the effect of the disk thickness. Figure 10 shows the maximum growth rate as a function of $\tau_{\mathrm{s}}$ and $\alpha$ for the same parameters as in Figure 8. Although the maximum growth rate is smaller than that in Figure 8, the extent of the unstable region is almost the same, and the maximum $\alpha$ for a certain stopping time $\tau_{\mathrm{s}}$ does not change more than a factor of two.

The dust disk is generally thinner than the gaseous disk (see Equation (49)). In reality, the gas above the dust disk does not interact with the dust grains through the friction, while in the above analysis on the infinitesimally thin disk the back reaction is assumed to be exerted on all of the gas. We need to exclude the gas located above the dust disk from our analysis (see also, Latter & Rosca, 2017). In this subsection, we simply assume that the vertical density profile is given by the following Gaussian function and investigate stability in the dust disk:

[TABLE]

[TABLE]

where $\rho_{0}$ and $\rho_{\mathrm{d},0}$ are the mass density of the gas and the dust, respectively. The vertical distance from the disk midplane is denoted by $z$. In the following analysis, we use the surface density obtained from the vertical integration of those gas and dust densities. There is a large uncertainty in the appropriate range of the vertical integration. Here, we integrate these densities in $-3H_{\mathrm{d}}\leq z\leq 3H_{\mathrm{d}}$. Although the dust density at $z=3H_{\mathrm{d}}$ is about 0.01 times smaller than that at the midplane, the dust-to-gas mass ratio is still large because the dust scale height is much smaller than the gas scale height and the dust concentrates in the midplane. In addition, the strength of the friction on the dust is independent from the dust density, which means that the dust even at a low density region can be unstable as a result of the friction. The vertical integration of the mass densities gives

[TABLE]

[TABLE]

where $\mathrm{erf}(x)$ is the error function. The relation between the mid-plane dust-to-gas mass ratio and $\Sigma_{\mathrm{d},0}^{\prime}/\Sigma_{0}^{\prime}$ is

[TABLE]

We perform the linear analysis by using these surface densities and the dust-to-gas mass ratio $\Sigma_{\mathrm{d},0}^{\prime}/\Sigma_{0}^{\prime}$. In this analysis, we introduce the modified Toomre parameter for the gas disk $\tilde{Q}\equiv c_{\mathrm{s}}\Omega/\pi G\Sigma_{0}^{\prime}$. We modify the self-gravitational potential by using $H_{\mathrm{d}}$ for both gas and dust:

[TABLE]

Figure 11 shows the maximum growth rate as a function of total dust-to-gas mass ratio $\Sigma_{\mathrm{d},0}/\Sigma_{0}$ and $Q=c_{\mathrm{s}}\Omega/\pi G\Sigma_{0}$. The normalized stopping time $\tau_{\mathrm{s}}$ and the strength of the turbulence $\alpha$ are set to be $0.1$ and $10^{-4}$, respectively. In this case, $H_{\mathrm{d}}/H$ is about 0.03, and, thus, $\tilde{Q}/Q$ is about 14, meaning that $\tilde{Q}$ is larger than 10 in the whole parameter space shown in Figure 11. The dust-to-gas mass ratio in the dust disk $\Sigma_{\mathrm{d},0}^{\prime}/\Sigma_{0}^{\prime}$ is also about ten times larger than $\Sigma_{\mathrm{d},0}/\Sigma_{0}$. The secular GI becomes the most unstable mode only in very small parameter space, which is because $\tilde{Q}$ is too large and the dust GI mode is unstable enough to grow faster than the secular GI. On the other hand, TVGI can grow in much larger parameter space. Since the self-gravity of the dust is important for the growth, TVGI grows even for large $\tilde{Q}$. Therefore, we expect that TVGI operates in the dust disk even if we consider the vertical structure. The quantitative results do not change even if we vary the range of the vertical integration. In Appendix C, we show results for cases where we integrate the densities in $-H_{\mathrm{d}}\leq z\leq H_{\mathrm{d}}$ and $-2H_{\mathrm{d}}\leq z\leq 2H_{\mathrm{d}}$. We note that gas above the dust disk, which is free from the frictional force, would affect the motion in the dust disk through the gravitational interaction. Thus, we might underestimate the self-gravity in this analysis. The parameter space unstable to TVGI and the secular GI might become larger if we evaluate the self-gravity more precisely. In order to examine the effect of the upper gas on the instabilities in the dust disk, we need to perform the multidimensional analysis, which is beyond the scope of this paper. Moreover, the dust is expected to be less diffusive in high dust-to-gas ratio regions. This was also reported by Schreiber & Klahr (2018), in which they performed numerical simulations on the streaming instability. The less radial diffusivity makes the disk more unstable so that the instabilities grow in the larger parameter space than that shown in Figure 11.

5.2 Ring Formation through TVGI

The secular GI has been proposed as one of the possible mechanisms to create multiple ring-like structures recently observed in some protoplanetary disks. Takahashi & Inutsuka (2016) performed the linear analysis by adopting physical values obtained from the observation of HL Tau (ALMA Partnership et al., 2015) as the unperturbed state and showed that the multiple rings observed in the HL Tau disk can form via the secular GI. TVGI also has potential to be the mechanism to form multiple rings because dust grains accumulate in the radial direction through the growth of TVGI. Here, we show results of the linear analysis based on one of the disk model used in Takahashi & Inutsuka (2016) and discuss where TVGI operates and creates rings in the HL Tau disk.

We use the following dust surface density profile and gas temperature profile that reflect the observation of HL Tau (Pinte et al., 2016; Kwon et al., 2015), referred to as the “exponentially cutoff disk model” in Takahashi & Inutsuka (2016):

[TABLE]

[TABLE]

where $R_{\mathrm{c}}=80.2$ au, $\gamma=-0.2$. We assume that the mass of the central star is $1M_{\odot}$, and the dust size and the internal density are 3 mm and 3 $\mathrm{g\;cm}^{-3}$, respectively. The dust-to-gas mass ratio is set to be 0.02 as in Takahashi & Inutsuka (2016) (see also, Kwon et al., 2015). The strength of the turbulence $\alpha$ is assumed to be $3\times 10^{-4}$. Such weak turbulence reflects the observed very thin dust disk (Pinte et al., 2016). Figure 12 shows the most unstable wavelength and its growth time $n^{-1}$. In order to compare the extent of the regions unstable to the secular GI and/or TVGI, we show the results of the linear analysis with and without the turbulent viscosity. Only the secular GI grows when we do not consider the viscosity, and it forms multiple rings in the region $80\;\mathrm{au}\lesssim r\lesssim 100\;\mathrm{au}$. When we include the viscosity in the equations, there exist three instabilities: the secular GI, TVGI and the viscous overstability (cf., Schmit & Tscharnuter, 1995). The region where the viscous overstability is the most unstable mode is $r\lesssim 50\;\mathrm{au}$. We do not discuss the viscous overstability since the most unstable wavelength of the viscous overstability is larger than the radius and its growth time is longer than the typical disk lifetime. The regions where TVGI is the most unstable mode are $50\;\mathrm{au}\lesssim r\lesssim 80\;\mathrm{au}$ and $r\gtrsim 100\;\mathrm{au}$. As shown in Figure 12, the most unstable wavelength is about 10 au, which is almost independent of the radius. This wavelength is comparable to the width of the rings observed in the HL Tau disk. As discussed in Takahashi & Inutsuka (2016), the growth timescale is required to be less than $10^{6}$ yr since HL Tau is thought to be young. We find that the requirement on the timescale is satisfied in $r\gtrsim 53\;\mathrm{au}$. Therefore, we conclude that TVGI can form some of the multiple rings observed in HL Tau. We also find that TVGI grows in the larger extent of the disk than the secular GI.

Recent observations with ALMA have found ring structures in other disks (e.g., Andrews et al., 2016; Isella et al., 2016; Sheehan & Eisner, 2017; Fedele et al., 2018; Andrews et al., 2018). The Disk Substructures at High Angular Resolution Project (DSHARP, Andrews et al., 2018) observed 20 disks with high spatial resolution ($\sim 5$ au) and found that most of them have rings and gaps although non-axisymmetric structures were also observed (see also, Huang et al., 2018). According to Dullemond et al. (2018), some of the observed rings have width comparable to the gas scale height or smaller. Such relatively narrow rings can form via TVGI. Equation (47) gives the most unstable wavelength $\lambda_{\mathrm{max}}$ normalized by the gas scale height:

[TABLE]

For $D\Omega c_{\mathrm{s}}^{-2}\sim\alpha\ll\tau_{\mathrm{s}}$, $\lambda_{\mathrm{max}}/H$ is roughly given by

[TABLE]

Equation (44) determines the upper limit on the turbulent strength $\alpha$ for the growth of TVGI as follows:

[TABLE]

where $\epsilon\ll 1$ is assumed. From Equations (59) and (60), we obtain the upper limit on the maximum wavelength independent from any parameters:

[TABLE]

The most unstable wavelength is observable as the distance between peaks of adjacent rings. We note that Dullemond et al. (2018) fitted ring structures by using a Gaussian intensity profile and defined ring widths as standard deviation, which is different from the length that $\lambda_{\mathrm{max}}$ represents. The standard deviation should be smaller than $\lambda_{\mathrm{max}}$ if the observed rings formed via TVGI, and this is consistent with the narrowness of the observed rings. We therefore conclude that TVGI can be the origin of the observed narrow ring structures.

6 Conclusions

In this work, we re-formulate the equations of gas and dust in a protoplanetary disk in which the turbulence of the gas diffuses dust grains, based on the Reynolds averaging. The equations guarantee the conservation of the total angular momentum of the disk, which is not conserved in the equations used in previous studies. By using the newly formulated equations, we investigate the linear growth of the secular GI. We find that the secular GI is a monotonically growing mode while the secular GI found in the previous work is overstable depending on parameters. The appearance of the overstability is ascribed to the non-conservation of the total angular momentum. We also find the new instability referred to as TVGI if the turbulent viscosity is considered. TVGI grows in regions where $\alpha$ is larger and $\tau_{\mathrm{s}}$ is smaller compared with the secular GI. Since the normalized stopping time $\tau_{\mathrm{s}}$ increases as a dust grain grows, we can expect that TVGI operates in protoplanetary disks earlier than the secular GI.

We investigate the ring formation in the HL Tau disk via TVGI using the result of the linear analysis. Assuming that the dust size and the strength of turbulence $\alpha$ is constant in the radial direction, we show TVGI can form the observed multiple rings at $50\;\mathrm{au}\lesssim r\lesssim 80\;\mathrm{au}$ and $r\gtrsim 100\;\mathrm{au}$. For $80\;\mathrm{au}\lesssim r\lesssim 100\;\mathrm{au}$, the secular GI grows faster than TVGI. We expect that the relatively narrow rings observed recently form through TVGI (see Equation (61)). TVGI is also a promising mechanism to form planetesimals since the dust accumulates during the growth of TVGI. It is important to investigate non-linear growth and to discuss the planetesimal formation via TVGI, which is our future work.

We thank Shugo Michikoshi and Hiroshi Kobayashi for fruitful discussion. We also thank the anonymous referee for insightful and constructive comments that helped to improve the manuscript. This work was supported by JSPS KAKENHI Grant Number JP18J20360, 16H02160, 18H05436, 18H05437 and 19K14764, and in part by the National Astronomical Observatory of Japan (NAOJ) Atacama Large Millimeter/submillimeter Array (ALMA) Scientific Research Grant Numbers 2016-02A (SZT) .

Appendix A Formulation of basic equations with turbulent diffusion

In this appendix, we summarize the derivation of the phenomenological equations for the dust. We use the Reynolds averaging in this work. The Reynolds averaging is one of the techniques to analyze turbulence by focusing on longer-term evolution of systems than the typical time scale of the turbulence. We decompose a physical variable $A$ into the time-averaged term $\left<A\right>$ and the short-term fluctuation due to the turbulence $\Delta A\equiv A-\left<A\right>$, where $\left<\Delta A\right>=0$. We average the following equations for a infinitesimally thin and axisymmetric dust disk:

[TABLE]

[TABLE]

[TABLE]

First, we average the continuity equation and obtain the following equation that describes the time evolution of the mean surface density:

[TABLE]

We model the term $\left<\Delta\Sigma_{\mathrm{d}}\Delta v_{i}\right>$, based on the gradient diffusion hypothesis (see also, Cuzzi et al., 1993):

[TABLE]

[TABLE]

The Reynolds-averaged continuity equation is given by

[TABLE]

This equation is equivalent to Equation (1). Similarly, we average Equations (A2) and (A3) and obtain

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

represent the so called Reynolds stress. By using a closure relation $\left<\Delta v_{r}^{2}\right>=\left<\Delta v_{\phi}^{2}\right>=c_{\mathrm{d}}^{2}$, we obtain the effective pressure gradient force as follows (Shariff & Cuzzi, 2011):

[TABLE]

[TABLE]

[TABLE]

We neglect the terms $\sigma_{rr}^{\prime},\sigma_{r\phi},\sigma_{\phi\phi}^{\prime}$ for simplicity since a closure relation on these terms is uncertain. Moreover, we only consider the case that the dust grains are small and the friction is strong enough and assume $\Delta v_{r}=\Delta u_{r}$, $\Delta v_{\phi}=\Delta u_{\phi}$. Equations (A8) and (A9) become

[TABLE]

[TABLE]

The fourth term on the right hand side of Equation (A16) stands for gravity from the density fluctuation generated by the turbulence. The self-gravity is determined by the volume integral of the density. Here, we assume that the volume-integrated density fluctuation is small and neglect the fourth term on the right hand side of Equation (A16). In this work, we follow Cuzzi et al. (1993) and neglect the sixth term on the right hand side by assuming the term is smaller than the time derivative of $\left<\Sigma_{\mathrm{d}}\right>\left<v_{r}\right>$. The seventh term represents the advection of the linear momentum $\left<\Delta\Sigma_{\mathrm{d}}\Delta v_{r}\right>$ with the mean velocity $\left<v_{r}\right>$. This term is the same order of the advection of $\left<\Sigma_{\mathrm{d}}\right>\left<v_{r}\right>$ along the diffusion flow that appears in the left hand side. As a consequence, we rearrange Equations (A16) and (A17), and obtain the equations describing the time evolution of the mean velocity:

[TABLE]

[TABLE]

From Equation (A19), we obtain

[TABLE]

which is equivalent to Equation (10) since the mean specific angular momentum is $r\left<v_{\phi}\right>$. We can also obtain the evolutionary equation for the angular momentum of the dust that is equivalent to Equation (11) by using Equations (A7) and (A20):

[TABLE]

If we use Equations (4) and (12) for gas, we can derive an equation equivalent to Equation (16). Therefore, the equations formulated above hold the total angular momentum conservation. We omit the brackets representing the averaged value in the main part of this paper for convenience.

Appendix B Static mode in self-gravitating disks: No friction, Diffusion or Viscosity

In this appendix, we briefly summarize the dispersion relation of the dusty-gas disk for the case without the friction, the dust diffusion and turbulent gas viscosity. The dust and the gas interact with each other only through the self-gravity. In this case, the linearized equations of motion are the followings:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From the above four equations and Equations (17), (20), (23) with $D=0$, we obtain the six-order equation for $n$:

[TABLE]

[TABLE]

Equation (B5) gives four density waves ($F_{\mathrm{DWs}}(n,k)=0$) and two static modes ($n=0$). In a case where we ignore the self-gravity of the gas, which corresponds to $\Sigma_{0}\to 0$ in the above equation, $F_{\mathrm{DWs}}(n,k)=0$ gives

[TABLE]

[TABLE]

The latter one is the dispersion relation of the density wave in the dust disk (see also, Youdin, 2011). We can also reproduce the dispersion relation of the density wave in the gas disk from $F_{\mathrm{DWs}}(n,k)=0$ in a case where we ignore the self-gravity of the dust, which means that the dispersion relation $F_{\mathrm{DWs}}(n,k)=0$ describes two dust density waves and two gas density waves. The density wave becomes unstable if the self-gravity is strong enough. This unstable density wave is the classical gravitational instability.

A static mode is a steady solution of the linearized equations. This mode is referred to as a “neutral mode” in Youdin (2011). For both static modes we obtain here, $\delta u_{x}=\delta v_{x}=0$ (Equations (B2) and (B4)) and the radial force balance holds (Equations (B1) and (B3)). For one static mode, the gas and the dust have different azimuthal velocity while they have the same azimuthal velocity for the other static mode. The former static mode is destabilized by the friction. The friction does not change the latter static mode since there is no relative motion between the gas and the dust. The latter mode is destabilized once we consider both friction and viscosity.

Appendix C Linear analysis in a thin dust disk

In Section 5.1, we discuss the stability in the thin dust disk by assuming the vertical density profiles. However, the range of the vertical integration to obtain the surface densities, $\Sigma_{0}^{\prime}$ and $\Sigma_{\mathrm{d},0}^{\prime}$, is uncertain. Thus, we here show results for cases where the range of the vertical integration is $-H_{\mathrm{d}}\leq z\leq H_{\mathrm{d}}$ and $-2H_{\mathrm{d}}\leq z\leq 2H_{\mathrm{d}}$.

Figure 13 shows the maximum growth rate as a function of total dust-to-gas mass ratio $\Sigma_{\mathrm{d},0}/\Sigma_{0}$ and $Q=c_{\mathrm{s}}\Omega/\pi G\Sigma_{0}$ for the case with the integral in $-H_{\mathrm{d}}\leq z\leq H_{\mathrm{d}}$. The normalized stopping time $\tau_{\mathrm{s}}$ and the strength of the turbulence $\alpha$ are set to be $0.1$ and $10^{-4}$, which is same with those in Figure 11. In this case, $\tilde{Q}/Q$ is about 42, and $\Sigma_{\mathrm{d},0}^{\prime}/\Sigma_{0}^{\prime}$ is about 28 times larger than the total dust-to-gas mass ratio $\Sigma_{\mathrm{d},0}/\Sigma_{0}$. TVGI grows fastest above the gray dotted line. The maximum value of $Q$ for each total dust-to-gas mass ratio is smaller than that seen in Figure 11, which is because $\tilde{Q}/Q$ is larger, and the gas disk is self-gravitationally more stable. For the same reason, the region where the secular GI is fastest growing mode does not appear, and the dust GI grows fastest below the gray dotted line.

Figure 14 is a result for the case with the integral in $-2H_{\mathrm{d}}\leq z\leq 2H_{\mathrm{d}}$. The ratio $\tilde{Q}/Q$ is about 21, and $\Sigma_{\mathrm{d},0}^{\prime}/\Sigma_{0}^{\prime}$ is about 20 times larger than $\Sigma_{\mathrm{d},0}/\Sigma_{0}$ in this case. The secular GI is the fastest growing mode in the region enclosed by the black dotted and dashed lines. TVGI grows fastest in most of the colored region, which is qualitatively same with Figures 11 and 13.