Filament Formation due to Diffusive Instabilities in Dusty Protoplanetary Disks

We report the finding of a new, local diffusion instability in a protoplanetary disk which can operate in a dust fluid, subject to mass diffusion, shear viscosity, and dust–gas drag, provided the diffusivity, viscosity, or both, decrease sufficiently rapidly with increasing dust surface mass density. We devise a vertically averaged, axisymmetric hydrodynamic model to describe a dense, midplane dust layer in a protoplanetary disk. The gas is modeled as a passive component, imposing an effective, diffusion-dependent pressure, mass diffusivity, and viscosity onto the otherwise collisionless dust fluid, via turbulence excited by the gas alone, or dust and gas in combination. In particular, we argue that such conditions are met when the dust–gas mixture generates small-scale turbulence through the streaming instability, as supported by recent measurements of dust mass diffusion slopes in simulations. We hypothesize that the newly discovered instability may be the origin of filamentary features, almost ubiquitously found in simulations of the streaming instability. In addition, our model allows for growing oscillatory modes, which operate in a similar fashion as the axisymmetric viscous overstability in dense planetary rings. However, it remains speculative if the required conditions for such modes can be met in protoplanetary disks.


INTRODUCTION
Protoplanetary disks are the birthplaces of planets.
One key stage within the core-accretion scenario for planet formation involves the conversion of small dust particles into km-sized planetesimals.The formation of planetesimals is associated with a multitude of challenges.Specifically, coagulational growth is thought to be inhibited at around meter sizes by both radial drift and fragmentation of small solids (Birnstiel et al. 2012;Blum 2018).
In its linear phase, the streaming instability utilizes the relative equilibrium velocity between dust and gas, in the classical picture induced by the background gas pressure gradient, to drive exponentially-growing modes (Youdin & Goodman 2005).The streaming instability saturates, after a few dynamic timescales, into a quasisteady-state characterized by turbulent particle density and velocity fluctuations (Johansen & Youdin 2007).Eventually, this system self-organizes into azimuthally elongated filaments which can drift inwards and merge (see e.g., Yang & Johansen 2014;Li et al. 2018;Li & Youdin 2021).The above three-step evolution has been readily observed in 3D shearing-box simulations of both vertically stratified and unstratified protoplanetary disks where drag and dust feedback is included.
The formation of planetesimals within the streaming instability framework requires the additional component of dust self-gravity, and albeit a priori not obvi-ous, is thought to occur during the streaming instability's non-linear phase, either before or after the emergence of the over-dense filaments.The non-linear phase has been investigated numerically on numerous occasions.Specifically, Schreiber & Klahr (2018) found in two-dimensional simulations that dust diffusivities tend to decrease with dust-to-gas ratio.This behavior was also seen in three-dimensional, stratified simulations by Gerbig & Li (2023), and is typically attributed to the particles carrying too much collective inertia to be effectively diffused away by residual gas turbulence, and conversely, the back-reaction of the particles' inertia onto the gas may lead to a decrease in diffusion with increasing dust-to-gas ratio.An alternative picture views particle diffusion similar to gas pressure (a model we aim to discuss thoroughly in this paper): a region of high diffusion expels particles towards regions of low diffusion.
Either way, the implication of particle diffusion decreasing with increasing dust density has hitherto not been investigated analytically in the context of the stability of dusty protoplanetary disks.Previous models by Chen & Lin (2020); Umurhan et al. (2020) had diffusion depend on stopping time and gas viscosity only, both of which are taken to be constant.In this paper, we perform an instability analysis of a sheet of particles, subject to dust-gas drag forces, and mass and momentum diffusion, where the diffusion coefficients are allowed to vary with particle density.The existence of such a dependence has been established in hydrodynamic simulations by Schreiber & Klahr (2018); Gerbig & Li (2023).Our treatment bears thus some similarity to hydrodynamic studies of the viscous instability (Lin & Bodenheimer 1981;Ward 1981;Salo & Schmidt 2010) and viscous overstability (Schmit & Tscharnuter 1995, 1999;Schmidt et al. 2001;Latter & Ogilvie 2009, 2010;Lehmann et al. 2017Lehmann et al. , 2019) ) in planetary rings.
This paper is structured as follows.We outline our hydrodynamical model, specifically focusing on the diffusion terms and their physical relevance, as well as perform a linear perturbation analysis in Sect. 2. Next, we discuss the arising non-oscillatory and overstable modes in Sects.3 and 4 respectively.We discuss and contextualize our results in Sect. 5. Lastly, Sect.6 concludes the paper with a summary of our findings.

Diffusion and viscosity in particle-laden protoplanetary disks
Dust diffusion has long been identified to be of immense importance for dust dynamics and consequently planetesimal formation (see e.g., Cuzzi et al. 1993).We deem it worth explicitly defining for this work the rele-vant properties terms and putting them into the context of previous studies related to particle diffusion.For a recent comprehensive discussion on turbulent diffusion in protoplanetary disks we refer to Binkert (2023).
Generally speaking, diffusion acts to minimize free energy.In this work, we describe dust as a fluid subject to diffusion of mass, driven by a gradient in dust concentration, and momentum, driven by pressure gradients and shear stresses.
Under typical conditions, dust particles in protoplanetary disks are not collisional, and therefore do not experience collisional pressure forces.Instead, their dynamics are influenced by their coupling to the gas, namely via their stopping time t s , which appears in the drag term of the momentum equation, i.e.F d ∝ (v − u)/t s , where v and u are particle and gas velocity respectively.If the gas turbulence were fully characterized by the gas velocities u, no additional diffusion terms would be needed in modeling dust diffusion in protoplanetary disks.Indeed, numerical simulations typically do not employ explicit diffusivity or viscosity, and instead compute diffusive effects indirectly via dust-gas interaction reflected in u (e.g., Yang et al. 2018;Riols et al. 2020).As such a treatment is often not practical for analytical progress, we employ a diffusion subgrid model and describe diffusion and viscosity due to the particles' coupling to the gas by using explicit terms in the hydrodynamical equations.

