text:
\section{Introduction}\label{sec:intro}
Over the past few years, observations obtained with the Atacama Large Millimeter Array (ALMA) have revealed small scale 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 astronomical units (au) \citep{almapartnetship2015, andrews2016ringed, isella2016ringed, andrews2018science, huang2018co, fedele2018alma, perez2018first, andrews2018disk, huang2018disk2, guzman2018disk, dullemond2018disk, isella2018disk}, while a smaller fraction of objects show prominent dust crescents \citep{fujiwara2006asymmetric,isella2010millimeter,isella2013azimuthal,van2013major,perez2014large,van2016vortices,boehler2017close,benisty2018shadows,cazzoletti2018evidence,perez2018disk,boehler2018complex} and spiral arms \citep{grady2012spiral,benisty2015asymmetric,akiyama2016spiral,boehler2018complex,dong2018eccentric,huang2018disk3,uyama2018subaru}.

The origin of these structures is debated but it is common understanding 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 \citep{ou2007disk,fung2014empty,zhu2014dust,dipierro2015planet,dong2015observational,jin2016modeling,liu2018new,van2018rings,jin2019new}. Numerical simulations indicate that azimuthally symmetric rings might result from the gravitational torque exerted by planets with masses below that of Jupiter \citep{bryden1999tidally,bryden2000interaction,bae2017formation,dong2017multiple,dong2018multiple}, while crescents might be caused by the onset of the Rossby Wave Instability at the unstable gap edge of more massive planets \citep[RWI,][]{lovelace1999rossby, li2000rossby, li2001rossby, li2005potential, fu2014long, lovelace2014rossby, zhu2014particle, huang2018identifying}. 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 \citep{regaly2011possible,regaly2013trapping,flock2015gaps,huang2019observability}. If the viscosity variation is steep, the rings might become Rossby Wave unstable and generate gas vortices capable of trapping dust particles \cite[][hereafter M17]{miranda2017long}. A third hypothesis is that rings form at the snow line of volatile molecules \citep{zhang2015evidence,zhang2016commonality,yen2016gas}. However this scenario does not seem to be consistent with the rather random distribution of rings radii \citep{huang2018co, long2018gaps}.

