Secular Gravitational Instability of Drifting Dust in Protoplanetary Disks: Formation of Dusty Rings without Significant Gas Substructures

We consider the time evolution of axisymmetric dust and gas disks. Because secular GI requires the dust–gas friction, the dust layer around the disk midplane seems to be the most important region. The gas above the dust layer will also contribute to the growth of secular GI through gravitational interaction while its frictional interaction with the dust is relatively weak. Thus, we consider a "lower layer" that includes those dust and gas driving secular GI. Hereafter, we do not consider the vertical extent of the lower layer as it is beyond the scope of this study. The following equations govern the evolution of the physical quantities vertically integrated within the lower layer orbiting around a central star with mass M_* in cylindrical coordinates (r, ϕ):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r{\rm{\Sigma }}{u}_{r}\right)=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }}\left(\displaystyle \frac{\partial {u}_{r}}{\partial t}+{u}_{r}\displaystyle \frac{\partial {u}_{r}}{\partial r}\right) & = & -\displaystyle \frac{\partial }{\partial r}\left({c}_{{\rm{s}}}^{2}{\rm{\Sigma }}\right)-{\rm{\Sigma }}\displaystyle \frac{\partial }{\partial r}\left({\rm{\Phi }}-\displaystyle \frac{{{GM}}_{* }}{r}\right)\\ & & +\,\displaystyle \frac{2}{r}\displaystyle \frac{\partial }{\partial r}\left(r{\rm{\Sigma }}\nu \displaystyle \frac{\partial {u}_{r}}{\partial r}\right)-\displaystyle \frac{2{\rm{\Sigma }}\nu {u}_{r}}{{r}^{2}}\\ & & -\,\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{2{\rm{\Sigma }}\nu }{3r}\displaystyle \frac{\partial {{ru}}_{r}}{\partial r}\right)\\ & & +\,{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{{v}_{r}-{u}_{r}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Sigma }}\left(\displaystyle \frac{\partial {j}_{{\rm{g}}}}{\partial t}+{u}_{r}\displaystyle \frac{\partial {j}_{{\rm{g}}}}{\partial r}\right) & = & \displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left[{r}^{3}{\rm{\Sigma }}\nu \displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{u}_{\phi }}{r}\right)\right]\\ & & +\,{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{{j}_{{\rm{d}}}-{j}_{{\rm{g}}}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left(r{{\rm{\Sigma }}}_{{\rm{d}}}{v}_{r}\right)=\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({rD}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r}\right),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Sigma }}}_{{\rm{d}}}\left[\displaystyle \frac{\partial {v}_{r}}{\partial t}+\left({v}_{r}-\displaystyle \frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r}\right)\displaystyle \frac{\partial {v}_{r}}{\partial r}\right]\\ \ =\,{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{{v}_{\phi }^{2}}{r}-\displaystyle \frac{\partial }{\partial r}\left({c}_{{\rm{d}}}^{2}{{\rm{\Sigma }}}_{{\rm{d}}}\right)-{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{\partial }{\partial r}\left({\rm{\Phi }}-\displaystyle \frac{{{GM}}_{* }}{r}\right)\\ \quad -\,{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{{v}_{r}-{u}_{r}}{{t}_{\mathrm{stop}}}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left({{rv}}_{r}D\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}\left[\displaystyle \frac{\partial {j}_{{\rm{d}}}}{\partial t}+\left({v}_{r}-\displaystyle \frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r}\right)\displaystyle \frac{\partial {j}_{{\rm{d}}}}{\partial r}\right]=-{{\rm{\Sigma }}}_{{\rm{d}}}\displaystyle \frac{{j}_{{\rm{d}}}-{j}_{{\rm{g}}}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}^{2}{\rm{\Phi }}=4\pi G\left({\rm{\Sigma }}+{{\rm{\Sigma }}}_{{\rm{d}}}\right)\delta (z),\end{eqnarray} \tag{ 7 }$$

where Σ and ${{\rm{\Sigma }}}_{{\rm{d}}}$ are the gas and dust surface densities of the lower layer, u_i and v_i are the ith component of the velocities, j_g and j_d denote the specific angular momenta of the gas and the dust, ${c}_{{\rm{s}}}$ is the sound speed, and ${c}_{{\rm{d}}}$ and D are the velocity dispersion and diffusion coefficient of the dust, respectively. The self-gravitational potential is denoted by Φ and the gravitational constant by G. We denote the coefficient of the turbulent viscosity by ν, measured from the dimensionless parameter $\alpha \equiv \nu {\rm{\Omega }}{c}_{{\rm{s}}}^{-2}$ (Shakura & Sunyaev 1973), where Ω is the angular velocity of the gas disk.

The terms including $D\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$ in Equations (5) and (6) were introduced by Tominaga et al. (2019). Those terms are derived from the mean-field approximation and represent the mean momentum flow caused by the diffusion. By including those terms, both total linear and angular momenta are conserved (see also Appendix A of Tominaga et al. 2019). The last term on the right-hand side of Equation (5) affects the dynamics if α ≳ 10⁻³.

Because the growth timescale of secular GI is much longer than the Keplerian periods, numerical errors due to time integration need to be minimized as much as possible. In addition to this, we should also minimize numerical viscosity, which unphysically smooths out advecting density fluctuations. As we will show in Section 4, secular GI originates from the "dust-streaming mode" where the dust surface density perturbation propagates with the radial drift velocity. Thus, we need to treat the advection of dust while keeping numerical dissipation as minimal as possible.

For these reasons, we use a numerical method that incorporates the symplectic integrator and the Lagrangian-cell method (Tominaga et al. 2018). The symplectic integrator allows us to avoid a monotonic increase of the error due to the time integration. The Lagrangian-cell method makes it possible to properly deal with the advection without numerical dissipation due to spatial discretization. In this scheme, density and pressure are defined at cell centers, and velocity and angular momentum are defined at cell boundaries. We use ${m}_{{\rm{d}},i}$, ${r}_{{\rm{d}},i}$, ${m}_{{\rm{g}},i}$, and ${r}_{{\rm{g}},i}$ to denote the masses and radial positions of the ith dust and gas cells (rings), respectively. Mass ${m}_{s,i+1/2}\,\equiv ({m}_{s,i+1}+{m}_{s,i})/2$ is also assigned at a cell boundary at ${r}_{s}\,={r}_{s,i+1/2}$ between the ith and (i + 1)th cells, where $s={\rm{d}},{\rm{g}}$. The gas and dust velocities and angular momenta are defined at the cell boundaries and denoted by ${u}_{r,i+1/2},{v}_{r,i+1/2},{j}_{{\rm{g}},i+1/2}$ and ${j}_{{\rm{d}},i+1/2}$. We note that the position of a dust cell does not necessarily coincide with that of a gas cell.

We adopted the operator-splitting technique with second-order accuracy in time (Inoue & Inutsuka 2008). The basic equations are split into four parts: (A) the frictional interaction part, (B) the dust diffusion part, which includes the last term on the right-hand side of Equation (5), (C) the gas viscosity part, and (D) the hydrodynamic part with the other physical processes for gas and dust.

Parts (A) and (D) follow Tominaga et al. (2018). In part (A), the frictional interaction between dust and gas is solved with the piecewise exact solution, which is free from the limitation on the time step Δt due to the small stopping time ${t}_{\mathrm{stop}}$. In this step, we use interpolation functions shown in Appendix 3 of Tominaga et al. (2018) for the radial velocities and the specific angular momenta to calculate the frictional interaction. In part (D), we calculate the forces (e.g., pressure gradient force, self-gravity) acting on the dust and gas cell boundaries and update their positions and velocities with the leap-frog integrator. Self-gravity is represented as a superposition of infinitesimally thin-ring gravity (see, Tominaga et al. 2018). We weaken the ring gravity with a softening length of half a cell's width (∼0.1 au) and rescale it when the cell width becomes larger in time.

We newly implement parts (B) and (C) based on the second-order Runge–Kutta integrator. In part (B), we first interpolate the dust surface density using a quadratic function to evaluate the mass flux $D\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$ at the cell boundaries. We then displace the dust-cell boundaries using velocity $-{{\rm{\Sigma }}}_{{\rm{d}}}^{-1}D\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$, which corresponds to the diffusion term in Equation (4). We do not change the angular momenta of the dust-cell boundaries when displacing the dust as the Lagrange derivative in the equation of motion includes the advection along the diffusion flow (Equation (6)). On the other hand, the radial linear momentum of the ith cell boundary, ${m}_{{\rm{d}},i+1/2}{v}_{r,i+1/2}$, changes according to the last term on the right-hand side of Equation (5). We evaluate $F(r)\equiv {{rv}}_{r}D\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$ at $r={r}_{{\rm{d}},i+1},{r}_{{\rm{d}},i}$ and update ${m}_{{\rm{d}},i+1/2}{v}_{r,i+1/2}$ using $F({r}_{{\rm{d}},i+1})\,-F({r}_{{\rm{d}},i})$.

In part (C), we interpolate the radial velocity u_r and the specific angular momentum j_g in order to obtain the radial and angular momentum fluxes due to the gas viscosity at $r={r}_{{\rm{g}},i+1},{r}_{{\rm{g}},i}$. From Equations (2) and (3), the radial and angular momentum changes are given by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial {\rm{\Sigma }}{u}_{r}}{\partial t}\right)}_{\mathrm{vis}}\equiv \displaystyle \frac{1}{r}\displaystyle \frac{\partial {F}_{\mathrm{vis},r}}{\partial r}-\displaystyle \frac{2{u}_{r}}{3{r}^{3}}\displaystyle \frac{\partial {r}^{2}{\rm{\Sigma }}\nu }{\partial r},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial {\rm{\Sigma }}{j}_{{\rm{g}}}}{\partial t}\right)}_{\mathrm{vis}}\equiv \displaystyle \frac{1}{r}\displaystyle \frac{\partial {F}_{\mathrm{vis},\phi }}{\partial r},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{vis},r}\equiv \displaystyle \frac{4}{3}r{\rm{\Sigma }}\nu \displaystyle \frac{\partial {u}_{r}}{\partial r},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{vis},\phi }\equiv {r}^{3}{\rm{\Sigma }}\nu \displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{u}_{\phi }}{r}\right).\end{eqnarray} \tag{ 11 }$$