Governing equations
Specifically, in this work, we will consider an isothermal, infinitesimally thin, axisymmetric particle disk in the absence of self-gravity, embedded in a gas that enters the system through diffusion, viscosity and drag.In polar coordinates (r, ϕ), the system is governed by the set of vertically averaged fluid equations (1) (3) Eqs.
(1) to (3) describe the dynamical evolution of the surface mass density Σ, the radial velocity v r and azimuthal velocity v ϕ , respectively, where Ω = GM * /r 3 is the Keplerian angular frequency, with stellar mass M * and gravitational constant G.The continuity equation (1) takes the form of an advection-diffusion equation, with mass diffusivity for dust D. The momentum equations incorporate advection both by v r and by the diffusion flux.This leads to the modified gradient advection terms on the left hand sides of the momentum equations, as well as the modified curvature-related advection term on the right hand side of Eq. ( 3).In addition, the fifth term on the right hand side of Eq. ( 2) incorporates v r -advection of the momentum carried in the diffusion flux itself.We refer to Tominaga et al. (2019) for a discussion on these additional advection terms associated with the diffusion flux and their implications when the full dust-gas mixture is considered, but note that they allow the total gas and dust angular momentum to be conserved.We provide a more rigorous justification for our set of hydrodynamical equations in Appendix A, using meanfield theory based on Reynolds averaging and the application of a set of plausible closure relations.Given this context, the fields Σ, v should be interpreted as mean fields separated in scale from the underlying small fluctuations that characterize turbulence.
Eq. ( 2) includes the vertically averaged, effective dust pressure P d = Σc 2 d , with velocity dispersion c 2 d of the dust fluid.Since the dust is assumed to be collisionless, this effective velocity dispersion is assumed to be generated solely by the particles' coupling to the turbulence with c 2 d ∝ D. Specifically, we follow Klahr & Schreiber (2021); Gerbig & Li (2023) and write which follows from a balance between diffusion and sedimentation.The latter approximation requires D ≪ t s c 2 s , which is the case in numerical simulations (e.g., Schreiber & Klahr 2018;Gerbig & Li 2023).We discuss this pressure model in Appendix B.
Finally, we include explicit momentum diffusion terms, modeled by Navier-Stokes stress terms F r and F ϕ .They can be calculated via with viscous stress tensor T , the components of which are where ν is the effective vertically averaged, kinematic, shear viscosity of the particle fluid.The inclusion of these shear stress terms differentiates Eq. ( 2) and Eq. ( 3) from the dust momentum equations used by Tominaga et al. (2019).Our model also distinguishes itself from other works concerning dusty protoplanetary disks on the scales of planetesimal formation as we consider diffusivity D and viscosity ν, and consequently via Eq.( 4) also the velocity dispersion c d and dust pressure, to depend on surface mass density Σ.Specifically, we assert power law dependencies of the form where β diff and β visc are dimensionless exponents.We also introduce the corresponding dimensional slopes At high dust-to-gas ratios, β D has been found to be negative (Schreiber & Klahr 2018;Gerbig & Li 2023).We are not aware of any numerical constraints on β ν within the context of turbulent diffusion in dusty protoplanetary disks.

Model applicability
Equations.
(1)-( 3) can be applied if the combined inertia of the disk is dominated by the particle fluid, for dust-to-gas volume mass density ratios ρ p /ρ g ≳ 1.In this case, the presence of the gas can be reduced to perturbations that evoke drag, mass diffusion, and momentum diffusion (aka viscosity).While the analytic model itself is agnostic to the source of diffusion and viscosity, in this paper, we specifically apply it to particle layers in the midplane of a protoplanetary disk, subject to non-linear streaming instability.
Given this context, our model is restricted to radial length scales exceeding the characteristic scale of the underlying turbulence, in this case, the characteristic scale of streaming instability, which for τ s ∼ 1, is of order l SI ∼ ηr (e.g., Youdin & Goodman 2005;Squire & Hopkins 2018;Gerbig et al. 2020), where η ∼ 0.01 characterizes the radial pressure gradient in the disk, and thus scales with the equilibrium relative velocity between dust and gas.For smaller stopping times, this restriction is relaxed, as the characteristic scale of the linear streaming instability decreases (e.g.Lin & Youdin 2017, Appendix D).Also note that in the vertical direction, this restriction is formally always satisfied as our model is vertically unstratified, implying a vanishing vertical wavenumber of all modes.Whether or not such modes are supported by an actual vertically stratified dust layer requires a stratified analysis and is thus subject to future work.
We further point out that a fluid description for particles in protoplanetary disks as applied here, is strictly only valid if t s < Ω −1 , since for decoupled grains the dynamical evolution of the stress tensor can, in general, not be ignored (see e.g., Garaud et al. 2004;Jacquet et al. 2011), and must be modeled using a kinetic approach.In this work, we instead assume that the 'external' turbulence is able to establish a simple Newtonian stress-strain relation for the particle fluid, characterized by shear viscosity ν and isotropic velocity dispersion c d , as discussed in Sect.2.2 and Appendix A. This is to some extent similar to the effect of mutual particle collisions in planetary rings, which indeed, if frequent enough, are known to establish a Newtonian stress-strain relation of the particle flow (e.g.Stewart et al. 1984;Shu & Stewart 1985).However, streaming instability turbulence, which is the main application of our model, is not expected to occur for large stopping times, which withdraws the physical justification for this assumption (at least given this context).In addition, mutual collisions, which are ignored in our model, may in principle become relevant if the velocity dispersion c d ∝ t −1/2 s becomes sufficiently small.
Despite these limitations, we will also present and discuss results assuming larger particles with t s > Ω −1 .In doing so, we retain a concrete connection to the viscous instability and overstability in planetary rings.Also, if large grains in protoplanetary disks do experience momentum and mass diffusion by some means (see discussion in Sect.5.5), our model may still provide useful insights, despite lacking the stress tensor evolution contained in a kinetic approach.