Regardless of their physical origin, it has been shown that dusty rings and crescents are places where the density of sub-millimeter dust particles is enhanced compared to that of the circumstellar gas \citep{isella2016ringed, boehler2017close, boehler2018complex, dullemond2018dust}. For example, in the case of HD~142527 \citep{boehler2017close}, the dust-to-gas ratio within the prominent crescent that characterizes this disk is so large ($\sim$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 \citep{youdin2005streaming,youdin2007protoplanetary,johansen2007protoplanetary} 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 dust-to-gas ratio, but they might nevertheless influence the overall dynamics of the disk. For example, the dust-to-gas ratio 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 dust-to-gas ratio is constrained to be as high as about 0.05 \citep{jin2016modeling,isella2016ringed}, i.e., 5 times higher than the typical dust-to-gas ratio in the interstellar medium.

The effect of a high dust-to-gas ratio on the evolution of a vortex generated by RWI induced by the planet-disk interaction has been studied by \cite{fu2014long}, 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\times10^5$ yr. Alternatively, M17 found that a crescent caused by a sharp radial decrease of the gas viscosity might have a life time 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 dust-to-gas ratio 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 ring causes the dust to continue accumulating in this region. Ultimately, as the dust to gas ratio approaches unity in this ring, 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~\ref{sec:setup}, we describe the setup of hydrodynamic model and the radiative transfer simulations. In Section~\ref{sec:hydro_results}, we present the hydrodynamical results of dust and gas, and discuss the dependence of instabilities on the adopted initial conditions.In Section~\ref{sec:images_results}, we present the synthetic ALMA images of our simulations. Discussion of the results in the framework of current and future observations are presented in Section~\ref{sec:discussion}, while the conclusions are summarized in Section~\ref{sec:summary}.

\section{Numerical Setup}\label{sec:setup}
\subsection{Hydrodynamic Model}
We calculate the temporal evolution of the circumstellar gas and dust using the multi-fluid hydrodynamic code LA-COMPASS \citep{li2005potential,li2008type,fu2014effects}. We assume a stellar mass $M_\star = \mrm{1\;M_\sun}$ and a normalization radius $R_0 = \mrm{50\;au}$, so that 1 orbit corresponds to a time interval of 353.5 years. The inner and outer boundaries of the numerical simulation are set to be $0.4\times R_0 = \mrm{20\;au}$ and $12\times R_0 = \mrm{600\;au}$, respectively. Initially, the gas disk surface density $\Sigma_\mrm{g}$ follows the viscous disk profile \citep{andrews2010protoplanetary}:
\begin{equation}
\label{eq:density}
\Sigma_\mrm{g} = 0.377 \left(\frac{R}{R_0}\right)^{-\gamma} \exp{\left(-\left(\frac{R}{R_\mrm{c}}\right)^{2-\gamma}\right)}
\end{equation}
where $\gamma = 0.8$, $R_\mrm{c} = 4\times R_0 = 200\;\mrm{au}$. The initial gas surface density at $R_0$ is $\Sigma_0 = 0.312\;\mrm{g\;cm^{-2}}$ and the total gas mass is 3 Jupiter masses (M$_\mrm{J}$). At the beginning of the simulation, the dust surface density ($\Sigma_\mrm{d}$) also follows Equation~\ref{eq:density} with a dust-to-gas ratio of 0.05, corresponding to a total dust mass of about 45 Earth masses (M$_\mrm{\earth}$). The adopted initial dust-to-gas ratio is larger than the typical value of 0.01 but is consistent with the observations of disks in Lupus \citep{ansdell2016alma}. Furthermore, as discussed in Section~\ref{sec:discussion}, we believe that the choice of the initial dust-to-gas ratio has no strong effects on the results of our analysis.

We calculate the disk evolution by neglecting the disk self-gravity and any effect due to magnetic fields. We treat dust and gas as fluids controlled by the following continuity and momentum equations \citep{fu2014effects}
\begin{equation}
\frac{\partial \Sigma_\mrm{d}}{\partial t} +\nabla \cdot \left(\Sigma_\mrm{d} \mbf{V}_\mrm{d}\right) = \nabla \cdot \left(\Sigma_\mrm{g} D_\mrm{d} \nabla \left(\frac{\Sigma_\mrm{d}}{\Sigma_\mrm{g}}\right)\right),
\label{eq:dust_con}
\end{equation}

\begin{equation}
\frac{\partial \Sigma_\mrm{g}}{\partial t} + \nabla \cdot \left(\Sigma_\mrm{g} \mbf{V}_\mrm{g}\right) =0,
\label{eq:gas_con}
\end{equation}

\begin{equation}
\frac{1}{\Sigma_\mrm{d}}\frac{\partial \Sigma_\mrm{d}\mbf{V}_\mrm{d}}{\partial t} +\frac{\nabla\left(\mbf{V}_\mrm{d}\cdot\Sigma_\mrm{d} \mbf{V}_\mrm{d}\right)}{\Sigma_\mrm{d}}=-\nabla\Psi_\mrm{G}+\mbf{F}_\mrm{drag},
\label{eq:dust_mom}
\end{equation}

\begin{equation}
\begin{aligned}
    \frac{1}{\Sigma_\mrm{g}}\frac{\partial \Sigma_\mrm{g}\mbf{V}_\mrm{g}}{\partial t} + \frac{\nabla\left(\mbf{V}_\mrm{g}\cdot\Sigma_\mrm{g} \mbf{V}_\mrm{g}\right)}{\Sigma_\mrm{g}}=&\\
-\nabla\Psi_\mrm{G}-\frac{\nabla P_\mrm{g}}{\Sigma_\mrm{g}}+&\mbf{F}_{\nu}-\frac{\Sigma_\mrm{d}}{\Sigma_\mrm{g}} \mbf{F}_\mrm{drag},
\end{aligned}
\label{eq:gas_mom}
\end{equation}
where $D_\mrm{d}$ is the dust mass diffusivity, $\Psi_\mrm{G}$ is the stellar gravitational potential, and $\mbf{V}_\mrm{d,g}$ are the dust and gas velocities. The viscous force is expressed by $\mbf{F}_\nu$, while $\mbf{F}_\mrm{drag}$ is the gas drag given by
\begin{equation}
\mbf{F}_\mrm{drag} = \frac{\Omega_\mrm{K}}{St} \left(\mbf{V}_\mrm{g} - \mbf{V}_\mrm{d}\right),
\label{eq:drag}
\end{equation}
where $St$ is the Stokes number of the dust particles in Epstein regime expressed as \citep{takeuchi2002radial}:
\begin{equation}
\label{eq:stokes}
    St = \frac{\pi \rho_\mrm{p} a_\mrm{p}}{2 \Sigma_\mrm{g}}~~,
\end{equation}
where $\rho_\mrm{p}$ and $a_\mrm{p}$ are the internal density and radius of dust particles, respectively. The last term of the right hand-side of Equation~\ref{eq:gas_mom} is the dust feedback term. The dust feedback term usually is weak in interstellar mediums because of the low dust-to-gas ratio. However, in protoplanetary disks, dust will be settled around mid-plane until it is balanced by turbulence stirring then form thin dust sub-layers~\citep{chiang2010forming,testi2014dust} and be trapped in pressure bumps to form dusty rings. Or gas can be removed by photoevaporation~\citep{alexander2006photoevaporationI,alexander2006photoevaporationII} and MHD disk wind~\citep{blandford1982hydromagnetic,bai2013wind1,bai2013wind2,bai2016magneto}. All these effects can increase the dust-to-gas ratio in protoplanetary disks then enhance the role of dust feedback.

In our analysis, we assume a single dust species characterized by $\rho_\mrm{p} = 2.0\;\mrm{g\;cm^{-3}}$ and $a_\mrm{p} = 0.25\;\mrm{mm}$. The initial $St$ at $R_0$ is about 0.25. The justification of this choice will be discussed in Section~\ref{subsubsec:stokes}. We evolve the system assuming a vertically isothermal equation of state, $P_\mrm{g} = c_\mrm{s}^2 \Sigma_\mrm{g}$, where $P_\mrm{g}$ and $c_\mrm{s}$ are the vertically integrated pressure and sound speed of gas. The radial profile of the sound speed is described as $c_\mrm{s}(R) = c_\mrm{s,0}(R/R_0)^{-1/4}$. 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 ($H = c_\mrm{s}/\Omega_\mrm{K}$, $\Omega_\mrm{K}$ is the Keplerian angular frequency). The disk temperature profile corresponding to the adopted sound speed is $T\left(R\right) = 12.4\left(R/R_0\right)^{-1/2}\;\mrm{K}$ and it agrees with the disk temperature calculated using the radiative transfer simulations discussed in the next section.

In our hydrodynamic simulations, the kinematic viscosity is parameterized as $\nu = \alpha c_\mrm{s} H$, where $\alpha$ is the dimensionless viscous parameter \citep{shakura1973black,pringle1981accretion}. We further assume that $\alpha$ 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 \citep[MRI,][]{balbus1998instability} dead-zone models \citep{gammie1996layered,armitage2011dynamics}. In this paper, we consider the Ohmic dead-zone viscosity transition model, which ignore the Hall effect~\citep{bai2014hall1,bai2014hall2} and the ambipolar diffusion~\citep{bai2011ambipolar}. Following \cite{regaly2011possible}, we assume that the radial profile of $\alpha$ in presence of a dead-zone is parameterized as
\begin{equation}
\label{eq:alpha}
    \alpha\left(R\right) = \alpha_0 - \left(\frac{\alpha_0-\alpha_{\mrm{DZ}}}{2}\right)\left( 1-\tanh\left(\frac{R-R_{\mrm{DZ}}}{\Delta_{\mrm{DZ}}}\right)\right)
\end{equation}
where $\alpha_0 = 10^{-2}$ and $\alpha_{\mrm{DZ}} = 10^{-5}$ are the viscosity parameters outside and within the dead-zone, respectively.

The transition between low and high viscosity happens at $R_\mrm{DZ} = 1.5 R_0 = 75$ au, and the width of the transition zone is controlled by the parameter $\Delta_{\mrm{DZ}}$. The transition from high to low viscosity leads to a corresponding formation of a gas over-density between $R_0$ and $R_\mrm{DZ}$. This over density leads to a local maximum in the gas pressure that is capable of forming a dusty ring by trapping dust particles as described in \cite{weidenschilling1977aerodynamics}.

Numerical simulations have shown that if $\Delta_{\mrm{DZ}} < 2H$, this ring becomes Rossby wave unstable and develop a large-scale vortex \citep[M17,][]{lyra2009planet, regaly2011possible}. Since the focus of this paper is on azimuthally symmetric rings, we avoid the formation of a large scale vortex by adopting $\Delta_{\mrm{DZ}} \geq 2H$. More precisely, we explore three cases in which $\Delta_{\mrm{DZ}}$ is equal to $2H, 3H$ and $4H$. The corresponding radial profiles of $\alpha$ are shown in Figure~\ref{fig:viscosity}.

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. At $R_0$, the disk pressure scale height is resolved by about 90 cells. 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 dust to gas ratio in the dusty rings approaches unity. A low level of noise is added to the 2D flow to seed non-axisymmetric instabilities. Performing the 1D simulation does not affect our results since no instabilities develop at the beginning of the simulation in these two models.

\begin{figure}[!t]
\centering
\includegraphics[width=0.50\textwidth]{Viscosity_profile.jpeg}
\caption{Radial profile of the viscosity parameter $\alpha$ 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 $\Delta_{\mrm{DZ}}$ is equal to $2H$, $3H$ and $4H$, respectively, where $H$ is the pressure scale height of the disk.}
\label{fig:viscosity}
\end{figure}

\subsection{Radiative Transfer and Synthetic Images}\label{subsec:radiative_images}
The outputs of LA-COMPASS 2D hydrodynamic simulations are post-processed using the radiative transfer code RADMC-3D\footnote{\url{http://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/}}~\citep{dullemond2012radmc} and the Common Astronomy Software Applications (CASA)\footnote{\url{https://casa.nrao.edu/}}~\citep{mcmullin2007casa} 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 \cite{jin2016modeling} and \cite{liu2018ngvla}. In brief, we convert the 2D gas surface density ($\Sigma_\mrm{g}$) into a 3D density distributions ($\rho_\mrm{g}$) assuming hydrostatic equilibrium for a vertically isothermal disk such that $\rho_\mrm{g}\left(R,\Phi,Z\right) = \frac{\Sigma_\mrm{g}\left(R,\Phi\right)}{\sqrt{2\pi}H\left(R\right)}\exp\left({\frac{-Z^2}{2H\left(R\right)^2}}\right)$ \cite[][Chapter 5]{frank2002accretion}.

We then assume that the disk temperature is controlled by sub-micron size dust particles which are well coupled to the gas, and therefore follow the gas density distribution scaled for a constant dust-to-gas ratio equal to 0.01. The disk temperature is calculated using RADMC-3D adopting a dust opacity proper of interstellar grains made of astronomical silicates \citep{weingartner2001dust}, organic carbonates \citep{zubko1998size}, and water ice with fractional abundances as in \cite{pollack1994composition} and a MNR grain size distribution ($n(a) \propto a^{-3.5}$) with a minimum grain size of 0.5 $\mu$m and a maximum grain size of 10 mm. The calculation of dust opacity is discussed in \cite{isella2009structure}. The disk temperature of mid-plane at $R_0$ calculated by RADMC-3D and that corresponding to the sound speed adopted in the hydrodynamic simulations differ 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 ray-tracing 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 \citep{pinte2015dust, isella2016ringed}. Note however that, since we present face on ($i = 0^\circ$) images of the dust emission, variations in the adopted dust scale height will not significantly affect our results. The dust opacity at (sub-)millimeter wavelength is assumed to be $2.3\;\mrm{cm^2\;g^{-1}}\times(\nu/230\;\mrm{GHz})^{0.4}$\citep{beckwith1990survey} and 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 observation that properly account for observational noise. We adopt an array configuration that deliver a synthetic beam FWHM of $0.039\arcsec \times 0.036\arcsec$ when a Briggs robust parameter of 0.5 is adopted in the image deconvolution. We integrate for 3 hours under typical weather conditions, resulting in an RMS noise of 6.1 $\mu$Jy beam$^{-1}$.

\section{Hydrodynamical simulation results} \label{sec:hydro_results}
\subsection{General properties}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Radial_profiles_hydrodynamical_runs_3columns.jpeg}
\caption{Panel a: Azimuthally averaged radial profile of the gas pressure at orbit 4000 normalized
by $P_0 = c_\mrm{s,0}^2 \Sigma_0$. The blue, green, and red curves correspond to the $2H$, $3H$, and $4H$, models. The vertical arrows mark the location of the local pressure maxima induces by the radial viscosity variations. Panel b: Azimuthally averaged radial profile of the dust-to-gas ratio at orbit 4000. The gray dashed line indicate a dust-to-gas ratio equal to unity. Dense dusty rings form around the location of the gas pressure maxima. Panel c: Maximum of the dust-to-gas ratio as a function of time.
}
\label{fig:hydro_summary}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Density_figures_for_specific_frames_2H_1.jpeg}
\caption{Top row: polar maps of the gas surface density for the model $2H$ at orbit 260 (left), 600 (center), and 1740 (right) normalized by initial gas surface density $\Sigma_0$. Bottom row: snapshots of dust-to-gas ratio at the same orbits. At orbit 260, an instability occurs at the outer edge of dusty ring. At orbit 600, small dusty vortices (dust clumps) form closed to the outer edge of the dusty ring while its inner edge becomes unstable. At orbit 1740, more than 50 vortices characterized by dust-to-gas ratio above unity are present along the dusty ring. We will discuss the properties of dust clumps in Section~\ref{subsubsec:GI} and Figure~\ref{fig:hist_clumps}. }
\label{fig:density_2H_1}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Density_figures_for_specific_frames_2H_2.jpeg}
\caption{Same as Figure~\ref{fig:density_2H_1}, but for orbit 2740, 3060 and 4000. At orbit 2740, the edges of the dusty ring show billow-shaped features but no dense clump are present. At orbit 3060 and 4000, the dusty ring shows both strong instabilities and several dense clumps.}
\label{fig:density_2H_2}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Density_figures_for_specific_frames_3H.jpeg}
\caption{Same as Figure~\ref{fig:density_2H_1}, but for the $3H$ model. Panels from left to right correspond to orbit 1000, 1060, 2000, 3000 and 4000. Billow-shaped features form at the outer edge of the dusty ring around orbit 1060. Differently from the $2H$ case, dense dust clumps never form in this model.}
\label{fig:density_3H}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Density_figures_for_specific_frames_4H.jpeg}
\caption{
Same as Figure~\ref{fig:density_2H_1}, but for the $4H$ model. Panels from left to right correspond to orbit 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.}
\label{fig:density_4H}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Mechanism.jpeg}
\caption{Schematic diagram of the evolution of gas velocity around the position of the gas pressure maximum without (left) and with dust feedback (right). 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 (grey 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 (grey 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 dust-to-gas ratio. 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 onset of the instability. This instability at the inner edge of dusty rings in $2H$ happens later than that at outer edge, because it needs to take some time to produce strong enough gradient at the inner edge.}
\label{fig:origin}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Vortisity_PV_2H.jpeg}
\caption{Panel a: Azimuthally averaged vorticity of gas. Panel b: Azimuthally averaged potential vorticity. Blue, green, and red solid curves correspond to time frames in which instability has not developed yet (orbit 200), in a linear regime (orbit 260), and in a non linear regime (orbit 600), respectively. The vorticity and potential vorticity are all normalized by code units and Keplerian components are subtracted. The x-axis indicates radial distances from the center of dusty ring divided normalized by ring radius. The vertical dashed-dotted lines mark the inner and outer edges of dusty rings at different orbits.}
\label{fig:Vorticity_PV_2H}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Vortisity_PV_3H_4H.jpeg}
\caption{Similar with the Figure~\ref{fig:Vorticity_PV_2H}, but for the $3H$ and $4H$ models. We chose the frames without instability (``$3H$ 1000'' and ``$4H$ 1700'') and the frames with instability (``$3H$ 1060'' and ``$4H$ 1780'').}
\label{fig:Vorticity_PV_3H_4H}
\end{figure*}

