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

Our previous linear analysis presents a new instability driven by dust coagulation in protoplanetary disks. The coagulation instability has the potential to concentrate dust grains into rings and assist dust coagulation and planetesimal formation. In this series of papers, we perform numerical simulations and investigate nonlinear outcome of coagulation instability. In this paper (Paper I), we first conduct local simulations to demonstrate the existence of coagulation instability. Linear growth observed in the simulations is in good agreement with the previous linear analysis. We next conduct radially global simulations to demonstrate that coagulation instability develops during the inside-out disk evolution due to dust growth. To isolate the various effects on dust concentration and growth, we neglect effects of backreaction to a gas disk and dust fragmentation in Paper I. This simplified simulation shows that either of backreaction or fragmentation is not prerequisite for local dust concentration via the instability. In most runs with weak turbulence, dust concentration via coagulation instability overcomes dust depletion due to radial drift, leading to the formation of multiple dust rings. The nonlinear development of coagulation instability also accelerates dust growth, and the dimensionless stopping time $\tau_{\mathrm{s}}$ reaches unity even at outer radii (>10 au). Therefore, coagulation instability is one promising process to retain dust grains and to accelerate dust growth beyond the drift barrier.


INTRODUCTION
Growth from (sub-)micron-sized dust grains to kmsized planetesimals is one long-standing problem in planet formation theory. Recently, exploration by the New Horizons spacecraft provided data that constrains planetesimal forming processes at outer radii in early solar nebula (Stern et al. 2019;McKinnon et al. 2020). The bi-lobe shape of one Kuiper Belt object, Arrokoth, indicates that the object might form via gravitational collapse of a dust clump in the presence of a gas disk. This scenario seems consistent with one proposed scenario of planetesimal formation (see also Nesvorný et al. 2019): streaming instability and subsequent gravitational instability (GI) (e.g., Youdin & Goodman 2005;Bai & Stone 2010;Simon et al. 2016;Abod et al. 2019;Carrera et al. 2021). However, recent studies show that the streaming instability is substantially stabilized by diffusion in gas turbulence (Chen & Lin 2020;Umurhan et al. 2020). McNally et al. (2021) investigated growth rates of the streaming instability in the presence of the turbulent diffusion with dust size distribution. They showed that even weak turbulence with α ∼ 10 −5 can quench the streaming instability (see Figure 8 therein), where α is the turbulence strength (Shakura & Sunyaev 1973). Therefore, it is important to investigate other processes potentially leading to outer planetesimal formation, including secular GI (e.g., Ward 2000;Youdin 2005aYoudin ,b, 2011Shariff & Cuzzi 2011;Michikoshi et al. 2012;Takahashi & Inutsuka 2014;Tominaga et al. 2018Tominaga et al. , 2019Tominaga et al. , 2020Pierens 2021).
Our previous work found a new instability driven by dust coagulation in a protoplantary disk (Tominaga et al. 2021). The instability that we named "coagulation instability" leads to clumping of drifting dust even in the presence of the turbulent diffusion with α ∼ 10 −4 . Another remarkable difference from the previous dust-gas instabilities is that coagulation instability can develop at tens of orbital periods even when dust-to-gas surface density ratio is so small (∼ 10 −3 ) that the previous instabilities hardly develop.
According to the previous studies on coagulation, dust grains are depleted as a result of the inside-out coagulation and radial drift, leading to a low dust-to-gas ratio of ∼ 10 −3 (e.g., Brauer et al. 2008;Birnstiel et al. 2009;Okuzumi et al. 2012). This indicates that coagulation instability will operate earlier than the other instabilities proposed in the previous studies such as streaming instability or secular GI. Tominaga et al. (2021) indicates that coagulation instability accelerates dust growth and may form planetesimals directly. At the same time, dust concentration and growth via coagulation instability can also set up preferable regions for the other instabilities to subsequently develop toward planetesimal formation. Since the instability potentially has a significant contribution to planetesimal formation, the viability of coagulation instability and the above hypothesis should be tested by numerical simulations.
In this series of papers, we perform numerical simulations and investigate the nonlinear outcome of coagulation instability. This paper (Paper I) presents the first results of numerical simulations. We conduct local simulations to prove the existence of coagulation instability as predicted in our linear analysis (Tominaga et al. 2021). We also conduct radially global simulations and show that coagulation instability does develop during the inside-out disk evolution via dust growth. The instability is found to cause radial dust concentration and to accelerate dust growth. In Paper I, we intentionally neglect frictional backreaction (deceleration of the drift) and dust fragmentation to isolate the contributions to dust concentration and growth from various effects. This simplification allows clear distinction between linear/nonlinear coagulation instability and selfinduced dust trap process due to a combination of coagulation and backreaction (Gonzalez et al. 2017). In other words, such simplified simulations demonstrate that either of backreaction or fragmentation is not prerequisite for dust concentration via the instability. Therefore, based on this simplification, we aim to obtain clear physical understanding of nonlinear outcome of coagulation instability in Paper I. Paper II (Tominaga et al. 2022b, submitted) presents simulations of coagulation instability under the influence of backreaction and fragmentation.
Paper I is organized as follows. We describe basic equations and numerical methods in Section 2. In Section 3, we first briefly review the linear analysis of coagulation instability (Section 3.1) and show results of local simulations that we compare with the linear analysis. We present radially global simulations of coagulation instability in Section 4 and discuss to what extent coagulation instability concentrates dust grains before dust falls onto a central star. In Section 5, we discuss further evolution of resulting dust rings and possible way to retain dust grains in a disk. We summarize in Section 6.