Linearized equations and dispersion relation
We adopt a local, co-rotating Cartesian reference frame at distance R from the star such that (x, y) = (r − R, R(ϕ − Ωt)) and v x = v r , v y = v ϕ − RΩ.We then perturb the system around a background state such that x , v y = −qΩx + v ′ y with Σ 0 = const., and linearize in perturbed quantities.Following Appendix B in Klahr & Schreiber (2021), we neglect the perturbed gas velocity u ′ such that the linearized drag terms are ∝ −v ′ /t s , which is justified if the mean field quantities derived in Appendix A are time-averaged over one turbulent correlation time.This assumption conveniently decouples dust and gas equations and allows us to isolate the effects of dust density-dependent turbulence alone.
For Keplerian shear where q = 3/2, Eqs. ( 1 This set of linearized equations is novel in that it includes a Navier-Stokes viscosity for the particle fluid, relates the particle pressure to diffusion and stopping time via Eq.( 4), which produces the second term on the r.h.s. of the radial momentum equation, and takes into account the dependence of diffusion and viscosity on the particle surface mass density, as motivated by simulations of Schreiber & Klahr (2018) and Gerbig & Li (2023).As a consequence, the radial and azimuthal momentum equations respectively contain the slope of diffusion and viscosity with respect to particle surface mass density.Depending on the slope, these terms can act both stabilizing or destabilizing on perturbations to the equilibrium state defined above, as we discuss below.
The mass diffusion term in the continuity equation, the terms ∝ ∂/∂x 2 (assuming ν > 0), and the term describing advection of background shear by the diffusion flux (second term on r.h.s of Eq. ( 3)) are always stabilizing.Note, that this diffusion flux term is the only term from the four angular momentum conserving terms in Eqs.
(2) and (3) added by Tominaga et al. (2019) that survives linearization.Lastly, the drag terms have a stabilizing effect.In our analysis, we include both radial and azimuthal drag terms, which is in contrast to Klahr & Schreiber (2021) who drop the azimuthal drag term.We proceed by introducing axisymmetric modes of the form with complex frequency n and (radial) wavenumber k.
We take k > 0 without loss of generality.Modes grow and decay for ℜ(n) > 0 and ℜ(n) < 0 respectively, and ℑ(n) corresponds to the oscillation frequency, the sign of which sets the wave travel direction.We get This system is solved by a cubic dispersion relation of the form with coefficients (20)

Dimensionless quantities
It is convenient to write the dispersion relation in terms of commonly used dimensionless quantities.The orbital frequency introduces a time unit, such that we write the dimensionless stopping time as which, in general, is distinct from the so-called Stokes number defined as the ratio of stopping time and turbulent correlation time (also known as eddy time or integral time) (Cuzzi et al. 1993;Youdin & Lithwick 2007).Particles with τ s ≪ 1 are well-coupled to the gas; while τ s ≫ 1 applies to loosely-coupled dust.
Our reference length unit is the gas pressure scale height, which is the ratio of the (gas) sound speed c s and orbital frequency, H = c s /Ω.We write the dimensionless wave number as K ≡ kH.Dimensionless versions for ground state diffusivity in Eq.( 7) and viscosity in Eq. ( 8) are introduced as where we use the same nomenclature as the well-known Shakura Sunyaev α-parameter (Shakura & Sunyaev 1973), albeit in our work, α describes the effective viscosity of the dust fluid, and not the viscosity of the gas that is often adopted to model angular momentum transport in protoplanetary disks.The corresponding power law slopes of diffusivity and viscosity were defined in Eqs. ( 7) and ( 8 2018) measured a shear viscosity of α ∼ 10 −4 in the gas.However, that value includes (albeit presumably small) contributions from Maxwell stresses, and constitutes an average value throughout the particle layer.We are not aware of any other applicable constraints on α in the particle layer, or any measurements of β visc .
We define the hydrodynamic Schmidt number as the ratio of viscosity and mass diffusion coefficient for our particle fluid, i.e.
which is analogous to the definition used for pure gas in protoplanetary disks by Johansen & Klahr (2005); Carballido et al. (2006).As pointed out in Youdin & Lithwick (2007), in the context of particles in protoplanetary disks, the Schmidt number suffers from being over-subscribed, as it is also used to describe the ratio of gas viscosity and particle diffusivity (Cuzzi et al. 1993;Schreiber & Klahr 2018;Binkert 2023).Here, we are primarily interested in the dust component, hence we assign Sc to characterize the relative importance of particle viscosity and particle diffusivity, both of which stem from the particles' coupling to gas turbulence.Lastly, we notice that for δ > 0, the dispersion relation includes a degeneracy in wave number and diffusivity, such that we can define Note, that since the herein utilized description for dustpressure is only appropriate for τ s ≫ δ, we equivalently require ξ ≪ τ s K 2 .As noted in Sect.2.3, our model is only appropriate for scales larger than the scale of the under-lying turbulence.We thus require, ξ ≲ 4π 2 δ/((H/r) 2 η 2 ).For typical values of H/r = 0.1, η ∼ 1%, and for fiducial diffusivities of δ ∼ 10 −5 − 10 −4 , this corresponds to maximum ξ of order unity.This is why, for the remainder of this study, we will mostly be concerned with the long-wavelength limit where ξ < 1.
The dimensionless complex frequency is denoted as where γ = ℜ(σ) is the growth rate and ω = ℑ(σ) is the oscillation frequency.The dimensionless version of the cubic dispersion relation in Eq. ( 18) is given by with where A 0 , A 1 , A 2 are real coefficients.This implies that A 0 < 0 is a sufficient condition for non-oscillatory instability, as then there is at least one real positive root of Eq. ( 28).The full, complex dispersion relation can also be cast into two equations for the growth rate and oscillation frequency, which need to hold independently:

DIFFUSIVE INSTABILITY
This section investigates the diffusive instability associated with the real roots of the dispersion relation in Eq. ( 28).Consider, however, first the case with zero diffusivity and viscosity, i.e.Sc = 0, ξ = 0.In this limit, the above dispersion relation reduces to which is solved by damped, epicyclic waves with oscillation frequency Ω and negative growth rate −t −1 s .The purely real root is the static null solution with σ = γ = ω = 0.It is this mode that will get destabilized by diffusion and/or viscosity, as we discuss in the following.As this section is only concerned with the non-oscillatory, purely real solution, we set ω = 0, and replace σ = γ.We will redirect our attention to oscillatory modes in Sect. 4.
Growth rates for the inviscid (Sc = 0) diffusive instability driven by the diffusion-dependent pressure for ξ = 0.1.The black dashed line corresponds to Eq. ( 36), below which the instability can operate, and above which, perturbations are exponentially damped.
Growth rates for the inviscid (Sc = 0) diffusive instability driven by the diffusion-dependent pressure for β diff = −3 and various stopping times.The hatched region corresponds to ξ > 1, for which our model seizes to be appropriate.

Inviscid case
We first investigate the purely diffusive case, where Sc = 0.In this limit the dispersion relation can be written as which leads to unconditional instability if This value is of order the power law slopes in Schreiber & Klahr (2018); Gerbig & Li (2023).For small growth rates γ ≪ 1, we find from Eq. ( 35), which for small stopping times equals On the other hand.for large stopping times with τ s ≫ 1, growth rates are uninterestingly small.Fig. 1 shows the growth rate for the inviscid diffusive instability for choices of slope β diff and stopping time τ s for modes with ξ = 0.1, which for a fiducial δ = 10 −5 corresponds to K = 100.Fig. 2 shows the growth rates of the inviscid diffusive instability for β diff = −3 for various stopping times.The system is unstable for all ξ although growth rates become uninterestingly small for large stopping times.In fact, we note that inviscid diffusive instability is not damped at small scales, but displays fastest growth rates at scales ξ > 1, which is where our model breaks down, as indicated by the hatched region in Fig. 2.
The physical origin of the inviscid diffusive instability lies in the dust pressure term in Eq. (2).For β diff < 0, this pressure term acts destabilizing and accelerates particles towards annuli of high density and low diffusivity.If β diff < −2, this pressure forcing can overcome the stabilizing mass diffusion term in Eq. ( 1) and the drag terms in the momentum equations.This increases the density, in the process further decreasing diffusion, resulting in positive feedback and instability.While in the large stopping time limit, the drag terms vanish, the destabilizing pressure term does too.The only remaining terms are the mass diffusion term in the continuity equation, and the diffusion flux term in Eq. ( 13), both of which act to repel particles away from density maxima.Hence, only for small stopping times are the growth rates of this instability appreciably fast.Next, we consider the case where D = ν, equivalently Sc = 1, such that the equilibrium value of momentum diffusion equals that of the mass diffusion, and also the power law slopes of diffusion and viscosity are identical, i.e. β visc = β diff .Under these assumptions, the dispersion relation is

Sc = 1 case
Pure instability is achieved if In the long wavelength limit of ξ < 1, we recover our finding for the inviscid case where a power law gradient β diff < −2 suffices for instability, if stopping times are sufficiently small: τ s ≪ 1.This is because for small stopping times and at large radial length scales, the destabilizing pressure term dominates over the viscosity terms.
For large stopping times of τ s ≫ 1, the explicit drag terms and the pressure term become neglible.In this case, the long wavelength limit ξ < 1 yields The diffusive instability behaves now analogously to the classical viscous instability (Ward 1981;Lin & Bodenheimer 1981), in that the instability is driven by the density slope of shear viscosity that appears in the azimuthal momentum equation ( 13).If the slope is sufficiently steep, there is a net flux towards density maxima that amplifies the linear perturbation.The criterion for the classical viscous instability is given by β visc < −1.
It is easy to see in Eq. ( 31), that if τ s ≪ 1, the slope of mass diffusion will dominate, whereas for τ s ≫ 1, the slope in momentum diffusion will dominate, thus leading to diffusive instability driven by the pressure term and viscosity terms respectively.For marginally coupled particles τ s ∼ 1, the instability can utilize both slopes.Because of this, on large scales ξ ≲ 1, if β diff is sufficiently negative, the system is unstable regardless of stopping time.
For small growth rates γ ≪ 1, and for small ξ, the real root of Eq. ( 39) behaves as which in the limits of well-coupled and decoupled particles equals respectively.We thus recover Eq. ( 38), in the limit of small stopping times.Fig. 3 shows the real root of the full cubic in Eq. ( 39) for ξ = 0.1.The black curve corresponds to Eq. ( 40).For a given β diff , growth rates are greater for large stopping times.Fig. 4 depicts the growth rate of the viscousdiffusive instability for β diff = β visc = −3 for various stopping times.Unlike in the invisicd case, where if Eq. ( 36) is satisfied, all modes ξ are unstable (see Fig. 2), the viscous diffusive instability is damped at small scales (ξ ≳ 1) by viscosity.This regularizes the system, by prohibiting growth on arbitrarily small scales.While our model itself is not applicable to large ξ, one does not expect growth rates to increase without bound at ever-decreasing scales, which is an advantage of including viscosity in the model.

General case
Given that the numerical constraints on the effective viscosity in high-density particle mid-plane of protoplanetary disks are sparse, we finally investigate the most general case our model allows, where we remain agnostic to the value of Sc and retain two independent power law slopes in δ and α.As before, pure instability is achieved for A 0 < 0, which for small ξ results in the condition In the general case, the growth rate of the diffusive instability, as well as the existence thereof in the first place, thus depends on five parameters: the stopping time τ s , Schmidt number Sc, diffusion slope β diff , viscosity slope β visc and dimensionless wave number ξ.The Schmidt number Sc acts as an amplification to the viscosity slope in the context of the diffusive instability.Like in the previous two cases, well-coupled particles with τ s ≪ 1 result in instability if β diff < −2, on account of the diffusive instability.On the other hand, looselycoupled particles with τ s ≫ 1 are unstable if We recover the classical criterion for viscous instability in planetary rings (Ward 1981;Lin & Bodenheimer 1981) .
For pure instability with small growth rates, γ ≪ 1, Eq. 28 yields γ = −A 0 /A 1 , which for ξ < 1 results in In the well-coupled limit τ s ≪ 1, this equals the inviscid case in Eq. ( 38).In the loosely-coupled limit with τ s ≫ 1, we get 4. OVERSTABILITY We now direct our attention to overstable modes with nonzero ω.Figure.5 shows growth rates of such oscillatory modes as obtained from the full dispersion relation, and demonstrates that large stopping times are required to achieve growing modes.We caution that in this regime, the fluid approximation for the dust grains that underlies our model starts to break down (also see Sects.2.3, 5.5 as well as Appendix A).We proceed with the analysis for completeness.
For oscillating solutions with ω ̸ = 0, Eq. ( 32) implies which is positive definite for ξ ≪ 1 as seen in Fig. 5.When inserted into (33), this results in a cubic for the growth rate of the wave, i.e.
For small growth rates γ ≪ 1, A 1 , A 2 , the root is For ξ ≪ 1, we have Assuming large stopping times τ s ≫ 1, this becomes which is plotted in Fig. 5 in comparison to the root of the full cubic.For overstability, i.e. growing oscillations, the growth rate must be positive, i.e. γ > 0, or equivalently, which equals 1/9 for Sc = 1.Since overstability is characterized by a restoring force which lacks in case of pure instability, this requirement on β visc is opposite in direction compared to the criteria for instability discussed in the previous sections (also compare to the classical viscous overstability e.g., in Latter & Ogilvie (2006a)).
We also note that the 1/Sc term in Eq. ( 53) originates from advection of angular momentum carried by background shear due to the diffusive flux, i.e. the second term on the right hand side of Eq. ( 13).For Sc ≫ 1, i.e. negligible diffusion compared to viscosity, this term vanishes and we recover the classical criterion for the viscous overstability in planetary rings (Schmit & Tscharnuter 1995).This term does not allow overstability in the invicsid case.That is, by setting Sc = 0, the growth rate resulting from Eq. ( 51) becomes which tends to γ = −ξ/2 for large stopping times τ s ≫ 1 (compare to Eq. ( 52)).
In Fig. 6, we depict growth rates of the diffusive overstability obtained from Eq. ( 51) as a function of the stopping time.The smallest scales have fastest growth rates, and, have the least stringent restrictions on stopping time.For the depicted set of parameters, τ s ≳ 10 would lead to overstability on the smallest scales.While this requirement can be further relaxed for larger Sc or more positive β visc , for τ s ≲ 1, all oscillating modes are damped.
For this reason, the overstable modes discussed in this section have little physical relevance in the context of particles that generate streaming instability turbulence, where the largest available dust grains are limited at typically τ s ≲ 1 (e.g., Birnstiel et al. 2016).We elaborate on this, as well as explore alternative situations where the diffusive overstability may find applicability in Sect.5.5.

DISCUSSION
In the previous sections, we showed that a dust fluid in protoplanetary disks, which is governed by Eqs. ( 1) -(3) can be unstable (Sect.3) and over-stable (Sect.4) when an otherwise constant background state is linearly perturbed, depending on the stopping time of particles and the steepness of the diffusion and viscosity slopes with respect to the dust surface mass density.Here, we discuss and contextualize our findings.

Physical picture
First, we reiterate the physical mechanisms that drive the newly found instabilities.Figures 7 and 8 show the  28).We set βvisc = β diff = 1/2.The dashed lines correspond to Eq. ( 52).The lines end once modes become non-oscillatory (and thus damped), which occurs as ξ approaches the hatched region where our model is not applicable.51), which assumes ξ ≪ 1.For the smallest scales our model applies to, and, given the chosen set of parameters, stopping times of τs ≳ 8 would allow for overstability.
growth rates of the full system including the inviscid (Sc = 0) and viscous (Sc = 1).Guided by the depicted parameter space we can broadly assign the unstable regions into three categories 1. Diffusive instability driven by diffusion-dependent pressure: For small stopping times, the utilized dust pressure prescription (with velocity dispersion c 2 d ∝ D ∝ Σ β diff ) can lead to linear instability if the diffusion slope β diff is sufficiently negative (see Eqs. ( 36) and ( 40)), and drives particles towards density maxima.The instability is only damped on small scales if viscosity is nonzero, and has fastest growth rates on the smallest unstable scales.This instability requires τ s ≲ 1, β diff < −2, and if Sc ≳ 1 also ξ ≲ 1, which for a fiducial δ ∼ 10 −5 corresponds to wavelengths of λ ≳ 0.01H.Associated growth rates are shown in the top two panels of Fig. 7 as well the left side of the top panel of Fig. 8 for small stopping times.

