© ESO, 2018

1. Introduction

High resolution, high contrast imaging at optical and sub-millimeter wavelengths has in recent years revealed that protoplanetary disks are highly structured. While some disks show complex non-axisymmetric structures, there are also numerous disks that display multiple nearly perfect concentric ringlike structures. The first and most prominent example was the disk around the star HL Tau (ALMA Partnership et al. 2015). Many other examples have since followed such as TW Hydra (Andrews et al. 2016; van Boekel et al. 2017), RX J1615.3 (de Boer et al. 2016), HD 97048 (Ginski et al. 2016), HD 163296 (Isella et al. 2016) and HD 169142 (Momose et al. 2015; Fedele et al. 2017). The rings are seen at submillimeter wavelength thermal dust emission as well as in optical/near-infrared scattered light. These rings and gaps have been interpreted as resulting from newly formed planets opening up gaps within the disk (e.g. Gonzalez et al. 2015; Kanagawa et al. 2015; Picogna & Kley 2015). This is an attractive scenario, because it would mean that we are indirectly “seeing” the young planets as they are being formed in their birth-disk. However, since we do not yet have direct indications of these planets, it is important to also investigate other explanations. Could these rings be caused by something entirely different altogether? And could this “something”, though unrelated to already existing planets, still teach us something about the formation of planets?

Various non-planet-related explanations for these rings have been proposed. In fact, their existence was predicted before their discovery, through a simple argument: It is known that dust aggregates of millimeter size and larger tend to radially drift toward the star on a very short time scale (Whipple 1972; Brauer et al. 2007). If the disk lives a few millions years, and the dust aggregates grow to the millimeter sizes inferred from millimeter-wave observations (e.g. Tazzari et al. 2016), then these disks should by the time they are observed already have lost most of their dust particles. By contrast, we see large quantities of dust in these disks, so something must be holding up the dust, preventing it from taking part in the rapid radial drift mechanism. It was suggested by Pinilla et al. (2012) that perhaps a multitude of local ring-shaped gas pressure bumps in the disk could trap the dust after it has grown to sizes of about a millimeter. The physical mechanism behind this dust trapping is the well-known, and unavoidable effect, that dust particles tend to drift toward regions of increased gas pressure (e.g. Whipple 1972; Adachi et al. 1976). Or in other words, that they drift in the direction of the gas pressure gradient ∇P. In smooth protoplanetary disks the gas pressure decreases with radius (∂_rP< 0), which is the reason for the inward drift of dust. But if the gas pressure has wiggles that are strong enough that they cause ∂_rP to flip sign, then there will exist local pressure maxima for which ∂_rP = 0, and for which the dust drift is converging. In the absense of gas turbulence, all dust grains sufficiently close by would get trapped in these traps. Turbulence can mix the small dust grains (that are most strongly coupled to the gas) out of the traps, so that they may, on average, continue to drift inward. The bigger dust particles, which are less coupled to the gas, may still remain trapped. It depends on the amplitude of the pressure bumps and the strength of the turbulence, which grain sizes remain trapped and which not.

While the origin of these gas pressure bumps was not addressed in Pinilla et al. (2012), it was shown that if they exist, and if they are strong enough, the radial drift paradoxon could be solved. As a direct consequence, however, it was shown that high-resolution ALMA observations should then be able to see these dust rings, and the predicted ALMA images bear striking resemblance to HL Tau and similar sources.

This scenario does not, however, explain why these gas pressure bumps form in the first place. Takahashi & Inutsuka (2014, 2016) propose an elegant scenario in which dust rings are formed through a secular gravitational instability (Ward 2000). The idea is that if the dust density is high enough for a dust-driven gravitational instability to occur, the gas drag will slow this process down. The slowness of this process allows the information about the gravitational contraction to shear out and spread along azimuth, so that grand-design rings are formed instead of gravitationally contracting clumps. Takahashi & Inutsuka (2016) argue that as these rings contract further, this eventually leads to planets being formed, which, in their turn, open up gaps in the dust distribution (Paardekooper & Mellema 2004). They suggest that the ring-like structures could therefore be witnesses of both the initial and the final stages of planet formation.

A completely different scenario was proposed by Zhang et al. 2015, who argue that the locations of the rings suggest their association with the snow lines of a series of different volatile molecular species. A physical mechanism by which snow lines could lead to rings was worked out by Okuzumi et al. (2016). The physics of ice sublimation and deposition near these snow lines is complex, because it interacts strongly with the coagulation and fragmentation of dust aggregates, as well as with the radial drift and turbulent mixing in the disk (e.g. Stammler et al. 2017).

Gonzalez et al. (2017) propose an alternative scenario of spontaneous ring formation. In their model the radial drift of dust particles coupled to the dust coagulation process tends to lead to ring-like regions of high dust concentration and fast growth (see also Dra¸żkowska et al. 2016). While this might explain transition disks with a single dust ring, it may be more difficult to explain the multi-ring structures discussed here.

Global magnetohydrodynamical disk simulations with dead-zones also tend to create rings-shaped structures (Flock et al. 2015), and zonal flows (Johansen et al. 2009).

In this paper we investigate whether the viscous disk evolution could lead to the spontaneous formation of rings. This idea is not new. For instance, Wünsch et al. (2005) show, using 2D radiation-hydrodynamics models, that a disk with an active surface layer, self-gravity and “dead” midplane region could become viscously unstable and lead to the formation of several concentric rings. Our proposal is, however, based on a different physical driving mechanism. A version of this mechanism was already studied in a local box model by Johansen et al. (2011).

The ring instability works as follows. If we perturb an otherwise smooth disk with an infinitesimal-amplitude wiggle in the gas pressure (of the concentric ring type), it will tend to cause a slight enhancement of the dust-to-gas ratio in these pressure enhancements. This is the same physical mechanism as for dust trapping, but we do not necessarily need a flip of sign of ∂_rP to cause this effect. It is just that the radial inward drift velocity | v_drift | of the dust is slightly increased on the outer side of the pressure enhancement and slightly reduced on the inner side, leading to a traffic-jam density enhancement effect. This effect works for big and small grains alike. It is known that dust has a negative influence on the disk viscosity, if it is caused by the magnetorotational instability (e.g. Sano et al. 2000; Ilgner & Nelson 2006; Okuzumi 2009; Dzyurkevich et al. 2013). A viscous disk reacts to the resulting wiggle in the viscosity ν(r) by adapting the surface density Σ_g(r) such that the steady-state ${\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}\mathit{\nu }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{const}\mathit{,}$(1)is restored. So, whereever ν is reduced, Σ_g is increased. This change in Σ_g(r) amplifies the initial perturbation, resulting in a positive feedback loop. We therefore expect an initial perturbation in the gas disk to be amplified by the combined effect of dust drift and the effect the dust has on the viscosity of the disk. This process is depicted in cartoon form in Fig. 1.

Fig. 1 Mechanism of the ring instability studied in this paper.

To test whether this mechanism indeed works requires a linear perturbation analysis. This is what we present in this paper.

In Sect. 2 we give the basic equations that stand at the basis of our model. These are the standard viscous disk equations coupled to a single dust component of a given Stokes number.

The linear perturbation analysis of the combined gas and dust system is rather cumbersome. So in Sect. 3 we first simplify the system of equations radically, so that the mechanism presents itself more clearly. In Sect. 4 we then tackle the full set of equations, with some mathematics moved to the appendix.

2. Basic disk equations and model assumptions