The general behaviour of the gas and dust evolution for disk models with different values of $\Delta_\textrm{DZ}$ is presented in Figure~\ref{fig:hydro_summary}. The imposed radial variation of the viscosity parameter $\alpha$ (Equation~\ref{eq:alpha}) 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 dust-to-gas ratio (panel b and c). At the end of the simulation, the azimuthally averaged dust-to-gas ratio at the position of the dusty rings approaches unity in all three models (panel b and c).

The temporal evolution of the dust-to-gas ratio strongly depends on the radial profile of $\alpha$. In the $2H$ model, the maximum value of the dust-to-gas ratio 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 (Figure~\ref{fig:density_2H_1} and \ref{fig:density_2H_2}). We will discuss these structures in detail in the next section. Conversely, the maximum of dust-to-gas ratio in model $3H$ remains below about 2 for the entire simulation and shows small temporal variation, while, in model $4H$, the maximum dust-to-gas ratio settles around a value of 0.7 after reaching unity
after about 500 orbits. The total mass of solids also evolves with time and it decreases from the initial value of 45~$\mrm{M_\earth}$ to about 35 $\mrm{M_\earth}$, 42 $\mrm{M_\earth}$ and 43 $\mrm{M_\earth}$ at the end of simulation in the $2H, 3H$ and $4H$ models, respectively. Note that all the dust particles in the outer disk are collected in the ring.

Figure~\ref{fig:density_2H_1} to \ref{fig:density_4H} 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 time scale for exciting the instability and their amplitude 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 on 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 \citep{takeuchi2002radial}. Without the dust feedback, the dusty ring would be very narrow and dense, and become gravitational unstable if self-gravity is included. The maxima of the dust-to-gas ratio would easily exceed unity. Conversely, if dust feedback is considered, the gas dynamics will be significantly affected by solid particles as soon as the dust-to-gas ratio approaches unity. This leads to a substantial change in the distribution and kinematics of the gas that might hinder building up of pressure bump \citep{taki2016dust}. And the widths of dusty rings become broader~\citep{kanagawa2018impacts} including dust feedback. Our simulations indicate that dust feedback leads to an unstable rings and might results in the formation of dense clumps. The details of these instabilities are discussed in the following sections.

\subsection{$2H$ model}
Figure~\ref{fig:density_2H_1} and \ref{fig:density_2H_2} present the temporal evolution of the gas and dust surface density in the $2H$ model. Around orbit 260, the dust-to-gas ratio 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 dusty ring became 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 dust-to-gas ratio as high as 50 are visible along the dusty ring. The dust clumps have masses between 0.1 and $0.5\;\mrm{M_\earth}$. They are roughly elliptical in shape and have a characteristic diameter of 0.6 au, or 1/5 of the disk scale height. Our simulations resolved the dust clumps with about 10 grid cells. As time proceeds, the strong dust feedback destroys the vortices \citep{fu2014effects} and, after about 2740 orbits, the dusty ring returns to be almost featureless with a dust-to-gas ratio around 2. After that, the disk undergoes a new phase of instabilities that lead to the generation and dissipation of new small vortices. In the last part of our simulation, between orbit 3000 and 4000, the dusty ring features dense dust clumps and feather-like features that extend over $180^\circ$ along the azimuthal direction.