Diffusive instability driven by viscosity slope:
For nonzero viscosity, the viscosity term dominates over the pressure term for sufficiently large stopping times, leading to a version of the viscous instability (Ward 1981;Lin & Bodenheimer 1981).This instability requires τ s ≳ 1, Sc ≳ 1, ξ ≲ 1, β visc ≲ −1, and is seen in the top two panels of Fig. 8 for large stopping times.

Diffusive overstability driven by viscosity slope:
Analogously to the classical viscous overstability (Schmit & Tscharnuter 1995), the repellent term that amplifies oscillations is provided by a viscosity slope.This overstability requires τ s ≳ 10, Sc ≳ 1, ξ ≲ 1, and β visc ≳ +1/9 (depending on the value of Sc).The associated growth rates are shown in the bottom panel of Fig. 8.The diffusive instability driven by diffusion-dependent pressure is active for all modes ξ if β diff < −2 (see Eq. ( 36)), which is the case in the upper two panels for small stopping times.For large stopping times, growth rates become uninterestingly small.Sc = 0, does not allow growth of oscillatory modes.
The linear theory describing the diffusive instabilities and overstabilities presented in this work is largely similar to that of the classical viscous instability (Ward 1981;Lin & Bodenheimer 1981) and axisymmetric viscous overstability (Schmit & Tscharnuter 1995) in planetary rings, at least under the neglect of self-gravity and thermal effects, and given the appropriate limits in our model, i.e. τ s ≫ 1, Sc ≫ 1.It should be noted though that the physical origin of viscosity and pressure in planetary rings lies in mutual particle collisions, in contrast  3).In the panel second from the top, Eq. ( 40) remains satisfied, while Eq. ( 36) is not, so instability is driven by the viscosity slope alone and thus restricted to larger τs.In the third panel from the top, neither instability is active.When Eq. ( 53) is satisfied (for Sc = 1 when βvisc exceeds +1/9, the (modified) viscous overstability driven by the viscosity slope can be active for large stopping times.
to the situation depicted in this work, where it results from self-generated turbulence, that is external insofar as the model is concerned.Also, a 'thin disk' version of the viscous overstability (i.e. on large radial length scales ≫ H g ), generated by a constant kinematic shear viscosity, can in principle also exist in gaseous protoplanetary disks (Latter & Ogilvie 2006a).
In addition, we mention the dust-driven viscous ringinstability pioneered by Dullemond & Penzlin (2018) as well as the related instability in Johansen et al. (2011), both of which, although operating on larger scales, are similar in spirit to this paper's diffusive instabilities.They consider a disk where the gas viscosity is set by turbulence generated from magnetorotational instability (Balbus & Hawley 1991), which weakens with increasing dust density through the ionization fraction.It is the viscous gas disk itself that can now be unstable to linear perturbations: an increase in gas density attracts dust grains as they tend to drift towards gas pressure maxima (e.g.Sano et al. 2000), turbulence weakens and the gas viscosity drops, thereby increasing the gas density further, which attracts more dust, and so on.In our model, we only consider the dust fluid but the decrease in its viscosity, diffusivity, and dust pressure as the dust surface density increases is likewise physically motivated by dust feedback lowering the local diffusive properties of the turbulence generated on small-scales.
Indeed, while we motivate our model with streaming instability turbulence and associated measurements of diffusion slopes (Schreiber & Klahr 2018;Gerbig & Li 2023), it may as well be applicable to other sources of turbulence.For example, the azimuthal streaming instability discovered by Hsu & Lin (2022) has likewise been observed to evolve into filaments once linear modes are saturated.Moreover, if one considers pure gas instabilities as sources of turbulence, the vertical shear instability (Urpin & Brandenburg 1998;Nelson et al. 2013;Lin & Youdin 2015;Barker & Latter 2015;Pfeil & Klahr 2021) or the convective overstability (Klahr & Hubbard 2014;Lyra 2014;Latter 2016) may generate environments suitable for secondary diffusive instabilities.The dependence of diffusivity and viscosity on disk parameters would need to be clarified with detailed simulations, but to zeroth order, one expects a drop in turbulence with increasing dust-loading, similar to that for the magnetorotational instability discussed above, because dust feedback tends to stabilize the vertical shear instability (Lin 2019;Lehmann & Lin 2022), and linear convective overstability (Lehmann & Lin 2023).

Filament formation through diffusive instability
Based on the analytic findings in Sects.3 and 4, we hypothesize, that the diffusive instability, driven by a sufficiently negative density dependence of the dust pressure, physically motivated by dust loading reducing diffusivity, may act to amplify density perturbations.These perturbations then could saturate to become marginally Diffusive instability driven by diffusiondependent pressure with same parameters as in Fig. 9. Shown is a map in x − y space of the particle surface density after a time of 3Ω −1 , with super-imposed re-scaled velocity vectors (v ′ x , v ′ y ).The fluid motion is towards density maxima, thus amplifying the perturbation and producing filaments.We show the y-coordinate for better visualization, even though the model itself is axisymmetric.
We visualize the growth of the linear perturbations into such filament-like overdensities in a space -time diagram of the perturbed density Σ ′ following Eq.( 14).Fig. 9 shows a perturbation subject to diffusive instability.The chosen set of parameters produces overdensities with a radial spacing of about 0.02H, which is consistent with the first emergent filaments found in streaming instability simulations (e.g., Li & Youdin 2021).
Indeed, the diffusive instability has the fastest growth rates on the smallest scale not limited by viscosity ξ max , which is of order ξ max ∼ 1.The corresponding fastest growing mode is thus expected to be around  55).On the other hand, the commonly used fiducial value for the filament feeding zone of 0.2H (Yang & Johansen 2014), cannot be directly compared to the scales of interest, as it already involves post-formation nonlinear dynamics, such as mergers and breakups.Fig. 11 visualizes the streaming motion that arises from the diffusive instability.Dust is moving towards density maxima which act as particle traps, in the process amplifying the perturbation.
For comparison, Fig. 10 shows traveling waves driven by the diffusive overstability.The shearing sheet is symmetric in x, so waves traveling towards positive x correspond to the complex conjugate of the depicted solution and are equally valid.Note, that the radial bulk velocity is entirely caused by to the overstability.If one were to include the background pressure gradient, the flow would pick up a to-first-order, constant background drift, which does not affect stability since we can shift into the co-moving frame.