The standard viscous disk equations for the gas together with a single dust component are: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}\mathit{r}{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{v}}_{\mathrm{rg}}\mathrm{)}& \mathrm{=}& \\ \frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{d}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{\mathit{r}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}\mathit{r}{\mathrm{\Sigma }}_{\mathrm{d}}{\mathit{v}}_{\mathrm{rd}}\mathrm{)}& \mathrm{=}& \end{array}$The gas radial velocity is given by the usual viscous disk equation: ${\mathit{v}}_{\mathrm{rg}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{{\mathrm{\Sigma }}_{\mathrm{g}}\sqrt{\mathit{r}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }{\sqrt{\mathit{r}}}^{\mathrm{)}}\mathit{,}$(4)with the turbulent viscosity defined by $\mathit{\nu }\mathrm{=}\mathit{\alpha }\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{,}$(5)where α is the usual alpha-turbulence parameter. The Kepler frequency is ${\mathrm{\Omega }}_{\mathrm{K}}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}^{\mathrm{3}}}}\mathit{,}$(6)with G the gravitational constant and M_∗ the stellar mass. The isothermal sound speed squared is ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{B}}\mathit{T}}{\mathit{\mu }{\mathit{m}}_{\mathrm{p}}}\mathit{,}$(7)with k_B the Boltzmann constant, μ = 2.3 the mean molecular weight, and m_p the proton mass. The vertical pressure scale height of the disk is ${\mathit{H}}_{\mathrm{p}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{.}$(8)The dust diffusion constant is ${\mathrm{𝒟}}_{\mathrm{d}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}{\mathrm{𝒟}}_{\mathrm{g}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\frac{\mathit{\nu }}{\mathrm{Sc}}\mathit{,}$(9)where St is the Stokes number of the dust and Sc is the Schmidt number of the gas. The radial velocity of the dust is ${\mathit{v}}_{\mathrm{rd}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}{\mathit{v}}_{\mathrm{rg}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{St}\mathrm{+}{\mathrm{St}}^{-1}}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{g}}{\mathrm{\Omega }}_{\mathrm{K}}}\frac{\mathit{\partial }{\mathit{P}}_{\mathrm{g}}}{\mathit{\partial r}}\mathit{,}$(10)where ρ_g is the midplane gas density and ${\mathit{P}}_{\mathrm{g}}\mathrm{\equiv }{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$ is the midplane gas pressure.

The α determines the strength of the turbulence, and hence the strength of the viscosity and the dust diffusivity. In our analysis we allow the dust surface density to affect the value of α: a higher dust concentration will lead to a lower α (e.g. Sano et al. 2000; Ilgner & Nelson 2006; Okuzumi 2009; Dzyurkevich et al. 2013). Since the physics of magnetorotational turbulence in non-ideal MHD is not yet fully understood, we parameterize this effect. We consider the following general prescription: $\mathit{\alpha }\mathrm{=}{\mathit{\alpha }}_{\mathrm{1}}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{d}\mathrm{1}}}\right)}^{{\mathit{\phi }}_{\mathrm{d}}}{\left(\frac{{\mathrm{\Sigma }}_{\mathrm{g}}}{{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{1}}}\right)}^{{\mathit{\phi }}_{\mathrm{g}}}\mathit{,}$(11)where α₁ is the unperturbed value of α, and likewise Σ_d1 and Σ_g1 are the unperturbed values of Σ_d and Σ_g respectively. The parameters φ_d and φ_g are powerlaw indices that parameterize how α depends on the change in dust and/or gas surface density. We focus on cases with φ_d< 0, meaning that an increase in Σ_d/ Σ_d1 leads to a decrease in α. The φ_g is used to allow for the following two cases of interest: $\begin{array}{ccc}{\mathit{\phi }}_{\mathrm{g}}& \mathrm{=}& \\ {\mathit{\phi }}_{\mathrm{g}}& \mathrm{=}& \end{array}$Case 1 could, at least in principle, allow for the expected instability even without dust drift relative to the gas, since without dust drift, Σ_d will necessarily increase if Σ_g does. In contrast, case 2 strictly requires dust drift for the instability to operate, since lack of dust drift keeps Σ_d/ Σ_g constant.

3. Simplified perturbation analysis

It is cumbersome to perform the linear stability analysis for the full set of equations of Sect. 2, because these equations contain factors $\sqrt{\mathit{r}}$ and the like. We postpone this full analysis to Sect. 4. For now, let us first simplify the equations.

3.1. Simplified equations

We simplify the equations to the following form: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{v}}_{\mathit{x}\mathrm{g}}\mathrm{)}& \mathrm{=}& \\ \frac{\mathit{\partial }{\mathrm{\Sigma }}_{\mathrm{d}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }}{\mathit{\partial x}}\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{d}}{\mathit{v}}_{\mathit{x}\mathrm{d}}\mathrm{)}& \mathrm{=}& \end{array}$where r₀ is the radius at which we wish to locally carry out the perturbation analysis, and x is the dimensionless radial coordinate starting from that location: $\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}\mathrm{)}\mathit{.}$(16)The gas radial velocity formula (Eq. (4)) is simplified as: ${\mathit{v}}_{\mathit{x}\mathrm{g}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{{\mathrm{\Sigma }}_{\mathrm{g}}{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }\mathrm{\left(}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }\mathrm{\right)}}{\mathit{\partial x}}\mathit{,}$(17)with the viscosity still defined by Eq. (5). However, for simplicity c_s and Ω_K are now assumed be be constant. On the other hand, α is allowed to depend on x, but only due to the perturbation. The radial velocity of the dust is given, in our simplified description, by ${\mathit{v}}_{\mathit{x}\mathrm{d}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}{\mathit{v}}_{\mathit{x}\mathrm{g}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{St}\mathrm{+}{\mathrm{St}}^{-1}}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathrm{\Omega }}_{\mathrm{K}}{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }\ln {\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathit{.}$(18)The stationary solution is Σ_g = constant, Σ_d =constant and α =constant. This also yields v_xd = v_xg = 0 for that stationary solution. This stationary solution is the backdrop of our perturbation analysis.

Now we introduce an infinitesimal perturbation: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}& \mathrm{=}& \\ {\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}& \mathrm{=}& \end{array}$where Σ_g1 and Σ_d1 are the stationary (constant) solutions for gas and dust, respectively. From here on the subscript 1 denotes this stationary solution. We will also omit the (x,t) notation, to not clutter the equations too much. The symbols σ_g and σ_d are the dimensionless infinitesimal perturbations on the gas and the dust respectively. For the α, according to Eq. (11), we obtain, to leading order: $\mathit{\alpha }\mathrm{=}{\mathit{\alpha }}_{\mathrm{1}}(\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{g}})\mathit{,}$(21)on account of the fact that | σ_d/g | ≪ 1.