\subsection{$3H$ and $4H$ models}
Figure~\ref{fig:density_3H} and \ref{fig:density_4H} present the temporal evolution of the gas and dust surface density for models $3H$ and $4H$, respectively. As shown in Figure~\ref{fig:hydro_summary}, shallower variations of $\alpha$ result in pressure bumps that are both shallower and located closer to the central star. In both model, the outer edge of the dusty ring becomes unstable. However the instability grows at a lower pace and to a lower amplitude compared to the $2H$ model.

In model $3H$, billow-shapes 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~\ref{fig:density_2H_1}. Furthermore, no instabilities develop at the inner edge of the dusty ring. In model $4H$, weak billow-shapes features appear around orbit 1780 and remain quasi-stable for the entire length of the simulation. At orbit 4000, the azimuthally averaged dust-to-gas ratio in the $3H$ and $4H$ model, is 0.8 and 0.6, respectively.

\subsection{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~\ref{fig:origin}. A circumstellar disk characterized by a monotonic decreasing profile of the gas surface density and temperature typically rotates at a slightly sub-Keplerian velocity \citep{adachi1976gas,weidenschilling1977aerodynamics}.  For the assumed initial disk parameters, the gas orbital velocity $V_\mrm{g,\Phi}$ at $R_0$ is about $1-O\left(\frac{H}{R}\right)^2$ of the Keplerian velocity $V_\mrm{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 $V_\mrm{g,\Phi}$, 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 dust-to-gas ratio to grow. Neglecting feedback of dust particles on the gas dynamics, and without any other process capable to halt the dust drift, solid particles will smoothly concentrate in an 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 would push the gas to rotate a Keplerian velocity creating inflection points in the gas velocity and vorticity profile. As the dust-to-gas ratio 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 dust-to-gas ratio 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.
The azimuthally averaged orbital vorticity ($\omega = \left[\nabla \times \left(\mbf{V}_\mrm{g} - \mbf{V}_\mrm{K}\right)\right]_z$) and potential vorticity ($PV$, vortensity) profiles where $PV = \omega/\Sigma_\mrm{g}$ for $2H$ are shown in Figure~\ref{fig:Vorticity_PV_2H}. Same profiles for $3H$ and $4H$ models are shown in Figure~\ref{fig:Vorticity_PV_3H_4H}, respectively. The three curves in Figure~\ref{fig:Vorticity_PV_2H} show the development of instability (orbit 200, 260) before it becomes completely nonlinear (orbit 600). Ultimately, the velocity knee at the outer edge the disk smooths out, eventually leading to a relative quiescent phase around orbit 2740. Similar behaviors are captured for $3H$ and $4H$ models as shown in Figure~\ref{fig:Vorticity_PV_3H_4H}.

Figures~\ref{fig:origin} - \ref{fig:Vorticity_PV_3H_4H} 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 Rossby wave instability \citep{lovelace1999rossby,li2000rossby,li2001rossby,ono2016parametric,ono2018parametric}. In fact, Figure~\ref{fig:density_2H_1} to \ref{fig:density_4H} show that the linear development of unstable modes have an azimuthal mode number between 3 and 5, similar to the predictions given in \cite{li2000rossby} for the most unstable modes. We plan to carry out 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 $\sim H$ to tens of $H$.

Density waves excited by planets~\citep{goldreich1978excitation,goldreich1979excitation,goldreich1980disk} carry angular momentum flux~\citep[AMF,][]{goodman2001planetary,rafikov2002planet,dong2011density1,dong2011density2,bae2018planet1,bae2018planet2,miranda2019gaps,miranda2019multiple} while propagating. When the density waves shock, AMF 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. \cite{pierens2019vortex} showed interactions between a low-mass planet and an inviscid pebble-rich disk. They found that a RWI-liked instability happens at gap edges of planet and a lot of dusty vortices emerged in their simulations. \cite{yang2019morphological} 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 both what they found are consistent with the instability in this paper.

\section{Synthetic ALMA images} \label{sec:images_results}
\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Synthetic_image_chosen_frames_2H.jpeg}
\caption{The synthetic ALMA images of specific orbits at 1.3 mm wavelength for $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 right ascension and declination. The synthetic beam is marked as a white ellipse. The FWHMs of synthetic beams are $0.039\arcsec$ and $0.036\arcsec$, i.e. the spatial resolution is $5.5\times5.0\;\mrm{au}$. The RMS of noise is about $6.1 \;\mrm{\mu Jy\; beam^{-1}}$.}
\label{fig:image_2H}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Synthetic_image_chosen_frames_3H_4H.jpeg}
\caption{The synthetic ALMA images of specific orbits at 1.3 mm wavelength for $3H$ and $4H$ models. The beam size and the noise level are same as the $2H$ model in Figure~\ref{fig:image_2H}. Because the instability is quite weak in $3H$ and $4H$, we just chose the frames with instability occurring and last frames of simulations. We can see that the rings in $3H$ and $4H$ are quite smooth compared with rings in $2H$.}
\label{fig:image_3H_4H}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Delta_intensity.jpeg}
\caption{Normalized intensity variation along the dusty rings for $2H$, $3H$ and $4H$ models. The black dashed lines mark the variations equaled to zero. The normalized intensity variation is defined as $\Delta I = (I - \overline{I})/\overline{I}$, where $\overline{I}$ is the azimuthally averaged intensity along dusty rings. The definition of the azimuthal angle $\Phi$ is started from right with anti-clockwise direction (from the west to the north), same as Figure~\ref{fig:density_2H_1} - \ref{fig:density_4H}, Figure~\ref{fig:image_2H} and \ref{fig:image_3H_4H}. The radial width of azimuthal cut is about 1 radial gaussian FWHM width for each dusty ring. For convenience, we calculate the variations of rings from intensity images directly, instead of calculating visibility on uv-plane. The amplitude of fluctuations caused by noise is an order of $1/\left(S/N\right)\sim 0.1$, as the black dotted lines show in each panel.}

\label{fig:delta_intensity}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Delta_tau.jpeg}
\caption{The optical depth $\tau$ of specific orbits for $2H$, $3H$ and $4H$ models. The black dashed lines mark the optical depth equaled to unity. The small fluctuation of optical depths is caused by noise. The optical depth of dusty rings for $2H$ present a large oscillation by strong instabilities. The optical depth of dusty rings for $3H$ are about 0.6 and 1.1 at 1060 and 4000 orbits, respectively. The optical depth of dusty rings for $4H$ are about 0.8 and 1.2 at 1780 and 4000 orbits, respectively.}
\label{fig:delta_tau}
\end{figure*}

\begin{figure}[!t]
\centering
\includegraphics[width=0.5\textwidth]{Stokes_number_4000_orbits.jpeg}
\caption{Azimuthally averaged Stokes number at 4000 orbits for $2H$, $3H$ and $4H$ models. The black line is the Stokes number of the initial frame and the dashed-dotted lines marked the boundaries of the dusty rings for each model. The Stokes numbers in the dusty rings are about 0.045, 0.055 and 0.060 for $2H$, $3H$ and $4H$, respectively.}
\label{fig:stokes_number}
\end{figure}

\begin{figure}[!t]
\centering
\includegraphics[width=0.50\textwidth]{Drag_Force.jpeg}
\caption{The radial and azimuthal components of drag force in the equation~\ref{eq:drag} varied with the Stokes number and dust size for different $\alpha$. The x axes are the Stokes number and dust size. The y axis is the strength of the drag force in arbitrary unit. The dashed lines and dashed-dotted lines mark the Stokes numbers and dust size when the drag force equal the mean drag force.}
\label{fig:drag_force}
\end{figure}

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.

\subsection{Intensity of the continuum emission}