Diffusive instability in the context of past numerical simulations
The aim of this study is to present a simple linear model for the emergence of the first filaments out of streaming instability turbulence, as seen in a number of simulations including Schreiber & Klahr (2018); Gerbig et al. (2020); Li & Youdin (2021); Gerbig & Li (2023).Within our model, the required linear growth rates and corresponding required stopping times that are expected to result in filament formation depend on the values of the diffusion and viscosity slopes β diff and β visc , respectively.Therefore, we aim to contextualize our findings with existing numerical studies and measurements of the aforementioned slopes.As such, our presented theory can most readily be compared to the numerical results presented in Schreiber & Klahr (2018), who performed two-dimensional, non-axisymmetric, non-selfgravitating shearing sheet simulations of the streaming instability.By conducting a number of simulations at different dustto-gas ratios and measuring the average particle diffusivity in each simulation, they were able to obtain the slope of diffusion with respect to the dust surface mass density, resulting in values −2 ≲ β diff ≲ −1, where the exact value depended on box size and particle stopping time used in their simulations.While some of the simulations of Schreiber & Klahr (2018) revealed the emergence of filaments with particle concentrations significantly enhanced relative to the ambient background, most did not.In the context of the linear diffusive instability discovered in our work, which requires β diff < −2 (see Fig. 3) for instability, the simulations in Schreiber & Klahr (2018) are hence expected to be marginally stable.More recently, Gerbig & Li (2023) measured the diffusion slope within a single vertically stratified shearing box simulation of the streaming instability.They also found β diff ∼ −2.As this value was obtained after the formation of filaments in their simulations, this is also consistent with their simulations being marginally unstable with respect to diffusive instability, provided the saturation of filament growth results in such a marginally-stable state.