Inserting these formulae into Eqs. (14), (15), and keeping in mind that both v_xg and v_xd are linear in the perturbations, and that we omit all terms of higher order in σ_g/d, we arrive, to leading order, at: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}\mathrm{g}}}{\mathit{\partial x}}& \mathrm{=}& \\ \frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial t}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}}\frac{\mathit{\partial }{\mathit{v}}_{\mathit{x}\mathrm{d}}}{\mathit{\partial x}}& \mathrm{=}& \end{array}$By inserting Eqs. (5), (19) together with Eqs. (21) into (17), and assuming that c_s and Ω_K are constant (see above), the radial velocity for the gas becomes ${\mathit{v}}_{\mathit{x}\mathrm{g}}\mathrm{=}\mathrm{-}\mathrm{3}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}}\left(\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }\mathrm{(}{\mathit{\phi }}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial x}}\right)\mathit{.}$(24)Similarly that of the dust becomes $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{v}}_{\mathit{x}\mathrm{d}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}}\{\mathrm{-}\mathrm{3}\left(\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }\mathrm{(}{\mathit{\phi }}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial x}}\right)\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}\frac{\mathrm{St}}{{\mathit{\alpha }}_{\mathrm{1}}}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\}\mathit{.}\end{array}& & \end{array}$(25)Inserting these into the continuity equations (Eqs. (22), (23)) yields $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial t}}& \mathrm{=}& \mathrm{3}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{(}{\mathit{\phi }}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\mathit{,}\\ \frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial t}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\{\mathrm{3}\left(\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{(}{\mathit{\phi }}_{\mathrm{d}}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\\ & & \mathrm{-}\frac{\mathrm{St}}{{\mathit{\alpha }}_{\mathrm{1}}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\}\mathrm{+}\frac{\mathrm{1}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{𝒟}}_{\mathrm{d}\mathrm{1}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{(}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathit{.}\end{array}$Now let us assume the linear perturbations to be plane waves in space x and time t. Since we seek modes for which the dust and the gas perturbations growth due to their mutual coupling, we can assume a single spatial frequency k and time frequency ω for both modes: $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{g}}& \mathrm{=}& \mathit{A}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathit{kx}}\mathit{,}\\ {\mathit{\sigma }}_{\mathrm{d}}& \mathrm{=}& \mathit{B}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathit{kx}}\mathit{.}\end{array}$The complex amplitudes A (for the gas) and B (for the dust) can be set independently. Inserting this mode into the above set of equations yields: $\begin{array}{ccc}\mathrm{i}\mathit{\omega A}& \mathrm{=}& \mathrm{-}\mathrm{3}{\mathit{k}}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}\mathit{A}\mathrm{+}{\mathit{\phi }}_{\mathrm{d}}\mathit{B}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{A}}^{\mathrm{)}}\mathit{,}\\ \mathrm{i}\mathit{\omega B}& \mathrm{=}& \frac{\mathrm{-}{\mathit{k}}^{\mathrm{2}}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\{\mathrm{3}\mathrm{(}\mathit{A}\mathrm{+}{\mathit{\phi }}_{\mathrm{d}}\mathit{B}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}{\mathit{A}}^{\mathrm{)}}\mathrm{-}\frac{\mathrm{St}}{{\mathit{\alpha }}_{\mathrm{1}}}\mathit{A}\\ & & \mathrm{+}\frac{\mathrm{1}}{\mathrm{Sc}}\mathrm{(}\mathit{B}\mathrm{-}\mathit{A}\mathrm{)}\}\mathit{,}\end{array}$where we made use of Eq. (9) to replace . This can be put into matrix form: $\mathrm{i}\mathit{\omega }\left(\begin{array}{c}\\ \mathit{A}\\ \mathit{B}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ {\mathit{M}}_{\mathit{aa}}& {\mathit{M}}_{\mathit{ab}}\\ {\mathit{M}}_{\mathit{ba}}& {\mathit{M}}_{\mathit{bb}}\end{array}\right)\left(\begin{array}{c}\\ \mathit{A}\\ \mathit{B}\end{array}\right)\mathit{,}$(32)with $\begin{array}{ccc}{\mathit{M}}_{\mathit{aa}}& \mathrm{=}& \mathrm{-}\mathrm{3}{\mathit{k}}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\mathrm{)}\mathit{,}\\ {\mathit{M}}_{\mathit{ab}}& \mathrm{=}& \mathrm{-}\mathrm{3}{\mathit{k}}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\phi }}_{\mathrm{d}}\mathit{,}\\ {\mathit{M}}_{\mathit{ba}}& \mathrm{=}& \mathrm{-}\mathrm{3}{\mathit{k}}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\left(\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\mathrm{)}\mathrm{-}\frac{\mathrm{St}}{\mathrm{3}{\mathit{\alpha }}_{\mathrm{1}}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{3}\mathrm{Sc}}\right)\mathit{,}\\ {\mathit{M}}_{\mathit{bb}}& \mathrm{=}& \mathrm{-}\mathrm{3}{\mathit{k}}^{\mathrm{2}}\frac{{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\left({\mathit{\phi }}_{\mathrm{d}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}\mathrm{Sc}}\right)\mathit{.}\end{array}$The eigenvalues of this matrix are found from: $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{\Gamma }}_{\mathrm{\pm }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{(}{\mathit{M}}_{\mathit{aa}}\mathrm{+}{\mathit{M}}_{\mathit{bb}}\mathrm{)}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{\pm }\sqrt{\mathrm{(}{\mathit{M}}_{\mathit{aa}}\mathrm{+}{\mathit{M}}_{\mathit{bb}}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{4}\mathrm{(}{\mathit{M}}_{\mathit{aa}}{\mathit{M}}_{\mathit{bb}}\mathrm{-}{\mathit{M}}_{\mathit{ba}}{\mathit{M}}_{\mathit{ab}}\mathrm{)}})\mathit{.}\end{array}& & \end{array}$(37)So we have $\mathrm{i}\mathit{\omega }\mathrm{=}{\mathrm{\Gamma }}_{\mathrm{\pm }}\mathit{.}$(38)The solution is stable if $\mathrm{Re}\mathrm{\left(}\mathrm{i}\mathit{\omega }\mathrm{\right)}\mathrm{\le }\mathrm{0}\mathit{,}$(39)for all possible values of k. One can see that if this is true/untrue for one value of k, it is true/untrue for all values of k, because k only enters as a multiplicative factor. The above stability condition requires that both M_aa + M_bb< 0 and M_aaM_bb>M_abM_ba. These two conditions simplify to: $\begin{array}{ccc}\mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\left({\mathit{\phi }}_{\mathrm{d}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{3}\mathrm{Sc}}\right)& \mathrm{\ge }& \\ \mathrm{1}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\mathrm{+}\left(\mathrm{1}\mathrm{+}\frac{\mathrm{St}}{{\mathit{\alpha }}_{\mathrm{1}}}\hspace{0.17em}\mathrm{Sc}\right){\mathit{\phi }}_{\mathrm{d}}& \mathrm{\ge }& \end{array}$Both conditions have to be fulfilled for the disk to be stable against the dust-driven viscous instability. Otherwise it is unstable for all k.

3.2. Results for the growth rates

We apply the model to the case of a young solar mass star (M_∗ = M_⊙) with still substantial luminosity (L_∗ = 10 L_⊙), to mimic the case of HL Tau (ALMA Partnership et al. 2015). We perform the analysis at a radius of r₀ = 60 au. The Schmidt number of the gas is set to Sc = 1. To compute the temperature of the midplane gas we assume a simple irradiated disk model, in which the irradiation angle is ϕ = 0.05. By setting the midplane temperature to the effective temperature of the disk assuming thermal equilibrium we obtain $\mathit{T}\mathrm{=}\mathrm{(}\mathit{\varphi }{\mathit{L}}_{\mathrm{\ast }}\mathit{/}\mathrm{(}\mathrm{4}\mathit{\pi }{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{SB}}\mathrm{)}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathrm{=}\mathrm{43}\hspace{0.17em}\mathrm{K}$, which is a very rough estimate of the disk temperature, but sufficient for the present purpose. The orbital time is 465 years. The vertical pressure scale height of the disk is H_p = c_s/ Ω_K = 6 au. We now apply the perturbation analysis for wave numbers k corresponding to dimensionless wavelengths λ = 2π/k in the range between the smallest possible wavelength λ = H_p/r₀ and the largest reasonable one λ = 1. The choice of H_p/r₀ is the smallest wavelength is based on the assumption that the turbulent viscosity in the disk cannot lead to radial structures that are narrower than about one vertical scale height. In the current example this means that the smallest wavelength we should consider is λ = 0.1.

We now compute the growth rate $\mathrm{\Gamma }\mathrm{=}\mathrm{Re}\mathrm{\left(}\mathrm{i}\mathit{\omega }\mathrm{\right)}\mathit{,}$(42)of each of these modes for a range of different dust particle sizes. We express the dust particle size in terms of its Stokes number St defined as St = t_stop/t_orbit, where t_stop is the stopping time of the dust particle defined as t_stop = f_fric/m_grain | v_grain−v_gas | where f_fric is the friction force between the gas and the dust particle, m_grain is the dust grain mass, and | v_grain−v_gas | is the absolute value of the velocity difference between the gas and the dust particle. The precise translation between particle mass m_grain and Stokes number is not trivial to express, because it depends on much detailed physics, such as the porosity or fractility of the dust aggregate, its size compared to the gas mean free path, the gas density and temperature etc. For typical disk parameters a Stokes number of unity at 60 AU would correspond to a compact silicate dust particle of about a centimeter or a decimeter, and is smaller for smaller particles. We refer to the literature for an in-depth discussion of the relation between particle size and Stokes number (see e.g. Birnstiel et al. 2010). For our analysis only the Stokes number is relevant. The results of the present analysis are shown in Fig. 2.

Fig. 2 Growth rate in units of the reciprocal orbital time (1 /torbit = ΩK/ 2π) of the dust-driven viscous instability, as a function of dust particle size (expressed as Stokes number St), according to the simplified analysis of Sect. 3. Dashed lines: case 1 (i.e. φg = 0), solid lines: case 2 (i.e. φg = −φd). The different lines show modes of different dimensionless wavelength λ in the dimensionless coordinate x, where λ = 2π/k. A dimensionless wavelength λ = 1 means a wavelength as large as r0. The parameters of the model shown here are φd = −1.0, α1 = 10-4, M∗ = 1 M⊙, L∗ = 10 L⊙, r0 = 60 au, T = 43,K, Sc = 1 (i.e. cs = 0.39 km s-1 and torbit = 465 yr).

The growth rate Γ = Re(iω) is expressed in terms of the reciprocal orbital time scale. This means that if this value is larger than unity, the perturbation grows faster than the gas can orbit around the star. This would not lead to “grand-design” rings, but instead to small arc-shaped clumps, because if a perturbation is triggered at some azimuthal position, the information about this event does not have time to propagate around the entire orbit before the perturbation has grown to much larger amplitude. In other words: to create global-scale rings we need a slow instability (a “secular instability”). It must be slow compared to the time it takes for a perturbation to shear out over 2π in azimuth. This time scale depends on the radial width of the perturbation, which is related to the dimensionless wavelength λ of the unstable mode through Δr = r₀λ. Through keplerian orbital dynamics we can then define this shear time scale as ${\mathit{t}}_{\mathrm{shear}}\mathrm{=}\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathrm{1}}{\mathit{\lambda }}\hspace{0.17em}{\mathit{t}}_{\mathrm{orbit}}\mathit{.}$(43)So instead of comparing Γ to the reciprocal orbital time, we should express Γ in terms of 1 /t_shear. If we do so, the curves for small λ will move up. The results are shown in Fig. 3.

Fig. 3 Same as Fig. 2, but now with the growth rate Γ expressed in units of the reciprocal shear time over a radial distance of Δr = r0λ.

Whereever the curve lies above unity, the growth is faster than the azimuthal communication. In that case small-scale arcs form instead of global scale rings. Wherever the curve lies sufficiently below unity but above zero (which in this log-representation means that the curve is visible in the plot), the perturbation may lead to large scale rings.

One can also see that, for case 2 (φ_g = −φ_d, solid lines in the figure) the instability does not operate for St ≤ 10^-4 (for this set of model parameters), consistent with Eqs. (40), (41). This can be understood because for very small grains (small St) the dust is so well-coupled to the gas that dust drift is virtually inhibited, meaning that Σ_d/ Σ_g remains constant.

For both case 1 and case 2, however, the instability does not occur for St → 0, i.e. for the case in which the dust does not drift at all. This changes if we set φ_d< −1. In Fig. 4 the results for φ_d = −2 are shown (both case 1 and case 2). Now, at least for case 1 (φ_g = 0, dashed lines), the instability even operates for St → 0, i.e. without dust drift. The reason is that the convergent flow of gas, dragging along the dust with it, increases the gas and dust density enough to set the instability in motion. We then recover the instability by Hasegawa & Takeuchi (2015) and others. We discuss this in Sect. 5.

Fig. 4 Same as Fig. 3, but now for φd = −2 instead of φd = −1.

As can be seen, however, the strongest growth occurs around Stokes numbers of unity, and for the shortest wavelength λ. In this regime the growth rate is faster than 1 /t_shear. In a disk with a grain size distribution spanning from tiny to large, this seems to suggest that the instability will be mainly driven by the comparetively large St ≃ 1 grains, which would then lead not to rings but to numerous small arcs. Perhaps these arcs can later merge into large scale rings is something that cannot be studied using this linear perturbation analysis.

In reality the situation is likely more subtle. In a disk with a dust size distribution it is typically the smallest grains that are affecting the α the most, because the smallest grains have the largest total surface area and can thus be most effective in removing free electrons and ions from the gas. We therefore speculate that in spite of the strong growth rate for St ≃ 1 particles that results from our analysis, it is mostly the smallest dust grains that drive the instability, if at all. If most/all of the dust has Stokes numbers below the cut-off for case 1, then if case 1 is applicable the instability would not operate at all.

3.3. A speculative scenario

Let us speculate about the following scenario: We assume that only relatively small dust affects the viscosity parameter α. From Fig. 3 for, say, St ≃ 3 × 10^-4 the instability is driven at a low enough rate, even for the smallest wavelengths, that a set of global rings can form. As these rings grow in strength, we will enter into the non-linear regime. The radial derivative of the gas pressure will start to display sign-changes, and thus form actual dust traps. Since the disk also has large grains (even though they did not participate in the instability), these large grains get trapped into the dust traps and form dense dust rings, possibly even dominating the local density over the gas density. At this point the frictional back-reaction of the dust onto the gas will have to be taken into account, and the streaming instability (Johansen & Youdin 2007) may set in within these dust rings. Also the self-gravity of the population of large grains may start to play a role. Perhaps a combination with the secular instability of Takahashi & Inutsuka (2016) could occur. We are aware that these are mere speculations, and more investigation (in particular: numerical modeling) is required.

So far we have only looked at the growth rates of the modes, not their spatial propagation. In other words, we looked at Re(iω) but not yet at Im(iω). The above speculative scenario is only possible if the initial ring-instability occurs more or less in situ, or in other words, that it is not a moving wave that is slowly amplifying but instead a standing wave that is growing in amplitude. To verify this we need to study the ratio Im(iω)/Re(iω) for all cases where Re(iω) > 0. For the simplified analysis of this section it turns out that Im(iω) = 0. The mode grows exactly in-situ, meaning that the above speculative scenario is plausible.

4. Full perturbation analysis

The perturbation analysis for the full system of equations of Sect. 2 is substantially more tedious than the simplified analysis of Sect. 3, but the results are overall consistent with each other. There are also some simplifications that we keep: we still assume that the Stokes number does not change with time and space, and the same holds for the Schmidt number. In reality, for a given particle size the Stokes number changes if the gas density changes. Such effects are not included.

4.1. Stationary powerlaw solution

Let us assume the following Ansatz for the stationary solution: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \\ {\mathrm{\Sigma }}_{\mathrm{d}\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \end{array}$where we deliberately took the same powerlaw index for both the dust and the gas component. For the temperature profile we also assume a powerlaw of the form $\mathit{T}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{T}}_{\mathrm{0}}{\left(\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{0}}}\right)}^{\mathit{q}}\mathit{.}$(46)The isothermal sound speed c_s follows from this temperature by Eq. (7). The dimensionless vertical scale height h is defined as h = H_p/r, where the scale height H_p is given by Eq. (8).