Basic equations
In this work, we consider dust evolution in a steady axisymmetric gas disk around a star of 1M ⊙ . Equation of continuity for dust in cylindrical coordinates (r, φ, z) is given by 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 K = GM ⊙ /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): where τ s ≡ t stop Ω is the stopping time of dust grains, t stop , normalized by the Keplerian angular velocity Ω ≡ v K /r = GM ⊙ /r 3 . The value of η for a given gas density ρ g and temperature T is determined by where c s = k B T /µm 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): We calculate the radial dust diffusion using the gradient ∂Σ d /∂r (e.g., Cuzzi et al. 1993) instead of Σ g ∂ (Σ d /Σ g ) /∂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, is small compared to the mean drift velocity in weakly turbulent disks considered in this work. For D ∼ αc 2 s /Ω and η ∼ (c s /v K ) 2 , one obtains The right hand side is much smaller than unity in weakly turbulent gas disks (α ∼ 10 −4 ) 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, Σ −1 g ∂Σ g /∂r is smaller than Σ −1 d ∂Σ d /∂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): where m p and a are the mass and the radius of the representative dust, ∆v pp is the collision velocity, and H d is the dust scale height. Equation (8) 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 ∆v pp = √ Cατ s c 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 : where ∆v B is the collision velocity due to Brownian motion, ∆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): Because our model is based on the vertically averaged equations, we adopt z = 2/πH d,1,2 to calculate the vertically averaged ∆v 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 : where ρ int = 1.4 g cm −3 is an internal mass density of dust grains (cf. Appendix D), 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 ρ g = Σ g / √ 2πH, we obtain (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.

Numerical methods
We use the Lagrangian-cell method developed in Tominaga et al. (2018) for both local and global simulations in the following sections. The advantage to use the Lagrangian scheme is the fact that the scheme is free from numerical diffusion associated with advection, which enable us to accurately describe the linear growth of coagulation instability. Using Equations (2) and (3), we update position of 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-dustrest 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 v r (r d,i ) 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 both time stepping, we use the second-order Runge-Kutta scheme 1 . 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 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. 2012Meyer et al. , 2014.

LOCAL SIMULATIONS
In this section, we first briefly review physical properties of coagulation instability from the previous lin-1 Although the Lagrangian-cell method developed in Tominaga et al. (2018) is used with a symplectic integrator in our previous studies, we adopt the Runge-Kutta method in this work because we assume the terminal velocity for dust instead of solving equations of motion. ear analyses (Tominaga et al. 2021). We next show results of local simulations of dust density and size evolution with initial perturbations. The simulations show exponential growth of the input perturbations. We compare the results with the previous linear analysis (Tominaga et al. 2021) and demonstrate the existence of coagulation instability.

Physical properties of coagulation instability
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 one-fluid mode (see Section 3 of Tominaga et al. 2021). We thus review 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 ∆v pp = √ Cατ s c 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 1965a), the equation is where t 0 ≡ (4/3 √ C)Ω −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: 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, the dimensionless stopping time in the unperturbed state is 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 0 (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.
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 where we set 1 ± τ 2 s,0 ≃ 1 assuming τ s,0 < 1. Equation (19) shows that the growth rate is determined by a geometric mean of the coagulation rate ε 0 /3t 0 and k| v x,0 |. The latter represents a rate of traffic jam induced by τ s -perturbations at a scale of k −1 . 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).
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).

