Nonlinear Outcome of Coagulation Instability in Protoplanetary Disks. I. First Numerical Study of Accelerated Dust Growth and Dust Concentration at Outer Radii

In this work, we consider dust evolution in a steady axisymmetric gas disk around a star of 1 M_⊙. The equation of continuity for dust in cylindrical coordinates (r, ϕ, z) is given by

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

$$\begin{eqnarray}&&{v}_{r}\equiv \left\langle {v}_{r}\right\rangle -\frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r},\end{eqnarray} \tag{ 2 }$$

where Σ_d is the dust surface density and D is the diffusion coefficient. The first term in Equation (2) denotes a mean drift velocity of dust grains. Assuming that the gas rotates at sub-Keplerian velocity (1 − η)v_K, where ${v}_{{\rm{K}}}=\sqrt{{{GM}}_{\odot }/r}$ is the Keplerian velocity and G is the gravitational constant, one obtains the radial terminal velocity in the test particle limit as follows (e.g., Adachi et al. 1976; Weidenschilling 1977; Nakagawa et al. 1986):

$$\begin{eqnarray}&&\left\langle {v}_{r}\right\rangle =-\frac{2{\tau }_{{\rm{s}}}}{1+{\tau }_{{\rm{s}}}^{2}}\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 3 }$$

where τ_s ≡ t_stopΩ is the stopping time of dust grains, t_stop, normalized by the Keplerian angular velocity ${\rm{\Omega }}\,\equiv {v}_{{\rm{K}}}/r\,=\sqrt{{{GM}}_{\odot }/{r}^{3}}$. The value of η for a given gas density ρ _g and temperature T is determined by

$$\begin{eqnarray}&&\eta \equiv -\frac{1}{2{\rho }_{g}r{{\rm{\Omega }}}^{2}}\frac{\partial {c}_{{\rm{s}}}^{2}{\rho }_{g}}{\partial r},\end{eqnarray} \tag{ 4 }$$

where ${c}_{{\rm{s}}}=\sqrt{{k}_{{\rm{B}}}T/\mu {m}_{{\rm{H}}}}$ is the sound speed, and k_B and m_H are the Boltzmann constant and the hydrogen mass, respectively. We adopt the mean molecular weight μ of 2.34. We use the analytical formula of the radial diffusion coefficient derived by Youdin & Lithwick (2007):

$$\begin{eqnarray}&&D=\frac{1+{\tau }_{{\rm{s}}}+4{\tau }_{{\rm{s}}}^{2}}{{\left(1+{\tau }_{{\rm{s}}}^{2}\right)}^{2}}\alpha {c}_{{\rm{s}}}H.\end{eqnarray} \tag{ 5 }$$

We calculate the radial dust diffusion using the gradient ∂Σ_d/∂r (e.g., Cuzzi et al. 1993) instead of ${{\rm{\Sigma }}}_{{\rm{g}}}\partial \left({{\rm{\Sigma }}}_{{\rm{d}}}/{{\rm{\Sigma }}}_{{\rm{g}}}\right)/\partial r$, where Σ_g is the gas surface density (e.g., Dubrulle et al. 1995). We note that coagulation instability is hardly affected by the choice of the gradients (see Section 4.1 of Tominaga et al. 2021). The difference in the diffusion flux,