Figures~\ref{fig:image_2H} and \ref{fig:image_3H_4H} show synthetic images of different snapshots for $2H$, $3H$ and $4H$ models in the 1.3 mm dust continuum emission calculated as discussed in Section~\ref{sec:setup}. The spatial resolution of the observations is about 5 au and is comparable to the radial extent of the ring. The RMS noise is $6.1 \;\mrm{\mu Jy\;beam^{-1}}$ and the emission is detected with a peak signal-to-noise ratio (S/N) of $20\sim25$. At this sensitivity and resolution, none of the billow-shaped features that characterize the edges of the dusty ring in $2H$ model are visible. Also, none of the dense dust clumps is 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. Only at orbit 1060 in the $3H$ model, the azimuthal variations in the intensity marginally trace the location of the billow-shapes structures present at the outer edge of the ring. At all 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~\ref{fig:delta_intensity} oscillates between $\pm20\%$ of the mean intensity in accordance with the position of the billow-shaped features shown in \ref{fig:density_2H_1}. 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\% which are caused by the non-uniform distribution in azimuth of the dense dust clumps. Indeed, at this orbit, the dust clumps are organized in two clusters centered around $\Phi=30^\circ$ and $\Phi=250^\circ$, 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 behaviour 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 ($I_\mrm{max}(\Phi)/I_\mrm{min}(\Phi) = (1+\Delta I_\mrm{max}(\Phi))/(1+\Delta I_\mrm{min}(\Phi)$) up to 4, as for orbit 3060.

\subsection{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 $\tau$ of our synthetic intensity maps from the relation:
\begin{equation}
\label{eq:tau}
    I_\nu = B_\nu\left(T_\mrm{d}\right)\left(1 -e^{-\tau_\nu}\right),
\end{equation}
where $T_\mrm{d} \simeq 12$~K is the dust temperature inside the dusty ring, and $B_\nu$ 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 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, $\tau$ 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~\ref{fig:delta_tau} shows the optical depth for the $2H, 3H$, and $4H$ models calculated using Equation~\ref{eq:tau}. 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, $\tau$ 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, and varies between $0.6\sim1.1$ and $0.8\sim1.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 ring. For example, the synthetic observations of model $2H$ provide a total dust mass between 11 and 18 $\mrm{M_\earth}$, while the true dust mass trapped within the ring is between 20 and 35 $\mrm{M_\earth}$ in $2H$ model. In the case of the $3H$ model, the dust mass inferred from the observations is between 18 and 27 $\mrm{M_\earth}$, to be compared to a true dust mass between 30 and 42 $\mrm{M_\earth}$. Finally, for model $4H$, the mass estimated from the observations is between 35 and 37 $\mrm{M_\earth}$, to be compared with a true mass between 37 and 42 $\mrm{M_\earth}$.

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.

\section{Discussion} \label{sec:discussion}

\subsection{Dependence on Model Parameters}\label{subsec:parameters}
The results discussed above 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 used 2D simulations whereas circumstellar disks are 3D objects. Here we discuss how some these assumptions might affect our results and present ideas for future follow-up studies.

\subsubsection{Dust Size and Stokes Number}\label{subsubsec:stokes}
Dust feedback is a necessary condition for the instability found in this study. The dust feedback depends on the coupling between gas and dust, which, in turn, is controlled by the dust Stokes number. Our simulations assume a single grain size $a_\mrm{p}=0.25$ mm, 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), assumes values between 0.05 and 0.15 across the dusty ring (Figure~\ref{fig:stokes_number}).

It is important to understand how different sizes of the particles contribute to the dust feedback. To do that, we adopt the approach of \cite{takeuchi2002radial} and express the difference between the azimuthal velocity of gas and dust as
\begin{equation}
\label{eq:dvphid}
V_{\mrm{d,\Phi}}-V_{\mrm{g,\Phi}}= -\frac{1}{2}V_{\mrm{d,R}}St,
\end{equation}
where $V_\mrm{d,R}$ is the radial velocity of the dust
\begin{equation}
\label{eq:vrd}
    V_{\mrm{d,R}} =-\frac{ \eta V_\mrm{K}}{St + St^{-1}},
\end{equation}
and
\begin{equation}
\eta=-\left( \frac{H}{R} \right)^2 \frac{\partial \ln P}{\partial \ln R}
\end{equation}
describes how much slower the gas rotates compared to the Keplerian velocity due to the radial pressure gradient. Combining these two equations with Equation~\ref{eq:drag}, we find
\begin{equation}
F_{\mrm{drag,\Phi}} = \frac{\eta \Omega_\mrm{K} V_\mrm{K}}{2(St+St^{-1})}
\end{equation}
Along the radial direction, the velocity difference between gas and dust is
\begin{equation}
\label{eq:dvrd}
    V_\mrm{d,R}-V_\mrm{g,R} = -\frac{\eta V_\mrm{K}}{St+St^{-1}} + \frac{\alpha c_\mrm{s} H}{R},
\end{equation}
which combined to Equation~\ref{eq:drag} gives
\begin{equation}
    F_\mrm{drag,R} = \frac{\eta V_\mrm{K} \Omega_\mrm{K}}{St(St+St^{-1})}- \frac{\alpha c_\mrm{s}^2}{R St}.
\end{equation}

If the dust population is composed of grains with different sizes, the average of the dust drag force per unit of mass is
\begin{equation}
    \overline{F}_\mrm{drag,R} = \frac{\int_{a_{min}}^{a_{max}} F_\mrm{drag,R} n(a) a^{3} da}{\int_{a_{min}}^{a_{max}}n(a) a^{3} da},
\end{equation}
\begin{equation}
    \overline{F}_\mrm{drag,\Phi} = \frac{\int_{a_{min}}^{a_{max}} F_\mrm{drag,\Phi} n(a) a^{3} da}{\int_{a_{min}}^{a_{max}}n(a) a^{3} da},
\end{equation}
where $n(a) \propto a^{-3.5}$ between $a_{min} =0.1$ $\mu$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 the latter equation and calculate the characteristic particles size by solving the relations $F_\mrm{drag,R}=\overline{F}_\mrm{drag,R}$ and $F_\mrm{drag,\Phi}=\overline{F}_\mrm{drag,\Phi}$. As Figure~\ref{fig:drag_force} shows, the drag force along 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\simeq 0.27$. For the assumed initial profiles of $\alpha$ and $c_\mrm{s}$, we find that the characteristic Stokes number of 0.27 along the $\Phi$ direction corresponding 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'' model in Figure 4 of M17 is similar with our $2H$ model. However, the dusty ring in ``ViscBroad'' just presents weak asymmetries and does not show a similar instability compared with $2H$ model. Because M17 used the dust particle with initial Stokes number equaled to unity, this is consistent with our analysis above. \cite{fu2014effects} found that dust will be concentrated effectively by Rossby Wave vortices and dust feedback will destroy vortices when the dust density is comparable and higher than the gas density. They used dust with Stokes number from $4\times 10^{-4}$ to 0.16 in their simulations. They found that the dust with 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 hundred orbits of our simulations.

In the paper, we just use one size dust component in our simulations. The simulations on full dust size distribution will be quite costly computationally. And we do not consider dust coagulation and fragmentation in the paper, which could play an important role in the dust evolution in the protoplanetary disk~\citep{drazkowska2019including}.

\subsubsection{2D approximation, Dust to Gas Ratio and Viscosity Profile}\label{subsubsec:ratio_viscosity}

We have used 2D approximations in our simulations. The initial instability wavelength is significantly longer than $H$ though the transition width is on the order of $H$. Some of the non-linear features generated in the $2H$ model can have sizes less than $H$. We used about 26 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 this study, we used three models with different viscosity transition widths to generate pressure bumps to trap dust. We find that the instability gets triggered when the dust-to-gas ratio approaches unity. We use a relatively high dust-to-gas ratio at the beginning ($= 0.05$), however, it does not affect the main conclusions. In our three models, the dust-to-gas ratio in the dusty rings all exceed or approach unity. We infer the critical dust-to-gas ratio is about 0.6 based on the weak instability at the outer edge of dusty rings in $4H$. In actual disks, the dust is likely to settle towards the mid-plane, producing higher dust to gas ratio. The processes that are studied here will likely to be applicable to the mid-plane region, though detailed 3D studies have to be carried out to test this scenario.

\subsubsection{Dust Opacity, Optical Depth and Dust Scattering}\label{subsubsec:opacity}
The dust opacity at (sub-)millimeter wavelengths in protoplanetary disk is commonly assumed as $\kappa \propto \nu^{\beta}$, where $0< \beta \lesssim 1$ \citep{beckwith1990survey} and we adopt $\beta = 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 $\sim 0.04\arcsec$ resolution for our models. Note that the width of dusty rings in $2H$ ($\sim 3$ au) is smaller than the beam size ($\sim 5$ au). The local optical depth of the dust clumps in $2H$ model are optically thick without the dilution of observational beam. The dusty rings will be optically thick at higher frequency and optically thin at lower frequency, respectively. When dusty rings become optically thick, the intensity variation (see Figure~\ref{fig:delta_intensity}) of dusty rings will become small and the dust clumps will be hard to distinguish.

The opacity of dust can be affected significantly by temperature, chemical composition, grain size, structure and topology of dust particle~\citep{semenov2003rosseland}. The dust component in protoplanetary disks is different with the interstellar grains~\citep{draine2006submillimeter}. The uncertainties of dust opacity also comes from the opacity models and the numerical computing~\citep{birnstiel2018disk}. There could be an order of variation on the dust opacity for different models~\citep{demyk2017low1,demyk2017low2}. For optically thin medium, the scattering can be ignored. In the situation with intermediate and high optical depth, the scattering opacity at (sub-)millimeter wavelength will create polarization emission of disks~\citep{kataoka2015millimeter} and extinct part of CO emission~\citep{isella2018disk}. \cite{liu2019anomalously} showed that the self-scattering of optically thick dust disk can produce low (sub-)millimeter spectral indices. \cite{zhu2019one} proposed that optically thick dust with scattering can lead to underestimating dust masses. The effects of scattering on interpreting our simulation images will be explored in the future.

\subsection{Self Gravity, Gravitational Instability and Streaming Instability}\label{subsubsec:GI}
When the thermal pressure and the rotation of disk can not support the self-gravity of disk, the gravitational instability (GI) will happen. For gas disk, this can be characterised by Toomre Q value~\citep{toomre1964gravitational,goldreich1965gravitational}:
\begin{equation}
Q \equiv \frac{\kappa c_\text{s}}{\pi G \Sigma}
\label{eq:toomre_Q}
\end{equation}
where $\kappa$ is the radially epicyclic frequency ($\kappa^2\equiv \frac{1}{R^3} \frac{d}{dR}\left(R^4 \Omega^2\right)$), and $G$ is gravitational constant. For pure dust disk (such as Saturn rings), the velocity dispersion $\sigma$ between dust particle replaces the sound speed $c_\text{s}$ in Equation~\ref{eq:toomre_Q}~\citep{armitage2007lecture}. When $Q$ is smaller than 1, GI will happen when the self-gravity is included.

\begin{figure*}[!t]
\centering
\includegraphics[width=1.0\textwidth]{Hist_masses_clumps.jpeg}
\caption{The histograms of properties of dust clumps in orbit 1740 for $2H$ shown in Figure~\ref{fig:density_2H_1}. Panel a: Radial extents of dust clumps. Panel b: Aspect ratios between $\Phi$ extent and $R$ extent of dust clumps. Panel c: Masses of dust clumps. We assume the dust clumps are with ellipse shape, then measure the peak density and densities at radial and azimuthal edges so as to estimate the masses in dust clumps. There are more than 50 clumps in the dusty ring in orbit 1740 for $2H$. The total masses of dust clumps are about $23\;\mrm{M_\oplus}$, compared with $34\;\mrm{M_\oplus}$ for the whole dusty ring. The medians of radial extents, aspect ratios and dust masses are 0.6 au, 3.6 and $0.4\;\mrm{M_\oplus}$.
}
\label{fig:hist_clumps}
\end{figure*}

Dust is trapped in the small vortices created by the instability where dust-to-gas ratio is high ($\Sigma_\text{d}/\Sigma_\text{g} \gtrsim 10$), such as ``$2H$ 1740'', ``$2H$ 3060'' and ``$2H$ 4000'' (see panel $c$ of Figure~\ref{fig:hydro_summary}).  Because dust are marginally coupled with gas ($St \lesssim 0.1$ in dusty rings, see Figure~\ref{fig:stokes_number}), we can still use the sound speed $c_\text{s}$ to calculate $Q$ as an approximation. Figure~\ref{fig:hist_clumps} describes the overall properties of the dust clumps developed at late time in the $2H$ model. We can see that the median values of radial extent, aspect ratio and masses of dust clumps are 0.6 au, 3.6 and 0.4 $\mrm{M_\earth}$, respectively. The Toomre $Q$ of the dust clumps is then about 5, so the dust clumps are marginally gravitational unstable. Since we did not include self-gravity in our simulations, it might be interesting to explore
the effects of GI in determining the properties of these dust clumps. Planetesimals might be formed in the dust clumps by the aid of GI.

The former studies \citep{goodman2000secular,youdin2005streaming,youdin2007protoplanetary,johansen2007protoplanetary} showed that dust can be concentrated under the dust-gas mutual aerodynamic drag, leading the local, linear Streaming Instability (SI). SI can generate high dust-density enhancement without self-gravity. The instability found in this paper is on 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 interplay between this instability and SI still need to be carried out in 3D simulations in the future.

\subsection{Implication of ALMA Observations}\label{subsec:implication}
Dusty rings are the most common features in the dust continuum of protoplanetary disks~\citep{andrews2018disk,long2018gaps}. Here we discuss some protoplanetary disks with significant ring-features observed by ALMA, in the context of the instability presented in this paper.

\subsubsection{HD~163296}\label{subsubsec:HD163296}
HD~163296 is a Herbig Ae star located at the distance of 101 pc~\citep{brown2018gaia}. It has extended dust and gas disks around it~\citep{isella2016ringed,isella2018disk}. HD~163296 disk is the third circumstellar disk which is revealed by ALMA with multiple rings in dust continuum observation (after HL~Tau \citep{almapartnetship2015} and TW~Hya \citep{andrews2016ringed}). \cite{isella2007millimeter} used SMA and VLA observations to show asymmetric dust continuum and CO isotopomers depletion in HD~163296 disk. \cite{isella2016ringed} 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 by three embedded Saturn-like planets. \cite{liu2018new} 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 ($\alpha < 10^{-4}$) in the inner part and large-viscosity ($\alpha \sim 7.5\times 10^{-3}$) in the outer part of disk. The viscosity transition is consistent with the traditional dead-zone model caused by Ohmic dissipation~\citep{gammie1996layered}.

\cite{isella2018disk} presented that there are three concentric dusty rings in the HD~163296 with the Gaussian radial widths 8.7, 6.6 and 5.8 au. They also found some asymmetric features in the dust emission of HD~163296 and a significant Arc-like feature in the inner side of the dusty ring at 67 au (Figure 1 of their paper). There is about $\pm 15\%$ amplitude intensity variation relative to the mean ring intensity in the dusty rings of HD~163296. \cite{zhang2018disk} run 0.15 $\mrm{M_J}$ planet-disk interaction model with 2D hydrodynamic simulation and radiative transfer post processing. They showed that the Arc-like feature can be explained by dust trapping in 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 planet on the gas kinematics in the disk.

From our simulations, we can see the 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~\ref{fig:image_2H}). The dust clumps and crescents are along the dusty rings in $2H$. They differ from the Arc-like feature of HD~163296, which is located outside of the dusty ring. The large intensity variations for $2H$ model might not be matched well with $\pm 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 non-linear outcome of this instability.