Integrating the above equations in the radial and azimuthal directions, we obtain the following equations for updating ${m}_{{\rm{g}},i+1/2}{u}_{r,i+1/2}$ and ${m}_{{\rm{g}},i+1/2}{j}_{{\rm{g}},i+1/2}$:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {m}_{{\rm{g}},i+1/2}{u}_{r,i+1/2}}{\partial t} & = & {\left[2\pi {F}_{\mathrm{vis},r}\right]}_{{r}_{{\rm{g}},i}}^{{r}_{{\rm{g}},i+1}}\\ & & -\,\displaystyle \frac{4\pi }{3}{\displaystyle \int }_{{r}_{{\rm{g}},i}}^{{r}_{{\rm{g}},i+1}}\displaystyle \frac{{u}_{r}}{{r}^{2}}\displaystyle \frac{\partial {r}^{2}{\rm{\Sigma }}\nu }{\partial r}{dr},\end{array}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {m}_{{\rm{g}},i+1/2}{j}_{{\rm{g}},i+1/2}}{\partial t}={\left[2\pi {F}_{\mathrm{vis},\phi }\right]}_{{r}_{{\rm{g}},i}}^{{r}_{{\rm{g}},i+1}},\end{eqnarray} \tag{ 13 }$$

where $[A(r){]}_{{r}_{1}}^{{r}_{2}}\equiv A({r}_{2})-A({r}_{1})$. We evaluate the last term on the right-hand side of Equation (12) using $-(4\pi {u}_{r,i+1/2}/3{r}_{{\rm{g}},i+1/2}^{2}){\left[{r}^{2}{\rm{\Sigma }}\nu \right]}_{{r}_{{\rm{g}},i}}^{{r}_{{\rm{g}},i+1}}$.

We determine the time step based on the following:

$$\begin{eqnarray}&&{\rm{\Delta }}t=\min \left({\rm{\Delta }}{t}_{{\rm{g}}},{\rm{\Delta }}{t}_{{\rm{d}}},{\rm{\Delta }}{t}_{\mathrm{diff}},{\rm{\Delta }}{t}_{\mathrm{vis}}\right),\end{eqnarray} \tag{ 14 }$$

where ${\rm{\Delta }}{t}_{{\rm{g}}}$ and ${\rm{\Delta }}{t}_{{\rm{d}}}$ are the time steps determined by the Courant–Friedrich–Levy condition for the gas and dust equations without dust diffusion or viscosity, and ${\rm{\Delta }}{t}_{\mathrm{diff}}$ and ${\rm{\Delta }}{t}_{\mathrm{vis}}$ are those determined by the diffusion term and the viscosity term: ${\rm{\Delta }}{t}_{\mathrm{diff}}=0.125\times \min ({({r}_{{\rm{d}},i+1/2}-{r}_{{\rm{d}},i-1/2})}^{2}/D)$ and ${\rm{\Delta }}{t}_{\mathrm{vis}}=0.125\times \min ({({r}_{{\rm{g}},i+1/2}-{r}_{{\rm{g}},i-1/2})}^{2}/\nu )$. As shown in the next section, the dust rings become spiky as secular GI grows, and then the dust diffusion tends to limit ${\rm{\Delta }}t$. Thus, we adopt super-time-stepping (STS) (Alexiades et al. 1996) in part (B) when the time step Δt is limited by the dust diffusion. The number of substeps and the stability parameter in STS are fixed to be 4 and 0.1, respectively. We just moderately accelerate time stepping as the physical diffusion timescale becomes short (≪1 Keplerian period), and the total time step in STS should not exceed it.

The number of cells N_r is 1024 for both dust and gas. The initial position of the inner boundaries is 10 au. We space the dust and gas domains so that each cell has the same mass and the outer boundaries are located at r ≃ 300 au. We fixed the gas and dust outer boundaries in simulations. The gas inner boundary moves so that the innermost gas surface density is constant in time. We have dust cells at $r\lt {r}_{{\rm{g}},1}$ move inward with the steady drift velocity (e.g., Nakagawa et al. 1986) estimated with the initial density profiles at $r={r}_{{\rm{g}},1}$.

We consider gas and dust disks around a $1{M}_{\odot }$ mass star, and those initial surface density profiles are given by the following power-law function:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{100}{\left(\displaystyle \frac{r}{100\mathrm{au}}\right)}^{-q}\exp \left(-\displaystyle \frac{r}{100\,\mathrm{au}}\right),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}(r)={{\rm{\Sigma }}}_{{\rm{d}},100}{\left(\displaystyle \frac{r}{100\mathrm{au}}\right)}^{-q}\exp \left(-\displaystyle \frac{r}{100\,\mathrm{au}}\right),\end{eqnarray} \tag{ 16 }$$

where Σ₁₀₀ and ${{\rm{\Sigma }}}_{{\rm{d}},100}$ are constants. In this study, we use Toomre's Q value of the gas, $Q={c}_{{\rm{s}}}{\rm{\Omega }}/\pi G{\rm{\Sigma }}$, to show how massive the lower layer is. One obtains Σ₁₀₀ from the Q value at r = 100 au (Table 1), Keplerian angular velocity, and the temperature shown below. The dust surface density at r = 100 au is determined by the assumed initial dust-to-gas ratio in the lower layer ${{\rm{\Sigma }}}_{{\rm{d}}}$/Σ. We note that ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ represents the dust-to-gas mass ratio averaged in the lower layer, which is different from ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$, where ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}$ and ${{\rm{\Sigma }}}_{\mathrm{tot}}$ are the total surface densities of the dust and gas disks including both upper and lower layers. In weakly turbulent gas disks, ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ easily becomes higher than ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$ by an order of magnitude. In this work, considering dust-rich disks ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}=0.05$, we assume ${{\rm{\Sigma }}}_{{\rm{d}}}$/Σ = 0.1. The power-law index is one of the important parameters that characterize the dust drift. Kitamura et al. (2002) observed 13 disks around T Tauri stars and found that the power-law index is 0–1 in most cases. Motivated by this study, we take the median value and fix q = 0.5 in this work.

Table 1.  Summary of Parameters and Results

Label	Q100a	αb	Resultsc	Time That Simulations Last for
Q4a10	4	1 × 10−3	thin dense dust rings	2.1 × 104 yr
Q4a20	4	$2\times {10}^{-3}$	transient low-contrast rings	5.9 × 104 yr
Q5a5	5	5 × 10−4	thin dense dust rings	1.9 × 104 yr
Q5a8	5	8 × 10−4	transient low-contrast dust rings	5.6 × 104 yr
Q5a8Ld	5	8 × 10−4	thin dense dust rings	2.5 × 104 yr
Q6a3	6	3 × 10−4	thin dense dust rings	2.7 × 104 yr
Q6a5	6	5 × 10−4	transient low-contrast dust rings	6.2 × 104 yr

Notes.

^aToomre's Q value for the lower layer of a gas disk at r = 100 au. ^bThe strength of turbulence. ^cResults of simulations: formation of "thin dense dust rings" or "transient low-contrast dust rings." ^dThe letter "L" means a run with six times larger perturbations (see Section 4.2).

Download table as:  ASCIITypeset image