Connection to planetesimal formation
Within the streaming instability paradigm of planetesimal formation, filaments are often thought to be a necessary precursor for planetesimal formation, as they provide the necessary conditions for subsequent selfgravitational fragmentation.It is therefore of interest to compare the parameter space for active diffusive instability within our model to that determined in nu- merical simulations of planetesimal formation, specifically clumping thresholds in streaming instability simulations.For this purpose, we can relate metallicity Z to the dust-to-gas ratio ϵ = ρ p /ρ g and diffusivity via (e.g., Lin 2021) where H d is the vertical scale height of dust.
Assuming τ s ≫ δ, we have δ ∼ τ s (Z/ϵ) 2 .Since our model, and therefore also the diffusive instability mechanism only implicitly depends on δ via ξ, a change in diffusivity due to metallicity increase (decrease) can be compensated by a decrease (increase) in wavenumber K towards larger (smaller) radial length scales.For example, the typically fast growing mode of ξ ∼ 1, would correspond to wave numbers of K ∼ ϵ/(Z √ τ s ).
We proceed by imposing a minimum scale of λ crit ∼ 0.01H, below which the model does not apply (cf.Sect.2.2), and use this to restrict the maximum allowed ξ, given some value of Z. Specifically, we have We then assume that the fastest growing mode admitted by our model is given by ξ(Z, τ s ) = min(ξ fgm , ξ max ), where ξ fgm is the mathematically fastest growing mode, which is obtained numerically from the full dispersion relation in Eq. ( 28).ξ fgm is typically of order unity (see Fig. 2), but may be larger, either in the absence of viscosity or for very small stopping times (see Figs. 2 and 4 respectively); or smaller, if the viscosity slope is only marginally more negative than the critically required slope (see top panel in Fig. 8).
Figures 12 and 13 show the growth rates associated with this preferred mode assuming ϵ = 1, and λ crit = 0.01H.We also take β diff = β visc = −2.2 in both Figures.Fig. 12 shows the inviscid case with Sc = 0, where only the diffusive instability driven by diffusiondependent pressure is possible given the set of parameters, and Fig. 13 shows growth rates associated with the same set of parameters, except that now Sc = 1.Over-plotted are the clumping thresholds obtained by Yang et al. (2017) and Li & Youdin (2021) (also see Carrera et al. (2015) for another version), below which, streaming instability clumping is unlikely to occur.
For either model, the entire parameter space, probed by the studies Yang et al. (2017) and Li & Youdin (2021) would be subject to diffusive instabilities (although with differences in expected growth rates), rendering the models consistent with the hypothesis that planetesimal formation is triggered by the emergence of filaments induced by diffusive instability.
We note, that the models displayed in Figs. 12 and 13 contain a number of parameters (λ crit , ϵ, β diff , β visc , Sc), which we chose relatively freely and assumed to be independent of τ s in order to get a prototype idea if there may be a connection to planetesimal formation.More detailed calculations, in concert with additional numerical constraints on diffusion and viscosity slopes, are required to further assess diffusion instability's role in planetesimal formation.