Setups of local simulations
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: δ 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: 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 ∆v pp = √ Cατ s c 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 Σ g,0 = Σ g,MMSN (20 au) where Σ g,MMSN (r) = 53.75 g cm −2 (r/10au) −1.5 is the gas surface density in the MMSN disk model. We assume uniform dust-to-gas ratio of 10 −3 in the local simulations. We initially adopt τ s = t stop Ω(R) = 10 −2 throughout the domain, and set the corresponding dust sizes using Equation (13). We adopt α = 10 −4 , and 10 −5 in the local simulations. The results for α = 10 −5 are presented in Appendix A. Figure 1 schematically shows the initial condition explained below. We simulates the time evolution of a si- The initial dust surface density perturbation is damped by the diffusion at a rate of Dk 2 (the black dotted line). The perturbation exponentially grows after the initial damping. The green dashed lines show the growth rates derived in the previous linear analyses (Equation (18)). We use physical values at time marked by the orange filled circles to derive the growth rate. The linear analyses are in good agreement with the simulation at each point.
nusoidal 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 1st 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 δξ: This gives a sinusoidal perturbation on the dust surface density with one percent 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 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Ω ε 0 , τ 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. Figure 3 shows 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 v x,0 (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)). Figure 4 compares the results of simulations for three different wavelengths: λ = 0.25 au, 1 au, and 4 au. Shorter wavelength perturbations experience more significant initial damping due to the diffusion (e.g., the blue solid line). At later time, each simulation shows exponential growth consistent with the previous linear analyses (see dashed lines), which in turn validates our numerical scheme.

Results and comparison with the linear analyses
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.

RADIALLY GLOBAL SIMULATIONS
In this section, we show results of radially global simulations. We investigate (1) whether the instability develops into nonlinear growth phase before dust grains fall onto a central star, and τ s locally increases up to unity, and (2) how much dust surface density and dust-to-gas ratio increase.  ) and δτs,0/τs,0 (the right panel) obtained from the local simulation with λ = 1 au and α = 10 −4 . Note that the horizontal axis is the radial position X in the background-dust-rest frame (Equation (22)). The phase shift between δΣ d and δτs and the radial propagation of the perturbations are consistent with the previous linear analyses (Tominaga et al. 2021).   (18)). We use physical values at time marked by the orange filled circles to derive the growth rate. The linear analyses have a good agreement with each simulation.

Disk models
We assume the following gas surface density and temperature profile: T (r) = T 100 r 100 au where Σ g,10 , q and T 100 are constant parameters. We assume initial dust-to-gas ratio of 10 −2 in all runs: 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 Σ g,10 = Σ g,MMSN (10 au). We adopt smaller Σ g,10 for q = 1 so that the gas disks are stable to nonaxisymmetric gravitational instability 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 (cf. Figure 1): where ξ 0 = 0.2(r out − r in )/N , r out is the outer boundary radius, and φ j is random number less than unity.
As described below, we adopt the radial domain size of 95 au in most runs with N = 2 14 = 16384. Thus, the 128th mode included in δξ has a wavelength less than 1 au. The above displacement gives initial Σ dperturbations of about 10-20 percents at most, but dust Note-"CI" on the top of 7th column stands for coagulation instability. The dust-to-gas ratio at the resulting rings (8th column) is derived from the last time step of each simulation.
diffusion decreases the perturbation amplitudes down to one percent or a few percents 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: where R 0 = 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 Ω −1 ∝ r 3/2 (see also Appendix B). Labels of runs, parameters, the adopted initial perturbation are listed in Table 1.

Numerical setups
The number of cells N is 2 14 = 16384. 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 struc-tures 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 to use 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 1st 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 100 = 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 . Time evolution of dust surface density (top panels) and dimensionless stopping time (bottom panels) obtained from the q2T20a1 run. The left two panels and the right two panels show evolution with and without initial perturbations, respectively. Coagulation instability locally increases dust-to-gas ratio up to ∼ 10 −1 . We stop the simulation at t = 9.0 × 10 4 yr since the dimensionless stopping time τs becomes larger than unity around the most collapsed dust ring. One can see that the ring structure is well resolved even though the radial width is less than 0.1 au. We find that radial gradient in vr around the left-side tail of Σ d profile (32.3 au r 32.35 au) while vr is almost uniform within the ring. Thus, the ring will sweep inner dust grains up, and the surface density increases through further evolution.
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 is 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 5 yr.