The gas disk in our simulations is locally isothermal, and the temperature profile T(r) is

$$\begin{eqnarray}&&T(r)=10\,{\rm{K}}{\left(\displaystyle \frac{r}{100\mathrm{au}}\right)}^{-1/2},\end{eqnarray} \tag{ 17 }$$

where we mimic a disk passively heated by stellar radiation (e.g., Chiang & Goldreich 1997). We calculate the sound speed ${c}_{{\rm{s}}}$ at each radius assuming the molecular weight to be 2.34. The gas scale height $H\equiv {c}_{{\rm{s}}}/{\rm{\Omega }}$ is $H\simeq 6.3\,\mathrm{au}{\left(r/100\mathrm{au}\right)}^{5/4}$. The initial azimuthal velocity is determined by the radial force balance without using the friction, the diffusion, or the viscosity. We put initially random perturbations to the cell positions and the velocities. Amplitudes of the position and velocity perturbations are 5% of the cell widths and ${c}_{{\rm{d}}}$, respectively.

The normalized stopping time ${t}_{\mathrm{stop}}{\rm{\Omega }}$ is an important parameter characterizing the growth rate and the unstable wavelengths of secular GI. The timescale of dust coagulation is shorter than that of secular GI when dust grains are small (∼100 orbits), and thus dust grains would grow up to the drift-limited size (see also Okuzumi et al. 2012). We mimic this situation by simply setting ${t}_{\mathrm{stop}}{\rm{\Omega }}$ constant. Assuming the total dust–gas mass ratio ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$, one can estimate the drift-limited stopping time. In Appendix A, we derive the drift-limited stopping time based on our disk model and show that for ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}=0.05$, the drift-limited ${t}_{\mathrm{stop}}{\rm{\Omega }}$ is about 0.6. We do not consider fragmentation-limited dust sizes because we are focusing on dynamics at outer regions in weakly turbulent gas disks (see also Birnstiel et al. 2009, 2012).

Although the diffusion coefficient D and the velocity dispersion ${c}_{{\rm{d}}}$ depend on ${t}_{\mathrm{stop}}{\rm{\Omega }}$ and α (Youdin & Lithwick 2007), those are well approximated by $D\simeq \alpha {c}_{{\rm{s}}}^{2}{{\rm{\Omega }}}^{-1}$ and ${c}_{{\rm{d}}}\simeq \sqrt{\alpha }{c}_{{\rm{s}}}$ for ${t}_{\mathrm{stop}}{\rm{\Omega }}\lt 1$. Thus, we use the relations $D=\alpha {c}_{{\rm{s}}}^{2}{{\rm{\Omega }}}^{-1}$ and ${c}_{{\rm{d}}}=\sqrt{\alpha }{c}_{{\rm{s}}}$ in the present simulations for simplicity. In this study, we perform numerical simulations with different Toomre's Q values and turbulent strength α. We summarize the parameters in Table 1.

Our choice of parameters for the density distribution to study the global and nonlinear outcomes of secular GI is optimistic. For example, disks are massive (see also Table 2). On the other hand, the disk masses and the dust-to-gas surface density ratio ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$ are not observationally well constrained because of uncertainty in the dust opacity. In addition, recent work shows that ignoring the scattering effect of the dust thermal emission underestimates the dust mass in a disk (Zhu et al. 2019). In the present study, we thus assume massive disks where the secular GI relatively easily grows.

Table 2.  List of the Evaluated Masses

	${R}_{{\rm{d}},\mathrm{out}}=120$ au			${R}_{{\rm{d}},\mathrm{out}}=200$ au
