A Double-logarithmic derivatives of the velocities

# Dust-driven viscous ring-instability in protoplanetary disks

###### Key Words.:
a

Protoplanetary disks often appear as multiple concentric rings in dust continuum emission maps and scattered light images. These features are often associated with possible young planets in these disks. Many non-planetary explanations have also been suggested, including snow lines, dead zones and secular gravitational instabilities in the dust. In this paper we suggest another potential origin. The presence of copious amounts of dust tends to strongly reduce the conductivity of the gas, thereby inhibiting the magneto-rotational instability, and thus reducing the turbulence in the disk. From viscous disk theory it is known that a disk tends to increase its surface density in regions where the viscosity (i.e. turbulence) is low. Local maxima in the gas pressure tend to attract dust through radial drift, increasing the dust content even more. We have investigated mathematically if this could potentially lead to a feedback loop in which a perturbation in the dust surface density could perturb the gas surface density, leading to increased dust drift and thus amplification of the dust perturbation and, as a consequence, the gas perturbation. We find that this is indeed possible, even for moderately small dust grain sizes, which drift less efficiently, but which are more likely to affect the gas ionization degree. We speculate that this instability could be triggered by the small dust population initially, and when the local pressure maxima are strong enough, the larger dust grains get trapped and lead to the familiar ring-like shapes. We also discuss the many uncertainties and limitations of this model.

ccretion, accretion disks – circumstellar matter – dust – stars: formation, pre-main-sequence – infrared: stars

## 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 . In smooth protoplanetary disks the gas pressure decreases with radius (), which is the reason for the inward drift of dust. But if the gas pressure has wiggles that are strong enough that they cause to flip sign, then there will exist local pressure maxima for which , 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 2-D 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 to cause this effect. It is just that the radial inward drift velocity 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 by adapting the surface density such that the steady-state

 Σg(r)ν(r)=const, (1)

is restored. So, whereever is reduced, is increased. This change in 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.

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

In Section 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 Section 3 we first simplify the system of equations radically, so that the mechanism presents itself more clearly. In Section 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:

 ∂Σg∂t+1r∂∂r(rΣgvrg) = 0 (2) ∂Σd∂t+1r∂∂r(rΣdvrd) = 1r∂∂r(rDdΣg∂∂r(ΣdΣg)). (3)

The gas radial velocity is given by the usual viscous disk equation:

 vrg=−3Σg√r∂∂r(Σgν√r), (4)

with the turbulent viscosity defined by

 ν=αc2sΩK, (5)

where is the usual alpha-turbulence parameter. The Kepler frequency is

 ΩK=√GM∗r3, (6)

with the gravitational constant and the stellar mass. The isothermal sound speed squared is

 c2s=kBTμmp, (7)

with the Boltzmann constant, the mean molecular weight, and the proton mass. The vertical pressure scale height of the disk is

 Hp=csΩK. (8)

The dust diffusion constant is

 Dd=11+St2Dg=11+St2νSc, (9)

where is the Stokes number of the dust and is the Schmidt number of the gas. The radial velocity of the dust is

 vrd=11+St2vrg+1St+St−11ρgΩK∂Pg∂r, (10)

where is the midplane gas density and 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:

 α=α1(ΣdΣd1)ϕd(ΣgΣg1)ϕg, (11)

where is the unperturbed value of , and likewise and are the unperturbed values of and respectively. The parameters and are powerlaw indices that parameterize how depends on the change in dust and/or gas surface density. We focus on cases with , meaning that an increase in leads to a decrease in . The is used to allow for the following two cases of interest:

 ϕg= 0 (Case 1: α depends on Σd) (12) ϕg= −ϕd (Case 2: α depends on Σd/Σg) (13)

Case 1 could, at least in principle, allow for the expected instability even without dust drift relative to the gas, since without dust drift, will necessarily increase if does. In contrast, case 2 strictly requires dust drift for the instability to operate, since lack of dust drift keeps constant.

## 3 Simplified perturbation analysis

It is cumbersome to perform the linear stability analysis for the full set of equations of Section 2, because these equations contain factors and the like. We postpone this full analysis to Section 4. For now, let us first simplify the equations.

### 3.1 Simplified equations

We simplify the equations to the following form:

 ∂Σg∂t+1r0∂∂x(Σgvxg) = 0, (14) ∂Σd∂t+1r0∂∂x(Σdvxd) = Missing or unrecognized delimiter for \right (15)

where is the radius at which we wish to locally carry out the perturbation analysis, and is the dimensionless radial coordinate starting from that location:

 r=r0(1+x). (16)