Results
The results are briefly summarized in 7th and 8th columns of Table 1. On the 8th 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 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 −3 . 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 . Time evolution of dust surface density (the left panel) and dimensionless stopping time (the right panel) obtained from the q2T20a5 run. As in the q2T20a1 run, we observe a local increase in dust-to-gas ratio by an order of magnitude. The dimensionless stopping time locally becomes larger than unity at t = 8.7 × 10 4 yr. the evolution in the absence of the initial perturbation on 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 −3 (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 −1 . The dimensionless stopping time increases above unity around the most collapsed dust ring at t = 9.0 × 10 4 yr.
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 that 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 width 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 with 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.
Figures 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 of this is smaller growth rates of coagulation instability in the q2T20a5 run because of stronger turbulent diffusion. Indeed, in smallerdomain-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 time evolution of cumulative mass profiles for the q2T20a1 and q2T20a5 runs. The cumulative mass also includes 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 black dashed line and each cumulative mass profiles. We also plot a line which 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 2 . This interpretation is consistent with our previous linear analysis that shows the phase velocity is roughly given by the dust drift velocity.
We investigate ηv K -dependence of the ring position and resulting dust-to-gas ratio by changing T 100 for a given q and α = 1 × 10 −4 . Figure 9 summarizes the results (solid lines with marks). On 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.
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 value of ηv K (see dark blue and orange lines). This trend originates from the radial profile 2 The background evolution in the q2T20a5 run is slightly faster than in the q2T20a1 run when dust grains are small and α > τs. We find that difference in the dust scale height leads to the difference in the dust growth time in such a case. Note that when dust grains become large (α ≫ τs) and we can approximate the dust scale height as H d ≃ α/τsH, the coagulation timescale is hardly dependent of α (see also Brauer et al. 2008;Okuzumi et al. 2012;Sato et al. 2016). In particular, the coagulation timescale is independent from α for the turbulence-driven coagulation (∆vpp = ∆vt ∝ √ αcs).
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. The resulting Σ d,ring /Σ g is almost similar in all run although we can see slightly larger value 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 on the left panel 3 . 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.

On possible ways to retain dust grains in a disk
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 Figure 9. Dependence of rring (left panel) and Σ d,ring /Σg (right panel) on ηvK. On both panels, solid lines with marks show the resulting rring and Σ d,ring /Σg. On the left panel, we plot the radial location of dust cells that are initially located in 49 au ≤ r ≤ 50 au with the shaded regions. We find that larger ηvK makes the ring radius smaller. This is because dust grains with a certain τs move faster for larger ηvK.  Figure 10. Radial profiles of τs in the q1T20a1 run (the left panel) and the q3T20a1 run (the right panel). In each run, τs reaches unity in the most collapsed ring at t = 9.6 × 10 4 yr and 6.1 × 10 4 yr, respectively. The background τs increases as r decreases in the q1T20a1 run and becomes larger than the background τs in the q3T20a1 run. This means that dust grains drift inward more efficiently in the former run, which explains the trend that rring is smaller for smaller q. 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.
The 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 the 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 where a ring and a drift are the sizes of the dust in a ring and the drifting dust, respectively, and m p,ring ≡ 4πρ int a 3 ring /3. We assume H d /H = α/τ 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 1 + (∆v φ /∆v r ) 2 ∆R ring /l and thus 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 parenthesis on the second row of Equation (30) is less than unity, and thus we can neglect it. The last factor comes from 1 + (∆v φ /∆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 the similar radial velocity, i.e., τ s,drift ≃ τ −1 s,ring . 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 on 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 the process (B) is inefficient in the present case.
As described above, the process (A) is more important to retain dust grains in a disk. Note that Equation (30) validates the hypothesis of the 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 twopopulation simulations that can treat drifting dust and large solids (e.g., Drażkowska & Alibert 2017). This will be addressed in our future studies.

Comments on the development to secular GI
Our previous work proposed a scenario that, 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 region 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., Youdin 2011;Shariff & Cuzzi 2011;Takahashi & Inutsuka 2014;Tominaga et al. 2019).
Although the present simulations demonstrate the possible path toward planetesimal formation, we have to treat the backreaction for rigorous discussions. The effect of the backreaction 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 −3 . On the other hand, the backreaction 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 backreaction but also fragmentation.

SUMMARY
Formation of planetesimals has been widely studied but still under debate. In the case of outer planetesimal formation, disk instabilities have been proposed as a possible mechanism. In particular, some studies indicate that planetesimal formation via gravitational instability subsequent to streaming instability may explain the origin of Kuiper Belt Objects (e.g., McKinnon et al. 2020). However, the viability of such instabilities is still under debate because even weak turbulence (α ∼ 10 −5 ) can suppress streaming instability (e.g., McNally et al. 2021).
Our previous study (Tominaga et al. 2021) proposed a new instability named coagulation instability, which is expected to trigger dust concentration and to assist planetesimal formation in a protoplanetary disk. Coagulation instability develops even in the presence of turbulence of α ∼ 10 −4 .
In this work (Paper I), we perform numerical simulations and prove that coagulation instability operates during the inside-out evolution via dust growth and even for the background dust-to-gas ratio of ∼ 10 −3 . The development of coagulation instability causes radial dust concentration into rings against the dust depletion due to the radial drift. Nonlinear growth also increases τ s up to unity in the rings even at outer radii (> 10 au), showing the possible path of dust growth beyond drift barrier. We intentionally exclude the backreaction and simplify the simulations to isolate the contribution to the dust concentration. The simplified simulations demonstrate that the backreaction is not prerequisite for the dust concentration via coagulation instability, which is in contrast to the self-induced dust trap (Gonzalez et al. 2017).
We measure the enhanced dust-to-gas ratio in the most collapsed ring at the time when the dimensionless stopping time τ s becomes unity via coagulation instability. We find that the dust-to-gas ratio increases by an order of magnitude or more and that the ratio can reach ∼ 0.1 in the present simulations without backreaction. The resulting dust-to-gas ratio Σ d,ring /Σ g at the final time step is only slightly higher in a disk with larger ηv K (Figure 9). This is because background evolution leads to smaller τ s (drift-limited), and coagulation instability increases dust density more until τ s = 1 is achieved.
We also find that stronger turbulence slows the growth of coagulation instability, and as a result rings form at inner radii as α increases (Table 1 and Figure 8). Too large α prevents coagulation instability from developing and creating dense rings before perturbations flow out of the numerical domain.
We discuss the further development of resulting rings and dust grains in Section 5. There are two possible ways to retain dust grains after the formation of the first generation rings: (A) coagulation instability of the remaining dust grains operates and form next generation rings, or (B) the first generation rings catch the remaining dust. We expect that the process (A) is more important for the dust retention. The process (B) is insignificant and dust grains coming from the outer region will pass though the ring without collisions once dust grains in the rings becomes large enough (τ s,ring ≫ 1; see Equation (30)). Tominaga et al. (2021) proposed a scenario that a combination of coagulation instability and secular GI will explain planetesimal formation. The results of the present simulations support the scenario. As briefly mentioned in Section 5.2, the radial concentration observed in the present simulations will be the key to trig-ger secular GI in a disk (see also Appendix C). Further discussion is given in Paper II.
In this paper, we neglect the effects of backreaction and dust fragmentation to isolate the contributions to dust growth and concentration. The simulations under this assumption are worth being presented since this simplification allows clear distinction between coagulation instability and the self-induced dust trap proposed in Gonzalez et al. (2017): the backreaction is not prerequisite for the local dust concentration via coagulation instability. We expect that a combined process of coagulation instability and the self-induced dust trap is one promising mechanism to save dust grains (see discussion in Paper II). Linear growth phase is insignificantly affected by the backreaction since coagulation instability starts to grow after dust-to-gas ratio decreases (∼ 10 −3 ). On the other hand, nonlinear growth of coagulation instability will be affected by the backreaction. The backreaction decelerates the radial drift and makes it easier for dust grains to overcome drift barrier in rings before Σ d /Σ g increases significantly. Besides, the linear growth rate is reduced by fragmentation according to the linear analysis in Tominaga et al. (2021). In Paper II, we present results of simulations where the effects of the backreaction (the drift deceleration) and the fragmentation are taken into account.

APPENDIX
A. LOCAL SIMULATIONS FOR WEAKER TURBULENCE Figure 11 shows results of the local simulations with α = 10 −5 . We perform three runs with three different wavelengths: λ = 0.25 au, 1 au, 4 au. As in the case of α = 10 −4 , we find a good agreement between the previous linear analyses and the local simulations. The perturbations grows faster than in the simulations with α = 10 −4 because the diffusion is ineffective. The ratio of the growth rates and the coagulation rate (see the gray dotted line) is also larger. Thus, coagulation instability more quickly develops at shorter spatial scales in more weakly turbulent disks than assumed in Section 4, leading to dust growth beyond the drift barrier. ηv K ≃ 54 m/s, α = 10 −5 λ = 0.25 au λ = 1.00 au λ = 4.00 au τ s, 0 (t) Figure 11. Results of the local simulation with α = 10 −5 for λ = 0.25 au, 1 au, and 4 au. The blue, green and yellow solid lines show the evolution of δΣ d /Σ d,0 for each wavelength. The gray dotted line shows time evolution of the background stopping time τs,0. The blue, green and yellow dashed lines show the growth rates derived in the previous linear analyses (Equation (18)). We use physical values at time marked by the orange filled circles to derive the growth rate. The linear analyses have a good agreement with each simulation.

B. CRITICAL RADIUS FOR COAGULATION INSTABILITY TO SAVE DUST GRAINS
In the present simulations, we input initial perturbations only in the outer half region and observe that the perturbations can grow before they flow out of the numerical domain in most runs. This process helps dust grains to avoid radial drift and depletion. In this appendix, we roughly estimate a "critical radius" that is a boundary dividing the radial regions where coagulation instability can save dust grains or not.
The necessary condition for the instability to save dust grains is that the radial drift timescale is longer than the growth timescale of coagulation instability. Similar discussion is given in Section 4.4 in Tominaga et al. (2021). We here focus on the critical radius at which the ratio of the above two timescales is unity while we focused on τ s -and εdependences of the ratio in the previous study. Following the notation in Tominaga et al. (2021), we denote the radial drift timescale by t travel , which is given by For simplicity, we assume |v r | ≃ 2τ s,0 ηv K . The necessary condition to save dust is t travel > n −1 , which yields Note that the growth rate has τ s,0 -and ηv K -dependences. Roughly speaking, we can expect the growth rate to be proportional to the square root of those physical quantities. Thus, as a net dependence, the critical radius given the right hand side of Equation (B2) depends roughly on the square root of τ s,0 and ηv K . We also note that the growth rate also depends on α and ε (see Tominaga et al. 2021). Equation (B2) states that dust grains initially located at r ≃ 10 au are potentially saved by coagulation instability. The adopted value of n in Equation (B2) is about the maximum growth rate for τ s,0 = 0.1, α = 1 × 10 −4 and ε = 1 × 10 −3 in the MMSN disk (e.g., see Figure 8 in Tominaga et al. (2021)). According to the results of the radial global simulations, coagulation instability already operates when ε becomes just a few times lower than the initial value (0.01). In such cases, the growth rate is larger than the above value according to the linear analyses in Tominaga et al. (2021), and the critical radius will decrease. On the other hand, the critical radius becomes larger if the radial diffusion is stronger since the growth rate decreases.
Note that perturbations that start growing at inner radii have larger growth rates since the growth rate is scaled by Ω ∝ r −3/2 , and thus the critical radius for such inner perturbations is smaller. For example, the critical radius is about 3.8 au for n = 3 × 10 −3 Ω(10 au), τ s,0 = 0.1 and ηv K = 54 m/s. In this case, inner dust grains will be saved as the inner perturbations grow via the instability. Finally, we also note that Equation (B2) is the necessary condition. If the initial amplitudes of perturbations are small, it takes multiple growth timescales for coagulation instability to concentrates dust grains significantly. In such a case, the critical radius is larger than shown in Equation (B2).

C. SECULAR GRAVITATIONAL INSTABILITY IN RESULTING DUST DENSE REGIONS
The radial dust concentration shown in the present simulations is the possible connection to planetesimal formation via secular GI (e.g., Ward 2000;Youdin 2005aYoudin ,b, 2011Michikoshi et al. 2012;Takahashi & Inutsuka 2014Tominaga et al. 2018Tominaga et al. , 2019Tominaga et al. , 2020Pierens 2021). Although the backreaction should be taken into account in the simulations for more rigorous discussion (Section 5.2), it is still worth noting whether secular GI is operational in the rings observed in the present simulations without the backreaction. In this appendix, we thus briefly examine the stability of the rings.
Growth efficiency of secular GI is determined by four parameters: dust-to-gas ratio, dimensionless stopping time, turbulence strength, and Toomre's Q value for gas. A condition for the growth of secular GI is given as follows (Tominaga et al. 2019) (1 + ε) τ s (ε + c 2 d /c 2 s ) +D(1 + ε) whereD ≡ Dc −2 s Ω and c d ∝ √ αc s is velocity dispersion of dust grains (see also Youdin 2005a;Youdin & Lithwick 2007). We briefly describe the physical insight into the growth condition of secular GI below. ForD ∼ c 2 d /c 2 s ∼ α ≪ τ s andD ≪ τ s ε, Equation (C3) is reduced as follows (see also Equation (21) in Takahashi & Inutsuka (2014)): This simplified equation clearly represents the key physical conditions for the onset of secular GI. The parameter Q gd is a measure of how massive the dusty gas disk is. The other parameter Q drag is a ratio of two timescales: (1) diffusion time over a spatial scale of H, and (2) a timescale for which dust moves over a distance H at the terminal velocity t stop πGΣ d . Secular GI can grow in self-gravitationally stable disk with Q gd > 1 when the dust diffusion is weak: Q drag < 1 is necessary. Note that the diffusion coefficient D in the linear analyses of secular GI is related to the radial diffusion. If the vertical diffusivity is also given by D, we obtain another form of the parameter Q drag as follows where ρ d = Σ d / √ 2πH d is the midplane dust density, and we also adopt H d = D/τ s Ω (e.g., Youdin & Lithwick 2007). If we neglect the factor difference, the numerator of the first term, Ω 2 /G, roughly corresponds to the Roche-like density that determines the stability of three dimensional mode in a self-gravitating disk (e.g., see Goldreich & Lynden-Bell 1965b;Sekiya 1983). We can see that the above necessary condition Q drag < 1 can be satisfied even for ρ d Ω 2 /G since the dust layer is much thin (H d /H ≪ 1) for grown dust grains. Therefore, secular GI does not require gravitationally unstable dust layer for its onset in contrast to the classical GI.
We include c d in the following discussion although our simulations does not include the effect of the velocity dispersion in the equation of motion for dust. Neglecting c d in Equation (C3) makes 1 − F slightly larger, and a ring appears more unstable to secular GI. Figure 12 shows dust surface density profiles at the final time steps from the six runs. We overplot filled circles whose color shows whether the condition for secular GI is satisfied at the dust cells or not. Secular GI is operational at red-colored regions where 1 − F is positive (Equation (C3)). In some runs, we find that secular GI becomes unstable at the dust rings.
The rings forming at outer radii can be sufficiently unstable to secular GI while inner rings tend to be stable or marginally unstable even when Σ d increases by an order of magnitude (e.g., see the panels of the q2T20a5 and q2T20a1 runs). This is partly due to difference in the turbulence strength α. Dust diffusion with larger α stabilize secular GI. Thus, less turbulent disks are still preferable even when dust surface density is increased by coagulation instability. Toomre's Q value of the background gas disk also significantly affect the stability of the rings. From the definition of Toomre's Q, we have for T ∝ r −1/2 , Ω ∝ r −3/2 and Σ g ∝ r −q/2 , meaning that Q increases as r decreases for q = 1, 2, and 3. As already mentioned in the previous studies, secular GI grows more easily for smaller Q (e.g., see Youdin 2005aYoudin , 2011Takahashi & Inutsuka 2014;Tominaga et al. 2019) although Q < 1 is not required. In other words, secular GI becomes operational more easily at outer radii, and thus in outer rings. Note that secular GI is stable at dust-depleted regions (blue regions in Figure 12). Thus, coagulation instability and resulting dust concentration are important for the onset of secular GI.

D. ON THE EFFECT OF POROSITY OF DUST AGGREGATES
Although we assume the spherical compact dust in the present simulations, dust collisions will produce porous aggregates (e.g., Dominik & Tielens 1997;Suyama et al. 2008;Wada et al. 2008). It is known that the coagulation