The radial gas velocity (Eq. (4)) becomes ${\mathit{v}}_{\mathrm{rg}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\nu }}{\mathit{r}}\frac{\mathit{\partial }\ln \mathrm{\left(}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }\sqrt{\mathit{r}}\mathrm{\right)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\nu }}{\mathit{r}}\left(\mathit{p}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathit{\partial }\ln \hspace{0.17em}\mathit{\nu }}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\right)\mathit{,}$(47)where ν is given by Eq. (5). It can be aposteriori verified that stationary powerlaw solutions only exist if the α coefficient is independent of radius, which we will, from here on, assume to be the case. With Eq. (5) we then obtain ${\mathit{v}}_{\mathrm{rg}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\nu }}{\mathit{r}}\left(\mathit{p}\mathrm{+}\mathit{q}\mathrm{+}\mathrm{2}\right)\mathit{.}$(48)Inserting this into the gas continuity equation (Eq. (2)), and setting the time-derivative to zero, yields Σ_gν = constant. This means (with Eqs. (5), (6), (7), (46), (44)) that $\mathit{p}\mathrm{+}\mathit{q}\mathrm{=}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\mathit{.}$(49)Inserting this into Eq. (48) yields ${\mathit{v}}_{\mathrm{rg}}\mathrm{=}\mathrm{-}\frac{\mathrm{3}\mathit{\nu }}{\mathrm{2}\mathit{r}}\mathit{.}$(50)Now let us do the dust, Eq. (3). Since by our Ansatz Σ_d(r)/Σ_g(r) is a constant (because both have the same powerlaw index), the right-hand-side of Eq. (3) is zero. This then immediately means that v_d(r) /v_g(r) must also be a constant. If we now look at the two terms in Eq. (10), and we assume that v_d(r) must have the same radial powerlaw dependency as v_g(r), then both terms in Eq. (10) must have the same radial powerlaw dependency. This means that c_sh must have the same radial powerlaw dependency on r as v_rg: $\frac{\mathrm{d}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}\mathit{h}\mathrm{\right)}}{\mathrm{d}\ln \hspace{0.17em}\mathit{r}}\mathrm{=}\frac{\mathrm{d}\ln \hspace{0.17em}\mathrm{\left|}{\mathit{v}}_{\mathrm{rg}}\mathrm{\right|}}{\mathrm{d}\ln \hspace{0.17em}\mathit{r}}\mathrm{=}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}$(51)where we used Eqs. (50), (5) and the definition of q (Eq. (46)) in the last step. However, given that ${\mathit{c}}_{\mathrm{s}}\mathit{h}\mathrm{=}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}{\mathrm{\Omega }}_{\mathrm{K}}\mathit{r}$ we already independently know that $\frac{\mathrm{d}\ln \mathrm{\left(}{\mathit{c}}_{\mathrm{s}}\mathit{h}\mathrm{\right)}}{\mathrm{d}\ln \hspace{0.17em}\mathit{r}}\mathrm{=}\mathit{q}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{,}$(52)which confirms that we indeed have a stationary powerlaw solution for both the dust and the gas. We can now calculate the ratio of the dust radial velocity to the gas radial velocity using Eq. (10) with Eq. (44): $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{v}}_{\mathrm{rd}}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{rg}}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{h}}{\mathrm{St}\mathrm{+}{\mathrm{St}}^{-1}}\left(\mathit{p}\mathrm{+}\frac{\mathit{q}\mathrm{-}\mathrm{3}}{\mathrm{2}}\right)\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\left[\mathrm{1}\mathrm{+}\mathit{L}\left(\mathit{p}\mathrm{+}\frac{\mathit{q}\mathrm{-}\mathrm{3}}{\mathrm{2}}\right)\right]{\mathit{v}}_{\mathrm{rg}}\mathit{,}\end{array}& & \end{array}$(53)where we define $\mathit{L}\mathrm{\equiv }\frac{\mathrm{St}\hspace{0.17em}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathrm{\Omega }}_{\mathrm{K}}\hspace{0.17em}\mathit{r}\hspace{0.17em}{\mathit{v}}_{\mathrm{rg}}}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}}\frac{\mathrm{St}}{\mathit{\alpha }}\mathit{,}$(54)where in the last step we used the stationary solution for v_rg (Eq. (50)) and the equation for ν (Eq. (5)).