Label and Time	${M}_{\mathrm{ring},\mathrm{tot}}$ [M⊕]	Md [M⊕]	${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$	${M}_{\mathrm{ring},\mathrm{tot}}$ [M⊕]	Md [M⊕]	${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$
Q4a10 (t = 2.0 × 104 yr)	4.7 × 102	1.2 × 103	0.38	8.5 × 102	1.8 × 103	0.46
Q5a5 (t = 1.6 × 104 yr)	4.5 × 102	9.8 × 102	0.46	7.4 × 102	1.4 × 103	0.51
Q5a8L (t = 2.3 × 104 yr)	3.2 × 102	9.8 × 102	0.33	6.9 × 102	1.4 × 103	0.47
Q6a3 (t = 1.8 × 104 yr)	3.5 × 102	8.1 × 102	0.42	7.3 × 102	1.2 × 103	0.60

Download table as:  ASCIITypeset image

Figure 1 shows the time evolution of the dust and gas surface densities from Q4a10 run. We also show the surface density evolution from a run in which we switch off self-gravity. We note that the gas disk is self-gravitationally stable and the dust GI is also stabilized because of dust diffusion. We observe the formation of multiple rings and gaps only when we switch on self-gravity. Thus, the ring-gap formation results from secular GI (see also Section 4.1). Secular GI grows at wavelengths $\sim {c}_{{\rm{s}}}/{\rm{\Omega }}$, resulting in the ring-gap formation in the dust disk (see the dashed line in Figure 1). The resultant rings become much thinner, and the dust surface density increases. We do not observe significant substructures in the gas disk, which we discuss in more detail in Section 4.

Zoom In Zoom Out Reset image size

On the left panel of Figure 2, we show the motion of dust cells that compose the resultant rings and gaps. The red dashed line shows the motions of the 525th cell as a reference cell comprising one dust ring. We note that we reduced the number of cells to plot the left figure of Figure 2. Because our numerical scheme is based on the Lagrangian-cell method, the motion of cells represents the actual motion of dust, and a cell-concentrating region corresponds to a high-density region. The dust initially moves inward with the steady drift velocity v_dri (Figure 3):

$$\begin{eqnarray}&&{v}_{\mathrm{dri}}\equiv -\displaystyle \frac{2{t}_{\mathrm{stop}}{\rm{\Omega }}}{{(1+{\epsilon }_{0})}^{2}+{\left({t}_{\mathrm{stop}}{\rm{\Omega }}\right)}^{2}}\eta ^{\prime} r{\rm{\Omega }},\end{eqnarray} \tag{ 18 }$$

where ₀ is the initial surface density ratio ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ and $\eta ^{\prime} \equiv -\left({c}_{{\rm{s}}}^{2}/2{r}^{2}{{\rm{\Omega }}}^{2}\right)\partial \mathrm{ln}\left({c}_{{\rm{s}}}^{2}{\rm{\Sigma }}\right)/\partial \mathrm{ln}r$, and finally concentrates in the rings. One can see that the dust density perturbations also move inward with the drift velocity. A significant dust concentration (t = 1.3 × 10⁵ yr) results from the self-gravitational collapse of the dust rings, which was also found by Tominaga et al. (2018). On the right panel of Figure 2, we plot the radial velocity at the dust-cell boundaries ($r(t)={r}_{{\rm{d}},i+1/2}$) that compose one dust ring (i = 520–530). The position and the radial velocity of the 525th dust-cell boundary are shown by the red solid line on the right panels. Those dust cells initially move inward with the drift velocity. As secular GI grows (r ≳ 77 au), the radial velocities deviate from the drift velocity and spread in the ${v}_{r}-r$ plane. When the cell width becomes smaller than the softening length as the dust surface density increases, the collapsing motion is suppressed (e.g., r ≃ 77 au for i = 526–530). As the dust-to-gas surface density ratio gradually increases, the radial velocities of those cell boundaries decrease.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the cumulative mass of the dust disks ${M}_{{\rm{d}}}(\lt r)$. We can roughly evaluate the ring mass from Figure 4 and find that large amounts (∼tens of ${M}_{\oplus }$) of dust grains reside in the individual ring. We define the radius of an individual dust ring as the position of the local maximum of the dust surface density and estimate the mass of the ith ring ${M}_{\mathrm{ring},i}$ by summing up the mass of the dust cells between adjacent local minima. We then compare the dust disk mass M_d and the total ring mass ${M}_{\mathrm{ring},\mathrm{tot}}\equiv {\sum }_{i=1}^{{N}_{\mathrm{ring}}}{M}_{\mathrm{ring},i}$, where N_ring is the number of dust rings whose amplitudes increase or are saturated at a certain time. In order to discuss the mass fraction of dust grains that are concentrated into rings, we set a certain radius denoted by ${R}_{{\rm{d}},\mathrm{out}}$ and use the dust cells whose initial radii are smaller than ${R}_{{\rm{d}},\mathrm{out}}$ to calculate M_d and ${M}_{\mathrm{ring},\mathrm{tot}}$. The ring masses are calculated at the time when the dust-to-gas mass ratio in one of the resultant dust rings exceeds unity.

Zoom In Zoom Out Reset image size

The results for ${R}_{{\rm{d}},\mathrm{out}}=120$ and 200 au are summarized in Table 2. We find that secular GI can convert tens of percent of the dust mass into ring mass. The ring mass ${M}_{\mathrm{ring},\mathrm{tot}}$ increases as ${R}_{{\rm{d}},\mathrm{out}}$ increases because dust grains drifting from the outer disk are concentrated into rings at the inner region (see Figure 2).

The mass M_d includes the mass of dust grains that are initially located at the secular-GI-stable region (see also Section 4.1) and never concentrated into rings. Thus, the mass fraction ${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$ increases if one excludes those dust grains from the above estimation. Table 3 summarizes ${M}_{\mathrm{ring},\mathrm{tot}},{M}_{{\rm{d}}}$ and ${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$ estimated for dust grains initially located at ${R}_{{\rm{d}},\mathrm{in}}\leqslant r\leqslant {R}_{{\rm{d}},\mathrm{out}}$, where we set ${R}_{{\rm{d}},\mathrm{in}}=60\,\mathrm{au}$ because secular GI can grow at r > 60 au in all runs. Over half of the dust masses are collected into the rings in most of the runs. In fact, 88% of the dust grains are saved in the dust rings in the case of the Q6a3 run with ${R}_{{\rm{d}},\mathrm{in}}=60\,\mathrm{au}$ and ${R}_{{\rm{d}},\mathrm{out}}=200\,\mathrm{au}$.

Table 3.  List of the Evaluated Masses

	${R}_{{\rm{d}},\mathrm{in}}=60\,\mathrm{au},{R}_{{\rm{d}},\mathrm{out}}=120\,\mathrm{au}$			${R}_{{\rm{d}},\mathrm{in}}=60\,\mathrm{au},{R}_{{\rm{d}},\mathrm{out}}=200\,\mathrm{au}$
Label and Time	${M}_{\mathrm{ring},\mathrm{tot}}$ [M⊕]	Md [M⊕]	${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$	${M}_{\mathrm{ring},\mathrm{tot}}$ [M⊕]	Md [M⊕]	${M}_{\mathrm{ring},\mathrm{tot}}/{M}_{{\rm{d}}}$
Q4a10 (t = 2.0 × 104 yr)	4.7 × 102	6.5 × 102	0.72	8.5 × 102	1.2 × 103	0.68
Q5a5 (t = 1.6 × 104 yr)	4.5 × 102	5.2 × 102	0.85	7.4 × 102	9.9 × 102	0.75
Q5a8L (t = 2.3 × 104 yr)	3.2 × 102	5.2 × 102	0.62	6.9 × 102	9.9 × 102	0.69
Q6a3 (t = 1.8 × 104 yr)	3.5 × 102	4.3 × 102	0.79	7.3 × 102	8.2 × 102	0.88

Download table as:  ASCIITypeset image

The dust-to-gas mass ratio becomes comparable to or higher than unity as a consequence of ring collapse, which makes the dust drift velocity smaller (see also the right panel of Figure 2). The increase in dust density is saturated because of the balance between diffusion and self-gravitational collapse. Figure 5 compares two velocities around one dust ring whose radius is $r\simeq 73.6\,\mathrm{au}$ at $t=2.1\times {10}^{4}\,\mathrm{yr}$: (1) the dust velocity with respect to the ring velocity v_ring = −1.5 × 10⁻³ au yr⁻¹ that we measured from the data and (2) the velocity due to the dust diffusion ${v}_{\mathrm{diff}}\equiv -D{{\rm{\Sigma }}}_{{\rm{d}}}^{-1}\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$. The red line shows their sum. We note that the red filled circles in Figure 5 are the data points, which show that the dusty ring is well resolved while the adjacent gaps are not. In the well-resolved ring, we find $| {v}_{r}-{v}_{\mathrm{ring}}| \simeq | {v}_{\mathrm{diff}}| \gt | {v}_{r}-{v}_{\mathrm{ring}}+{v}_{\mathrm{diff}}|$, which means that the dust diffusion quenches further self-gravitational collapse of the spiky dust rings. Specifically, we can see the equilibrium at the inner half of the well-resolved ring. At the outer half of the ring, the red line shows an increasing trend, which means that the radial convergence speed is lower than the velocity due to diffusion. This slow converging flow results from the deceleration due to the last term on the right-hand side of Equation (5), i.e., ${r}^{-1}\partial F(r)/\partial r$, where $F(r)={{rv}}_{r}D\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$. Figure 6 compares four forces per mass acting on the dust: the self-gravity, the pressure gradient force (${{\rm{\Sigma }}}_{{\rm{d}}}^{-1}{c}_{{\rm{d}}}^{2}\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$), the sum of the curvature term and the stellar gravity, and ${{\rm{\Sigma }}}_{{\rm{d}}}^{-1}{r}^{-1}\partial F(r)/\partial r$. The magnitude of the force ${{\rm{\Sigma }}}_{{\rm{d}}}^{-1}{r}^{-1}\partial F(r)/\partial r$ is comparable to that of the self-gravity, meaning that the term decelerates the radial speed of the dust coming from the outer gap toward the ring.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because we weaken the self-gravity with a constant softening length for the rings, the final dust surface density might be underestimated.⁴ To check whether or not the dust density is underestimated, we evaluate a resultant dust surface density ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}$, which we would obtain if we ignore the softening. The diffusion timescale becomes comparable to the self-gravitational-collapse timescale after nonlinear growth:

$$\begin{eqnarray}&&{D}^{2}{k}_{{\rm{c}}}^{4}\sim 2\pi G{{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}{k}_{{\rm{c}}},\end{eqnarray} \tag{ 19 }$$

where ${k}_{{\rm{c}}}^{-1}$ represents the length scale of the spiky ring. This gives ${k}_{{\rm{c}}}\sim {\left(2\pi G{{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}/{D}^{2}\right)}^{1/3}$. Assuming that the ring mass does not change throughout the linear and nonlinear growth, one has a relation between the surface density and the wavenumber: ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}/{k}_{{\rm{c}}}={{\rm{\Sigma }}}_{{\rm{d}},0}/{k}_{0}$, where ${{\rm{\Sigma }}}_{{\rm{d}},0}$ and k₀ denote the unperturbed dust surface density and the wavenumber at which secular GI grows in the linear phase, respectively. Using this relation and Equation (19), one obtains

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}} & = & {\left(\displaystyle \frac{2\pi G{{\rm{\Sigma }}}_{{\rm{d}},0}}{{D}^{2}}\right)}^{\tfrac{1}{2}}{k}_{0}^{-\tfrac{3}{2}}\\ & \simeq & 9.3{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},0}/{{\rm{\Sigma }}}_{0}}{0.1}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{\alpha }{1\times {10}^{-3}}\right)}^{-1}{\left(\displaystyle \frac{Q}{4.5}\right)}^{-\tfrac{1}{2}}{\left(\displaystyle \frac{{k}_{0}H}{8}\right)}^{-\tfrac{3}{2}},\end{array}\end{eqnarray} \tag{ 20 }$$

where Σ₀ is the unperturbed gas surface density. This estimation is in good agreement with the resultant dust density of the ring at r ≃ 73.6 au in the Q4a10 run (see Figure 1). When the assumed α is smaller, numerical simulations with a constant softening length would underestimate ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}$.

Figure 7 shows the surface density and radial velocity profiles from the Q5a8 run. We can see the formation of dust rings and gaps at t = 3.2 × 10⁴ yr, but these structures decay as they move inward, which is in contrast to the Q4a10 run (Figure 1). The mean radial velocity of dust is in good agreement with the steady drift velocity as in the other runs. The relative velocity with respect to the drift motion is so low that significant dust concentration does not occur. As a result, the dust rings drift further without significant deceleration due to an increase of the dust-to-gas mass ratio.

Zoom In Zoom Out Reset image size

The decay of the substructures in the dust surface density profile can be clearly seen in Figure 8, where we show the deviation of ${{\rm{\Sigma }}}_{{\rm{d}}}$ from the initial profile ${{\rm{\Sigma }}}_{{\rm{d}}}(t=0)$ as a function of radius and time. Perturbations with a wavelength of ∼5 au grow while moving inward from r ≳ 100 au. The amplitudes decrease after they cross the radius r ≃ 60–70 au. We also observe the formation of those transient low-contrast rings in the Q4a20 and Q6a5 runs. In these cases, the unstable region is so small that the dust moves across the region without a significant increase in the surface density. We discuss this point in detail in the next section.

Zoom In Zoom Out Reset image size

In this subsection, we derive the dispersion relation of secular GI in the presence of drift motion in order to understand the physics captured in the present simulations. We also show what kind of mode becomes the secular GI in the dust-drifting system. Although we consider both dust and gas in the present simulations, we start with the linear analyses of a one-fluid system to readily understand the mode properties of secular GI. This simplified analysis seems valid even for comparison with two-component simulations because the gas disk does not significantly evolve in our simulations.⁵ The two-fluid analysis is described later in this section.