The gas radial velocity formula (Eq. 4) is simplified as:

 vxg=−3Σgr0∂(Σgν)∂x, (17)

with the viscosity still defined by Eq. (5). However, for simplicity and are now assumed be be constant. On the other hand, is allowed to depend on , but only due to the perturbation. The radial velocity of the dust is given, in our simplified description, by

 vxd=11+St2vxg+1St+St−1c2sΩKr0∂lnΣg∂x. (18)

The stationary solution is constant, constant and constant. This also yields for that stationary solution. This stationary solution is the backdrop of our perturbation analysis.

Now we introduce an infinitesimal perturbation:

 Σg(x,t) = Σg1(1+σg(x,t)), (19) Σd(x,t) = Σd1(1+σd(x,t)), (20)

where and are the stationary (constant) solutions for gas and dust, respectively. From here on the subscript denotes this stationary solution. We will also omit the notation, to not clutter the equations too much. The symbols and are the dimensionless infinitesimal perturbations on the gas and the dust respectively. For the , according to Eq. (11), we obtain, to leading order:

 α=α1(1+ϕdσd+ϕgσg), (21)

on account of the fact that .

Inserting these formulae into Eqs. (14, 15), and keeping in mind that both and are linear in the perturbations, and that we omit all terms of higher order in , we arrive, to leading order, at:

 ∂σg∂t+1r0∂vxg∂x = 0, (22) ∂σd∂t+1r0∂vxd∂x = 1r20Dd∂2(σd−σg)∂x2. (23)

By inserting Eqs. (5, 19) together with Eq. (21) into Eq. (17), and assuming that and are constant (see above), the radial velocity for the gas becomes

 vxg=−3ν1r0(∂σg∂x+∂(ϕdσd+ϕgσg)∂x). (24)

Similarly that of the dust becomes

 vxd=11+St2ν1r0{−3(∂σg∂x+∂(ϕdσd+ϕgσg)∂x)+Stα1∂σg∂x}. (25)

Inserting these into the continuity equations (Eqs. 22, 23) yields

 ∂σg∂t = 3ν1r20(∂2σg∂x2+∂2(ϕdσd+ϕgσg)∂x2), (26) ∂σd∂t = 11+St2ν1r20{3(∂2σg∂x2+∂2(ϕdσd+ϕgσg)∂x2) (27) −Stα1∂2σg∂x2}+1r20Dd1∂2(σd−σg)∂x2.

Now let us assume the linear perturbations to be plane waves in space and time . 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 and time frequency for both modes:

 σg = Aeiωt−ikx, (28) σd = Beiωt−ikx. (29)

The complex amplitudes (for the gas) and (for the dust) can be set independently. Inserting this mode into the above set of equations yields:

 iωA = Missing or unrecognized delimiter for \right (30) iωB = Missing or unrecognized delimiter for \bigg (31) +1Sc(B−A)},

where we made use of Eq. (9) to replace . This can be put into matrix form:

 iω(AB)=(MaaMabMbaMbb)(AB), (32)

with

 Maa = −3k2ν1r20(1+ϕg), (33) Mab = −3k2ν1r20ϕd, (34) Mba = −3k2ν1r2011+St2((1+ϕg)−St3α1−13Sc), (35) Mbb = −3k2ν1r2011+St2(ϕd+13Sc). (36)

The eigenvalues of this matrix are found from:

 Missing or unrecognized delimiter for \bigg (37)

So we have

 iω=Γ±. (38)

The solution is stable if

 Re(iω)≤0, (39)

for all possible values of . One can see that if this is true/untrue for one value of , it is true/untrue for all values of , because only enters as a multiplicative factor. The above stability condition requires that both and . These two conditions simplify to:

 1+ϕg+11+St2(ϕd+13Sc) ≥ 0, (40) 1+ϕg+(1+Stα1Sc)ϕd ≥ 0. (41)

Both conditions have to be fulfilled for the disk to be stable against the dust-driven viscous instability. Otherwise it is unstable for all .

### 3.2 Results for the growth rates

We apply the model to the case of a young solar mass star () with still substantial luminosity (), to mimic the case of HL Tau (ALMA Partnership et al., 2015). We perform the analysis at a radius of . The Schmidt number of the gas is set to . To compute the temperature of the midplane gas we assume a simple irradiated disk model, in which the irradiation angle is . By setting the midplane temperature to the effective temperature of the disk assuming thermal equilibrium we obtain , 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 . We now apply the perturbation analysis for wave numbers corresponding to dimensionless wavelengths in the range between the smallest possible wavelength and the largest reasonable one . The choice of 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 .