If we insert a standard example, p = −1, q = −0.5, then this becomes: ${\mathit{v}}_{\mathrm{rd}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\left[\mathrm{1}\mathrm{+}\mathrm{1.833}\hspace{0.17em}\frac{\mathrm{St}}{\mathit{\alpha }}\right]{\mathit{v}}_{\mathrm{rg}}\mathit{.}$(55)

4.2. Linearization

Now we impose a perturbation on the stationary solution. We introduce the coordinate x: $\mathit{r}\mathrm{=}{\mathit{r}}_{\mathrm{0}}{\mathit{e}}^{\mathit{x}}\mathrm{\simeq }{\mathit{r}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}\mathrm{)}\mathit{,}$(56)and the perturbations: $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \\ {\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \end{array}$similar to Sect. 3. Again, the subscript 1 denotes the stationary solution. We keep the temperature and sound speed stationary (i.e. static in time, but varying in space). The perturbations σ_d(x,t) and σ_g(x,t) are allowed to affect α(x,t). We use the same recipe for α as before (Eq. (11)). To first order in the perturbations we can write: $\begin{array}{ccc}\mathrm{d}\ln \hspace{0.17em}\mathit{r}& \mathrm{=}& \\ \mathrm{d}\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \mathrm{d}\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathrm{d}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}\\ \mathrm{d}\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \mathrm{d}\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{d}\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}\mathrm{d}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}\\ \mathrm{d}\ln \left(\frac{{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}{{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{r,t}\mathrm{\right)}}\right)& \mathrm{=}& \mathrm{d}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathrm{-}\mathrm{d}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}\\ \mathrm{d}\ln \hspace{0.17em}\mathit{\alpha }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& {\mathit{\phi }}_{\mathrm{d}}\hspace{0.17em}\mathrm{d}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\hspace{0.17em}\mathrm{d}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\mathit{,}\\ \mathrm{d}\ln \hspace{0.17em}\mathit{\nu }\mathrm{\left(}\mathit{r,t}\mathrm{\right)}& \mathrm{=}& \mathrm{d}\ln \hspace{0.17em}{\mathit{\nu }}_{\mathrm{1}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{+}{\mathit{\phi }}_{\mathrm{d}}\hspace{0.17em}\mathrm{d}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{\left(}\mathit{x,t}\mathrm{\right)}\\ & & \end{array}$Or specifically for the double-logarithmic derivative with respect to r we obtain (omitting the (r,t) for notational convenience): $\begin{array}{ccc}\frac{\mathit{\partial }\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \\ \frac{\mathit{\partial }\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{d}}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \\ \frac{\mathit{\partial }\ln \mathrm{(}{\mathrm{\Sigma }}_{\mathrm{d}}\mathit{/}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{-}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathit{,}\\ \frac{\mathit{\partial }\ln \hspace{0.17em}\mathit{\alpha }}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& {\mathit{\phi }}_{\mathrm{d}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}{\mathit{\phi }}_{\mathrm{g}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathit{,}\\ \frac{\mathit{\partial }\ln \hspace{0.17em}\mathit{\nu }}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \end{array}$The gas velocity v_rg (Eq. (4)) is now: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{rg}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}\mathit{\nu }}{\mathit{r}}\frac{\mathit{\partial }\ln \mathrm{\left(}{\mathrm{\Sigma }}_{\mathrm{g}}\mathit{\nu }\sqrt{\mathit{r}}\mathrm{\right)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\\ & \mathrm{=}& \end{array}$where we used Eqs. (69), (49). The dust velocity v_rd (Eq. (10)) becomes, using Eqs. (6)–(8), (44), (65) and the identities h = H_p/r and $\mathit{P}\mathrm{=}{\mathit{\rho }}_{\mathrm{g}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}$: $\begin{array}{ccc}{\mathit{v}}_{\mathrm{rd}}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rg}}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\mathrm{+}\frac{{\mathit{c}}_{\mathrm{s}}\mathit{h}}{\mathrm{St}\mathrm{+}{\mathrm{St}}^{-1}}\left(\frac{\mathit{\partial }\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\mathrm{+}\frac{\mathit{q}\mathrm{-}\mathrm{3}}{\mathrm{2}}\right)\\ & \mathrm{=}& \end{array}$To be able to use v_rg and v_rd in the viscous disk equations Eqs. (2), (3), will be forced to compute their radial derivatives, which is where the cumbersome math comes in. To keep things as orderly as possible, we rewrite Eqs. (2), (3) into double-logarithmic form: $\begin{array}{ccc}\frac{\mathit{\partial }\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{g}}}{\mathit{\partial t}}\mathrm{+}\frac{{\mathit{v}}_{\mathrm{rg}}}{\mathit{r}}\frac{\mathit{\partial }\ln \mathrm{\left(}\mathit{r}{\mathrm{\Sigma }}_{\mathrm{g}}\mathrm{\right|}{\mathit{v}}_{\mathrm{rg}}\mathrm{\left|}\mathrm{\right)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \\ \frac{\mathit{\partial }\ln \hspace{0.17em}{\mathrm{\Sigma }}_{\mathrm{d}}}{\mathit{\partial t}}\mathrm{+}\frac{{\mathit{v}}_{\mathrm{rd}}}{\mathit{r}}\frac{\mathit{\partial }\ln \mathrm{\left(}\mathit{r}{\mathrm{\Sigma }}_{\mathrm{d}}\mathrm{\right|}{\mathit{v}}_{\mathrm{rd}}\mathrm{\left|}\mathrm{\right)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \\ \frac{\mathrm{1}}{\mathit{r}{\mathrm{\Sigma }}_{\mathrm{d}}}\frac{\mathit{\partial }}{\mathit{\partial r}}(\mathit{r}{\mathrm{𝒟}}_{\mathrm{d}}& {\mathrm{\Sigma }}_{\mathrm{d}}& \end{array}$The double-logarithmic derivatives of rΣ_g/dv_{r g/d} can, using Eqs. (65), (66), be written out as $\frac{\mathit{\partial }\ln \mathrm{\left(}\mathit{r}{\mathrm{\Sigma }}_{\mathrm{g}\mathit{/}\mathrm{d}}\mathrm{\right|}{\mathit{v}}_{\mathrm{rg}\mathit{/}\mathrm{d}}\mathrm{\left|}\mathrm{\right)}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\mathrm{=}\mathrm{1}\mathrm{+}\mathit{p}\mathrm{+}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}\mathit{/}\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\frac{\mathit{\partial }\ln \hspace{0.17em}\mathrm{\left|}{\mathit{v}}_{\mathrm{rg}\mathit{/}\mathrm{d}}\mathrm{\right|}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}\mathit{.}$(76)The right-hand-side of Eq. (75) can also be written into the derivatives of the perturbations. To first order in σ_g and σ_d we get: $\frac{\mathrm{1}}{\mathit{r}{\mathrm{\Sigma }}_{\mathrm{d}}}\frac{\mathit{\partial }}{\mathit{\partial r}}(\mathit{r}{\mathrm{𝒟}}_{\mathrm{d}}{\mathrm{\Sigma }}_{\mathrm{d}}\frac{\mathit{\partial }}{\mathit{\partial r}}\ln (\frac{{\mathrm{\Sigma }}_{\mathrm{d}}}{{\mathrm{\Sigma }}_{\mathrm{g}}}))\mathrm{=}\frac{{\mathrm{𝒟}}_{\mathrm{d}}}{{\mathit{r}}^{\mathrm{2}}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{(}{\mathit{\sigma }}_{\mathrm{d}}\mathrm{-}{\mathit{\sigma }}_{\mathrm{g}}\mathrm{)}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathit{,}$(77)where we made use of the fact that ∂(σ_d−σ_g) /∂x is already first order in σ_g and σ_d, and that is, to first order in σ_g and σ_d, constant. This is because the stationary solution (Sect. 4.1) obeys , which is constant.

What remains to be done is to derive expressions for the double-logarithmic derivatives of the gas and dust velocities (Eqs. (71) and (73), respectively) used in Eq. (76). This is somewhat tedious algebra, which we defer to Appendix A. After inserting the resulting expressions (Eqs. (A.3), (A.4)), into Eq. (76), we see that for both the gas and the dust version the constant term 1 + p drops out, and the expression of Eq. (76) becomes linear in the perturbations. This cancellation of the constant 1 + p is not surprising, because it follows from the fact that we start our perturbation analysis from the stationary solutions of Sect. 4.1. Inserting the resulting gas- and dust-versions of Eq. (76), together with Eq. (77), into the continuity equations Eqs. (74), (75), we find the following set of equations: $\begin{array}{ccc}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial t}}& \mathrm{+}& [\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{g}}\end{array}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial x}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gg}}\end{array}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gd}}\end{array}\hspace{0.17em}{\mathit{\phi }}_{\mathrm{d}}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\\ & & \\ \frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial t}}& \mathrm{+}& [\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{d}}\end{array}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dg}}\end{array}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dd}}\end{array}\hspace{0.17em}{\mathit{\phi }}_{\mathrm{d}}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\\ & & \end{array}$where the tilde-symbols are defined as: $\begin{array}{ccc}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{g}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rg}\mathrm{1}}}{\mathit{r}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{d}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rd}\mathrm{1}}}{\mathit{r}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gg}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rg}\mathrm{1}}}{\mathit{r}}{\mathit{C}}_{\mathrm{gg}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gd}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rg}\mathrm{1}}}{\mathit{r}}{\mathit{C}}_{\mathrm{gd}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dg}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rd}\mathrm{1}}}{\mathit{r}}{\mathit{C}}_{\mathrm{dg}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dd}}\end{array}& \mathrm{=}& \frac{{\mathit{v}}_{\mathrm{rd}\mathrm{1}}}{\mathit{r}}{\mathit{C}}_{\mathrm{dd}}\mathit{,}\\ \begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒟}}_{\mathrm{d}}\end{array}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{2}}}{\mathrm{𝒟}}_{\mathrm{d}}\mathit{,}\end{array}$where the symbols C_gg, C_gd, C_dg and C_dd are defined in Eqs. (A.5)–(A.8). We again insert trial functions $\begin{array}{ccc}{\mathit{\sigma }}_{\mathrm{g}}& \mathrm{=}& \mathit{A}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathit{kx}}\mathit{,}\\ {\mathit{\sigma }}_{\mathrm{d}}& \mathrm{=}& \mathit{B}{\mathrm{e}}^{\mathrm{i}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathit{kx}}\mathit{,}\end{array}$and obtain the matrix equation $\mathrm{i}\mathit{\omega }\left(\begin{array}{c}\\ \mathit{A}\\ \mathit{B}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ {\mathit{M}}_{\mathit{aa}}& {\mathit{M}}_{\mathit{ab}}\\ {\mathit{M}}_{\mathit{ba}}& {\mathit{M}}_{\mathit{bb}}\end{array}\right)\left(\begin{array}{c}\\ \mathit{A}\\ \mathit{B}\end{array}\right)\mathit{,}$(32)with $\begin{array}{ccc}{\mathit{M}}_{\mathit{aa}}& \mathrm{=}& \mathrm{i}\mathit{k}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{g}}\end{array}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gg}}\end{array}\mathrm{+}\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{i}\mathit{k}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gd}}\end{array}{\mathit{\phi }}_{\mathrm{g}}\mathit{,}\\ {\mathit{M}}_{\mathit{ab}}& \mathrm{=}& \mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{i}\mathit{k}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{gd}}\end{array}{\mathit{\phi }}_{\mathrm{d}}\mathit{,}\\ {\mathit{M}}_{\mathit{ba}}& \mathrm{=}& {\mathit{k}}^{\mathrm{2}}\mathrm{(}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dg}}\end{array}\mathrm{+}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒟}}_{\mathrm{d}}\end{array}\mathrm{)}\mathrm{+}\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{i}\mathit{k}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dd}}\end{array}{\mathit{\phi }}_{\mathrm{g}}\mathit{,}\\ {\mathit{M}}_{\mathit{bb}}& \mathrm{=}& \mathrm{i}\mathit{k}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{d}}\end{array}\mathrm{+}\mathrm{(}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{i}\mathit{k}\mathrm{+}{\mathit{k}}^{\mathrm{2}}\mathrm{)}\begin{array}{c}\mathrm{˜}\\ {\mathit{C}}_{\mathrm{dd}}\end{array}{\mathit{\phi }}_{\mathrm{d}}\mathrm{-}{\mathit{k}}^{\mathrm{2}}\begin{array}{c}\mathrm{˜}\\ {\mathrm{𝒟}}_{\mathrm{d}}\end{array}\mathit{.}\end{array}$The eigenvalues, and thereby the growth rates of the modes, follow from Eqs. (37), (38).