We use the following equations that govern the dust motion in the local coordinates (x, y) rotating around the central star with the angular velocity Ω₀ = Ω(r = R):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial t}+\displaystyle \frac{\partial }{\partial x}\left({{\rm{\Sigma }}}_{{\rm{d}}}{v}_{x}\right)=D\displaystyle \frac{{\partial }^{2}{{\rm{\Sigma }}}_{{\rm{d}}}}{\partial {x}^{2}},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {v}_{x}}{\partial t}+\left({v}_{x}-\displaystyle \frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial x}\right)\displaystyle \frac{\partial {v}_{x}}{\partial x} & = & 3{{\rm{\Omega }}}_{0}^{2}x+2{{\rm{\Omega }}}_{0}{v}_{y}\\ & - & \displaystyle \frac{{c}_{{\rm{d}}}^{2}}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial x}-\displaystyle \frac{\partial {\rm{\Phi }}}{\partial x}-\displaystyle \frac{{v}_{x}}{{t}_{\mathrm{stop}}}\\ & + & \displaystyle \frac{1}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial }{\partial x}\left({v}_{x}D\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial x}\right),\end{array}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {v}_{y}}{\partial t}+\left({v}_{x}-\displaystyle \frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial x}\right)\displaystyle \frac{\partial {v}_{y}}{\partial x} & = & -2{{\rm{\Omega }}}_{0}\left({v}_{x}-\displaystyle \frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial x}\right)\\ & - & \displaystyle \frac{{v}_{y}-{U}_{{\rm{g}},0}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 23 }$$

where ${U}_{{\rm{g}},0}\equiv -3{{\rm{\Omega }}}_{0}x/2-{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0}$ is the steady azimuthal velocity of the background gas disk, and ${\eta }_{0}^{{\prime} }\equiv \eta ^{\prime} (r=R)$. We include the dust diffusion so that the momentum is conserved in the absence of drag force (see Section 2.1 and Tominaga et al. 2019). An unperturbed state is the steady drift solution, ${v}_{x,0},{v}_{y,0}$, with uniform surface density ${{\rm{\Sigma }}}_{{\rm{d}}}={{\rm{\Sigma }}}_{{\rm{d}},0}$:

$$\begin{eqnarray}&&{v}_{x,0}=-\displaystyle \frac{2{t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}}{1+{\left({t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\right)}^{2}}{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{v}_{y,0}={U}_{{\rm{g}},0}+\displaystyle \frac{{\left({t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\right)}^{2}}{1+{\left({t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\right)}^{2}}{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0}.\end{eqnarray} \tag{ 25 }$$

Assuming that the perturbations $\delta {{\rm{\Sigma }}}_{{\rm{d}}},\delta {v}_{x},\delta {v}_{y}$ are proportional to exp[ikx + nt], one can derive the following linearized equations:

$$\begin{eqnarray}&&(n+{{ikv}}_{x,0})\delta {{\rm{\Sigma }}}_{{\rm{d}}}+{ik}{{\rm{\Sigma }}}_{{\rm{d}},0}\delta {v}_{x}=-{{Dk}}^{2}\delta {{\rm{\Sigma }}}_{{\rm{d}}},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}\begin{array}{rcl}(n+{{ikv}}_{x,0})\delta {v}_{x} & = & 2{{\rm{\Omega }}}_{0}\delta {v}_{y}-{{ikc}}_{{\rm{d}}}^{2}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}}-{ik}\delta {\rm{\Phi }}\\ & & -\,\displaystyle \frac{\delta {v}_{x}}{{t}_{\mathrm{stop}}}-{v}_{x,0}{{Dk}}^{2}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}},\end{array}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&(n+{{ikv}}_{x,0})\delta {v}_{y}=-\displaystyle \frac{{{\rm{\Omega }}}_{0}}{2}\left(\delta {v}_{x}-{ikD}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}}\right)-\displaystyle \frac{\delta {v}_{y}}{{t}_{\mathrm{stop}}}.\end{eqnarray} \tag{ 28 }$$

We also assume that the background gas is not perturbed and $\delta {\rm{\Phi }}=-2\pi G\delta {{\rm{\Sigma }}}_{{\rm{d}}}/k$. From the above linearized equations, we obtain the following dispersion relation for $\gamma \equiv n+{{ikv}}_{x,0}$:

$$\begin{eqnarray}&&\begin{array}{l}\left(\gamma +\displaystyle \frac{1}{{t}_{\mathrm{stop}}}\right)\left[{F}_{\mathrm{DW}}(\gamma ,k)+\gamma \left(\displaystyle \frac{1}{{t}_{\mathrm{stop}}}+{{Dk}}^{2}\right)\right.\\ \quad +\left.\,{{Dk}}^{2}\left(-{{ikv}}_{x,0}+\displaystyle \frac{1}{{t}_{\mathrm{stop}}}\right)\right]=\displaystyle \frac{{{\rm{\Omega }}}_{0}^{2}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{F}_{\mathrm{DW}}(\gamma ,k)\equiv {\gamma }^{2}+{{\rm{\Omega }}}_{0}^{2}+{c}_{{\rm{d}}}^{2}{k}^{2}-2\pi G{{\rm{\Sigma }}}_{{\rm{d}},0}k.\end{eqnarray} \tag{ 30 }$$

In the absence of dust diffusion (D = 0), Equation (29) is equivalent to Equation (22) of Youdin (2011) except for the softening term of the self-gravity due to the finite disk thickness. The softening term does not essentially change the mode properties, and thus we do not include it here. One of the solutions of Equation (29) with D = 0 corresponds to the one-component secular GI studied in Youdin (2005a, 2011), which is denoted by γ_SGI in the following, and we obtain $n=-{{ikv}}_{x,0}+{\gamma }_{\mathrm{SGI}}$. We note that Equation (29) with D = 0 is a cubic equation of γ with real coefficients, and $\mathrm{Im}[{\gamma }_{\mathrm{SGI}}]=0$.⁶ This is mathematically and physically expected because we can remove the drift motion via the Galilean transformation (see also Youdin 2005a). The phase velocity $\mathrm{Re}[n/(-{ik})]$ is the drift velocity ${v}_{x,0}$. This is consistent with our numerical simulations (see Figures 3 and 7). As mentioned in Youdin (2005a, 2011), secular GI originates from a neutral mode, referred to as a "static mode" in Tominaga et al. (2019), i.e., ${\gamma }_{\mathrm{SGI}}\to 0$ for ${t}_{\mathrm{stop}}{\rm{\Omega }}\to \infty$. Our dispersion relation also shows n → 0 for ${t}_{\mathrm{stop}}{\rm{\Omega }}\to \infty$ because ${v}_{x,0}$ also becomes zero.

The above statement does not qualitatively change in the presence of dust diffusion due to the weak gas turbulence. The dust diffusion just suppresses the growth of secular GI especially at short wavelengths. Although strong dust diffusion results in a nonzero Im[γ_SGI] at short wavelengths, the secular GI does not grow (Re[γ_SGI] < 0) with those wavelengths. Equation (29) is different from the dispersion relation derived by Youdin (2011) because the diffusion modeling is different. As mentioned in Tominaga et al. (2019), the diffusion modeling adopted in Youdin (2011) unphysically changes the dust angular momentum while ours does not. Hence, our dispersion relation (Equation (29)) describes the mode properties more precisely. For example, Youdin (2011) found mode coupling between the neutral mode and the GI mode (unstable density waves) due to the diffusion while Equation (29) does not show such a mode coupling (see also Tominaga et al. 2019).

We also perform two-fluid analyses with drift motion included. In contrast to the above one-fluid analyses without dust diffusion, we cannot remove the drift motion because the dust and gas have different drift speeds. We use the following equations for the gas, which are similar to those used to analyze the streaming instability (e.g., Youdin & Goodman 2005; Youdin & Johansen 2007; Jacquet et al. 2011; Chen & Lin 2020; Umurhan et al. 2020):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial t}+\displaystyle \frac{\partial }{\partial x}\left({\rm{\Sigma }}{u}_{x}\right)=0,\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {u}_{x}}{\partial t}+{u}_{x}\displaystyle \frac{\partial {u}_{x}}{\partial x} & = & 3{{\rm{\Omega }}}_{0}^{2}x+2{{\rm{\Omega }}}_{0}{u}_{y}+2{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0}^{2}-\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\rm{\Sigma }}}\displaystyle \frac{\partial {\rm{\Sigma }}}{\partial x}\\ & - & \displaystyle \frac{\partial {\rm{\Phi }}}{\partial x}+\displaystyle \frac{1}{{\rm{\Sigma }}}\displaystyle \frac{\partial }{\partial x}\left(\displaystyle \frac{4}{3}{\rm{\Sigma }}\nu \displaystyle \frac{\partial {u}_{x}}{\partial x}\right)+\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{\rm{\Sigma }}}\displaystyle \frac{{v}_{x}-{u}_{x}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial {u}_{y}}{\partial t}+{u}_{x}\displaystyle \frac{\partial {u}_{y}}{\partial x} & = & -2{{\rm{\Omega }}}_{0}{u}_{x}+\displaystyle \frac{1}{{\rm{\Sigma }}}\displaystyle \frac{\partial }{\partial x}\left({\rm{\Sigma }}\nu \displaystyle \frac{\partial {u}_{y}}{\partial x}\right)\\ & & +\,\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{\rm{\Sigma }}}\displaystyle \frac{{v}_{y}-{u}_{y}}{{t}_{\mathrm{stop}}}.\end{array}\end{eqnarray} \tag{ 33 }$$

