Meso-scale Instability Triggered by Dust Feedback in Dusty Rings: Origin and Observational Implications

High spatial resolution observations of protoplanetary disks by ALMA have revealed many substructures that are providing interesting constraints on disk physics as well as dust dynamics, both of which are essential for understanding planet formation. We carry out high-resolution, 2D global hydrodynamic simulations, including the effects of dust feedback, to study the stability of dusty rings. When the ring edges are relatively sharp and the dust surface density becomes comparable to the gas surface density, we find that dust feedback enhances the radial gradients of both the azimuthal velocity profile and the potential vorticity profile at the ring edges. This eventually leads to instabilities on meso-scales (spatial scales of several disk scale heights), causing dusty rings to be populated with many compact regions with highly concentrated dust densities. We also produce synthetic dust emission images using our simulation results and discuss the comparison between simulations and observations.


Introduction
Over the past few years, observations obtained with the Atacama Large Millimeter Array (ALMA) have revealed smallscale structures in the distribution of (sub)millimeter dust grains around many circumstellar disks surrounding nearby young stars. The majority of the millimeter bright disks observed so far feature multiple rings and gaps characterized by radii spanning from a few to a few hundreds of au (ALMA Partnership et al. 2015;Andrews et al. 2016Andrews et al. , 2018aAndrews et al. , 2018bIsella et al. 2016Isella et al. , 2018Fedele et al. 2018;Guzmán et al. 2018;Huang et al. 2018aHuang et al. , 2018bPérez et al. 2018a), while a smaller fraction of objects show prominent dust crescents (Fujiwara et al. 2006;Isella et al. 2010Isella et al. , 2013van der Marel et al. 2013;Pérez et al. 2014Pérez et al. , 2018bvan der Marel et al. 2016;Boehler et al. 2017Boehler et al. , 2018Benisty et al. 2018;Cazzoletti et al. 2018) and spiral arms (Grady et al. 2012;Benisty et al. 2015;Akiyama et al. 2016;Boehler et al. 2018;Dong et al. 2018b;Huang et al. 2018c;Uyama et al. 2018).
The origin of these structures is debated, but it is commonly understood that they might carry key information about the formation of planets. One hypothesis is that these features originate from the interaction between newly formed planets and the circumstellar material (Ou et al. 2007;Fung et al. 2014;Dipierro et al. 2015;Dong et al. 2015;Jin et al. 2016;Liu et al. 2018b;van der Marel et al. 2018;Jin et al. 2019). Numerical simulations indicate that azimuthally symmetric rings might result from the gravitational torque exerted by planets with masses below that of Jupiter (Bryden et al. 1999(Bryden et al. , 2000Bae et al. 2017;Dong et al. 2017Dong et al. , 2018a, while crescents might be caused by the onset of the Rossby Wave Instability (RWI) at the unstable gap edge of more massive planets (Lovelace et al. 1999;Li et al. 2000Li et al. , 2001Li et al. , 2005Fu et al. 2014a;Lovelace & Romanova 2014;Huang et al. 2018d). In an alternative scenario, rings and crescents form from radial variations in the disk viscosity and, consequently, in the gas radial drift velocity, which, in turn, leads to the accumulation of dust and gas in circular rings (Regály et al. 2011(Regály et al. , 2013Flock et al. 2015;Huang et al. 2019). If the viscosity variation is steep, the rings might become Rossby Wave unstable and generate gas vortices capable of trapping dust particles (Miranda et al. 2017, hereafter M17). A third hypothesis is that rings form at the snowline of volatile molecules (Zhang et al. 2015(Zhang et al. , 2016Yen et al. 2016). However, this scenario does not seem to be consistent with the rather random distribution of the rings' radii (Huang et al. 2018a;Long et al. 2018). Takahashi & Inutsuka (2014) proposed that the mutual friction and self gravity of dust and gas can lead to a secular gravitational instability (SGI), then form dusty rings in protoplanetary disks.
Regardless of their physical origin, it has been shown that dusty rings and crescents are places where the density of submillimeter-sized dust particles is enhanced compared to that of the circumstellar gas Boehler et al. 2017Boehler et al. , 2018Dullemond & Penzlin 2018). For example, in the case of HD142527 (Boehler et al. 2017), the dust-to-gas ratio (hereafter D2G) within the prominent crescent that characterizes this disk is so large (∼1.5) that we expect the dynamics of the system to be influenced by the dust component. When this happens, instabilities such as the streaming instability (SI; Youdin & Goodman 2005;) are likely to grow and lead to the formation of dense dust clumps, and, perhaps, planetary embryos. Other systems featuring crescents and rings show more moderate variations in the D2G, but they might nevertheless influence the overall dynamics of the disk. For example, the D2G in the crescents observed around MWC758 reaches a maximum of 0.2, while in the multiple ring systems observed around HLTau and HD163296, the D2G is constrained to be as high as about 0.05 (Jin et al. 2016;Isella et al. 2016), i.e., five times higher than the typical D2G in the interstellar medium. Pinte et al. (2016) showed that HLTau has a local D2G larger than 0.1 at the mid-plane.
The effect of a high D2G on the evolution of a vortex generated by RWI induced by the planet-disk interaction has been studied by Fu et al. (2014a), who found that dust feedback might lead to the destruction of the vortex itself, and, consequently, to the corresponding dust crescent, in about 10 3 -10 4 orbits of the planet responsible for the perturbation. For example, a dust crescent generated by a massive planet orbiting at 10 au around a solar-mass star might live for less than about 3×10 5 yr. Alternatively, M17 found that a crescent caused by a sharp radial decrease of the gas viscosity might have a lifetime comparable to that of the circumstellar disk itself, thanks to the continuous replenishing of gas accreting inward from the outer disk regions. Despite the differences in the predicted vortex lifetime, these studies indicate that a high D2G might lead to the fragmentation of dusty rings into dense clumps.
In this work, we extend the analysis of the effect of dust feedback on meso-scale disk structures by investigating the evolution of dusty rings. Similar to M17, we generate gas and dusty rings by imposing radial variations of the disk viscosity with a low viscosity region mimicking the dead-zone. Different from M17, however, we set the viscosity gradient so that the resulting radial density profile after a few thousand orbits is not steep enough to trigger RWI for the gas component. The presence of the gas surface density bump causes the dust to continue accumulating in this region, and then form a dusty ring. Ultimately, as the D2G approaches unity in the dusty rings, we find that instabilities get excited due to the dust feedback on the gas, and in certain cases, a large number of dust clumps are formed. In this work, we investigate the conditions for such an instability to develop, study the properties and evolution of the resulting clumps, and discuss the possibility of observing such structures using ALMA.
The outline of this paper is as follows. In Section 2, we describe the setup of the hydrodynamic model and the radiative transfer simulations. In Section 3, we present the hydrodynamical results of the dust and gas, and we discuss the dependence of instabilities on the adopted initial conditions. In Section 4, we present the synthetic ALMA images of our simulations. A discussion of the results in the framework of current and future observations is presented in Section 5, and the conclusions are summarized in Section 6.

Hydrodynamic Model
We calculate the temporal evolution of the circumstellar gas and dust using the multi-fluid hydrodynamic code LA-COM-PASS (Li et al. 2005(Li et al. , 2008Fu et al. 2014b , respectively. Initially, the gas disk surface density S g follows the viscous disk profile (Andrews et al. 2010): 200 au c 0 . The initial gas surface density at R 0 is S = -0.312 g cm 0 2 and the total gas mass is 3 Jupiter masses (M J ). At the beginning of the simulation, the dust surface density (S d ) also follows Equation (1) with a D2G of 0.05, corresponding to a total dust mass of about 45 Earth masses ( Å M ). The adopted initial D2G is larger than the typical value of 0.01 but is consistent with the observations of disks in Lupus (Ansdell et al. 2016). Furthermore, as discussed in Section 5 and the Appendix, we believe that the choice of the initial D2G has no strong effect on the results of our analysis.
We calculate the disk evolution by neglecting the disk self gravity and any effects due to magnetic fields. We treat the dust and gas as fluids controlled by the following continuity and momentum equations (Fu et al. 2014b (Frank et al. 2002, Chapter 4), while F drag is the gas drag given by where r p and a p are the internal density and radius of the dust particles, respectively. The last term of the right-hand side of Equation (5) is the dust feedback term. The dust feedback term usually is weak in interstellar mediums because of the low D2G ( = S S D2G d g in 2D and r r = D2G d g in 3D, where r d,g are the densities of the dust and gas, respectively). However, in protoplanetary disks, dust will be settled around the mid-plane until it is balanced by turbulence stirring, and then it will form thin dust sublayers (Chiang & Youdin 2010;Testi et al. 2014). Dust will be trapped in pressure bumps to form dusty rings. Or gas can be removed by photoevaporation (Alexander et al. 2006a(Alexander et al. , 2006b and MHD disk wind (Blandford & Payne 1982;Bai 2013;Bai & Stone 2013;Bai et al. 2016). All of these effects can increase the D2G in protoplanetary disks, enhancing the role of dust feedback.
In our analysis, we assume a single dust species characterized by r = -2.0 g cm p 3 and = a 0.25 mm p . The initial St at R 0 is about 0.25. The justification of this choice will be discussed in Section 5.1.1. We evolve the system assuming a vertically isothermal equation of state, = S P c g s 2 g , where P g and c s are the vertically integrated pressure and sound speed of gas. The radial profile of the sound speed is described as . The assumed sound speed corresponds to a disk aspect ratio H 0 /R 0 =0.05, where H 0 is the pressure scale height of disk ( = W H c s K , W K is the Keplerian angular frequency). The disk temperature profile corresponding to the adopted sound speed is = . The features of instabilities are scaled by the scale height of disks. In order to be compared with the results of M17, we use the same cold disk ( = R 50 au 0 and H 0 /R 0 =0.05) in this study. The luminosity of star is about 0.1 L e of the star in RADMC simulations, and it agrees with the disk temperature calculated above.
In our hydrodynamic simulations, the kinematic viscosity is parameterized as n a = c H s , where α is the dimensionless viscous parameter (Shakura & Sunyaev 1973;Pringle 1981). We further assume that α depends on the radius following a prescription designed to mimic the presence of a low viscosity region in the disk similar to that predicted by magnetorotational instability (MRI; Balbus & Hawley 1998) dead-zone models (Gammie 1996;Armitage 2011). In this paper, we consider the ohmic dead-zone viscosity transition model, for which the Hall effect (Bai 2014a(Bai , 2014b and ambipolar diffusion (Bai & Stone 2011) are ignored. Following Regály et al. (2011), we assume that the radial profile of α in the presence of a deadzone is parameterized as where α 0 =10 −2 and a = -10 DZ 5 are the viscosity parameters outside and within the dead-zone, respectively.
The transition between low and high viscosity happens at = = R R 1.5 75 au DZ 0 , and the width of the transition zone is controlled by the parameter D DZ . The transition from high to low viscosity leads to a corresponding formation of a gas over-density between R 0 and R DZ . This over-density leads to a local maximum in the gas pressure bump that is capable of forming a dusty ring by trapping dust particles as described in Weidenschilling (1977).
Numerical simulations have shown that if D < H 2 DZ , this gas bump becomes Rossby wave unstable and develops a largescale vortex (M17, Lyra et al. 2009;Regály et al. 2011). Since the focus of this paper is on azimuthally symmetric rings, we avoid the formation of a large-scale vortex by adopting D  H 2 DZ . More precisely, we explore three cases in which D DZ is equal to 2H, 3H, and 4H. The corresponding radial profiles of α are shown in Figure 1.
The hydrodynamic simulations are performed on a 2D numerical grid characterized by 6144 cells both in the radial and azimuthal directions, and open boundary conditions. We use log-scale grids along the radial direction. The aspect ratio d dF R R ( ) is about 0.5 at R 0 . The disk pressure scale height is resolved by about 90 cells along radial direction at R 0 and =d W R c St 1 5 1 20 K s at R 0 . In order to speed up the hydrodynamic simulations, we run the first 1000 orbits of the 3H model and the first 1500 orbits of the 4H model using the 1D version of LA-COMPASS, and switch to 2D when the D2G in the dusty rings approaches unity. A low level of noise is added to the 2D flow to seed nonaxisymmetric instabilities. Performing the 1D simulation does not affect our results since no instabilities develop at the beginning of the simulation in these two models.

Radiative Transfer and Synthetic Images
The outputs of LA-COMPASS 2D hydrodynamic simulations are post-processed using the radiative transfer code RADMC-3D 8 (Dullemond 2012) and the Common Astronomy Software Applications (CASA) 9  (McMullin et al. 2007) to generate synthetic images of the dust continuum emission at the wavelength of 1.3 mm (230 GHz, ALMA Band 6). The adopted procedure is similar to that described in Jin et al. (2016) and Liu et al. (2018a). In brief, we convert the 2D gas surface Figure 1. Radial profile of the viscosity parameter α for the three models studied in this paper. The blue, green, and red lines correspond to models in which the transition between the low viscosity inner disk and the high viscosity outer disk D DZ is equal to 2H, 3H, and 4H, respectively, where H is the pressure scale height of the disk. density (S g ) into a 3D density distribution (r g ) assuming hydrostatic equilibrium for a vertically isothermal disk such that r (Frank et al. 2002, Chapter 5).
We then assume that the disk temperature is controlled by submicron-sized dust particles that are well coupled to the gas and, therefore, follow the gas density distribution scaled for a constant D2G equal to 0.01. The disk temperature is calculated using RADMC-3D adopting a dust opacity proper of interstellar grains made of astronomical silicates (Weingartner & Draine 2001), organic carbonates (Zubko et al. 1998), and water ice with fractional abundances as in Pollack et al. (1994) and an MNR grain size distribution (n(a)∝a −3.5 ) with a minimum grain size of 0.5 μm and a maximum grain size of 10 mm. The calculation of dust opacity is discussed in Isella et al. (2009). The disk temperature of the mid-plane at R 0 calculated by RADMC-3D and that corresponding to the sound speed adopted in the hydrodynamic simulations differ by about 5%, and we believe that the difference in temperature is sufficiently small that it should not affect our results.
After the temperature has been calculated, we use the raytracing method of RADMC-3D to generate synthetic images of the dust continuum emission at the wavelength of 1.3 mm. To do this, we convert the surface density of 0.25 mm grains generated by the hydrodynamic simulations into a 3D density distribution as above, but we now assume that the pressure scale height of the dust disk is 10% of that of the gas (Pinte et al. 2015;Isella et al. 2016). We calculate the dust scale height H d varied with viscosity α (Birnstiel et al. 2010), then we get thatH H 0.01 0.05 d g  in our models. Because the locally isothermal approximation is applicable to the mid-plane and the face-on (i=0°) images are considered here, the effect of α on dust scale height can be neglected. In a more realistic 3D situation, the viscosity might be anisotropic between the vertical and radial directions . The dust absorption opacity at (sub)millimeter wavelengths is assumed to be ń -2.3 cm g 230 GHz 2 1 0.4 ( ) (Beckwith et al. 1990), and no dust scattering is included in this study. The disk is assumed to be at a distance of 140 pc from Earth.
Finally, synthetic disk images are processed using the "simobserve" task of CASA to simulate ALMA observations that properly account for observational noise. We adopt the ALMA configuration 24 that delivers a synthetic beam FWHM of 0 039×0 036 when a Briggs robust parameter of 0.5 is adopted in the image deconvolution. We integrate for 3 hr under typical weather conditions, resulting in an rms noise of 6.1 μJy beam −1 .

General Properties
The general behavior of the gas and dust evolution for disk models with different values of viscosity transition widthΔ DZ is presented in Figure 2. The imposed radial variation of the viscosity parameter α (Equation (8)) leads to the formation of local maxima in the gas pressure between R/R 0 =0.8 and 1.1 (panel (a)). As time passes, dust particles drifting inward from the outer disk regions are trapped within the gas pressure maxima, resulting in a strong radial and temporal variation of the D2G (panels (b) and (c)). At the end of the simulation, the azimuthally averaged D2G at the position of the dusty rings approaches unity in all three models (panels (b) and (c)).
The temporal evolution of the D2G strongly depends on the radial profile of α. In the 2H (D = H 2 DZ ) model, the maximum value of the D2G increases above 1 after about 300 orbits and shows large temporal oscillations caused by the generation and dissipation of small dense clumps of dust within the gas ring (Figures 3 and 4). We will discuss these structures in detail in the next section. Conversely, the maximum of D2G in model 3H remains below about 2 for the entire simulation and shows small temporal variation, while, in model 4H, the maximum D2G settles around a value of 0.7 after reaching unity after about 500 orbits. The total mass of the solids also evolves with time, and it decreases from the initial value of 45 Å M to about 35 Å M , 42 Å M , and 43 Å M at the end of simulation in the 2H, 3H, and 4H models, respectively. Note that all of the dust particles in the outer disk are collected in the ring. Figures 3-6 show the temporal evolution of the gas and dust surface density in the 2H, 3H, and 4H models. A common feature among models is the formation of ripples at the outer edge of the dusty rings, which, as discussed below, we attribute to the emergence of a meso-scale instability. Both the timescales for exciting the instability and their amplitudes strongly vary from model to model based on the steepness of the viscosity transition. The origin of these features is intrinsically related to the feedback that dust particles induce in the gas dynamics.
If dust feedback is not accounted for, dust grains would be passively transported by the gas and would concentrate within the gas pressure bump until balanced by dust diffusion (Takeuchi & Lin 2002). Without the dust feedback, the dusty ring would be very narrow and dense and would become gravitationally unstable if self gravity is included. The maxima of the D2G would easily exceed unity. Conversely, if dust feedback is considered, the gas dynamics will be significantly affected by solid particles as soon as the D2G approaches unity. This leads to a substantial change in the distribution and kinematics of the gas that might hinder building up of pressure bump (Taki et al. 2016). And the widths of the dusty rings become broader (Kanagawa et al. 2018), including dust feedback. Our simulations indicate that dust feedback leads to unstable rings and might result in the formation of dense clumps. The details of these instabilities are discussed in the following sections.

2H Model
Figures 3 and 4 present the temporal evolution of the gas and dust surface density in the 2H model. Around orbit 260, the D2G at the outer boundary of the dusty ring reaches unity, and five billow-shaped features emerge in the dust distribution. Around orbit 600, both the inner and outer boundaries of the , and 1740 (right panel) normalized by initial gas surface density Σ 0 . Bottom panels: snapshots of D2G at the same orbits. At orbit 260, an instability occurs at the outer edge of the dusty ring. At orbit 600, small dusty vortices (dust clumps) form nearby the outer edge of the dusty ring while its inner edge becomes unstable. At orbit 1740, more than 50 vortices characterized by D2G above unity are present along the dusty ring. The scale between the radial and azimuthal directions is about 1:10. We will discuss the properties of the dust clumps in Section 5.2 and Figure 16. dusty ring become unstable and feature the formation of small vortices that grow in number as the dust trapping becomes stronger. At orbit 1740, more than 50 dust clumps with D2Gs as high as 50 are visible along the dusty ring. The dust clumps have masses between 0.1 and . They are roughly elliptical in shape and have a characteristic radial diameter of 0.6 au (as discussed below) or one-fifth of the disk scale height. Our simulations resolve the dust clumps with about 20 and 40 cells along radial and azimuthal directions, respectively. As time proceeds, the strong dust feedback destroys the vortices (Fu et al. 2014b;Crnkovic-Rubsamen et al. 2015), and, after about 2740 orbits, the dusty ring returns to being almost featureless with a D2G around 2. After that, the disk undergoes a new phase of instabilities leading to the generation and dissipation of new small vortices. In the last part of our simulation, between orbits 3000 and 4000, the dusty ring features dense dust clumps and feather-like structures that extend over 180°along the azimuthal direction.

3H and 4H Models
Figures 5 and 6 present the temporal evolution of the gas and dust surface density for models 3H and 4H, respectively. As shown in Figure 2, shallower variations of α result in pressure bumps that are both shallower and located closer to the central star. In both models, the outer edge of the dusty ring becomes unstable. However, the instability grows at a slower pace and to a lower amplitude compared to the 2H model.
In model 3H, billow-shaped features appear at the outer edge of the ring after about 1780 orbits; however, the instability never leads to the formation of small-scale vortices as those observed in Figure 3. Furthermore, no instabilities develop at the inner edge of the dusty ring. In model 4H, weak billow-shaped features appear around orbit 1780 and remain quasi-stable for the entire length of the simulation. At orbit 4000, the azimuthally averaged D2Gs in the 3H and 4H models are 0.8 and 0.6, respectively.

Origin of Instability
To understand the instability revealed in the simulation results described above, we studied the gas and dust velocity profiles in more detail, which is shown in Figure 7. A circumstellar disk characterized by a monotonically decreasing profile of the gas surface density and temperature typically rotates at a slightly sub-Keplerian velocity (Adachi et al. 1976;Weidenschilling 1977). For the assumed initial disk parameters, the gas orbital velocity of the Keplerian velocity V K . As a radial maximum in the gas pressure forms as the result of the assumed viscosity profile, the gas located inside the pressure maximum is forced to rotate faster than F V g, , while the gas outside the pressure bump slows down. This variation in gas velocity forces dust particles to drift toward the pressure maximum and the D2G to grow. Neglecting the feedback of dust particles on the gas dynamics, and without any other process capable of halting the dust drift,   Figure 3, but for the 4H model. The panels from left to right correspond to orbits 1000, 1780, 2000, 3000, and 4000. Weak billow-shaped features form at the outer edge of the dusty ring around orbit 1780. As in the 3H case, no dense dust clumps form in this model. solid particles will smoothly concentrate in a very narrow dusty ring without affecting the rotational velocity of the gas. However, when the dust feedback is included, as the dust density increases, solid particles push the gas to rotate at Keplerian velocity V K , creating inflection points in the gas velocity and vorticity profile. As the D2G within the dusty ring approaches unity, the gas in the ring will be "forced" to rotate almost at the Keplerian velocity, while the gas just outside the ring, where the D2G is much smaller, will still move at sub-Keplerian velocity. This leads to a steep velocity shear that makes the outer edge of the dusty ring unstable. A similar process happens at the inner edges of the dusty rings of 2H, while the gas rotates as super-Keplerian and Keplerian.
The potential vorticity (PV; vortensity) profiles, where , for 2H are shown in Figure 8. The same profiles for the 3H and 4H models are shown in Figure 9, respectively. Figure 8 presents the development of the instability (orbit 240, 260, 2740) before it becomes completely nonlinear (orbit 600, 1740, 3060). There are local vortensity extremums at the edges of the dusty rings. Ultimately, the velocity knees at the edges are smoothed out, eventually leading to a relative quiescent phase around orbit 2740. Similar behaviors are captured for the 3H and 4H models as shown in Figure 9. Figures 7-9 indicate that the dust feedback is essential in producing the sharp features in the PV profiles near the edges of the dusty rings in all cases. We speculate that such profiles will render the flow unstable to the RWI (Lovelace et al. 1999;Li et al. 2000Li et al. , 2001Ono et al. 2016Ono et al. , 2018. The existence of a local extremum in the PV profile is a necessary condition for RWI. In fact, Figures 3-6 show that the linear development of unstable modes have an azimuthal mode number between 3 and 5, similar to the predictions given in Li et al. (2000) for the most unstable modes. We plan to carry out a more detailed linear theory analysis of this instability in a future study. For now, we will keep referring to the instability as a "meso-scale" instability because the features generated by the instability have length-scales ranging from ∼H to tens of H.
Density waves excited by planets (Goldreich & Tremaine 1978, 1979, 1980 carry angular momentum flux (Goodman & Rafikov 2001;Rafikov 2002aRafikov , 2002bDong et al. 2011aDong et al. , 2011bBae & Zhu 2018aLi et al. 2019a;Miranda & Rafikov 2019a while propagating. When the density waves steepen into shock waves, the angular momentum flux is transferred into the disk, and gaps can be opened in an inviscid disk, eventually. Because gap edges are also pressure bumps, similar instabilities could happen at dusty rings induced by planets, including dust feedback. Pierens et al. (2019) showed interactions between a low-mass planet and an inviscid pebblerich disk. They found that an RWI-like instability occurs at gap edges of the planet, generating dusty vortices. Yang & Zhu (2019) studied planet-disk interactions including dust feedback. They showed that large-scale lopsided vortices on either gap edge of planets would be broken into numerous small dusty vortices. When the small dusty vortices dissipate, the large-scale vortices reappear again. We believe that the instabilities found by them are the same as the meso-scale instability described here.

Synthetic ALMA Images
In this section, we investigate whether ALMA high angular resolution observations of nearby circumstellar disks could detect some of the small-scale structures resulting from the instability. Figures 10 and 11 show synthetic images of different snapshots for the 2H, 3H, and 4H models in the 1.3 mm dust continuum emission calculated as discussed in Section 2. The spatial resolution of the observations is about 5 au and is comparable to the radial extent of the ring. The rms noise is m -6.1 Jy beam 1 , and the emission is detected with a peak . The x-axis indicates the orbital radius normalized to the radius of the dusty ring forming due to the converging migration of solid particles toward the gas pressure maximum. The y-axis indicates the azimuthal velocity of gas normalized to the local Keplerian velocity (gray dashed lines). At the beginning of our simulation, the gas pressure decreases with the orbital radius. This results in a slightly sub-Keplerian gas rotational velocity (black). As the pressure bump develops, dust is trapped in a radially narrow dusty ring (gray dashed-dotted lines). The gas inside and outside the dusty ring rotates slightly faster (red) and slower (blue) than Keplerian, respectively. Without dust feedback, the velocity of the gas across the dusty ring is not affected by the increasing D2G. When the dust feedback is included, solid particles trapped inside the dusty ring push the gas to move at Keplerian velocity (green). The steep velocity gradient at the outer edge of the ring leads to the onset of the instability. This instability at the inner edge of the dusty rings in 2H happens later than that at outer edge, because it needs to take some time to produce a strong enough gradient at the inner edge.

Intensity of the Continuum Emission
signal-to-noise ratio (S/N) of 20∼25. At this sensitivity and resolution, none of the billow-shaped and feather-like features that characterize the edges of the dusty ring in the 2H model are visible. Also, none of the dense dust clumps are singularly resolved. However, the synthetic images show azimuthal variations in the intensity profile related to the instability in the 2H model. Models 3H and 4H are characterized by a lower level of instability, which lead to more azimuthally symmetric dusty rings. At orbit 1060 in the 3H model, the azimuthal variations in the intensity marginally trace the location of the billow-shaped structures present at the outer edge of the ring. For all of the other orbits in both models, the intensity variations are consistent with the noise of the observations.
More detailed analysis of the 2H model reveals the following features: At orbit 260, the azimuthal intensity profile shown in Figure 12 oscillates between ±20% of the mean intensity, in accordance with the position of the billow-shaped features shown in Figure 3. Essentially, at the azimuthal positions where these features occur, the dusty ring is more extended radially, and its radially integrated flux is larger. We argue that a power spectrum analysis of the azimuthal intensity profile of known dusty rings might reveal the presence of the instability. At orbit 600, the azimuthal intensity profile shows variations as large as 50% that are caused by the nonuniform distribution in azimuth of the dense dust clumps. Indeed, at this orbit, the dust clumps are organized in two clusters centered around Φ=30°F igure 8. Potential Vorticity (PV) of 2H. PV is normalized by code units, and the Keplerian component is subtracted. The x-axis indicates radial distances from the center of the dusty ring divided normalized by ring radius. We choose the frames without instability ("2H 240") in the linear regime ("2H 260" and "2H 2740"), and in the nonlinear regime ("2H 600," "2H 1740" and "2H 3060"). There are local vortensity minimums at the edges of the dusty rings, and then instabilities happen at the edges. When instabilities enter nonlinear, lots of small anticyclonic (negative PV compared with background flows) vortices emerge. and Φ=250°, containing about 10 dust clumps each. At orbit 1740, the dust surface density is characterized by a large number of clumps, which are more uniformly distributed along the ring compared to the previous case, but nevertheless lead to intensity variations of about 30% of the mean intensity. A similar behavior of the intensity profile is observed at orbits 2740, 3060, and 4000, in which amplitude variations correlate to the amount of instabilities in the disk and can lead to a contrast of intensity along the ring ( ( ( )) up to 4, as for orbit 3060.

Optical Depth
We calculate the fraction of dust mass trapped inside regions that are optically thick in the dust continuum emission at 1.3 mm. Following the usual procedure adopted to derive the optical depth of the continuum emission from ALMA observations, we calculate the optical depth τ of our synthetic intensity maps from the relation = - K is the dust temperature inside the dusty ring, and B ν is the Planck function. A clear advantage of working with simulated data is that we know the true value of the dust temperature, while in the case of real observations, the dust temperature is generally estimated based on radiative transfer models for dust or observations of optically thick molecular lines. Note that the optical depth calculated using the previous equation is an average across the synthesized beam. In the presence of structures smaller than the beam, τ would therefore differ from the "true" optical depth of the emission as calculated by multiplying the dust surface density map by the dust opacity. The comparison between the optical depth calculated from the synthetic images and that derived from the surface density map provides a measurement of the beam dilution. Figure 13 shows the optical depth for the 2H, 3H, and 4H models calculated using Equation (9). In the orbit 260 frame of model 2H, the optical depth of the dusty ring varies between 0.4 and 0.7. As dense dust clumps form as a result of the instability, τ increases to a maximum of 1.3 at orbit 600. This value should be compared with the optical depth at the same position calculated from the surface density, which, in this specific case, is 20. This implies that the limited angular resolution of the observations leads to a dilution of intensity of about a factor of 15, or that the structures that emit most of the emission at that specific location in the disk have a size about 15 times smaller than the angular resolution of the observations. The 3H and 4H models have lower dust surface densities than the 2H model. However, because the dust is more uniformly distributed, the optical depth inferred form the simulated observations is larger than in the 2H model. The optical depths of the 3H and 4H models are about 0.6∼1.1 and 0.8∼1.2, respectively. The maximum optical depth of the 3H model as calculated from the dust surface density map is about 2.2, implying that beam dilution plays a much smaller role than that in the 2H models. This is consistent with the lack of dense clumps in this model. Even smaller is the effect of beam dilution in the 4H model, where the maximum "true" dust optical depth peaks at 1.8.
A direct consequence of the presence of optically thick and spatially unresolved dust clumps is that the dust mass estimated from the observed intensity will underestimate the true amount of dust trapped in the dusty rings. For example, the synthetic observations of model 2H provide a total dust mass between 11 and 18 Å M , while the true dust mass trapped within the ring is between 20 and 35 Å M in the 2H model. In the case of the 3H model, the dust mass inferred from the observations is between 18 and 27 Å M , to be compared to a true dust mass between 30 and 42 Å M . Finally, for model 4H, the mass estimated from the observations is between 35 and 37 Å M , to be compared with a true mass between 37 and 42 Å M . These results stress the importance of imaging dusty rings at the highest angular resolution possible to mitigate the effect of beam dilution, and at wavelengths longer than 1 mm to reduce the limiting effect of the dust optical depth.

Dependence on Model Parameters
The asymmetries created by this meso-scale instability rely on several assumptions concerning the origin of the radial dust trap, the initial dust and gas densities, the dust grain size and opacity, the equation of state, the lack of self gravity, and the fact that we use 2D simulations whereas circumstellar disks are Figure 9. Similar to Figure 8, but for the 3H and 4H models. We choose the frames without instability ("3H 1040" and "4H 1760") and those with instability ("3H 1060" and "4H 1780"). This instability only happens at the outer edges of the dusty rings.
3D objects. Here, we discuss how some of these assumptions might affect our results and present ideas for future follow-up studies.

Dust Size and Stokes Number
Dust feedback is a necessary condition for the instability found in this study. The dust feedback depends on the coupling between the gas and dust, which, in turn, is controlled by the dust Stokes number. Our simulations assume a single grain size = a 0.25 mm p , corresponding to a Stokes number of about 0.25 at R 0 . As a result of the radial evolution of the gas density, the Stokes number varies with time and radius and, at the end of the simulations (at orbit 4000), reaches values between 0.05 and 0.15 across the dusty ring ( Figure 14).
It is important to understand how different sizes of the particles contribute to the dust feedback. To do that, we adopt the approach of Takeuchi & Lin (2002) and express the difference between the azimuthal velocity of the gas and dust as where V d,R is the radial velocity of the dust Figure 10. Synthetic ALMA images of specific orbits at 1.3 mm wavelength for the 2H model. The intensity is normalized by signal-to-noise ratios (S/N). The x-and y-axes present the relative location from the central star along R.A. and declination, respectively. The synthetic beam is marked as a white ellipse. The FWHMs of the synthetic beam are 0 039 and 0 036, i.e., the spatial resolution is5.5 5.0 au.  where n(a)∝a −3.5 between a min =0.1 μm (St=10 −4 at R 0 ) and a max =1 cm (St=10 at R 0 ) is a typical grain size distribution for circumstellar disks. We numerically integrate Equations (13) and (15), then we calculate the characteristic particle size by solving the relations drag, . As Figure 15 shows, the drag force along the radial direction is dominated by the small Stokes number (small dust size), and the drag force along the azimuthal direction is dominated by the dust with St;0.27. For the assumed initial profiles of α and c s , we find that the characteristic Stokes number of 0.27 along the Φ direction corresponds to a grain size of 0.27 mm, which is very close to the one adopted in our simulations. Initializing the simulations with a larger grain size (larger than 0.27) would have reduced the dust feedback and might have led to weaker instabilities.
Our findings are also consistent with some previous studies. The "ViscBroad" (Δ DZ =2H) model in Figure 4 of M17 is similar to our 2H model. However, the dusty ring in "ViscBroad" just presents weak asymmetries and does not show a similar instability compared with the 2H model. The initial Stokes number of dust in M17 is 1. This leads to a weaker back reaction of dust on gas, so they just get a weak Figure 11. Synthetic ALMA images of specific orbits at 1.3 mm wavelength for the 3H and 4H models. The beam size and the noise level are same as the 2H model in Figure 10. Because the instability is quite weak in the 3H and 4H models, we just choose the frames with instability occurring as well as the last frames of simulations. We can see that the rings in 3H and 4H are quite smooth compared with rings in 2H.
instability at the edges of dusty rings. This is consistent with our analysis above. Fu et al. (2014b) found that the dust will be concentrated effectively by Rossby Wave vortices, and dust feedback will destroy the vortices when the dust density is comparable and higher than the gas density. They used dust with a Stokes number from´-4 10 4 to 0.16 in their simulations. They found that the dust with a large Stokes number has a large back reaction on the vortices. This is consistent with the fact that the small vortices emerge and dissipate in 100 orbits of our simulations.
In this paper, we just use one size dust component in our simulations. The simulations on the full dust size distribution will be quite costly computationally. We do not consider dust coagulation and fragmentation in this paper, which could play an important role in the dust evolution in the protoplanetary disk. Tamfal et al. (2018) investigated the role of dust coagulation on the two-fluids hydrodynamical evolutions. They found that dust has different distributions in fixed/variable dust-size simulations. Drażkowska et al. (2019) ran full-sized dust coagulation and hydrodynamical simulations, and they found that the back reaction does not change the gas structure significantly. Li et al. (2019b) showed that the ring structures help the dust growth and produce a low spectral index. The instability in this paper might also generate some regions with low spectral indexes along dusty rings when we combine the dust coagulation model with the hydrodynamical model.

2D Approximation, Dust-to-gas Ratio, and Viscosity Profile
We have used 2D approximations in our simulations. The initial instability wavelength is significantly longer than H, though the transition width is of the order of H. Some of the nonlinear features generated in the 2H model can have sizes less than H. We use about 20 radial cells to capture these small vortices. On the other hand, features with size <H should be taken as indicative only because 2D approximations will not be fully consistent with such features. In our 2D simulations, the strong dust feedback can destroy vortices. This phenomenon is similar to previous studies in 2D. However, Lyra et al. (2018) showed that gas vortices can still survive under the dust disturbance in 3D, due to different vertical structures of the dust and gas.
In the 3D case, dust settles toward the mid-plane to form a dust layer. Gas beyond the dust layer rotates as sub-Keplerian, and the dust feedback lets gas in the dust layer rotate as Keplerian. The vertical shear rate of gas will be unstable and then trigger the Kelvin-Helmholtz Instability(KHI; Sekiya 1998; Youdin & Shu 2002). To combine this mesoscale instability and KHI in 3D simulations will be interesting in future studies.
In this study, we use three models with different viscosity transition widths to generate pressure bumps to trap dust. We find that the instability gets triggered when the D2G approaches unity. If a higher a DZ or lower α 0 is adopted, the Figure 12. Normalized intensity variation along the dusty rings for the 2H, 3H, and 4H models. The black dashed lines mark the variations equal to zero. The normalized intensity variation is defined as D = -I I I I ( ) , where I is the azimuthally averaged intensity along the dusty rings. The definition of the azimuthal angle Φ starts from the right in the counterclockwise direction (from the west to the north), the same as Figures 3-6, 10, and 11. The radial width of the azimuthal cut is about 1 radially Gaussian FWHM for each dusty ring. For convenience, we calculate the variations of rings from intensity images directly, instead of calculating visibility on the uv-plane. The amplitude of the fluctuations caused by noise is of the order of1 S N 0.1 ( ) , as the black dotted lines show in each panel.
viscous diffusion will let the pressure bump appear later. Because we use density exponential truncated models, the dust is barely supplied from the outer boundaries. We set a large viscosity transition and a high initial D2G; thus, the bumps are built up quickly before the dust drifts into the star. However, the relative viscosity and initial D2G of 0.05 (see the Appendix) do not change the main conclusions of this study. The D2Gs in the dusty rings in all three models exceed or approach unity. We infer that the critical D2G is about 0.5 based on the weak instability at the outer edge of the dusty rings in 4H. In actual disks, the settling dust can produce higher D2Gs. The processes that are studied here will likely be   applicable to the mid-plane region, though detailed 3D studies have to be carried out to test this scenario.

Dust Opacity, Optical Depth, and Dust Scattering
The dust absorption opacity at (sub)millimeter wavelengths in protoplanetary disks is commonly assumed as κ∝ν β , where 0<β1 (Beckwith et al. 1990), and we adopt β=0.4 in this paper. We calculate the optical depth of the rings in our synthetic images and find that the optical depth is between 0.4 and 1.4 at 1.3 mm wavelength with ∼0 04 resolution for our models. Note that the width of the dusty rings in 2H (∼3 au) is smaller than the beam size (∼5 au). The local optical depth of the dust clumps in the 2H model are optically thick without the dilution of the observational beam. The dusty rings will be optically thick at higher frequencies and optically thin at lower frequencies, respectively. When the dusty rings become optically thick, the intensity variation (see Figure 12) of the dusty rings will become small, and the dust clumps will be hard to distinguish.
The opacity of the dust can be affected significantly by the temperature, chemical compositions, grain sizes, structures, and topology of the dust particle (Semenov et al. 2003). The dust components in protoplanetary disks are different than the interstellar grains (Draine 2006). The uncertainties of dust opacity also comes from the opacity models and the numerical computing ). There could be an order of variation on the dust opacity for different models (Demyk et al. 2017a(Demyk et al. , 2017b. For an optically thin medium, the scattering can be ignored. In the situation with intermediate and high optical depth of dust, the scattering opacity at (sub)millimeter wavelengths can create a polarized emission of dust continuum (Kataoka et al. 2015) and an extinct part of the CO emission . Liu (2019) showed that the selfscattering of an optically thick dust disk can produce low (sub) millimeter spectral indices. Zhu et al. (2019) proposed that optically thick dust with scattering can lead to an underestimation of dust masses. The effects of scattering on interpreting our simulation images will be explored in the future.

Self Gravity, Gravitational Instability, and Streaming Instability
When the thermal pressure and the rotation cannot support the self gravity of the disk, gravitational instability (GI) will happen. For the gas disk, this can be characterized by the Toomre Q value (Toomre 1964;Goldreich & Lynden-Bell 1965): where κ is the radially epicyclic frequency (k º W R R d dR 2 1 4 2 3 ( )), and G is the gravitational constant. For a pure dust disk (such as Saturn's rings), the velocity dispersion σ between dust particles replaces the sound speed c s in Equation (18)  (Armitage 2007). When Q is smaller than 1, GI will happen if the self gravity is included.
Dust is trapped in the small vortices created by the instability where D2G is high (Σ d /Σ g 10), such as "2H 1740," "2H 3060," and "2H 4000" (see panel (c) of Figures 2-4). Because dust is tightly coupled with gas (St0.1 at dusty rings, see Figure 14), we can still use the sound speed c s to calculate Q as an approximation. Figure 16 describes the overall properties of the dust clumps developed at late time in the 2H model. The radial extents range from 0.3 au (∼0.1 H) to 1.1 au (∼0.6 H). The aspect ratio of the dust clumps is about 4. Since dust is tightly coupled with gas, the dust clumps have similar sizes to the gas vortices. The aspect ratio ∼4 is consistent with previous theoretical prediction (Lesur & Papaloizou 2009). These clumps contain a lot of dust between about Å M 0.1 and Å M 1.5 . The Toomre Q calculated by the medians of the dust clumps is then about 5, which is consistent with the calculations based on S = S + S d g ( ) (see Figure 17). The minimum of Q can be below unity, so the dust clumps might be gravitationally unstable, and planetesimals might be formed in dust clumps with the aid of GI. We did not include self gravity in our simulations; it will be interesting to explore the effects of GI to determine the properties of these dust clumps. Former studies (Goodman & Pindor 2000;Youdin & Goodman 2005; showed that dust can be concentrated under the dustgas mutual aerodynamic drag, leading to the local, linear Streaming Instability (SI). SI can generate high dust-density enhancement without self gravity. The instability found in this paper is on a much larger scale (at least initially), so we believe that this instability is different from the traditional SI. The nonlinear outcome seen in the 2H model, however, could interplay with the SI. In the simulations presented here, we do not have enough radial resolution per H to capture the SI, so presumably, the traditional SI is not excited here. The vertical motion of dust is also important for SI, so the interplay between this instability and SI still needs to be carried out in 3D simulations in the future.

Implication of ALMA Observations
Dusty rings are the most common features in the dust continuum of protoplanetary disks Long et al. 2018). Here, we discuss some protoplanetary disks with significant ring features observed by ALMA, in the context of the instability presented in this paper.

HD163296
HD163296 is a Herbig Ae star located at a distance of 101 pc (Gaia Collaboration et al. 2018). It has extended dust and gas disks around it . Isella et al. (2007) used SMA and VLA observations to show an asymmetric dust continuum and CO isotopomers depletion in ) for dusty rings of 2H. Q min presents the minimum of Q, and they can be below unity at the dusty rings. The self gravity will play an important role in the nonlinear regime of this instability. The Toomre Q values for the dusty rings of 3H and 4H are all above 40, so the self gravity can be ignored in 3H and 4H. the HD163296 disk. Isella et al. (2016) showed that there are three deleted dust gaps and reduced CO gas density in the dust gaps compared with the surroundings regions. The dust gaps and reduced CO density can be explained by the dynamical clearing of three embedded Saturn-like planets. Liu et al. (2018b) used three half-Jovian planets located at the dust gaps to give a reasonable fit to the dust continuum and CO gas emission. They also claimed that HD163296 presents a low viscosity (α<10 −4 ) in the inner part and a large viscosity (α∼7.5×10 −3 ) in the outer part of the disk. The viscosity transition is consistent with the traditional dead-zone model caused by ohmic dissipation (Gammie 1996). Isella et al. (2018) reported that there are three concentric dusty rings in HD163296, and they also found some asymmetric features in the dust emission of HD163296. There is a significant arc-like feature at the inner side of the dusty ring of 67 au (Figure 1 of their paper). There is about ±15% amplitude intensity variation relative to the mean ring intensity in the dusty rings of HD163296. Zhang et al. (2018) run a 0.15 M J planet-disk interaction model with a 2D hydrodynamic simulation and radiative transfer post-processing. They showed that the arc-like feature can be explained by dust trapping at the Lagrangian points along the orbit of planet. They also argued that the intensity variation of the dusty rings of HD163296 is caused by the perturbation of the planet on the gas kinematics in the disk.
From our simulations, we can see some transient bright clumps and crescent features along the dusty rings in our dust emission images of the 2H model (see "2H 1740," "2H 3060," and "2H 4000" in Figure 10). They differ from the arc-like feature of HD163296, which is located outside of the dusty ring. The large intensity variations for the 2H model might not be matched well with ±15% intensity variation in HD163296. Even though the purpose of this paper is not to give a fit on HD163296, we infer that the asymmetric substructure of HD163296 is probably not caused by the nonlinear outcome of this instability.

MWC758
MWC758 is a young (3.5 Myr; Meeus et al. 2012) Herbig Ae star with a circumstellar disk, located 160 pc away (Gaia Collaboration et al. 2018). Marino et al. (2015) showed that there are two bright bumps in the submillimeter continuum. Boehler et al. (2018) showed that MWC758 presents a complex morphology with the ALMA Band 7 dust continuum and 12 CO/ 13 CO gas line emission with an angular resolution 0 1∼0 2. Dong et al. (2018b) showed that there are three dusty rings, two bright bumps, and one eccentric cavity based on ALMA observations at 0.87 mm with 40 mas resolution. The contrasts of the southern and northern bright clumps are 4.4 and 10, respectively, compared with the averaged intensity at their radii.
Two bright bumps in MWC758 are located at different dusty rings. The contrast magnitude of the south bump is similar to the contrast from our simulations (∼4 in "2H 3060," See Section 2.2). There are also some asymmetries along the inner ring (see Figures 2 and 3 in Dong et al. 2018b). The northern bumps with high contrast (∼10) might come from dust trapping in a larger-scale vortex generated by RWI. Baruteau et al. (2019) showed that the bright bumps in MWC758 are the large-scale vortices of RWI triggered by two massive planets in gaps. They also claimed that if the inner vortex decays faster than the outer one, then this also can explain the formation of the eccentric inner ring. In our simulations, the dusty ring is also distorted in "2H 3060." The instability and its nonlinear outcome might provide some explanations for the complex features seen in MWC758.

HLTau
HLTau is a young TTauri Star located at about 140 pc from Earth (Gaia Collaboration et al. 2018). ALMA Partnership et al. (2015) used ALMA observations to present multiple bright and dark rings in HLTau. Their observations also showed the low spectral index in the rings, which implies that the dust might have evolved compared with the interstellar medium grains. Jin et al. (2016) calculated the optical depth at 1.0 mm of the dusty rings of HLTau, and they found that the bright rings (such as B1 and B2) are optically thick. Carrasco-González et al. (2016) used the VLA data on HLTau, and they showed that the bright rings of HLTau at the 7.0 mm are optically thin. There are some clump candidates at the B1 ring of HLTau (Figure 2 in their paper). The clump candidates of HLTau at 7.0 mm could be explained by this instability.
As discussed in Section 5.1.3, dusty rings will be smooth when they are optically thick. The bright rings in HLTau observed by ALMA are optically thick, so they are smooth and ordered at (sub)millimeter wavelengths. It will be difficult to assess whether there are dense clumps formed inside the dusty rings of protoplanetary disks at a single wavelength. Observations at multiple wavelengths are essential to confirm the asymmetries caused by this instability.

Summary
We have studied a meso-scale instability caused by dust-gas interaction at the edges of dusty rings by carrying out highresolution 2D two-fluid global hydrodynamic simulations, combined with the radiative transfer calculation and ALMA synthetic image processing. Our main conclusions are as follows: 1. If the conditions are favorable, there is an instability at the edges of the dusty rings produced by pressure bumps. Dust feedback alters the azimuthal gas velocity profile, producing a sharpened potential vorticity profile at the edges of the dusty rings. 2. This instability can develop when the local D2G approaches unity. A ring that is being continuously fed from the outer disk and dust settled toward the mid-plane can eventually accumulate enough dust, at least transiently. 3. The nonlinear outcome of such an instability can lead to many small vortices at dusty rings. These small vortices can contain a large amount of dust, reaching at least ∼10% of Earth mass within each dust clump. 4. These dust clumps are optically thick, and they might not be fully resolved by current observations at (sub) millimeters. However, these dust clumps can still lead to both nonaxisymmetric structures and intensity variations along dusty rings at specific wavelengths. This instability can provide some observational implications and explanations on the asymmetries in protoplanetary disks.
Future studies in 3D, including interactions between this meso-scale instability with other instabilities, such as SI and GI, will be quite interesting. In the future, the Next Generation Very Large Array (ngVLA) will be capable of 5 mas spatial resolution at 3 mm wavelength to detect substructures of protoplanetary disks in the closest star-forming regions (∼140 pc; Ricci et al. 2018). The asymmetries generated by this instability can be resolved by ngVLA. They will help to address whether dusty rings are suitable for planetesimal formation and eventually planet formation. in our runs, which is larger than the value (S S 0.01 d g  ) in the interstellar medium. To quantify the effect of the initial D2Gs on this meso-scale instability, we set another hydrodynamic runs with S S = d g init ( ) 0.01. The setup of this model is same as the 2H model above, except for the initial D2G. We run the simulation for 1000 orbits in 1D until the dust is trapped in the bump, then we turn it into 2D for another 200 orbits. Figure A1 shows the maximum D2G and one snapshot of the density profile. The maximum D2G varied with time has a similar pattern to that of 2H. The instability happens with an initial D2G equal to 0.01. The dusty ring has similar radial width compared with those of 2H. The initial D2G=0.05 does not change the main conclusions of this study, because there is dust trapping in pressure bumps in our simulations. This instability happened with a high D2G (0.5), which is larger than both 0.05 and 0.01. When the initial D2G is smaller, it just takes longer to build up the D2G at the ring. Then, the maximum of the D2G at the ring is determined by the competition between the build-up process and the diffusion. If the initial D2G is extremely low (for example, 10 −5 ), then the timescale to build up the ring via dust drift and the viscous diffusion to reduce the peak D2G at the ring will become comparable; so if the final D2G at the ring is quite low, then the 2D instability will never get triggered. Figure A1. The left panel shows the maximum of D2G varied with orbits of "2H" and "2H, D2G=0.01." The right panel shows one snapshot of the density profile of "2H, D2G=0.01."