\subsubsection{MWC~758}\label{subsubsec:MWC758}
MWC~758 is a young \citep[3.5 millions years,][]{meeus2012observations} Herbig Ae star with circumstellar disk, located at 160 pc away~\citep{brown2018gaia}. \cite{marino2015compact} showed that there are two bright bumps at sub-millimeter continuum. \cite{boehler2018complex} showed that MWC~758 presents complex morphology with ALMA Band 7 dust continuum and $^{12}$CO/$^{13}$CO gas line emission with angular resolution $0.1\arcsec \sim 0.2\arcsec$. \cite{dong2018eccentric} showed that there are three dusty rings, two bright bumps, and one eccentric cavity based on ALMA observation at 0.87 mm with 40 mas resolution. The contrast of the south and north bright clumps is 4.4 and 10, 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 with the contrast from our simulations ($\sim 4$ in ``$2H$ 3060'', See Section~\ref{subsec:radiative_images}). There are also some asymmetries along the inner ring (See Figure 2 and 3 in \cite{dong2018eccentric}). The north bumps with high contrast ($\sim 10$) might come from dust trapping in a larger scale vortex generated by RWI. \cite{baruteau2019dust} 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 the inner vortex decays faster than the outer one, then it 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 explanation for the complex features seen in MWC~758.