$$\begin{eqnarray}&&\frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r}-\frac{D{{\rm{\Sigma }}}_{{\rm{g}}}}{{{\rm{\Sigma }}}_{{\rm{d}}}}\frac{\partial }{\partial r}\left(\frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\right)=\frac{D}{{{\rm{\Sigma }}}_{{\rm{g}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}}{\partial r},\end{eqnarray} \tag{ 6 }$$

is small compared to the mean drift velocity in weakly turbulent disks considered in this work. For $D\sim \alpha {c}_{{\rm{s}}}^{2}/{\rm{\Omega }}$ and $\eta \sim {({c}_{{\rm{s}}}/{v}_{{\rm{K}}})}^{2}$, one obtains

$$\begin{eqnarray}&&\left|\frac{D}{{{\rm{\Sigma }}}_{{\rm{g}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}}{\partial r}\times {\left\langle {v}_{r}\right\rangle }^{-1}\right|\sim \frac{\alpha }{{\tau }_{{\rm{s}}}}.\end{eqnarray} \tag{ 7 }$$

The right-hand side is much smaller than unity in weakly turbulent gas disks (α ∼ 10⁻⁴) because the drifting dust has τ_s ∼ 0.1–1. As noted by Tominaga et al. (2021), coagulation instability creates Σ_d perturbations at a radial scale of∼ H ≡ c_s/Ω, where H is the gas scale height. In such a case, ${{\rm{\Sigma }}}_{{\rm{g}}}^{-1}\partial {{\rm{\Sigma }}}_{{\rm{g}}}/\partial r$ is smaller than ${{\rm{\Sigma }}}_{{\rm{d}}}^{-1}\partial {{\rm{\Sigma }}}_{{\rm{d}}}/\partial r$ by a factor of H/r. Therefore, the difference in the diffusion flux is minor (see also Laibe et al. 2020).

We describe dust size evolution via coagulation adopting the moment method (Estrada & Cuzzi 2008; Ormel & Spaans 2008; Sato et al. 2016). The evolutionary equation of mass-dominating dust sizes is derived in Sato et al. (2016):

$$\begin{eqnarray}&&\frac{\partial {m}_{{\rm{p}}}}{\partial t}+{v}_{r}\frac{\partial {m}_{{\rm{p}}}}{\partial r}=\frac{2\sqrt{\pi }{a}^{2}{\rm{\Delta }}{v}_{\mathrm{pp}}}{{H}_{{\rm{d}}}}{{\rm{\Sigma }}}_{{\rm{d}}},\end{eqnarray} \tag{ 8 }$$

where m_p and a are the mass and the radius of the representative dust, respectively, Δv_pp is the collision velocity, and H_d is the dust scale height. Equation (8) can be derived from the Smoluchowski equation for column number density for dust, that is, a vertically integrated equation. The dust scale height is a function of τ_s and α (Youdin & Lithwick 2007):

$$\begin{eqnarray}&&{H}_{{\rm{d}}}({\tau }_{{\rm{s}}})=H{\left(1+\frac{{\tau }_{{\rm{s}}}}{\alpha }\left(\frac{1+2{\tau }_{{\rm{s}}}}{1+{\tau }_{{\rm{s}}}}\right)\right)}^{-1/2}.\end{eqnarray} \tag{ 9 }$$

We use turbulence-driven collision velocity Δv_pp = Δv_t given by Ormel & Cuzzi (2007) in Section 3 for comparison with the previous linear analyses. Especially, we assume ${\rm{\Delta }}{v}_{\mathrm{pp}}=\sqrt{C\alpha {\tau }_{{\rm{s}}}}{c}_{{\rm{s}}}$ with a numerical constant C in Section 3.1 for simplicity, which is valid for intermediate dust sizes (see Section 3.4.2 in Ormel & Cuzzi 2007). In radially global simulations in Section 4, we include Brownian motion and differential drift and settling velocities in Δv_pp:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{{\rm{p}}{\rm{p}}}=\sqrt{{\left({\rm{\Delta }}{v}_{{\rm{t}}}\right)}^{2}+{\left({\rm{\Delta }}{v}_{{\rm{B}}}\right)}^{2}+{\left({\rm{\Delta }}{v}_{r}\right)}^{2}+{\left({\rm{\Delta }}{v}_{\phi }\right)}^{2}+{\left({\rm{\Delta }}{v}_{z}\right)}^{2}},\end{eqnarray} \tag{ 10 }$$

where Δv_B is the collision velocity due to Brownian motion and Δv_r , Δv_ϕ , and Δv_z are the relative velocities between dust grains of τ_s = τ_s,1 and τ_s,2. In the test particle limit, the azimuthal and vertical drift velocities of one dust particle are as follows (e.g., Nakagawa et al. 1986):

$$\begin{eqnarray}&&{v}_{\phi }({\tau }_{{\rm{s}}})={v}_{{\rm{K}}}-\frac{1}{1+{\tau }_{{\rm{s}}}^{2}}\eta {v}_{{\rm{K}}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{v}_{z}({\tau }_{{\rm{s}}})=-\frac{{\tau }_{{\rm{s}}}}{1+{\tau }_{{\rm{s}}}}z{\rm{\Omega }}.\end{eqnarray} \tag{ 12 }$$

Because our model is based on the vertically averaged equations, we adopt $z=\sqrt{2/\pi }{H}_{{\rm{d}},\mathrm{1,2}}$ to calculate the vertically averaged Δv_z , where ${H}_{{\rm{d}},\mathrm{1,2}}\equiv {({H}_{{\rm{d}}}{\left({\tau }_{{\rm{s}},1}\right)}^{-2}+{H}_{{\rm{d}}}{({\tau }_{{\rm{s}},2})}^{-2})}^{-1/2}$ (see Taki et al. 2021).

We focus on only evolution of Σ_d and m_p and ignore the higher moments. This simplification is valid until the onset of the runaway growth because the evolution of size distribution is well described by the evolution of the peak mass m_p (Kobayashi et al. 2016). Sato et al. (2016) formulated a closure by comparing full-size simulations and single-moment simulations. They found that using τ_s,1 = 0.5τ_s,2 for the turbulent and differential collision velocities well reproduces the m_p evolution from the full-size simulations. We thus adopt their formalism (for the details, see Sato et al. 2016).

We use the Epstein and Stokes laws to calculate t_stop:

$$\begin{eqnarray}&&\displaystyle {t}_{\mathrm{stop}}=\left\{\begin{array}{l}\sqrt{\frac{\pi }{8}}\frac{{\rho }_{\mathrm{int}}a}{{\rho }_{g}{c}_{{\rm{s}}}},\left(\frac{a}{{l}_{\mathrm{mfp}}}\leqslant \frac{9}{4}\right)\\ \sqrt{\frac{\pi }{8}}\frac{4{\rho }_{\mathrm{int}}{a}^{2}}{9{\rho }_{g}{c}_{{\rm{s}}}{l}_{\mathrm{mfp}}},\left(\frac{a}{{l}_{\mathrm{mfp}}}\gt \frac{9}{4}\right)\end{array}\right.,\end{eqnarray} \tag{ 13 }$$

where ρ_int = 1.4 g cm⁻³ is an internal mass density of dust grains (see Appendix D) and l_mfp is the mean free path of gas. In this work, we use the midplane gas density to calculate t_stop since most of the dust grains reside around the midplane (H_d < H) when the radial drift motion becomes significant. Assuming the Gaussian gas density profile in the vertical direction and ${\rho }_{g}={{\rm{\Sigma }}}_{{\rm{g}}}/\sqrt{2\pi }H$, we obtain

$$\begin{eqnarray}&&{\tau }_{{\rm{s}}}=\frac{\pi }{2}\frac{{\rho }_{\mathrm{int}}a}{{{\rm{\Sigma }}}_{{\rm{g}}}}\max \left(1,\frac{4a}{9{l}_{\mathrm{mfp}}}\right),\end{eqnarray} \tag{ 14 }$$

(see also Sato et al. 2016). The Epstein law is applicable in the present work because we focus on r ≳ 5 au where the gas density is low and l_mfp is large enough.

We use the Lagrangian-cell method developed in Tominaga et al. (2018) for both local and global simulations in the following sections. The advantage of using the Lagrangian scheme is the fact that the scheme is free from numerical diffusion associated with advection, which enables us to accurately describe the linear growth of coagulation instability. Using Equations (2) and (3), we update the position of the ith cell boundary, r_d,i+1/2 for i = 1, 2,..., N − 1, where N denotes the number of cells. Dust densities Σ_d,i are assigned at cell centers r = r_d,i while dust sizes a_i+1/2 are assigned at cell boundaries. Assuming cell masses to be constant in time, we update dust densities with time-variable cell widths. In the radially global simulations (Section 4 and Paper II), the cell positions correspond to the dust positions in the laboratory frame. On the other hand, we regard the cell positions as the coordinates in the unperturbed-dust-rest frame in the local simulations, which is explained in detail in the next section.

We adopt the operator-splitting method for time integration (e.g., Inoue & Inutsuka 2008). We first update the cell boundaries' positions r_d,i+1/2 with $\left\langle {v}_{r}({r}_{{\rm{d}},i})\right\rangle$ and Σ_d,i by a half time step, Δt/2. We again update r_d,i+1/2 and Σ_d,i with the diffusion flow velocity −DΣ_d,i ⁻¹∂Σ_d,i /∂r by a half time step (see Tominaga et al. 2020). For time integration, we use the second-order Runge–Kutta scheme. ³ We then update the dust sizes by one time step using Equation (8). Finally, we update r_d,i+1/2 and Σ_d,i with the mean radial velocity and the diffusion flow velocity in reverse order of the above time stepping to ensure the second-order accuracy (Inoue & Inutsuka 2008).

The time step Δt is determined by

$$\begin{eqnarray}&&{\rm{\Delta }}t=\min ({\rm{\Delta }}{t}_{{\rm{d}}},{\rm{\Delta }}{t}_{\mathrm{diff}},{\rm{\Delta }}{t}_{\mathrm{coag}}),\end{eqnarray} \tag{ 15 }$$

where Δt_d, Δt_diff and Δt_coag are the time steps limited by the Courant−Friedrich−Levy condition for the dust drift, dust diffusion, and coagulation. When the time step becomes too short because of the diffusion, we adopt the super-time-stepping scheme to accelerate the time integration (Alexiades et al. 1996; Meyer et al. 2012, 2014).

Our previous study (Tominaga et al. 2021) derived a dispersion relation of coagulation instability based on one-fluid and two-fluid equations. The previous study shows that coagulation instability is essentially a one-fluid mode (see Section 3 of Tominaga et al. 2021). We thus review the one-fluid dispersion relation below. We suggest referring to Tominaga et al. (2021) for detailed derivation and more comprehensive discussions.

The previous analyses assume the turbulence-driven coagulation and ${\rm{\Delta }}{v}_{\mathrm{pp}}=\sqrt{C\alpha {\tau }_{{\rm{s}}}}{c}_{{\rm{s}}}$ with C ≃ 2.3. Tominaga et al. (2021) reduced Equation (8) into an evolutionary equation for τ_s based on the Epstein law. In the local shearing sheet coordinates (x, y) = (r − R, R(ϕ − Ωt)) around a radius R (Goldreich & Lynden-Bell 1965b), the equation is

$$\begin{eqnarray}&&\frac{\partial {\tau }_{{\rm{s}}}}{\partial t}+{v}_{x}\frac{\partial {\tau }_{{\rm{s}}}}{\partial x}=\frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\frac{{\tau }_{{\rm{s}}}}{3{t}_{0}}+\frac{{\tau }_{{\rm{s}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}{u}_{x}}{\partial x}-\frac{{\tau }_{{\rm{s}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}{v}_{x}\frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}}{\partial x},\end{eqnarray} \tag{ 16 }$$

where ${t}_{0}\equiv (4/3\sqrt{C}){{\rm{\Omega }}}^{-1}$ and an axisymmetric disk is assumed. The first term on the right-hand side is a coagulation term that increases τ_s in time. The second term represents the change in τ_s caused by the enhanced Σ_g through gas compression. This term can be neglected in the one-fluid analyses with the steady gas disk. The third term represents a process where τ_s decreases as dust moves into an inner gas-dense region if dust size a is fixed. Tominaga et al. (2021) found that the third term reduces growth rates of the instability by a factor of a few. In this section, we assume uniform Σ_g = Σ_g,0 for simplicity and use the following equation:

$$\begin{eqnarray}&&\frac{\partial {\tau }_{{\rm{s}}}}{\partial t}+{v}_{x}\frac{\partial {\tau }_{{\rm{s}}}}{\partial x}=\frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\frac{{\tau }_{{\rm{s}}}}{3{t}_{0}};\end{eqnarray} \tag{ 17 }$$

note that the neglected third term is consistently taken into account in radially global simulations in the next section.

The unperturbed dust surface density is assumed to be uniform: Σ_d = Σ_d,0. Besides, only in the linear analysis is the dimensionless stopping time in the unperturbed state uniform and constant in time: τ_s = τ_s,0. We note that the background dust size and thus τ_s,0 should increase in time at the background coagulation rate Σ_d,0/3Σ_g,0 t₀ (see Equation (17)). However, coagulation instability can develop faster at short wavelengths than the background coagulation (see Section 2 of Tominaga et al. 2021). Therefore, the above assumption to derive the growth rates is applicable as the lowest-order approximation. Note that we consistently treat the background evolution in dust size, and thus τ_s, in the local and global simulations.

Using Equations (1), (2), and (17), one can derive a complex growth rate n as a function of wavenumber k:

$$\begin{eqnarray}\begin{array}{rcl}n & = & -{ik}\left\langle {v}_{x,0}\right\rangle +\frac{1}{2}\left(\frac{{\varepsilon }_{0}}{3{t}_{0}}-{{Dk}}^{2}\right)\\ & & +\frac{1}{2}\sqrt{{\left(\frac{{\varepsilon }_{0}}{3{t}_{0}}+{{Dk}}^{2}\right)}^{2}-\frac{4{ik}{\varepsilon }_{0}}{3{t}_{0}}\left\langle {v}_{x,0}\right\rangle \frac{1-{\tau }_{{\rm{s}},0}^{2}}{1+{\tau }_{{\rm{s}},0}^{2}}},\end{array}\end{eqnarray} \tag{ 18 }$$

where ε₀ = Σ_d,0/Σ_g,0 and $\left\langle {v}_{x,0}\right\rangle =-2{\tau }_{{\rm{s}},0}\eta R{\rm{\Omega }}/(1+{\tau }_{{\rm{s}},0}^{2})$. Equation (18) shows $\mathrm{Re}[n]\gt 0$ for intermediate wavenumbers.

The physics behind coagulation instability can be clearly understood in the absence of diffusion. Setting D = 0 in Equation (18), we obtain the growth rate at sufficiently high wavenumbers as follows:

$$\begin{eqnarray}&&\mathrm{Re}[n]\simeq \sqrt{\displaystyle \frac{1}{2}\displaystyle \frac{{\varepsilon }_{0}}{3{t}_{0}}k| \left\langle {v}_{x,0}\right\rangle | },\end{eqnarray} \tag{ 19 }$$

where we set $1\pm {\tau }_{{\rm{s}},0}^{2}\simeq 1$ assuming τ_s,0 < 1. Equation (19) shows that the growth rate is determined by a geometric mean of the coagulation rate ε₀/3t₀ and $k| \left\langle {v}_{x,0}\right\rangle |$. The latter represents a rate of traffic jam induced by τ_s-perturbations at a scale of k⁻¹. Therefore, coagulation instability can be regarded as a positive feedback process driven by a combination of coagulation and traffic jam (see Figure 1 and Section 2 of Tominaga et al. 2021).

Zoom In Zoom Out Reset image size

The dust diffusion neglected in Equation (19) damps short-wavelength perturbations. As a result, coagulation instability develops the most efficiently at intermediate wavelengths (see Section 4.1 of Tominaga et al. 2021).

Next, we explain setups and the initial condition of local simulations with a radial domain around r = R. In local simulations in this section, we update the cell positions according to the following perturbed drift velocity:

$$\begin{eqnarray}&&\delta {v}_{x}\equiv \delta \left\langle {v}_{x}\right\rangle -\frac{D}{{{\rm{\Sigma }}}_{{\rm{d}}}}\frac{\partial {{\rm{\Sigma }}}_{{\rm{d}}}}{\partial r},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\delta \left\langle {v}_{x}\right\rangle \equiv -\frac{2{\tau }_{{\rm{s}}}(x,t)}{1+{\tau }_{{\rm{s}}}{\left(x,t\right)}^{2}}\eta R{\rm{\Omega }}+\frac{2{\tau }_{{\rm{s}},0}(t)}{1+{\tau }_{{\rm{s}},0}{\left(t\right)}^{2}}\eta R{\rm{\Omega }},\end{eqnarray} \tag{ 21 }$$

where we calculate the dimensionless stopping time as τ_s = t_stopΩ(R). This frame is the unperturbed-dust-rest frame, X, and thus we can clearly observe propagation of waves relative to the dust. The radial coordinates x (dust positions in the laboratory frame) and X have the following relation:

$$\begin{eqnarray}&&X=x+{\int }_{0}^{t}\frac{2{\tau }_{{\rm{s}},0}(t^{\prime} )}{1+{\tau }_{{\rm{s}},0}{\left(t^{\prime} \right)}^{2}}\eta R{\rm{\Omega }}{dt}^{\prime} .\end{eqnarray} \tag{ 22 }$$

Note that we let dust evolve in time according to Equation (8), and thus the background stopping time τ_s,0(t) increases in time. For comparison with the linear analyses, we also assume ${\rm{\Delta }}{v}_{\mathrm{pp}}=\sqrt{C\alpha {\tau }_{{\rm{s}}}}{c}_{{\rm{s}}}$ in this subsection.

There are four parameters: (1) η RΩ, (2) background surface densities Σ_{g ,0}, Σ_{d ,0}, (3) initial dust size or dimensionless stopping time τ_s,0, and (4) turbulence strength α. We adopt η RΩ ≃ 54 m s⁻¹. This choice is based on the minimum-mass solar nebula (MMSN) disk model (Hayashi 1981). The gas surface density is uniformly set as ${{\rm{\Sigma }}}_{{\rm{g}},0}={{\rm{\Sigma }}}_{{\rm{g}},\mathrm{MMSN}}(20\,\mathrm{au})$, where ${{\rm{\Sigma }}}_{{\rm{g}},\mathrm{MMSN}}{(r)=53.75\,{\rm{g}}\,{\mathrm{cm}}^{-2}(r/10\,\mathrm{au})}^{-1.5}$ is the gas surface density in the MMSN disk model. We assume a uniform dust-to-gas ratio of 10⁻³ in the local simulations. We initially adopt τ_s = t_stopΩ(R) = 10⁻² throughout the domain and set the corresponding dust sizes using Equation (13). We adopt α = 10⁻⁴ and 10⁻⁵ in the local simulations. The results for α = 10⁻⁵ are presented in Appendix A.

Figure 1 schematically shows the initial condition explained below. We simulate the time evolution of a sinusoidal perturbation. The radial width of the domain is set to be one wavelength λ and divided into N = 128 cells. We adopt the periodic boundary condition: the boundaries of the domain move so that Σ_d becomes periodic. We stop simulations once the cell boundary of the first cell, r_d,3/2, crosses the initial position of the inner domain boundary r_in.

Since our code is based on the Lagrangian scheme, we can set up density perturbations by displacing the cell locations (see the bottom part of Figure 1). We adopt the following displacement δ ξ:

$$\begin{eqnarray}&&\delta \xi \equiv {10}^{-2}\displaystyle \frac{\lambda }{2\pi }\sin \left(\displaystyle \frac{2\pi }{\lambda }(r-{r}_{\mathrm{in}})\right).\end{eqnarray} \tag{ 23 }$$

This gives a sinusoidal perturbation on the dust surface density with 1% amplitudes. We input only dust surface density perturbations at t = 0. The perturbations in Σ_d naturally induce perturbations in τ_s in the subsequent time evolution.

Figure 2 shows the time evolution of amplitude of the dust surface density perturbation for λ = 1 au (the green solid line). The evolution of the background stopping time τ_s,0(t) is also plotted with the gray dotted line. The dust surface density perturbation is first damped by the diffusion. After the initial damping, the perturbation exponentially grows at a larger rate than the coagulation rate (see the gray dotted line). The growth rate seen in the simulation slightly changes in time. Using Equation (18), we estimate linear growth rates of coagulation instability at time marked by the orange filled circles in Figure 2. The green dashed lines show the estimated growth rates. These are in good agreement with the exponential growth rate observed in the simulation. This result thus demonstrates the existence of coagulation instability predicted by our previous analysis (Tominaga et al. 2021). Linear growth rates are determined by four parameters η RΩ, ε₀, τ_s,0, and α, and only the stopping time τ_s,0 is time dependent. We thus attribute the observed time-dependent growth rate to the background evolution of dust sizes.

Zoom In Zoom Out Reset image size

Figure 3 shows the time evolution of δΣ_d/Σ_{d ,0} and δτ_s,0/τ_s,0 as a function of the radial position and time. As predicted by Tominaga et al. (2021), the dust size perturbation is shifted outward relative to the surface density perturbation (see Section 2.2 of Tominaga et al. 2021). Note that the dust cells move at the perturbed velocity δ v_x , and we extracted the background drift motion $\left\langle {v}_{x,0}\right\rangle$ (Equations (20) and (21)). Therefore, the inward radial propagation of the perturbations in Figure 3 means that the phase speed is larger than the background drift speed, which is also consistent with the results in Tominaga et al. (2021) (see also Equation (18)).

Zoom In Zoom Out Reset image size

Figure 4 compares the results of simulations for three different wavelengths: λ = 0.25, 1, and 4 au. Shorter-wavelength perturbations experience more significant initial damping due to the diffusion (e.g., the blue solid line). At a later time, each simulation shows exponential growth consistent with the previous linear analyses (see dashed lines), which in turn validates our numerical scheme.

Zoom In Zoom Out Reset image size

We note that the linear analyses ignored the background evolution of the stopping time. Regardless of this assumption, the growth rates from the linear analyses well reproduce the exponential growth rates observed in the local simulations. This in turn validates the assumption in the linear analyses.

We assume the following gas surface density and temperature profile:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{g}}}(r)={{\rm{\Sigma }}}_{{\rm{g}},10}{\left(\frac{r}{10\,\mathrm{au}}\right)}^{-q/2},\end{eqnarray} \tag{ 24 }$$

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

where Σ_{g ,10}, q, and T₁₀₀ are constant parameters. We assume an initial dust-to-gas ratio of 10⁻² in all runs:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{d}}}(r)={10}^{-2}{{\rm{\Sigma }}}_{{\rm{g}},10}{\left(\frac{r}{10\,\mathrm{au}}\right)}^{-q/2}.\end{eqnarray} \tag{ 26 }$$

Initial dust size is assumed to be 10 μm. The choice of the initial dust sizes insignificantly affects background evolution unless the dimensionless stopping time is much less than unity because dust grains initially just grow in size at the initial location.

The parameter q on the exponent of the surface density profile takes three different values: q = 1, 2, and 3. In most runs, we set ${{\rm{\Sigma }}}_{{\rm{g}},10}={{\rm{\Sigma }}}_{{\rm{g}},\mathrm{MMSN}}(10\,\mathrm{au})$. We adopt smaller Σ_{g ,10} for q = 1 so that the gas disks are stable to nonaxisymmetric GI in the domain: Q ≡ c_sΩ/πGΣ_g > 2.

We next set up initial perturbations. As in the local simulations (Section 3), we displace the cells and input perturbations only in dust surface density at t = 0. The displacement δ ξ adopted in the radially global simulations is the following (see Figure 1):

$$\begin{eqnarray}&&\delta \xi ={\xi }_{0}\sum _{j=0}^{128}\cos \left(\displaystyle \frac{2\pi {jr}}{{r}_{\mathrm{out}}-{r}_{\mathrm{in}}}+2\pi {\phi }_{j}\right),\end{eqnarray} \tag{ 27 }$$

where ξ₀ = 0.2(r_out − r_in)/N, r_out is the outer boundary radius, and ϕ_j is a random number less than unity. As described below, we adopt the radial domain size of 95 au in most runs with N = 2¹⁴ = 16,384. Thus, the 128th mode included in δ ξ has a wavelength less than 1 au. The above displacement gives initial Σ_d-perturbations of about 10%–20% at most, but dust diffusion decreases the perturbation amplitudes down to 1% or a few percent before the onset of coagulation instability. We note that the adopted δ ξ in the present simulations is small enough that cell crossing by the initial displacement does not occur.

To see wave propagation during the development of the instability, we input initial perturbations only in the outer half region, reducing the displacement δ ξ in all runs as follows:

$$\begin{eqnarray}&&\delta \xi \to \delta \xi \times \displaystyle \frac{1}{2}\left[1+\tanh \left(\displaystyle \frac{r-{R}_{0}}{{\rm{\Delta }}r}\right)\right],\end{eqnarray} \tag{ 28 }$$

where R₀ = 50 au and Δr = 1 au. If we input perturbations without the reduction, we observe fast growth of coagulation instability at inner radii and dust growth up to τ_s = 1 well before outer dust grains start drifting because the timescales of the instability are proportional to Ω⁻¹ ∝ r^3/2 (see also Appendix B).

Labels of runs, parameters, the adopted initial perturbation are listed in Table 1.

Table 1. Summary of Parameters and Results

Runs	q	$\tfrac{{{\rm{\Sigma }}}_{{\rm{g}},10}}{{{\rm{\Sigma }}}_{{\rm{g}},\mathrm{MMSN}}(10\,\mathrm{au})}$	T100 (K)	α	rin (au)	τs → 1 via CI?	(rring, Σd ,ring/Σg)
q3T28a1	3	1	28	1 × 10−4	5	Yes	(34.9 au, 2.80 × 10−1)
q3T20a1	3	1	20	1 × 10−4	5	Yes	(37.7 au, 2.96 × 10−1)
q3T10a1	3	1	10	1 × 10−4	5	Yes	(40.1 au, 1.65 × 10−1)
q2T28a1	2	1	28	1 × 10−4	5	Yes	(31.5 au, 1.67 × 10−1)
q2T20a1	2	1	20	1 × 10−4	5	Yes	(32.4 au, 1.50 × 10−1)
q2T10a1	2	1	10	1 × 10−4	5	Yes	(39.7 au, 1.23 × 10−1)
q1T28a1*	1	0.5	28	1 × 10−4	20	Yes	(30.3 au, 1.48 × 10−1)
q1T20a1*	1	0.5	20	1 × 10−4	20	Yes	(33.9 au, 1.40 × 10−1)
q1T10a1*	1	0.5	10	1 × 10−4	20	Yes	(33.6 au, 5.42 × 10−2)
q3T20a3	3	1	20	3 × 10−4	5	Yes	(13.6 au, 1.00 × 10−1)
q3T20a5	3	1	20	5 × 10−4	5	Yes	(6.74 au, 7.25 × 10−2)
q2T20a3	2	1	20	3 × 10−4	5	Yes	(16.4 au, 7.50 × 10−2)
q2T20a5	2	1	20	5 × 10−4	5	Yes	(9.22 au, 6.24 × 10−2)
q1T20a1	1	0.5	20	1 × 10−4	5	Yes	(27.3 au, 9.70 × 10−2)
q1T20a3	1	0.5	20	3 × 10−4	5	Yes	(8.72 au, 1.22 × 10−2)
q1T20a5	1	0.5	20	5 × 10−4	5	Yes	(7.15 au, 8.99 × 10−3)
q1T20a3*	1	0.5	20	3 × 10−4	20	No	⋯
q1T20a5*	1	0.5	20	5 × 10−4	20	No	⋯

Note. "CI" at the top of the seventh column stands for coagulation instability. The dust-to-gas ratio at the resulting rings (eighth column) is derived from the last time step of each simulation.

Download table as:  ASCIITypeset image

The number of cells N is 2¹⁴ = 16,384. The domain is uniformly spaced by N cells at t = 0. The main reason to adopt this large number of cells is that the nonlinear evolution of the instability results in very narrow structures whose width is less than 0.1 au as shown in the next subsection. Note that the Lagrangian scheme automatically increases the spatial resolutions of high-density regions. However, coagulation instability develops after dust grains are depleted as a result of the radial drift, meaning that the resolution of the background state first becomes lower than the initial resolution. Therefore, to describe nonlinear coagulation instability, we need the relatively large number of cells even when we utilize the Lagrangian scheme. We note again that the Lagrangian scheme is free from numerical diffusion regardless of the number of cells, and thus this is still the advantage of using the Lagrangian scheme.

The inner and outer boundaries are set initially at r_in = 5 au and r_out = 100 au in most runs. As a disk evolves and dust grows, the cells move inward and flow out of the domain. We thus let the inner boundary move with the first cell boundary while the outer boundary radius is fixed. We discuss disk evolution in r_in(t = 0) ≤ r ≤ r_out.

We adopt r_in = 20 au at t = 0 for a run with T₁₀₀ = 10 K and q = 1 because even background dust evolution results in dust sizes of τ_s = 1 at inner radii (q1T10a1* run). We also adopt r_in = 20 au in other runs to investigate parameter dependence under the same configuration of the domain (the q1T28a1* and q1T20a1* runs) and to check the r_in-dependence of the results (the q1T20a3* and q1T20a5* runs).

The moment method we utilize in this work might be valid for τ_s ≤ 1 because dust grains of τ_s = 1 are the most significantly depleted because of the fastest drift speed, and dust size distribution will become bimodal if there are larger solids of τ_s > 1. Therefore, we stop simulations once the dimensionless stopping time becomes larger than unity as a result of coagulation instability. The simulations otherwise last for 2 × 10⁵ yr.

The results are briefly summarized in the seventh and eighth columns of Table 1. In the eighth column, we list radial locations of dust surface density maximum, r_ring, and resulting dust-to-gas ratios, Σ_{d ,ring}/Σ_g, at the time when τ_s = 1 is achieved. We find that in most runs coagulation instability develops into the nonlinear growth phase before perturbations drift into the outside of the domain, locally increasing τ_s up to unity. We also find that dust-to-gas surface density ratio can increase at least by an order of magnitude after dust depletion down to Σ_d/Σ_g ∼ 10⁻³.

Figure 5 shows the results of the q2T20a1 run. The left two panels show the time evolution of dust surface density and dimensionless stopping time. We also show the evolution in the absence of the initial perturbation in the right two panels. As mentioned in the previous studies, inside-out dust coagulation and radial drift result in dust depletion and dust-to-gas ratio of order 10⁻³ (e.g., Brauer et al. 2008). In the presence of the initial perturbations, they exponentially grow via coagulation instability, and dust-to-gas ratio locally increases up to ∼10⁻¹. The dimensionless stopping time increases above unity around the most collapsed dust ring at t = 9.0 × 10⁴ yr.

Zoom In Zoom Out Reset image size

Regardless of the presence of the radial dust diffusion, we find that the nonlinear coagulation instability creates narrow rings. This is because the radial concentration is driven by the radial drift velocity ∼ τ_s η v_K, which is much larger than D/r by a factor of τ_s/α. The difference in the drift velocity can increase the dust surface density significantly until the density gradient ∂Σ_d/∂r becomes large enough that the dust diffusion balances with the drift motion.

Figure 6 shows dust surface density, dust sizes, and radial velocity profiles around the most collapsed ring in the q2T20a1 run. Dust sizes are 1–10 cm within the ring in this run. We note that very narrow ring structures with widths of ∼0.05 au are well resolved in our simulations. The bottom panel of Figure 6 shows that there is a radial gradient in v_r in 32.30 au ≲ r ≲ 32.35 au, while v_r is almost uniform around the Σ_d-peak radius (the dashed line). Similar velocity profiles are found in other runs. This indicates that the ring will sweep inner dust grains up while moving inward, and the dust surface density in the ring will increase through further evolution until significant dust growth results in τ_s ≫ 1 and quenches the dust drift. Figure 5 shows that there are also some rings that are still growing at r ≳ 37 au. Therefore, it is expected that further evolution will lead to dust retention at multiple locations with τ_s ≥ 1.

Zoom In Zoom Out Reset image size

Figure 7 shows results of the q2T20a5 run. Diffusion with larger α quenches short-wavelength modes that will exponentially grow for smaller α. Nevertheless, coagulation instability locally increases τ_s up to unity as well as in the q2T20a1 run. The most collapsed ring in the q2T20a5 run is located at smaller radii than in the q2T20a1 run. The reason for this is smaller growth rates of coagulation instability in the q2T20a5 run because of stronger turbulent diffusion. Indeed, in smaller-domain-size simulations (q1T20a3* and q1T20a5*) perturbations flow out across r = r_in(t = 0) before a significant increase in their amplitudes. One can see that dust grains initially located at the perturbed region drift into smaller radii in Figure 8. Figure 8 shows the time evolution of cumulative mass profiles for the q2T20a1 and q2T20a5 runs. The cumulative mass also includes the mass of cells flowing out of the domain. Because of the mass conservation, one can see the motion of cells from Figure 8. For instance, a dust cell initially at r = 50 au is located at intersections of the black dashed line and each cumulative mass profile. We also plot a line that shows the cumulative mass of M(<49 au) at t = 0 yr (gray dashed line) as a reference. The positions of the dust rings in each run correspond to positions of "cliffs" on the cumulative mass profiles. The rings in both runs consist of dust cells initially located around r ≃ 49–50 au. We thus attribute the difference in r_ring to the difference in radial distance over which dust drifts within the growth time of coagulation instability. ⁴ This interpretation is consistent with our previous linear analysis that shows that the phase velocity is roughly given by the dust drift velocity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We investigate η v_K-dependence of the ring position and resulting dust-to-gas ratio by changing T₁₀₀ for a given q and α = 1 × 10⁻⁴. Figure 9 summarizes the results (solid lines with marks). In the left panel, the shaded regions show the final positions of dust cells that are initially located at 49 au ≤ r ≤ 50 au. All marks are in the shaded regions. This means that perturbations initially located around the inner edge of the perturbed region (r ≃ 50 au) have enough time to develop via coagulation instability toward the nonlinear phase; the ring would consist of dust cells initially located at r ≫ 50 au if the growth timescale was much longer than the drift timescale. We find that when q is fixed the radial position r_ring decreases as η v_K. Since the shaded regions show the similar trend, we attribute the η v_K-dependence seen for a given q to the difference in the drift velocity: dust grains with a certain τ_s have faster drift speed for larger η v_K and reach inner regions.

Zoom In Zoom Out Reset image size

The left panel of Figure 9 also shows that in most cases r_ring is smaller for smaller q, although smaller-q runs show lower values of η v_K (see dark-blue and orange lines). This trend originates from the radial profile of the background τ_s. Figure 10 shows the radial profile of τ_s in the q1T20a1 and q3T20a1 runs. At the final time step (dark-blue lines), background τ_s reaches the drift-limited value at inner radii. The q1T20a1 run shows r-dependent background τ_s, while the q3T20a1 run shows almost uniform background τ_s. The background τ_s thus becomes larger in the q1T20a1 run than in the q3T20a1 run, meaning that the drift velocity proportional to τ_s η v_K can become larger for q = 1 even though η v_K is smaller. Therefore, dust grains and perturbations move inward more efficiently for smaller q, which explains q-dependence in Figure 9.

Zoom In Zoom Out Reset image size

The resulting Σ_{d ,ring}/Σ_g is almost similar in all runs, although we can see slightly larger values for larger η v_K and larger q. This slight dependence is also explained based on the difference in the τ_s profile as in the q-dependence in the left panel. ⁵ Larger η v_K reduces the background drift-limited τ_s because of the higher drift speed. In such a case, coagulation instability locally increases dust surface density more before τ_s in rings reaches unity via the instability and we stop the simulations. If q is small, the background τ_s becomes larger especially at inner radii as mentioned above, and there is little space for τ_s perturbations to develop up to unity, leading to smaller dust surface density.

As shown above, we observe significant dust concentration into a ring via coagulation instability. There are some rings that are growing toward the nonlinear phase (see Figure 5). The radial drift of the rings halts as a result of further dust growth beyond τ_s = 1. On the other hand, remaining dust grains in gaps and other less dense regions still suffer from the radial drift toward the central star. For retaining such dust grains in the disk, either of the following two processes will be necessary: (A) coagulation instability of the remaining dust grains operates and triggers concentration into second-generation rings, or (B) the first-generation rings catch the remaining dust that radially drifts through the ring.

Process (A) will be expected since any disturbances on a disk give rise to seed perturbations; in the present simulations, we input seed perturbations only at t = 0 yr. If dust grains in the first-generation rings significantly grow in size (τ_s ≫ 1), coagulation instability of the remaining dust grains might be insensitive to the rings. Therefore, successive development of coagulation instability will play the key role in retaining dust grains in a disk.

Next, we discuss the efficiency of process (B). We evaluate collisional "optical depth" regarding a collision between a drifting dust grain and dust grains in the first-generation ring. The associated mean free path of a drifting dust at the midplane in a ring, l, is given by

$$\begin{eqnarray}&&{l}^{-1}=\frac{{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{ring}}}{\sqrt{2\pi }{H}_{{\rm{d}}}({a}_{\mathrm{ring}}){m}_{{\rm{p}},\mathrm{ring}}}\pi {\left({a}_{\mathrm{ring}}+{a}_{\mathrm{drift}}\right)}^{2},\end{eqnarray} \tag{ 29 }$$

where a_ring and a_drift are the sizes of the dust in a ring and the drifting dust, respectively, and ${m}_{{\rm{p}},\mathrm{ring}}\equiv 4\pi {\rho }_{\mathrm{int}}{a}_{\mathrm{ring}}^{3}/3$. We assume ${H}_{{\rm{d}}}/H=\sqrt{\alpha /{\tau }_{{\rm{s}}}}$ and the Epstein law, and thus τ_s,ring/τ_s,drift = a_ring/a_drift in the ring. For a ring width of ΔR_ring, collisional "optical depth," τ_coll, is $\sqrt{1+{({\rm{\Delta }}{v}_{\phi }/{\rm{\Delta }}{v}_{r})}^{2}}{\rm{\Delta }}{R}_{\mathrm{ring}}/l$, and thus

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}{\tau }_{\mathrm{coll}} & \simeq & 0.09\left(\frac{{\rm{\Delta }}{R}_{\mathrm{ring}}/H}{0.02}\right)\left(\frac{{{\rm{\Sigma }}}_{{\rm{d}},\mathrm{ring}}/{{\rm{\Sigma }}}_{{\rm{g}}}}{0.1}\right)\\ & & \times {\left(\frac{\alpha }{{10}^{-4}}\right)}^{-\tfrac{1}{2}}{\left(\frac{{\tau }_{{\rm{s}},\mathrm{ring}}}{1}\right)}^{-\tfrac{1}{2}}{\left(1+\frac{{\tau }_{{\rm{s}},\mathrm{drift}}}{{\tau }_{{\rm{s}},\mathrm{ring}}}\right)}^{2}\\ & & \times \sqrt{1+\frac{1}{4}{\left(\frac{1+{\tau }_{{\rm{s}},\mathrm{drift}}{\tau }_{{\rm{s}},\mathrm{ring}}^{-1}}{{\tau }_{{\rm{s}},\mathrm{ring}}^{-1}-{\tau }_{{\rm{s}},\mathrm{drift}}}\right)}^{2}}.\end{array}\end{eqnarray} \tag{ 30 }$$