The equations for the dust are the almost same as Equations (21)–(23) except for ${U}_{{\rm{g}},0}$ being replaced by u_y in Equation (23).

As in previous one-fluid analyses, the unperturbed state is the steady drift solution where the unperturbed surface densities, ${{\rm{\Sigma }}}_{{\rm{d}},0}$ and Σ₀, are spatially constant and the unperturbed velocities, ${v}_{x,0},{v}_{y,0},{u}_{x,0}$ and ${u}_{y,0}$, are the following (Nakagawa et al. 1986):

$$\begin{eqnarray}&&{v}_{x,0}=-\displaystyle \frac{2{t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}}{{(1+\epsilon )}^{2}+{\left({t}_{\mathrm{stop}}{\rm{\Omega }}\right)}^{2}}{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{v}_{y,0}=-\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{0}x-\left[1-\displaystyle \frac{{\left({t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\right)}^{2}}{{(1+\epsilon )}^{2}+{\left({t}_{\mathrm{stop}}{\rm{\Omega }}\right)}^{2}}\right]\displaystyle \frac{{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0}}{1+\epsilon },\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{u}_{x,0}=\displaystyle \frac{2{t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\epsilon }{{(1+\epsilon )}^{2}+{\left({t}_{\mathrm{stop}}{\rm{\Omega }}\right)}^{2}}{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0},\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{u}_{y,0}=-\displaystyle \frac{3}{2}{{\rm{\Omega }}}_{0}x-\left[1+\displaystyle \frac{{\left({t}_{\mathrm{stop}}{{\rm{\Omega }}}_{0}\right)}^{2}\epsilon }{{(1+\epsilon )}^{2}+{\left({t}_{\mathrm{stop}}{\rm{\Omega }}\right)}^{2}}\right]\displaystyle \frac{{\eta }_{0}^{{\prime} }R{{\rm{\Omega }}}_{0}}{1+\epsilon },\end{eqnarray} \tag{ 37 }$$

where $\epsilon \equiv {{\rm{\Sigma }}}_{{\rm{d}},0}/{{\rm{\Sigma }}}_{0}$. We fix ${t}_{\mathrm{stop}}$ to compare with the numerical simulations.⁷ The linearized continuity equation for the dust is the same as Equation (26). The other linearized equations are as follows:

$$\begin{eqnarray}&&(n+{{iku}}_{x,0})\delta {\rm{\Sigma }}+{ik}{{\rm{\Sigma }}}_{0}\delta {u}_{x}=0,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}\begin{array}{rcl}(n+{{iku}}_{x,0})\delta {u}_{x} & = & 2{{\rm{\Omega }}}_{0}\delta {u}_{y}-{{ikc}}_{{\rm{s}}}^{2}\displaystyle \frac{\delta {\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}-{ik}\delta {\rm{\Phi }}-\displaystyle \frac{4}{3}\nu {k}^{2}\delta {u}_{x}\\ & & +\,\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{0}}\displaystyle \frac{{v}_{x,0}-{u}_{x,0}}{{t}_{\mathrm{stop}}}-\epsilon \displaystyle \frac{\delta {\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\displaystyle \frac{{v}_{x,0}-{u}_{x,0}}{{t}_{\mathrm{stop}}}\\ & & +\,\epsilon \displaystyle \frac{\delta {v}_{x}-\delta {u}_{x}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}\begin{array}{rcl}(n+{{iku}}_{x,0})\delta {u}_{y} & = & -\displaystyle \frac{{{\rm{\Omega }}}_{0}}{2}\delta {u}_{x}-\displaystyle \frac{3}{2}{\rm{\Omega }}{ik}\nu \displaystyle \frac{\delta {\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}-\nu {k}^{2}\delta {u}_{y}\\ & & +\,\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{0}}\displaystyle \frac{{v}_{y,0}-{u}_{y,0}}{{t}_{\mathrm{stop}}}-\epsilon \displaystyle \frac{\delta {\rm{\Sigma }}}{{{\rm{\Sigma }}}_{0}}\displaystyle \frac{{v}_{y,0}-{u}_{y,0}}{{t}_{\mathrm{stop}}}\\ & & +\,\epsilon \displaystyle \frac{\delta {v}_{y}-\delta {u}_{y}}{{t}_{\mathrm{stop}}},\end{array}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}\begin{array}{rcl}(n+{{ikv}}_{x,0})\delta {v}_{x} & = & 2{{\rm{\Omega }}}_{0}\delta {v}_{y}-{{ikc}}_{{\rm{d}}}^{2}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}}-{ik}\delta {\rm{\Phi }}\\ & & -\,\displaystyle \frac{\delta {v}_{x}-\delta {u}_{x}}{{t}_{\mathrm{stop}}}-{v}_{x,0}{{Dk}}^{2}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}},\end{array}\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&(n+{{ikv}}_{x,0})\delta {v}_{y}=-\displaystyle \frac{{{\rm{\Omega }}}_{0}}{2}\left(\delta {v}_{x}-{ikD}\displaystyle \frac{\delta {{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}}\right)-\displaystyle \frac{\delta {v}_{y}-\delta {u}_{y}}{{t}_{\mathrm{stop}}},\end{eqnarray} \tag{ 42 }$$

where $\delta {\rm{\Phi }}=-2\pi G(\delta {\rm{\Sigma }}+\delta {{\rm{\Sigma }}}_{{\rm{d}}})/k$.

In Figure 9, we show the dispersion relation of secular GI for the Q4a10 run (cross marks). We set the unperturbed state using the physical values at r = 75 au. The qualitative properties are the same as in previous studies that ignored drift motion. The long-wavelength perturbations are stabilized by the backreaction from dust to gas. The dust diffusion limits the shortest unstable wavelength. Thus, secular GI is operational only at intermediate wavelengths.

Zoom In Zoom Out Reset image size

Drift motion does not modify the growth rate of secular GI significantly, which is similar to the above one-fluid analyses. To demonstrate this, we also show the growth rate without turbulent viscosity (filled circles) with the approximated one ${n}_{\mathrm{ap}}$ which does not include the turbulent viscosity or drift motion (gray line). The analytic expression of n_ap is given by Tominaga et al. (2019; see Equations (26)–(28) therein). Turbulent viscosity makes the growth rate larger (cross symbols) as also seen in the cases without drift motion (see Tominaga et al. 2019). The frequency $-\mathrm{Im}[n]$ is shown on the right panel of Figure 9. The frequency ${v}_{x,0}k$ (black line) expresses well the exact frequencies of secular GI (filled circles and cross symbols) as in the one-fluid analyses (see also Appendix B). At higher wavenumbers, the frequency deviates from ${v}_{x,0}k$ because of the dust diffusion term.

In the present analyses, we do not find an unstable mode corresponding to TVGI discovered by Tominaga et al. (2019). In the previous analyses of Tominaga et al. (2019), TVGI originates from a static mode in which the dust velocity perturbations are in phase with the gas velocity perturbations. This situation is not realized if the drift velocity is significant, which seems to be the reason for the stabilization of TVGI. We address and analyze in more detail the effects of the drift motion on TVGI in our future work.

The condition for the formation of thin dense rings is that radial drift motion is almost stopped because of the high dust-to-gas mass ratio in dust rings. Our numerical simulations indicate that the critical dust-to-gas mass ratio is about unity. In the Q4a10 run, the maximum dust-to-gas mass ratio is about unity, and the thin dense rings form. On the other hand, the maximum dust-to-gas ratio is smaller than unity in the Q5a8 run, and the transient low-contrast dust rings form (see Sections 3.1 and 3.2). We can understand the origin of these difference based on the linear stability analyses. Figure 10 shows the growth timescales of secular GI in the cases of the Q4a10 and Q5a8 runs. The wavelengths λ are 3, 5, and 10 au for each line. The minimum radii of the unstable regions for these wavelengths are larger for the Q5a8 run than those for the Q4a10 run. As shown in the previous section, the perturbations with λ ≃ 3–5 au grow and form the substructures in both cases. In the Q4a10 case, those perturbations can grow at r ≳ 50 au, which is consistent with Figure 1. The drift velocity decelerates as secular GI increases the dust-to-gas ratio, and the resultant rings tend to keep their radii once the dust-to-gas ratio approaches unity. As a result, the thin dense rings are not transient but sustain their dust-concentrating structure. On the other hand, in the Q5a8 case, the growth timescale at λ = 5 au starts increasing significantly at r ≃ 60–70 au. In other words, the growth of the perturbations decelerates, and the dust concentration gradually stops. Those perturbations keep drifting inward and eventually enter the stable region (r ≲ 58 au for λ = 5 au). In the stable region, the perturbations decay and decrease their amplitudes as shown in Figure 8.

Zoom In Zoom Out Reset image size

The maximum dust-to-gas mass ratio in the rings does not only depend on the background disk parameter Q and α but also depends on the amplitudes of the initial perturbations. We performed a simulation in which the parameters are the same as in the Q5a8 run, but the amplitude of the initial random perturbations is six times larger. We do not initially perturb a region of 10 au < r < 20 au in order to avoid sudden dust concentrations due to the initial large perturbations near the inner boundary. Figure 11 shows the time evolution of the surface density. Rings with larger amplitudes form before they enter the stable region, resulting in thin dense rings with ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}\simeq 1$ (see also Figure 12). The amplitude of the initial perturbations in a real disk is unknown, and we do not discuss its origin in the present article.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the present simulations, we fix the power-law index q of the surface densities. As shown in the previous subsection and Appendix B, q insignificantly affects the growth timescale of secular GI. In disks with a larger q, the inner region is more massive than smaller-q disks, and thus the unstable region shifts inward. If $\alpha ,{t}_{\mathrm{stop}}{\rm{\Omega }},Q$ and ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ are the same in the shifted region, the growth timescale becomes shorter as n/Ω is constant for the same parameters and ${n}^{-1}\propto {{\rm{\Omega }}}^{-1}\propto {r}^{3/2}$. The drift timescale also becomes shorter. The drift speed is faster for larger q. If one ignores the exponential cutoff term, $\eta ^{\prime} r{\rm{\Omega }}$ is independent of the radial distance. Thus, the drift timescale is proportional to r and becomes shorter in the shifted region. However, the r-dependence of the drift timescale is weaker than that of the growth timescale. Therefore, the secular GI can grow and more easily create thin dense rings in steeper disks. We note that the unstable wavelength becomes shorter as the gas scale height decreases inward. The widths of the resultant rings are thus smaller than what we present in Section 3.