\subsubsection{HL~Tau}\label{subsubsec:HL_Tau}
HL~Tau is a young T~Tauri Star located at about 140 pc from Earth~\citep{loinard2007vlba}. \cite{almapartnetship2015} used ALMA observation to present multiple bright and dark rings in HL~Tau. Their observation also showed the low spectral index in the rings, which implies that the dust might have evolved compared with the interstellar medium grain. \cite{jin2016modeling} used three Saturn mass planets to model the rings structures in the dust continuum of HL~Tau. They also calculated the optical depth at 1.0 mm of the dusty rings and they found that the bright rings (such as B1 and B2) are optically thick, and the dark rings are optically thin. \cite{carrasco2016vla} used the VLA data on HL~Tau, and they found that the brights rings of HL~Tau at the 7.0 mm are optically thin. There are some clump candidates on the B1 ring of HL~Tau (Figure 2 in their paper).

As discussed in Section~\ref{subsubsec:opacity}, the dusty rings will be smooth in the optically thick wavelength. The bright rings in HL~Tau at the wavelengths observed by ALMA are optically thick, so they are smooth and ordered. It will be difficult to assess whether there are dense clumps formed inside these rings. Observations at multiple wavelengths are essential to confirm this instability.

\section{Summary}\label{sec:summary}
We have studied a meso-scale instability caused by dust-gas interaction at the edges of dusty rings by carrying out high resolution 2D two fluid global hydrodynamic simulation, combined with the radiative transfer calculation and ALMA synthetic image processing. Here are our main conclusions:

\begin{enumerate}
    \item If the conditions are favorable, there is an instability at the edges of dusty rings produced by pressure bumps. Dust feedback alters the azimuthal gas velocity profile, producing a sharpened potential vorticity profile at the edges of dusty rings.
    \item The instability can develop when the local dust-to-gas ratio approaches unity. A ring that is being continuously ``fed" from outer disk can eventually reach such a condition, at least transiently.
    \item The non-linear outcome of such an instability can lead to many small vortices in dusty rings. These small vortices can contain large amount of dust, reaching at least $\sim 10\%$ of Earth mass within each clump.
    \item These dust clumps are highly optically thick and they might not be able to resolved by current observations at (sub-)millimeters. However, these dust clumps still can lead to both non-axisymmetric structures and intensity variations in dusty rings. This instability can be an origin of some asymmetries observed in protoplanetary disks.

\end{enumerate}

Future studies in 3D, including interaction between this meso-scale instability with other instabilities such as SI and GI, will be quite interesting. They will help to address whether dusty rings are suitable for planetesimals formation and eventually planet formation.

\acknowledgments
We thank Richard Lovelace, Douglas Lin, Andrew Youdin, Xuening Bai, Ruobing Dong, Shangfei Liu, Cong Yu for helpful discussions. This work is supported by the National Natural Science Foundation of China (grant Nos. 11773081, 11661161013, 11633009 and 11873097), the CAS Interdisciplinary Innovation Team, the B-type Strategic Priority Program of the Chinese Academy of Sciences, Grant No. XDB41000000 and the Foundation of Minor Planets of Purple Mountain Observatory. We also acknowledge the support by a LANL/CSES project. This work was partially performed at the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607611.

\software{LA-COMPASS \citep{li2005potential,li2008type}, RADMC-3D v0.41 \citep{dullemond2012radmc}, CASA ALMA pipeline 5.0.0 \citep{mcmullin2007casa}, Astropy \citep{robitaille2013astropy}, Numpy \citep{walt2011numpy}, Scipy \citep{jones2016scipy}}

\bibliographystyle{aasjournal}

\subsubsection*{#1}}

\pagestyle{headings}
\markright{Reference sheet: \texttt{natbib}}

\usepackage{shortvrb}
\MakeShortVerb{\|}

\begin{document}
\thispagestyle{plain}

\newcommand{\textsc{Bib}\TeX}{\textsc{Bib}\TeX}
\newcommand{\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}}{\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}}
\begin{center}{\bfseries\Large
    Reference sheet for \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ usage}\\
    \large(Describing version \fileversion\ from \filedate)
\end{center}

\begin{quote}\slshape
For a more detailed description of the \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package, \LaTeX\ the
source file \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\texttt{.dtx}.
\end{quote}
\head{Overview}
The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package is a reimplementation of the \LaTeX\ |\cite| command,
to work with both author--year and numerical citations. It is compatible with
the standard bibliographic style files, such as \texttt{plain.bst}, as well as
with those for \texttt{harvard}, \texttt{apalike}, \texttt{chicago},
\texttt{astron}, \texttt{authordate}, and of course \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}.

\head{Loading}
Load with |\usepackage[|\emph{options}|]{|\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}|}|. See list of
\emph{options} at the end.

\head{Replacement bibliography styles}
I provide three new \texttt{.bst} files to replace the standard \LaTeX\
numerical ones:
\begin{quote}\ttfamily
 plainnat.bst \qquad abbrvnat.bst \qquad unsrtnat.bst
\end{quote}
\head{Basic commands}
The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package has two basic citation commands, |\citet| and
|\citep| for \emph{textual} and \emph{parenthetical} citations, respectively.
There also exist the starred versions |\citet*| and |\citep*| that print
the full author list, and not just the abbreviated one.
All of these may take one or two optional arguments to add some text before
and after the citation.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90}| & Jones et al. (1990)\\
  |\citet[chap.~2]{jon90}| & Jones et al. (1990, chap.~2)\\[0.5ex]
  |\citep{jon90}| & (Jones et al., 1990)\\
  |\citep[chap.~2]{jon90}| & (Jones et al., 1990, chap.~2)\\
  |\citep[see][]{jon90}| & (see Jones et al., 1990)\\
  |\citep[see][chap.~2]{jon90}| & (see Jones et al., 1990, chap.~2)\\[0.5ex]
  |\citet*{jon90}| & Jones, Baker, and Williams (1990)\\
  |\citep*{jon90}| & (Jones, Baker, and Williams, 1990)
\end{tabular}
\end{quote}
\head{Multiple citations}
Multiple citations may be made by including more than one
citation key in the |\cite| command argument.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90,jam91}| & Jones et al. (1990); James et al. (1991)\\
  |\citep{jon90,jam91}| & (Jones et al., 1990; James et al. 1991)\\
  |\citep{jon90,jon91}| & (Jones et al., 1990, 1991)\\
  |\citep{jon90a,jon90b}| & (Jones et al., 1990a,b)
\end{tabular}
\end{quote}

\head{Numerical mode}
These examples are for author--year citation mode. In numerical mode, the
results are different.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90}| & Jones et al. [21]\\
  |\citet[chap.~2]{jon90}| & Jones et al. [21, chap.~2]\\[0.5ex]
  |\citep{jon90}| & [21]\\
  |\citep[chap.~2]{jon90}| & [21, chap.~2]\\
  |\citep[see][]{jon90}| & [see 21]\\
  |\citep[see][chap.~2]{jon90}| & [see 21, chap.~2]\\[0.5ex]
  |\citep{jon90a,jon90b}| & [21, 32]