4.3. Results

For the same parameters as in Sect. 3.2 we plot the resulting growth rates for the full perturbation analysis. The result is shown in Fig. 5 for case 2.

Fig. 5 Results for the full perturbation analysis of Sect. 4 (solid lines) compared to the results of the simplified perturbation analysis of Sect. 3 (dashed lines). Here case 2 is shown (φg = −φd) with φd = −1. For the rest the figure is the same as Fig. 3. If we would have plotted case 1 instead of case 2, the difference would only be that the curves would not be cut off for St ≲ 10-4, but would continue down to St → 0 in the same fashion as the dashed lines in Fig. 3.

It shows that for small enough particles and small enough wavelength λ the simple perturbation analysis of Sect. 3 agrees well with the full perturbation analysis. The same is true if we would have plotted this diagram for case 1, the only difference to Fig. 5 being, that the curves would continue down to St → 0 as in Fig. 3. However, for larger particles and/or larger wavelength, the full perturbation analysis yields substantially weaker growth rates. But we see that for λ ≲ 0.25 the growth rates are nevertheless everywhere positive where the simplified analysis predicts positive growth rates. Since we need the instability to be slow to obtain large scale rings, this is, in fact, advantageous for the model.

As we did for the simplified analysis of Sect. 3.3 we have to verify if the initial ring-instability occurs more or less in situ. To this end we plot Im(iω)/Re(iω) (Fig. 6). We see that, in contrast to the simplified analysis, the imaginary component of iω is not zero. However, we see in the plot that for most Stokes numbers of interest Im(iω) is sufficiently much smaller than Re(iω). That means that the growth can be considered to be sufficiently well in situ for the speculative scenario of Sect. 3.3 to remain plausible.

Fig. 6 Ratio of Im(iω)/Re(iω) for the full perturbation analysis. Dashed lines: case 1 (φd = −1, φg = 0), solid lines: case 2 (φd = −1, φg = 1).

5. Discussion

5.1. Limitations and speculations

The linear stability analysis of this paper shows that, at least in principle, the combination of viscous disk theory, radial drift of dust, and the negative feedback of the dust on the viscosity, could lead to ring-shaped patterns in a protoplanetary disk.