The nonlinear growth of secular GI results in dusty rings with ${{\rm{\Sigma }}}_{{\rm{d}}}\gtrsim {\rm{\Sigma }}$. In such massive dusty rings, dust-growth acceleration and planetesimal formation via gravitational collapse are expected as described below. Because the rings are massive with a mass of $\gtrsim 10{M}_{\oplus }$ and self-gravitational, they will finally fragment in the azimuthal direction into solid bodies (see also Section 4.5). The timescale of this ring fragmentation can be measured with the freefall time ${t}_{\mathrm{ff}}$. The mean dust surface density in a dust ring ${\tilde{{\rm{\Sigma }}}}_{{\rm{d}}}$ is given by

$$\begin{eqnarray}&&{\tilde{{\rm{\Sigma }}}}_{{\rm{d}}}=\displaystyle \frac{{M}_{\mathrm{ring}}}{2\pi {R}_{\mathrm{ring}}{w}_{\mathrm{ring}}},\end{eqnarray} \tag{ 43 }$$

where ${M}_{{\rm{ring}}},{R}_{{\rm{ring}}},\,{\rm{and}}\,{w}_{{\rm{ring}}}$ are the mass, radius, and width of a dust ring, respectively. We assume the mean mass density ${\tilde{\rho }}_{{\rm{d}}}$ to be ${\tilde{{\rm{\Sigma }}}}_{{\rm{d}}}/\sqrt{2\pi }{H}_{{\rm{d}}}$. The freefall time is given by ${t}_{\mathrm{ff}}\,=\sqrt{3\pi /32G{\tilde{\rho }}_{{\rm{d}}}}$. Taking the ring properties from the Q4a10 run and assuming ${M}_{\mathrm{ring}}=49.7{M}_{\oplus }$, ${R}_{\mathrm{ring}}=73.6\,\mathrm{au}$, and ${w}_{\mathrm{ring}}=0.3\,\mathrm{au}$, we obtain ${t}_{\mathrm{ff}}\simeq 55\,\mathrm{yr}$, which is much shorter than one Keplerian period at $r=73.6\,\mathrm{au}$ (≃631 yr).

Both the total dust-to-gas surface density ratio ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$ and the dust-to-gas ratio ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ in the lower disk increase with the growth of the dusty rings. Because the collisional growth timescale of dust grains ${t}_{\mathrm{grow}}$ is inversely proportional to the total dust-to-gas surface density ratio ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$ in the Epstein regime (see Equation (A6)), the dust growth should proceed faster in dusty rings. When the height of the lower layer is much larger than the dust scale height, one has ${{\rm{\Sigma }}}_{{\rm{d}}}\simeq {{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}$. Because Σ insignificantly changes during the growth of the secular GI, ${{\rm{\Sigma }}}_{\mathrm{tot}}$ will be almost constant in time. Thus, the enhancement factor of ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{\mathrm{tot}}$ will be similar to that of ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$. In the case of the Q4a10 run, if the dust-to-gas ratio ${{\rm{\Sigma }}}_{{\rm{d}},\mathrm{tot}}/{{\rm{\Sigma }}}_{{\rm{g}},\mathrm{tot}}$ is initially 0.05 and increases as well as ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ by a factor of 9.3 as expected from the estimation of ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}/{{\rm{\Sigma }}}_{{\rm{d}},0}$ (see Equation (20)), Equation (A6) gives ${t}_{\mathrm{grow}}\simeq 199\,\mathrm{yr}$ at r = 73.6 au, which is about 3.1 times smaller than one orbital period although it is still larger than t_ff.

In the case of the Q5a8 run with six times larger initial perturbations, ${{\rm{\Sigma }}}_{{\rm{d}}}/{\rm{\Sigma }}$ becomes an order of magnitude larger than the initial value at r ≃ 80 au (Figure 11). The expected dust-growth timescale in the dust ring is about 209 yr, which is much shorter than the radial drift timescale. Thus, the dust grains in the ring will grow and decouple from the gas before the dust ring drifts into the stable region, which will also prevent the ring from being transient. Because we only include the turbulence-driven relative velocity in Equation (A6), the dust-growth timescale is overestimated. Thus, considering both ring fragmentation and dust growth is important, and in both cases, we can expect the subsequent planetesimal formation if significant dust concentration occurs via secular GI. More quantitative analyses require multidimensional simulations with dust growth, which will be our future work.

As mentioned in Section 3, secular GI results in significant substructures only in a dust disk. This property is in contrast to the planet-based scenario in which unseen Jupiter-mass planets carve gaps in both dust and gas disks (e.g., Gonzalez et al. 2015; Kanagawa et al. 2015; Zhang et al. 2018). Therefore, observational determination of a midplane gas density profile would enable us to discriminate the ring formation mechanisms that require a gap-like profile in the gaseous component (i.e., the hidden high-mass planet hypothesis) and show a relatively smooth profile in gas (i.e., the present mechanism).

We should note that if the nonlinear growth of secular GI results in massive rocky objects, those objects would create some substructures in the gas disk. If the gas gap is as large as the dust gap, it might be difficult to clearly distinguish between the secular-GI-based mechanism and the planet-based mechanism. Even in this case, however, secular GI might explain the origin of the hypothetical planets in the gaps.

Even if a gas gap is not observed around a dust gap, it is possible that a low-mass planet carves the dust gap. A low-mass planet first creates prominent structures in a dust disk and takes a long time to carve a gas gap because the dust scale height is smaller than the gas scale height (e.g., Yang & Zhu 2020). Thus, there is a possibility that observations only identify the resolvable dust gaps but cannot resolve low-contrast gas gaps. Thus, the low-mass-planet-based scenario and the secular-GI-based scenario can be degenerate. Nevertheless, those scenarios depend on the background density and pressure profiles, the dust drift, and sizes in different ways. Therefore, it is important to accurately measure the gas density profile near the midplane as well as the dust surface density.