\end{tabular}
\end{quote}
\head{Suppressed parentheses}
As an alternative form of citation, |\citealt| is the same as |\citet| but
\emph{without parentheses}. Similarly, |\citealp| is |\citep| without
parentheses. Multiple references, notes, and the starred variants
also exist.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citealt{jon90}| & Jones et al.\ 1990\\
  |\citealt*{jon90}| & Jones, Baker, and Williams 1990\\
  |\citealp{jon90}| & Jones et al., 1990\\
  |\citealp*{jon90}| & Jones, Baker, and Williams, 1990\\
  |\citealp{jon90,jam91}| & Jones et al., 1990; James et al., 1991\\
  |\citealp[pg.~32]{jon90}| & Jones et al., 1990, pg.~32\\
  |\citetext{priv.\ comm.}| & (priv.\ comm.)
\end{tabular}
\end{quote}
The |\citetext| command
allows arbitrary text to be placed in the current citation parentheses.
This may be used in combination with |\citealp|.
\head{Partial citations}
In author--year schemes, it is sometimes desirable to be able to refer to
the authors without the year, or vice versa. This is provided with the
extra commands
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citeauthor{jon90}| & Jones et al.\\
  |\citeauthor*{jon90}| & Jones, Baker, and Williams\\
  |\citeyear{jon90}|   & 1990\\
  |\citeyearpar{jon90}| & (1990)
\end{tabular}
\end{quote}
\head{Forcing upper cased names}
If the first author's name contains a \textsl{von} part, such as ``della
Robbia'', then |\citet{dRob98}| produces ``della Robbia (1998)'', even at the
beginning of a sentence. One can force the first letter to be in upper case
with the command |\Citet| instead. Other upper case commands also exist.
\begin{quote}
\begin{tabular}{rl@{\quad$\Rightarrow$\quad}l}
  when & |\citet{dRob98}| & della Robbia (1998) \\
  then & |\Citet{dRob98}| & Della Robbia (1998) \\
     &   |\Citep{dRob98}| & (Della Robbia, 1998) \\
     &   |\Citealt{dRob98}| & Della Robbia 1998 \\
     &   |\Citealp{dRob98}| & Della Robbia, 1998 \\
     &   |\Citeauthor{dRob98}| & Della Robbia
\end{tabular}
\end{quote}
These commands also exist in starred versions for full author names.

\head{Citation aliasing}
Sometimes one wants to refer to a reference with a special designation,
rather than by the authors, i.e. as Paper~I, Paper~II. Such aliases can be
defined and used, textual and/or parenthetical with:
\begin{quote}
\begin{tabular}{lcl}
  |\defcitealias{jon90}{Paper~I}|\\
  |\citetalias{jon90}| & $\Rightarrow$ & Paper~I\\
  |\citepalias{jon90}| & $\Rightarrow$ & (Paper~I)
\end{tabular}
\end{quote}
These citation commands function much like |\citet| and |\citep|: they may
take multiple keys in the argument, may contain notes, and are marked as
hyperlinks.
\head{Selecting citation style and punctuation}
Use the command |\bibpunct| with one optional and 6 mandatory arguments:
\begin{enumerate}
\item the opening bracket symbol, default = (
\item the closing bracket symbol, default = )
\item the punctuation between multiple citations, default = ;
\item the letter `n' for numerical style, or `s' for numerical superscript
      style, any other letter for
      author--year, default = author--year;
\item the punctuation that comes between the author names and the year
\item the punctuation that comes between years or numbers when common author
      lists are suppressed (default = ,);
\end{enumerate}

The optional argument is the character preceding a post-note, default is a
comma plus space. In redefining this character, one must include a space if
one is wanted.

Example~1, |\bibpunct{[}{]}{,}{a}{}{;}| changes the output of
\begin{quote}
  |\citep{jon90,jon91,jam92}|
\end{quote}
into [Jones et al. 1990; 1991, James et al. 1992].

Example~2, |\bibpunct[; ]{(}{)}{,}{a}{}{;}| changes the output of
\begin{quote}
  |\citep[and references therein]{jon90}|
\end{quote}
into (Jones et al. 1990; and references therein).

\head{Other formatting options}
Redefine |\bibsection| to the desired sectioning command for introducing
the list of references. This is normally |\section*| or |\chapter*|.

Define |\bibpreamble| to be any text that is to be printed after the heading but
before the actual list of references.

Define |\bibfont| to be a font declaration, e.g.\ |\small| to apply to
the list of references.

Define |\citenumfont| to be a font declaration or command like |\itshape|
or |\textit|.

Redefine |\bibnumfmt| as a command with an argument to format the numbers in
the list of references. The default definition is |[#1]|.

The indentation after the first line of each reference is given by
|\bibhang|; change this with the |\setlength| command.

The vertical spacing between references is set by |\bibsep|; change this with
the |\setlength| command.

\head{Automatic indexing of citations}
If one wishes to have the citations entered in the \texttt{.idx} indexing
file, it is only necessary to issue |\citeindextrue| at any point in the
document. All following |\cite| commands, of all variations, then insert
the corresponding entry to that file. With |\citeindexfalse|, these
entries will no longer be made.

\head{Use with \texttt{chapterbib} package}

The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package is compatible with the \texttt{chapterbib} package
which makes it possible to have several bibliographies in one document.

The package makes use of the |\include| command, and each |\include|d file
has its own bibliography.

The order in which the \texttt{chapterbib} and \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ packages are loaded
is unimportant.

The \texttt{chapterbib} package provides an option \texttt{sectionbib}
that puts the bibliography in a |\section*| instead of |\chapter*|,
something that makes sense if there is a bibliography in each chapter.
This option will not work when \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ is also loaded; instead, add
the option to \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}.

Every |\include|d file must contain its own
|\bibliography| command where the bibliography is to appear. The database
files listed as arguments to this command can be different in each file,
of course. However, what is not so obvious, is that each file must also
contain a |\bibliographystyle| command, \emph{preferably with the same
style  argument}.
\head{Sorting and compressing citations}
Do not use the \texttt{cite} package with \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}; rather use one of the
options \texttt{sort} or \texttt{sort\&compress}.

These also work with author--year citations, making multiple citations appear
in their order in the reference list.

\head{Long author list on first citation}
Use option \texttt{longnamesfirst} to have first citation automatically give
the full list of authors.

Suppress this for certain citations with |\shortcites{|\emph{key-list}|}|,
given before the first citation.

\head{Local configuration}
Any local recoding or definitions can be put in \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\texttt{.cfg} which
is read in after the main package file.

\head{Options that can be added to \texttt{\char`\\ usepackage}}
\begin{description}
\item[\ttfamily round] (default) for round parentheses;
\item[\ttfamily square] for square brackets;
\item[\ttfamily curly] for curly braces;
\item[\ttfamily angle] for angle brackets;
\item[\ttfamily colon] (default) to separate multiple citations with
     colons;
\item[\ttfamily comma] to use commas as separaters;
\item[\ttfamily authoryear] (default) for author--year citations;
\item[\ttfamily numbers] for numerical citations;
\item[\ttfamily super] for superscripted numerical citations, as in
     \textsl{Nature};
\item[\ttfamily sort] orders multiple citations into the sequence in
     which they appear in the list of references;
\item[\ttfamily sort\&compress] as \texttt{sort} but in addition multiple
     numerical citations are compressed if possible (as 3--6, 15);
\item[\ttfamily longnamesfirst] makes the first citation of any reference
     the equivalent of the starred variant (full author list) and subsequent
     citations normal (abbreviated list);
\item[\ttfamily sectionbib] redefines |\thebibliography| to issue
     |\section*| instead of |\chapter*|; valid only for classes with a
     |\chapter| command; to be used with the \texttt{chapterbib} package;
\item[\ttfamily nonamebreak] keeps all the authors' names in a citation on
     one line; causes overfull hboxes but helps with some \texttt{hyperref}
     problems.
\end{description}
\end{document}