The linear growth rate depends on the size of the dust grains, or more precisely: on their Stokes number. Depending on the prescription of the feedback, the instability is inhibited for the very small Stokes numbers. But for Stokes numbers beyond a critical value, the growth rate increases with increasing St. Around St ≃ 1 the instability is suppressed again. Figure 5 shows the growth rates as a function of St for given wavelengths of the mode.

Not surprisingly, the instability grows the quickest for the shortest wavelengths. The shortest wavelength is expected to be the disk pressure scale height, which is therefore expected to be the dominant mode.

However, the feedback of the dust onto the viscosity of the gas is a surface-area effect (dust grains removing free electrons and ions from the gas), so one should expect the feedback to be the strongest for the smallest grains. In our model we did not include this effect: we took the same feedback recipe (Eq. (11)) independent on grain size. We speculate that if a disk contains a size distribution of dust, the smallest grains with Stokes numbers still beyond the critical one, will be the ones that drive the instability. Since the growth rate of the instability for such grain sizes is substantially slower than the time a blob would be sheared out into a ring (see Fig. 5), the information about the growth of the instability can be communicated over the full 2π azimuth of the disk, so that a global ring is formed instead of a set of independent arcs.

With our model we cannot study what happens if the instability becomes non-linear. Assuming we have a size distribution of dust grains, then, although only the smaller grains drive the instability, also the larger grains will undergo density enhancements: they will do so even stronger than the instability-driving smaller grains. Once the mode becomes so strong that rings of positive pressure gradient are produced, then the larger grains will get trapped. Turbulent mixing always leaks a few of these grains out of the traps, but on the whole, large grains would be trapped and produce large-amplitude rings made from millimeter to centimeter size grains, similar to what is seen with ALMA in many sources. Our perturbation analysis suggests that the shortest wavelengths grow the quickest. However, viscous disk theory works on scales equal to or larger than the pressure scale height. The spacing between the rings in the gas structure of the disk will therefore be at least a pressure scale height. The large dust grains, however, could conceivably get trapped in rings that are thinner than that.

In this two-stage scenario (small grains triggering the rings, large grains getting trapped) the rings have to be sustained. If the small grains coagulate and become large, their ability to affect α reduces, and the rings may dissipate. Maybe this can be prevented through a bit of fragmentation of the pebbles, producing fresh fine grained dust. Even for low levels of turbulence the collision velocities of the pebbles may be relatively high. The typical turbulent eddy velocity at the top of the Kolmogorov cascade is ${\mathit{v}}_{\mathrm{eddy}}\mathrm{=}\sqrt{\mathit{\alpha }}{\mathit{c}}_{\mathrm{s}}$. With a temperature of, say, 100 K one has c_s = 0.6 km s^-1. If the fragmentation velocity is 1 m/s (e.g. Güttler et al. 2010) one would need α ≲ 3 × 10^-6 to prevent fragmentation. Anything above that would lead to the production of small grains. It is therefore feasible that a sufficient amount of small grains are continuously regenerated to keep the rings in place. On the other hand, aggregates made up of icy grains are thought to be more robust, and may fragment only at collision velocities of ~ 10 m/s, or even up to ~50 m/s, depending on the size of the monomers of which they are made (Wada et al. 2009)

Dust aggregates are most likely porous or fractal. Large “pebbles” may therefore still have rather large surface-to-mass ratios, and thus still strongly affect the ionization degree of the disk. We would, however, also observe them as if they were much smaller than they are, since also the optical properties of such dust aggregates depend a lot on the surface-to-mass ratio (Kataoka et al. 2014). The “big grains” that we identify as “big” due to their opacity slope at millimeter wavelengths therefore presumably also have a relatively low surface-to-mass ratio, and thus have little influence on the ionization degree of the disk. To have such big grains arranged in rings, we really need this two-stage scenario, as the big grains cannot (according to our scenario) generate the rings themselves.

There is, however, a big uncertainty with the model: the role of the vertical structure. Small grains can be turbulently stirred to several pressure scale heights above the midplane, even for relatively low turbulent α. Is it therefore still justified to assume that the wavelength of one pressure scale height to be the strongest growing mode? Could it be that this would lead to sufficient radial smearing that larger wavelength modes dominate? Okuzumi & Hirose (2011) study the effect of the vertical structure on the viscosity of the disk. They conclude that the vertically averaged magneto-turbulent viscosity only depends on the resistivity profile (and thereby on the vertically averaged dust abundance) through three critical heights, and is largely insensitive to the details of the resistivity profile itself. This may mean that the vertically averaged dust abundance has, in most parts of the disk, only limited influence (S. Okuzumi, priv. comm.).

Multi-wavelength observations of the same sources can help answer these questions. For instance, for TW Hydra both millimeter (Andrews et al. 2016) and H-band observations (van Boekel et al. 2017) exist. van Boekel et al. (2017) study how the rings in the millimeter compare to the large scale rings in the H-band. Some correspondence is found, but overall the structures appear to be uncorrelated, which appears to argue against our model, at least for this source.

5.2. On rings and the slowness of the instability