The value of ΔR_ring/H is roughly given by the value of the most collapsed ring in the q2T20a1 run: ΔR_ring ∼ 0.05 au and H ≃ 2.17 au at r ≃ 32.3 au. The value of τ_s,drift/τ_s,ring in the last parentheses in the second row of Equation (30) is less than unity, and thus we can neglect it. The last factor comes from $\sqrt{1+{({\rm{\Delta }}{v}_{\phi }/{\rm{\Delta }}{v}_{r})}^{2}}$ and is less than 10 for τ_s,drift = 0.1 ≪ τ_s,ring. We then obtain τ_coll < 1. The collisional optical depth can be much larger than unity only in a special case where the ring and the drifting dust have a similar radial velocity, i.e., ${\tau }_{{\rm{s}},\mathrm{drift}}\simeq {\tau }_{{\rm{s}},\mathrm{ring}}^{-1}$. This is not our main focus of this discussion since we consider dust grains that have radial relative velocity and drift with respect to the ring. In the limit of τ_s,ring ≫ 1, the last factor converges to ∼ 0.5/τ_s,drift ∼ 5 for τ_s,drift = 0.1. In such a limit, we obtain τ_coll ≪ 1 because of the τ_s,ring-dependence appearing in the second row of Equation (30). In this way, Equation (30) indicates that the dust in the ring and the drifting dust are collisionless, i.e., the drifting dust grains can pass through the ring, especially after dust growth proceeds and τ_s,ring becomes larger than unity. Therefore, we expect that process (B) is inefficient in the present case.