We now compute the growth rate

 Γ=Re(iω), (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 defined as , where is the stopping time of the dust particle defined as where is the friction force between the gas and the dust particle, is the dust grain mass, and is the absolute value of the velocity difference between the gas and the dust particle. The precise translation between particle mass 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).

The growth rate 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 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 . Through keplerian orbital dynamics we can then define this shear time scale as

 tshear=231λtorbit. (43)

So instead of comparing to the reciprocal orbital time, we should express in terms of . If we do so, the curves for small will move up. The results are shown in Fig. 3.

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 (, solid lines in the figure) the instability does not operate for (for this set of model parameters), consistent with Eqs. (40, 41). This can be understood because for very small grains (small ) the dust is so well-coupled to the gas that dust drift is virtually inhibited, meaning that remains constant.

For both case 1 and case 2, however, the instability does not occur for , i.e. for the case in which the dust does not drift at all. This changes if we set . In Fig. 4 the results for are shown (both case 1 and case 2). Now, at least for case 1 (, dashed lines), the instability even operates for , 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 Section 5.

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 . 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 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 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, 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 but not yet at . 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 for all cases where . For the simplified analysis of this section it turns out that . 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 Section 2 is substantially more tedious than the simplified analysis of Section 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:

 Σg1(r) = Σg0(rr0)p, (44) Σd1(r) = Σd0(rr0)p, (45)

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

 T(r)=T0(rr0)q. (46)

The isothermal sound speed follows from this temperature by Eq. (7). The dimensionless vertical scale height is defined as , where the scale height is given by Eq. (8).