Instability and overstability at large stopping times?
The local axisymmetric viscous overstability of a thin astrophysical disk has extensively been studied in the context of Saturn's dense rings, employing hydrodynamic models (Schmit & Tscharnuter 1995, 1999;Spahn et al. 2000;Schmidt et al. 2001;Schmidt & Salo 2003;Latter & Ogilvie 2009, 2010;Lehmann et al. 2017Lehmann et al. , 2019)), kinetic models (Latter & Ogilvie 2006b, 2008), and Nbody simulations (Salo et al. 2001;Rein & Latter 2013;Ballouz et al. 2017;Lehmann et al. 2017;Mondino-Llermanos & Salo 2023).Based on results from hydrodynamic and N-body simulations, overstability in Saturn's rings typically saturates in the form of nonlinear traveling wave trains that could in principle carry appreciable amounts of angular momentum.Wave trains are often interspersed by defect structures, which may act as sources or sinks of the former.Indeed, viscous over-stability is the most promising mechanism to explain the occurrence of periodic fine structures on a ∼ 100m scale in parts of Saturn's A and B rings, that has directly been observed (Thomson et al. 2007;Colwell et al. 2007;Hedman et al. 2014).It is thus of interest to explore the conditions under which we expect the diffusive overstability to operate in over-dense particle layers in protoplanetary disks.
Since this requires τ s ≳ 1, we preface further discussion, by reiterating that the hydrodynamical model is strictly not applicable for large stopping times, and instead a kinetic model should be used (see Sect. 2.3).
As outlined in Sect.5.1, in addition to the diffusive instability driven by the pressure term, there are two types of instabilities that can arise for large stopping times, i.e. instability and overstability driven by the viscosity slope with their respective specific requirements on on β visc .
The questions are (a) to what extent the invisicd, hydrodynamic model can still appropriately describe the system for τ s ≳ 1, and whether or not particles with sufficiently large stopping times (b) can exist and (c) are in their turbulent behavior still well characterized by the diffusive flux model that shapes the continuity equation, the pressure model and the angular momentum conserving terms in the momentum equations.
Indeed, in the classical picture of particle growth in protoplanetary disks, stopping times are limited at around τ s ∼ 1 (e.g., Birnstiel et al. 2012Birnstiel et al. , 2016)), and as such, individual particles that qualify for diffusive overstability would already be considered planetesimal-sized objects.Another possible pathway of getting objects with large stopping times was suggested by Johansen & Youdin (2007).In an effort to explain unexpected drift rates of particle clumps found in their streaming instability simulations, they hypothesize that a clump may collectively have an increased stopping time relative to the individual grains, due to shielding each other from the gas stream and providing an order-ofmagnitude scaling with R clump being the clump's radius, and ∆v its velocity relative to the gas.While such clumps are far too few in number density to be appropriately modeled by a fluid approach, we adapt this hypothesized collective shielding effect to our situation.Specifically, applying Johansen & Youdin's argument to the quasi-steady, turbulent dust layer that we envision as the equilibrium, we set R clump = H d , the dust scale height.Then, using H d = ZH/ϵ from Eq. ( 56), we find τ eff s ∼ Z/Π, where Π ≡ η/(H/r) is the reduced radial pressure gradient parameter, again taking ∆v ∼ ηrΩ.
The diffusive overstability's requirements of τ eff s ≳ 1 translates to Z ≳ Π, and the same applies for the diffusive instability driven by viscosity slope.Coincidentally, Bai & Stone (2010) find that clumping via the streaming instability becomes easier with smaller Π.Similarly, at fixed τ s (the stopping time for individual grains), Sekiya & Onishi (2018) find filament formation in their streaming instability simulations if √ 2πZ/Π ≳ 1.We can interpret these results within our model as filaments only form if the effective stopping time, realized through Z/Π, exceeds unity in order to trigger the diffusive overstability (or the diffusive instability driven by viscosity slope).
The streaming instability formally still operates for τ s ≳ 1 as growth rates only decrease slowly with stopping time (Pan 2020).However, due to the lack of numerical constraints on diffusivity in this regime, it is unclear to what extent instabilities and overstabilities would develop on scales within our model's validity; even if diffusion and viscosity slopes were to remain unchanged.Since, in order to achieve instability on relevant scales, diffusivities cannot be arbitrarily small, we consider the possibility of diffusion to not be selfgenerated.Instead, diffusion may stem directly from a turbulent gas, the diffusivity and viscosity of which we denote as δ g and α g respectively.
In the preceding sections, we were exclusively concerned with particle diffusivity δ, which we treated as wholly independent from τ s .While this is mathematically self-consistent, it does not necessarily reflect physical conditions.Indeed, an increased particle response time to turbulence diminishes the diffusion experienced by the particles as (Youdin 2011) which reduces to δ ∼ δ g for small stopping times, but becomes δ ∼ δ g /τ 2 s large stopping times.That is, large grains feel a much reduced turbulence than that in the gas.
Consider, for example, gas with a fiducial diffusivity of δ g ∼ α g ∼ 10 −3 , which is a typical value in numerical simulations of vertical shear instability or magnetorotational instability (e.g., Flock et al. 2017).For τ s ∼ 10, this would lead to δ ∼ 10 −5 , that is comparable to the diffusivities generated by streaming instability with smaller particles.The fastest growing scale per Eq. ( 55) remains unchanged at ∼ 0.02H.Indeed, since the diffusive instabilities discussed in this paper depend on ξ = δK 2 only, any decrease in diffusivity due to particle response to turbulence would only shift the fastest growing mode down to smaller scales, but not prohibit the mechanism itself from operating.
We thus argue that the large stopping time instability and overstabilities may find applicability, if a big enough collection of large stopping time particles are available in protoplanetary disks, or if the dust layer has a collective stopping time that exceeds unity.Note, that this estimation ignored the effect of the particle layer on gas turbulence, i.e. δ g itself.For example, turbulence driven by the vertical shear instability is damped by dust feedback, even when the dust-to-gas ratio is less than unity (e.g., Lin 2019).

Caveats and additional considerations
Our vertically averaged model neglects vertical motions.Therefore, inertial waves are discarded, such that the classical, axisymmetric streaming instability is not captured by our model (Squire & Hopkins 2018), even if we were to include the gas equations with a radial background pressure gradient.While we attribute the underlying turbulence, characterized by diffusion and viscosity of the dust fluid, to the streaming instability (or an equivalent mechanism), filament formation in our model is not a direct result of this underlying, small-scale instability.Instead, it originates from one or more intrinsic large-scale instabilities of the mid-plane dust layer, described by our model.In order to neglect gas perturbations, our model is time-averaged over one turbulent correlation time.As a result, linear modes with frequency smaller than one inverse correlation time should be considered with caution.Future work should include the dynamical equations describing the gas, as well as the vertical dimension, to investigate filament formation via streaming-type instabilities with variable viscosity and diffusivity in a more rigorous manner.
In our model, the diffusion and viscosity slope is assumed to depend on dust surface density only, primarily because of the available numerical constraints from Schreiber & Klahr (2018); Gerbig & Li (2023) and due to the analogy to isothermal hydrodynamic models for viscous instabilities in planetary rings.However, if diffusion and viscosity indeed arise from small-scale streaming instability, other dependencies are conceivable, namely relative dust-gas streaming velocity which powers the streaming instability in the first place.
We also neglect particle self-gravity as filaments already form before self-gravity is turned on in simulations (e.g., Schreiber & Klahr 2018;Gerbig et al. 2020;Gerbig & Li 2023).Nonetheless, self-gravity may have importance, modifying instability criteria and growth rates.
Lastly, we also ignored the possibility of a polydisperse dust fluid with a distribution of stopping times, and instead asserted a single stopping time, which is reasonable for a top-heavy size distribution (e.g.Birnstiel et al. 2012).Still, the damping effect of particle size distributions on streaming instability (Krapp et al. 2019;Paardekooper et al. 2020;Zhu & Yang 2021;Yang & Zhu 2021) underscores that this is a relevant point to keep in mind in future investigations of diffusive instabilities.