As argued in this paper, for an instability to lead to ring-like structures rather than patchy/clumpy structures in a disk, the instability has to be slower than the shear. This is the case for the dust-driven viscous instability discussed in this paper. In hindsight this is not surprising, since the viscous time scale of the disk can be quite long. One can quantify this by considering a perturbation with radial width W (i.e. in our dimensionless form this is W = λr). If we express W in units of the pressure scale height, which is the narrowest viscous structures we can expect in the disk, we get $\mathit{W}\mathrm{=}\mathit{w}{\mathit{H}}_{\mathrm{p}}\mathrm{=}\mathit{w}\frac{{\mathit{c}}_{\mathrm{s}}}{{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{,}$(93)The viscous time scale for this perturbation is ${\mathit{t}}_{\mathrm{visc}}\mathrm{=}\frac{{\mathit{W}}^{\mathrm{2}}}{\mathit{\nu }}\mathrm{=}\frac{{\mathit{w}}^{\mathrm{2}}}{\mathit{\alpha }{\mathrm{\Omega }}_{\mathrm{K}}}\mathit{.}$(94)The shear time (from Eq. (43)) is ${\mathit{t}}_{\mathrm{shear}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\frac{\mathit{r}}{\mathit{w}{\mathit{c}}_{\mathrm{s}}}\mathit{.}$(95)The condition for the instability to be slow enough then becomes $\begin{array}{ccc}\frac{{\mathit{t}}_{\mathrm{visc}}}{{\mathit{t}}_{\mathrm{shear}}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{w}}^{\mathrm{2}}}{\mathit{\alpha }}\frac{{\mathit{H}}_{\mathrm{p}}}{\mathit{r}}\mathrm{\gg }\mathrm{1.}& & \end{array}$(96)The smallest possible value for w is 1, which is also the fastest growing mode. This shows that if α is much smaller than the disk’s dimensionless thickness (aspect ratio), then the instability (if it exists) is slow enough to create rings instead of patches. For a typical disk H_p/r ≃ 0.05···0.1 this means that α ≪ 10^-2 for the instability to be slow enough. Therefore we can conclude that, if the rings seen in numerous disks are due to any form of viscous instability, the viscosity of the disk must be substantially lower than the canonical value of 10^-2, or the instability must be slowed down by an even slower process such as dust drift of small enough dust grains.

5.3. Comparison to earlier work

The ring-instability we have investigated in this paper appears to have a relation to the ring-instability found by Wünsch et al. (2005). In their model the disk had an active surface layer and a passive (“dead”) midplane layer. If gas would accumulate at some radius, the surface density of the active layer stays the same, but that of the dead layer increases. In a vertically averaged sense the viscosity thus gets reduced. This is mathematically identical to our case of φ_d = −1 and φ_g = 0 (prescription 1, Eqs. (11), (12)), with St = 0. In our model, however, we do not find an instability for these parameters. However, Wünsch et al. (2005) include the effect of the perturbation on the disk midplane temperature, which we do not. In our case we indeed get the instability if φ_d< −1.

Hasegawa & Takeuchi (2015) also studied the behavior of the two-layered disk model. They find that if the effective α of the two-layered disk is simply an average of the active and dead layers (weighted by their respective surface densities), then the disk remains stable. Our model with φ_d = −1 and φ_g = 0 (prescription 1) and St = 0 (no dust drift) confirms this. When they apply a more sophisticated effective α recipe, based on Okuzumi & Hirose (2011), they find that the disk becomes unstable near the dead zone outer edge, because that is where the dependence of α on Σ is the steepest. Our model confirms this, because the more sophisticated α recipe has a steeper dependence of α on Σ, which would amount, in our work, to φ_d< −1 (again taking φ_g = 0 and St = 0), which we confirm to lead to instability. Our model is thus consistent with the earlier work by Hasegawa & Takeuchi (2015). But by including dust drift our model is more general.

Flock et al. (2015) perform 3D full disk non-ideal MHD models and find ringlike structures, too, although rather wide ones and only two of them. But like Hasegawa & Takeuchi (2015), this is unrelated to dust drift.

Our dust-drift induced viscous instability is, however, very similar to the instability found by Johansen et al. (2011). While the analysis in that paper is locally more detailed (including radial and azimuthal motions), our linear stability analysis includes the radial gradients and the global cylindrical geometry terms an is thus not just a local analysis. Furthermore our analysis includes the radial drift, as well as the turbulent diffusion, both of which play a key role in the mechanism. Moreover, we argue that the slowness of the instability is critical in getting grand-design rings rather than chaotic arc-shaped structures.

The analysis in this paper is by no means a proof of the feasibility of this scenario. It will require detailed 2D/3D viscous hydrodynamic disk modeling, and comparisons to observations, to test this scenario.

6. Conclusion

In this paper we show that it is conceivable that the combination of viscous disk theory, radial drift of dust, and the negative feedback of the dust on the viscosity, can produce (or at least trigger the formation of) ring-shaped patterns in protoplanetary disks similar to what is seen in protoplanetary disks at millimeter and optical/near-infrared wavelengths. From our present analysis we conclude:

• 1.

  Even without dust drift, if ∂ln α/∂ln Σ < −1, the disk is prone to the viscous ringinstability. This is a conclusion in agreement with work byHasegawa & Takeuchi (2015).

• 2.

  When dust grains are large enough to start drifting, yet small enough to have a substantial influence on the viscosity of the disk (through their ability to capture free electrons and ions from the gas), dust drift tends to cause a feedback loop on the disk viscosity, leading to dust-rich regions of low viscosity and dust-poor regions of high viscosity. In this way dust drift can trigger the viscous instability even when the disk would be stable otherwise (∂ln α/∂ln Σ ≥ −1). These findings are consistent with Johansen et al. (2011), and generalize them to global disk accretion with generalized viscosity description.

• 3.

  The radial drift due to the global pressure gradient in the disk does not suppress the instability for small grains, but does so for grains with Stokes number near unity. The ring instability (viscous instability) must therefore be driven by small enough grains. This also agrees with the issue that small grains more easily affect the viscosity of the disk, because due to their larger surface-to-mass ratio, they more easily capture free electrons and ions.

• 4.

  For grand-design rings to form, such as those observed in real protoplanetary disks, rather than pseudo-random patchy structures, the growth rate of the instability must be slower than the shear-out time scale. Our analysis shows that this is generally the case, if the instability is driven by small enough grains and/or if the viscous α ≪ 10^-2.

• 5.

  If the grains are too small, they hardly drift. Whether the disk is then stable or not depends strongly on the viscosity recipe. We have identified two cases: case 1 in which the absolute value of the dust density determines the viscous α and case 2 in which the ratio of dust-to-gas density determines the viscous α. For very small grains, the disk becomes stable for case 2, but may still become unstable for case 1, if the dependence of α on 1/Σ is steep enough.

• 6.

  Observations seem to show that the rings seen at millimeter wavelengths are populated by relatively large rains. Since the viscous ring instability seems to be driven by small grains, this may seem inconsistent at first. We suggest that the initial rings are created by the viscous instability, which then produces strong enough dust traps that the larger grains get trapped and produce the ALMA images observed.

Acknowledgments

This research was supported by the Munich Institute for Astro- and Particle Physics (MIAPP) of the DFG cluster of excellence “Origin and Structure of the Universe”. We thank Satoshi Okuzumi for interesting discussions regarding the effect of the vertical structure. We thank the anonymous referee for insightful comments which helped improve the manuscript.

References

Appendix A: Double-logarithmic derivatives of the velocities

The double-logarithmic derivative of the gas velocity with respect to r is computed from Eq. (71) by initially taking the single-logarithmic-derivative ∂v_rg/∂ln r, working out ∂(ν/r)∂ln r, and dividing again by Eq. (71), where one regularly makes use of the fact that | σ_g/d | ≪ 1 and that we expand only to first order in σ_g/d.

For the double-logarithmic derivative of the dust velocity with respect to r we start from Eq. (73), and make use of the result we already obtained for the double-logarithmic derivative of v_rg. Again we regularly make use of the fact that | σ_g/d | ≪ 1 and that we expand only to first order in σ_g/d. To make the algebra more convenient we define $\begin{array}{ccc}\mathit{L}& \mathrm{=}& \\ \mathit{K}& \mathrm{=}& \end{array}$where Eq. (A.1) is, in fact, the same as Eq. (54).

After substantial algebra we find: $\begin{array}{ccc}\frac{\mathit{\partial }\ln \mathrm{\left|}{\mathit{v}}_{\mathrm{rg}}\mathrm{\right|}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathrm{gg}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{gd}}{\mathit{\phi }}_{\mathrm{d}}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\\ & & \\ \frac{\mathit{\partial }\ln \mathrm{\left|}{\mathit{v}}_{\mathrm{rd}}\mathrm{\right|}}{\mathit{\partial }\ln \hspace{0.17em}\mathit{r}}& \mathrm{=}& \mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{C}}_{\mathrm{dg}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\\ & & \mathrm{+}{\mathit{C}}_{\mathrm{dd}}{\mathit{\phi }}_{\mathrm{d}}\left(\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial x}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{d}}}{\mathit{\partial }{\mathit{x}}^{\mathrm{2}}}\right)\\ & & \end{array}$with $\begin{array}{ccc}{\mathit{C}}_{\mathrm{gg}}& \mathrm{=}& \\ {\mathit{C}}_{\mathrm{gd}}& \mathrm{=}& \\ {\mathit{C}}_{\mathrm{dg}}& \mathrm{=}& \\ {\mathit{C}}_{\mathrm{dd}}& \mathrm{=}& \end{array}$The second identities in Eqs. (A.7), (A.8) are for the case p = −1, q = −1/2.

All Figures

	Fig. 1 Mechanism of the ring instability studied in this paper.
In the text

Fig. 2 Growth rate in units of the reciprocal orbital time (1 /torbit = ΩK/ 2π) of the dust-driven viscous instability, as a function of dust particle size (expressed as Stokes number St), according to the simplified analysis of Sect. 3. Dashed lines: case 1 (i.e. φg = 0), solid lines: case 2 (i.e. φg = −φd). The different lines show modes of different dimensionless wavelength λ in the dimensionless coordinate x, where λ = 2π/k. A dimensionless wavelength λ = 1 means a wavelength as large as r0. The parameters of the model shown here are φd = −1.0, α1 = 10-4, M∗ = 1 M⊙, L∗ = 10 L⊙, r0 = 60 au, T = 43,K, Sc = 1 (i.e. cs = 0.39 km s-1 and torbit = 465 yr).

In the text

	Fig. 3 Same as Fig. 2, but now with the growth rate Γ expressed in units of the reciprocal shear time over a radial distance of Δr = r0λ.
In the text

	Fig. 4 Same as Fig. 3, but now for φd = −2 instead of φd = −1.
In the text

Fig. 5 Results for the full perturbation analysis of Sect. 4 (solid lines) compared to the results of the simplified perturbation analysis of Sect. 3 (dashed lines). Here case 2 is shown (φg = −φd) with φd = −1. For the rest the figure is the same as Fig. 3. If we would have plotted case 1 instead of case 2, the difference would only be that the curves would not be cut off for St ≲ 10-4, but would continue down to St → 0 in the same fashion as the dashed lines in Fig. 3.

In the text

	Fig. 6 Ratio of Im(iω)/Re(iω) for the full perturbation analysis. Dashed lines: case 1 (φd = −1, φg = 0), solid lines: case 2 (φd = −1, φg = 1).
In the text