The separations and widths of rings/gaps depend on which wavelength perturbations are large and grow fast via secular GI. In our simulations with the initial random perturbations, the wavelength is comparable to the gas scale height at each radius, which is about the most unstable wavelength. The rings migrate inward while keeping separations of the adjacent rings almost constant (Figure 8), because the radial variation of the drift velocity is small within a spatial scale of the unstable wavelength. This is another feature of the ring-forming process via secular GI.

Thin dense rings forming via secular GI may be observed as marginally optically thin rings although the solid surface density in those rings can be much higher. Stammler et al. (2019) claimed that planetesimal formation via the streaming instability in a pressure bump can explain the marginally optically thick rings (see also Dullemond et al. 2018; Huang et al. 2018). The streaming instability in a bump converts dust grains into planetesimals and stalls itself as a result of the dust depletion. Because the resultant planetesimals are not observed at submillimeter emission, this self-regulating process can limit the optical depth in the dust ring. Similar effects are expected in the ring formation via secular GI.

As discussed above, dust growth to larger solid bodies can occur in the resultant dusty rings. If those solid bodies are large enough, those rings would be dark at millimeter wavelengths, and multiple spiky rings and adjacent gaps would be observed as a wide gap structure. The collisional fragmentation of the resultant larger bodies would supply small grains. The remaining and resupplied small grains take part in the growth of the secular GI and contribute to submillimeter emission observed by ALMA. To examine what kinds of substructures are expected in the intensity profile, we have to explicitly include dust growth and fragmentation in our simulations. Because self-gravitational fragmentation of the rings could occur along with dust growth at comparable timescales, non-axisymmetic analyses and simulations are important to explore observational signatures. Those will be the scope of our future studies.

In the presence of dust drift motion, streaming instability can operate in a disk (Youdin & Goodman 2005; Youdin & Johansen 2007; Jacquet et al. 2011). Nonlinear properties of the streaming instability have been extensively studied based on local-shearing-box simulations (e.g., Johansen & Youdin 2007; Johansen et al. 2007, 2009; Bai & Stone 2010a, 2010b, 2010c; Yang & Johansen 2014; Simon et al. 2016; Yang et al. 2017; Schreiber & Klahr 2018). These numerical studies showed that the streaming instability creates dust clumps and triggers GI if those clumps are massive enough. This self-gravitational clump collapse results in planetesimals. Thus, self-gravity and GI are required to trigger planetesimal formation even in the scenario based on the streaming instability.

In the presence of self-gravity, secular GI can also operate as we show in Section 4.1 and thus would be important for forming dust clumps and planetesimals. Abod et al. (2019) studied streaming instability based on local simulations with self-gravity and showed dust clumping even in the absence of a global pressure gradient. The clumping they found is caused by classical GI mediated by the dust–gas friction although they referred to the process as secular GI. Secular GI is essentially different from friction-mediated GI (Youdin & Goodman 2005; Tominaga et al. 2019), and thus the relation between secular GI and streaming instability is still unclear.

Our linear analyses and nonlinear simulations show that secular GI grows at wavelengths ∼H for α ∼ 10⁻³ and ${t}_{\mathrm{stop}}{\rm{\Omega }}=0.6$. Streaming instability grows at shorter wavelengths ≪H (e.g., see Youdin & Goodman 2005) and creates smaller structures. Therefore, secular GI and streaming instability would grow independently at each scale at least in the linear growth phase.

In the presence of dust diffusion as in our study, the growth of streaming instability would be suppressed (Chen & Lin 2020; Umurhan et al. 2020), and dust diffusion limits spatial scales at which subsequent gravitational collapse can proceed (Gerbig et al. 2020; Gole et al. 2020; Klahr & Schreiber 2020). The linear or nonlinear growth of secular GI will potentially provide a dust-rich environment that is necessary for streaming instability to grow against the diffusion. Therefore, the combination of dust concentration at large scales (∼H) due to secular GI and further concentration at small scales (≪H) due to streaming instability is one possible path for forming planetesimals.

Depending on the driving mechanism of gas turbulence, radial diffusion is inefficient in contrast to vertical diffusion (Yang et al. 2018; Schäfer et al. 2020). In this case, both secular GI and streaming instability can relatively easily develop against the diffusion. To explore more detailed mode properties, nonlinear coupling, and evolutionary paths toward planetesimal formation, multidimensional linear analyses and simulations with the radial extent of ∼H are necessary. These will also be the scope of our future work.

In our simulations, the dust diffusivity D is given by a constant parameter α. However, it will decrease when a large dust-to-gas mass ratio is realized in the rings formed by secular GI (see Figure 1). For example, Schreiber & Klahr (2018) reported that the turbulent diffusion caused by the streaming instability weakens when the dust-to-gas mass ratio is larger than unity. If the turbulent diffusion in dust rings is driven by streaming instability, the saturated surface density of the thin dense rings will be larger than that obtained from our simulations or estimated in Equation (20). It is also possible that inefficient diffusion cannot support the thin dense ring, and the ring just radially collapses. In either case, the planetesimal formation timescale via gravitational collapse or dust growth will be shorter than that estimated in the previous section.

We can estimate the saturated surface density at dense thin rings as follows. We assume that the diffusion coefficient of dust is given by

$$\begin{eqnarray}&&D={\alpha }_{1}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\rm{\Omega }}}{\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}}{{{\rm{\Sigma }}}_{0}}\right)}^{-\beta },\end{eqnarray} \tag{ 44 }$$

where α₁ is a dimensionless diffusion parameter when ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}/{{\rm{\Sigma }}}_{0}=1$. In the same way we derive Equation (20), we can evaluate the dust surface density when it achieves saturation as follows:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}}{{{\rm{\Sigma }}}_{{\rm{d}},0}} & = & f(\beta ){\left(\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}},0}/{{\rm{\Sigma }}}_{0}}{0.1}\right)}^{\tfrac{1+2\beta }{2-2\beta }}{\left(\displaystyle \frac{{\alpha }_{1}}{1\times {10}^{-3}}\right)}^{\tfrac{1}{\beta -1}}\\ & & \times \,{\left(\displaystyle \frac{Q}{4.5}\right)}^{\tfrac{1}{2\beta -2}}{\left(\displaystyle \frac{{k}_{0}H}{8}\right)}^{\tfrac{3}{2\beta -2}},\end{array}\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&f{(\beta )\simeq (8.7\times {10}^{1-2\beta })}^{\tfrac{1}{2-2\beta }}.\end{eqnarray} \tag{ 46 }$$

The equation has no solution and f(β) diverges if β = 1. This indicates that dust diffusion cannot sustain the gravitational collapse of the dust ring if the dust-to-gas mass ratio dependence of D is steeper than the case of β = 1. Although the surface density has a solution for β > 1, it is not achieved during the growth of the perturbation. In the case of β > 1, the diffusion timescale is smaller than the collapse timescale when the dust ring surface density is smaller than ${{\rm{\Sigma }}}_{{\rm{d}},{\rm{f}}}$.

The velocity dispersion of dust grains ${c}_{{\rm{d}}}$ is expected to decrease as the dust-to-gas ratio increases beyond unity. If the decrease of ${c}_{{\rm{d}}}$ is realized, the Coriolis force will be a unique repulsive process against self-gravitational collapse. The characteristic length scale λ_crit of such a system is

$$\begin{eqnarray}&&{\lambda }_{\mathrm{crit}}\equiv \displaystyle \frac{4{\pi }^{2}G{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Omega }}}^{2}}.\end{eqnarray} \tag{ 47 }$$

Dust clumps smaller than λ_crit will self-gravitationally collapse. Taking the ring properties from the Q4a10 run and adopting ${{\rm{\Sigma }}}_{{\rm{d}}}$ = 9.5 g cm⁻² at r = 73.6 au, one obtains λ_crit ≃ 17.7 au, which is much larger than the resultant ring width (see Figure 1). The circumference of a ring with the radius of 73.6 au is about 462 au, which is six times longer than λ_crit. Thus, non-axisymmetric modes with an azimuthal wavenumber m larger than 6 might grow and cause azimuthal fragmentation. Because higher m modes have larger growth rates, the sizes of the resultant dust clumps can be much smaller. If the dust ring in the Q4a10 run whose mass is M_ring = 49.7 ${M}_{\oplus }$ and radius is r = ≃73.6 au fragments, the mass of the dust clumps M_clump resulting from azimuthal fragmentation is

$$\begin{eqnarray}&&{M}_{\mathrm{clump}}=\displaystyle \frac{{M}_{\mathrm{ring}}}{m}\simeq 8.2{M}_{\oplus }{\left(\displaystyle \frac{m}{6}\right)}^{-1}.\end{eqnarray} \tag{ 48 }$$

Because non-axisymmetric modes with $m\lt 6$ are expected to be stable, the above clump mass for m = 6 gives an upper limit. For further studies, non-axisymmetric simulations are necessary.