The radial gas velocity (Eq. 4) becomes

 Missing or unrecognized delimiter for \right (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

 vrg=−3νr(p+q+2). (48)

Inserting this into the gas continuity equation (Eq. 2), and setting the time-derivative to zero, yields constant. This means (with Eqs. 5, 6, 7, 46, 44) that

 p+q=−32. (49)

Inserting this into Eq. (48) yields

 vrg=−3ν2r. (50)

Now let us do the dust, Eq. (3). Since by our Ansatz is a constant (because both have the same powerlaw index), the right-hand-side of Eq. (3) is zero. This then immediately means that must also be a constant. If we now look at the two terms in Eq. (10), and we assume that must have the same radial powerlaw dependency as , then both terms in Eq. (10) must have the same radial powerlaw dependency. This means that must have the same radial powerlaw dependency on as :

 dln(csh)dlnr=dln|vrg|dlnr=q+12, (51)

where we used Eqs. (50, 5) and the definition of (Eq. 46) in the last step. However, given that we already independently know that

 dln(csh)dlnr=q+12, (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):

 Missing or unrecognized delimiter for \right (53)

where we define

 L≡Stc2sΩKrvrg=−23Stα, (54)

where in the last step we used the stationary solution for (Eq. 50) and the equation for (Eq. 5).

If we insert a standard example, , , then this becomes:

 vrd=11+St2[1+1.833Stα]vrg. (55)

### 4.2 Linearization

Now we impose a perturbation on the stationary solution. We introduce the coordinate :

 r=r0ex≃r0(1+x), (56)

and the perturbations:

 Σg(r,t) = Σg1(r)(1+σg(x,t)), (57) Σd(r,t) = Σd1(r)(1+σd(x,t)), (58)

similar to Section 3. Again, the subscript denotes the stationary solution. We keep the temperature and sound speed stationary (i.e. static in time, but varying in space). The perturbations and are allowed to affect . We use the same recipe for as before (Eq. 11). To first order in the perturbations we can write:

 dlnr = dx, (59) dlnΣg(r,t) = dlnΣg1(r)+dσg(x,t), (60) dlnΣd(r,t) = dlnΣd1(r)+dσd(x,t), (61) dln(Σd(r,t)Σg(r,t)) = dσd(x,t)−dσg(x,t), (62) dlnα(r,t) = ϕddσd(x,t)+ϕgdσg(x,t), (63) dlnν(r,t) = dlnν1(r)+ϕddσd(x,t) (64) +ϕgdσg(x,t).

Or specifically for the double-logarithmic derivative with respect to we obtain (omitting the for notational convenience):

 ∂lnΣg∂lnr = p+∂σg∂x, (65) ∂lnΣd∂lnr = p+∂σd∂x, (66) ∂ln(Σd/Σg)∂lnr = ∂σd∂x−∂σg∂x, (67) ∂lnα∂lnr = ϕd∂σd∂x+ϕg∂σg∂x, (68) ∂lnν∂lnr = q+32+ϕd∂σd∂x+ϕg∂σg∂x. (69)

The gas velocity (Eq. 4) is now:

 vrg = −3νr∂ln(Σgν√r)∂lnr (70) = −3ν2r(1+2(1+ϕg)∂σg∂x+2ϕd∂σd∂x), (71)

where we used Eqs. (69, 49). The dust velocity (Eq. 10) becomes, using Eqs.(6, 7, 8, 44, 65) and the identities and :

 vrd = vrg1+St2+cshSt+St−1(∂lnΣg∂lnr+q−32) (72) = vrg1+St2[1+Stcshvrg(p+∂σg∂x+q−32)]. (73)

To be able to use and 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:

 ∂lnΣg∂t+vrgr∂ln(rΣg|vrg|)∂lnr = 0, (74) ∂lnΣd∂t+vrdr∂ln(rΣd|vrd|)∂lnr = 1rΣd∂∂r(rDd Σd ∂∂rln(ΣdΣg)). (75)

The double-logarithmic derivatives of can, using Eqs. (65, 66), be written out as

 ∂ln(rΣg/d|vrg/d|)∂lnr=1+p+∂σg/d∂x+∂ln|vrg/d|∂lnr. (76)

The right-hand-side of Eq. (75) can also be written into the derivatives of the perturbations. To first order in and we get:

 1rΣd∂∂r(rDdΣd∂∂rln(ΣdΣg))=Ddr2∂2(σd−σg)∂x2, (77)

where we made use of the fact that is already first order in and , and that is, to first order in and , constant. This is because the stationary solution (Section 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. 100, 101), into Eq. (76), we see that for both the gas and the dust version the constant term drops out, and the expression of Eq. (76) becomes linear in the perturbations. This cancellation of the constant is not surprising, because it follows from the fact that we start our perturbation analysis from the stationary solutions of Section 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:

 ∂σg∂t + [~Cg∂σg∂x+~Cgg∂2σg∂x2+~Cgdϕd(12∂σd∂x+∂2σd∂x2) (78) +~Cgdϕg(12∂σg∂x+∂2σg∂x2)]=0, ∂σd∂t + [~Cd∂σd∂x+~Cdg∂2σg∂x2+~Cddϕd(12∂σd∂x+∂2σd∂x2) (79) +~Cddϕg(12∂σg∂x+∂2σg∂x2)]=~Dd∂2(σd−σg)∂x2,

where the tilde-symbols are defined as:

 ~Cg = vrg1r, (80) ~Cd = vrd1r, (81) ~Cgg = vrg1rCgg, (82) ~Cgd = vrg1rCgd, (83) ~Cdg = vrd1rCdg, (84) ~Cdd = vrd1rCdd, (85) ~Dd = 1r2Dd, (86)

where the symbols , , and are defined in Eqs. (102, 103, 104, 105). We again insert trial functions

 σg = Aeiωt−ikx, (87) σd = Beiωt−ikx, (88)

and obtain the matrix equation

 iω(AB)=(MaaMabMbaMbb)(AB), (89)

with

 Maa = ik~Cg+k2~Cgg+(12ik+k2)~Cgdϕg, (90) Mab = (12ik+k2)~Cgdϕd, (91) Mba = k2(~Cdg+~Dd)+(12ik+k2)~Cddϕg, (92) Mbb = ik~Cd+(12ik+k2)~Cddϕd−k2~Dd. (93)

The eigenvalues, and thereby the growth rates of the modes, follow from Eqs. (37, 38).

### 4.3 Results

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

It shows that for small enough particles and small enough wavelength the simple perturbation analysis of Section 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 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 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 Section 3.3 we have to verify if the initial ring-instability occurs more or less in situ. To this end we plot (Fig. 6). We see that, in contrast to the simplified analysis, the imaginary component of is not zero. However, we see in the plot that for most Stokes numbers of interest is sufficiently much smaller than . That means that the growth can be considered to be sufficiently well in-situ for the speculative scenario of Section 3.3 to remain plausible.

## 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 . Around the instability is suppressed again. Fig. 5 shows the growth rates as a function of 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 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 . With a temperature of, say, 100 K one has . If the fragmentation velocity is 1 m/s (e.g. Güttler et al., 2010) one would need 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 , or even up to , 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 (i.e. in our dimensionless form this is ). If we express in units of the pressure scale height, which is the narrowest viscous structures we can expect in the disk, we get

 W=wHp=wcsΩK, (94)

The viscous time scale for this perturbation is

 tvisc=W2ν=w2αΩK. (95)

The shear time (from Eq. 43) is

 tshear=4π3rwcs. (96)

The condition for the instability to be slow enough then becomes

 tvisctshear=34πw2αHpr≫1. (97)

The smallest possible value for 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 this means that 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 , 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