As described above, process (A) is more important to retain dust grains in a disk. Note that Equation (30) validates the hypothesis of process (A) since it will occur when remaining dust grains are insensitive to large grains in the rings. For the further discussion, we need simulations with full-size distributions or at least two-population simulations that can treat drifting dust and large solids (e.g., Dra̧żkowska & Alibert2017). This will be addressed in our future studies.

Our previous work proposed a scenario in which if coagulation instability operates in a disk, the instability is expected to set up preferable locations for secular GI to develop toward planetesimal formation. In this work, we demonstrate that coagulation instability operates indeed during the inside-out dust evolution. Those results support our scenario. As shown in Appendix C, the rings observed in the present simulations can be unstable to secular GI if they form in outer regions where Toomre's Q is relatively small, for example, Q ≲ 10. Weaker turbulence is also required as shown in the previous studies on secular GI (e.g., Shariff & Cuzzi 2011; Youdin 2011; Takahashi & Inutsuka 2014; Tominaga et al.2019).

Although the present simulations demonstrate the possible path toward planetesimal formation, we have to treat the back-reaction for rigorous discussions. The effect of the back-reaction is expected to be insignificant in the linear growth phase. This is because coagulation instability starts to grow after dust density decreases and dust-to-gas surface density ratio becomes on the order of 10⁻³. On the other hand, the back-reaction will become significant in the nonlinear growth phase, which will limit the maximum dust density in rings. We give further discussion on the development to secular GI in Paper II, where we present simulations with not only the back-reaction but also fragmentation.