SUMMARY
In this work we present a novel, vertically averaged axisymmetric hydrodynamic model for a dense particle layer embedded in a gaseous protoplanetary disk.The dust being dominant, we model the effect of gas as a perturbation on the dust dynamics by evoking drag forces, mass diffusion, viscosity, and pressure of the dust.The pressure is assumed to depend on the dust diffusivity following a sedimentation-diffusion ansatz, and diffusivity and viscosity are allowed to depend on the dust-surface mass density.We find that our model supports a variety of linear diffusion and oscillatory instabilities.
The diffusion instabilities arise if the dust particles' diffusion and/or viscosity decrease sufficiently fast with increasing particle surface mass density, which is motivated by the results of past simulations.Specifically, for well-coupled particles with τ s ≪ 1, the diffusiondependent pressure can destabilize the particle flow if the mass diffusion slope with respect to dust surface mass density is sufficiently negative.On the other hand, for decoupled particles τ s ≳ 1, instability is driven by the viscosity slope, similar to the viscous instability in planetary rings.
The main application of our model is a dense midplane particle layer subject to turbulence generated by the streaming instability on small-scales.Indeed, the diffusivities associated with streaming instability turbulence measured in past simulations are found to be sufficient for the diffusive instabilities in our model to possess appreciable growth rates on the order of the dynamical time scale, and radial length scales that are characteristic of over-dense particle filaments seen in numerical simulations of the streaming instability.Based on these findings we argue that diffusive instabilities as captured by our model may play a role in filament formation within dusty protoplanetary disks, which is a key step within the streaming instability paradigm of planetesimal formation.
In addition, our model can also give rise to growing oscillatory modes in the presence of a sufficiently flat, or positive, viscosity slope that provides the necessary repellent acceleration.This overstability behaves similar to the axisymmetric viscous overstability in planetary rings.Whether or not it has applicability in protoplanetary disks is unclear, as it relies strongly on the particles possessing large stopping times τ s > 1, as well as their interaction with the gas turbulence.
More detailed analytical investigations including vertical motions and an explicit inclusion of gas within a two-fluid formalism, accompanied by additional numerical constraints on diffusivity, viscosity, as well as their slopes, will be crucial to further pinpoint the relevance of diffusive instabilities in dusty protoplanetary disk, filament formation therein, and planetesimal formation.
requires consideration of the gas also, leading to a dust layer thickness of H d = δ/(δ + τ s )H (see e.g., Lin 2021).We define the effective particle velocity dispersion via As also noted by Klahr & Schreiber (2021), the sedimentation-diffusion velocity dispersion c d and thus our closure relation, is in general different to the root mean square (RMS) velocity of the particles.The Hinze-Tchen formalism for turbulent transport neglecting orbital oscillations or other external forces (Tchen 1947;Hinze 1959), gives the RMS velocity as (also see e.g., Fan & Zhu 1998;Youdin & Lithwick 2007;Binkert 2023) with gas velocity dispersion u rms , and correlation time of the turbulence t c , which connect to the gas diffusivity D g via D g = t c u 2 rms .Only if D g ∼ D (see Eq. ( 59)) and for t s ≫ t c , the RMS velocity equals the velocity dispersion following from the sedimentation diffusion ansatz.While the former is fairly well grounded in numerical simulations (e.g., Schreiber & Klahr 2018), the latter is somewhat more ambigous.For Kolmogorov turbulence, the correlation time equals the turnover time of the largest eddies, which in protoplanetary disks equals Ω −1 (Youdin & Lithwick 2007).On the other hand, Schreiber & Klahr (2018) find for simulations values with active streaming instability values of t c ∼ 0.1Ω −1 .Indeed, ignoring orbital oscillations is only a good model if t s , t c ≪ 1.For larger particles, epicylic oscillations are important, particles can decouple from the turbulence, and the velocity dispersion and thus diffusion needs to be modified (see Youdin & Lithwick 2007;Youdin 2011, and also Eq. ( 59)).We neglect this effect in this work.While our pressure term does vanishes for large stopping times like the prescriptions used by e.g., Youdin (2011); Umurhan et al. (2020), the diffusivity itself does not, and is instead treated as an independent parameter.This is important to be kept in mind when evaluating our model in particular in the large stopping time limit.Klahr & Schreiber (2021) call c d the pseudo sound speed, which is appropriate as long as c d is constant.In our work, we allow the diffusivity to depend on density, and thus c 2 d ∝ D ∝ Σ β diff .As a result, the dust pressure takes the form of a polytropic equation of state, i.e.P d ∝ Σ 1+β diff .One can can now formally define a dust sound speed a d via If β D < −1, this effective squared sound speed is negative and as a result also the associated pressure perturbations.Indeed, such a negative pressure perturbation is a necessary (but not sufficient) requirement for the diffusion-dependent pressure driven diffusive instability that may operate for τ s ≲ 1, as discussed in this paper.

Figure 3 .Figure 4 .
Figure 3. Like Fig. 1.Growth rates for the viscous-diffusive instability for ξ = 0.1, assuming equal mass and momentum diffusion with Sc = 1 and β diff = βvisc.The black curve corresponds to Eq. (40), below which the instability can operate, and above which perturbations are exponentially damped.Note that while parts of the depicted parameter space allow for viscous overstability, this figure is only showing the purely real solution.See Sect. 4.

Figure 5 .
Figure5.Growth rates for diffusive overstability driven by viscosity slope (Sc = 1, left panel) compared to the damping rates in the absence of a viscosity (Sc = 0, right panel) for various stopping times obtained from the full dispersion relation in Eq. (28).We set βvisc = β diff = 1/2.The dashed lines correspond to Eq. (52).The lines end once modes become non-oscillatory (and thus damped), which occurs as ξ approaches the hatched region where our model is not applicable.

Figure 6 .
Figure6.Growth rates for diffusive overstability driven by viscosity slope (Sc = 1) vs stopping time for various values of ξ; setting βvisc = β diff = 1/2 following Eq.(51), which assumes ξ ≪ 1.For the smallest scales our model applies to, and, given the chosen set of parameters, stopping times of τs ≳ 8 would allow for overstability.

Figure 7 .
Figure7.Growth rates for inviscid (Sc = 0) diffusive instability.The diffusion slope β diff increases from top to bottom.The diffusive instability driven by diffusion-dependent pressure is active for all modes ξ if β diff < −2 (see Eq. (36)), which is the case in the upper two panels for small stopping times.For large stopping times, growth rates become uninterestingly small.Sc = 0, does not allow growth of oscillatory modes.

Figure 8 .
Figure8.Like Fig.7but for the viscous-diffusive instability and (modified) viscous overstability with Sc = 1.The top panel shows combination of diffusive instability driven by diffusion-dependent pressure and viscous instability driven by the negative viscosity slope for small and large stopping times respectively (compare to Fig.3).In the panel second from the top, Eq. (40) remains satisfied, while Eq.(36) is not, so instability is driven by the viscosity slope alone and thus restricted to larger τs.In the third panel from the top, neither instability is active.When Eq. (53) is satisfied (for Sc = 1 when βvisc exceeds +1/9, the (modified) viscous overstability driven by the viscosity slope can be active for large stopping times.

Figure 9 .Figure 10 .
Figure 9. Filament formation shown in a space time diagram of a linear perturbation subject to diffusive instability driven by pressure.For this plot, we chose Sc = 0, τs = 0.1, β diff = −3.We also chose a value of δ = 10 −4 , and a fast growing mode of ξ = 0.1, this corresponds to a physical wave number of k = 100/H.The eigenvector is scaled with Σ = 0.01Σ0.

Figure 13 .
Figure 13.Like Fig. 12 but now for Sc = 1, which allows for instability driven by the viscosity slope at large stopping times in addition to the instability driven by the diffusiondependent pressure term.Other parameter choices are identical to those in Fig. 12.