Stratified and vertically-shearing streaming instabilities in protoplanetary disks
DDraft version November 26, 2020
Typeset using L A TEX twocolumn style in AASTeX63
Stratified and vertically-shearing streaming instabilities in protoplanetary disks
Min-Kai Lin Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617, Taiwan
ABSTRACTUnder the right conditions, the streaming instability between imperfectly coupled dust and gas isa powerful mechanism for planetesimal formation as it can concentrate dust grains to the point ofgravitational collapse. In its simplest form, the streaming instability can be captured by analyzingthe linear stability of unstratified disk models, which represent the midplane of protoplanetary disks.We extend such studies by carrying out vertically-global linear stability analyses of dust layers inprotoplanetary disks. We find the dominant form of instability in stratified dust layers is one drivenby the vertical gradient in the rotation velocity of the dust-gas mixture, but also requires partial dust-gas coupling. These vertically-shearing streaming instabilities grow on orbital timescales and occuron radial length scales ∼ − H g , where H g is the local pressure scale height. The classic streaminginstability, associated with the relative radial drift between dust and gas, occur on radial lengthscales ∼ − H g , but have much smaller growth rates than vertically-shearing streaming instabilities.Including gas viscosity is strongly stabilizing and leads to vertically-elongated disturbances. We brieflydiscuss the potential effects of vertically-shearing streaming instabilities on planetesimal formation. INTRODUCTIONThe formation of planetesimals from mm–cm sizeddust grains or pebbles in protoplanetary disks (PPDs)is a key stage in planet formation (Birnstiel et al. 2016).Neither pair-wise collisions nor gravitational forces leadto effective growth on pebble scales (Chiang & Youdin2010; Blum 2018). However, if a swarm of solids can bemade sufficiently dense relative to the ambient gas, thenit can undergo direct self-gravitational collapse into kmor larger-sized planetesimals (Goldreich & Ward 1973).The critical dust-to-gas ratio for collapse is (cid:29) (Shi &Chiang 2013). This should be compared to the typ-ical value of ∼ expected uniformly throughout anewly-born PPD (Testi et al. 2014). Thus, an efficientmechanism is needed to first concentrate dust grains.These include dust settling, particle, trapping by pres-sure bumps, and dust-gas instabilities (Johansen et al.2014).The streaming instability (SI, Youdin & Goodman2005; Youdin & Lithwick 2007) is one such candidate.The SI is linear instability in rotating flows of partially-coupled dust and gas that mutually interact throughfrictional drag – conditions natural in PPDs — whichcan amplify dust-to-gas ratios by orders of magnitudeand trigger gravitational collapse (Johansen et al. 2009), [email protected] although the SI itself does not require self-gravity. In theoften considered, idealized case of a laminar disk with amonodisperse dust population, the SI is a robust processand has thus received considerable attention as the defacto mechanism for planetesimal formation.The SI has undergone intense studies through nu-merical simulations (Johansen & Youdin 2007; Bai &Stone 2010a; Yang & Johansen 2014; Yang et al. 2017).Modern simulations have generalized the SI to considermagnetic fields (Balsara et al. 2009; Tilley et al. 2010;Yang et al. 2018); various geometries (Kowalik et al.2013; Schreiber & Klahr 2018); multiple grain sizes (Bai& Stone 2010b,c; Schaffer et al. 2018; Benítez-Llambayet al. 2019; Krapp et al. 2019; Zhu & Yang 2020); turbu-lence (Schäfer et al. 2020; Gole et al. 2020); self-gravity(Simon et al. 2016; Schäfer et al. 2017; Li et al. 2019),pressure bumps (Carrera et al. 2020), etc. These effortsare necessary to understand planetesimal formation inrealistic PPDs. Indeed, sophisticated simulations showthat planetesimals formed through the SI have proper-ties consistent with that in the solar system (Nesvornýet al. 2019).On the other hand, a physical understanding of the SIthrough analytical studies has progressed more slowly.Jacquet et al. (2011) showed that the SI is driven bya process of runaway dust-trapping by pressure bumps,while Lin & Youdin (2017) gave a thermodynamic in-terpretation of the SI in which partial dust-gas couplingleads to ‘PdV’ work that acts to amplify oscillations. a r X i v : . [ a s t r o - ph . E P ] N ov M.-K. Lin
Recently, Squire & Hopkins (2020) presented detailedmodels of Youdin & Goodman (2005)’s classic SI, inwhich the mutual interaction between epicyclic motionsand the relative radial dust-gas drift in a PPD leadsto growing perturbations. In fact, for small dust-to-gas ratios, the classic SI belong to a broader class of‘resonant drag instabilities’ (Squire & Hopkins 2018a,b;Zhuravlev 2019) generic to dusty-gas in which the rela-tive dust-gas motions resonate with a wave mode in thegas. The classic SI is thus expected in PPDs since dustand gas naturally exhibit a relative radial drift, as thegas rotation is partially supported by a (negative) radialpressure gradient (Whipple 1972; Weidenschilling 1977).The classic SI can be captured in relatively simple diskmodels, such as that employed by Youdin & Goodman,who considered a small region near the disk midplane.In this limit, the vertical gravity from the central starcan be neglected, which produces a uniform vertical diskstructure. These ‘unstratified’ disk models allow signif-icant simplifications for analyses of the linear SI andgeneralizations thereof, which include the effect of pres-sure bumps, multiple species, and turbulence (Jaupart& Laibe 2020; Auffinger & Laibe 2018; Krapp et al. 2019;Paardekooper et al. 2020; Chen & Lin 2020; Umurhanet al. 2020; Pan 2020).However, PPDs do have a vertical structure. Dust-settling, which leads to the formation of a dense particlelayer about the disk midplane (Dubrulle et al. 1995), isoften considered as a prerequisite to trigger the SI, asit requires dust-to-gas ratios of order unity or above tooperate efficiently (Youdin & Goodman 2005). Dust-settling naturally produces a stratified dust layer. Infact, stratified simulations are now common, but thelinear SI has not been examined in stratified disks. Fill-ing this gap will be helpful in understanding how the SIoperates in realistic PPDs.The purpose of this work is to generalize previousstudies of the linear SI to account for the vertical struc-ture of dust layers in PPDs. To this end, we analyzethe stability of vertically-global, radially-local models ofdusty disks. We employ the standard, two-fluid descrip-tion of dusty-gas, as well as a simplified ‘one-fluid’ model(Lin & Youdin 2017) to verify some of our calculations.Our main result is that in stratified disks the verti-cal variation in the azimuthal velocity of the dusty-gas(or vertical shear) provides a significant source of ‘freeenergy’ that can be accessed by partial dust-gas cou-pling, which results in instability. We typically findthese vertically-shearing streaming instabilities (VSSIs)dominate over classic SIs, so the former should be thefirst to develop in settled dust layers. Our study con-firms and expands upon an earlier work from Ishitsu et al. (2009), who used direct simulations to study theeffect of vertical dust density gradients on the evolutionof dust layers.This paper is organized as follows. By way of introduc-ing notation, we first provide order-of-magnitude moti-vations to examine the SI in stratified disks in §2. Wethen describe our framework and disk models in §3. Re-sults from our linear stability analyses are presented in§5 for three examples: a high dust density layer, a lowdust density layer, and a viscous disk. We discuss impli-cations of our findings in §6 and conclude in §7. A listof frequently used symbols is summarized in AppendixA. PHYSICAL MOTIVATION2.1.
Geometric considerations
The classic SI of Youdin & Goodman (2005), discov-ered in unstratified disk models, is driven by the relativeradial drift between dust and gas, v drift = − (cid:15) ) ηr ΩSt + (1 + (cid:15) ) (1)(Nakagawa et al. 1986), where the Stokes number St isa dimensionless inverse measure of the strength of dust-gas coupling (and is proportional to the grain size), (cid:15) isthe dust-to-gas volume density ratio, r is the cylindricaldistance from the star, Ω = (cid:112) GM ∗ /r is the Keplerianfrequency ( M ∗ and G being the stellar mass and gravi-tational constant, respectively), and η is a dimensionlessmeasure of the global radial pressure gradient defined as η ≡ − r Ω ρ g ∂P∂r , (2)where ρ g is the gas density and P is the gas pressure.PPDs have η ∼ O ( h ) , where h g ≡ H g /r is the gas diskaspect-ratio and H g is the pressure scale-height, with h g ∼ . .We thus expect the SI to have characteristic length-scales ∼ ηr , which is significantly shorter than the globalradial lengthscales of typical PPD disk models ( ∼ r ).The SI can thus be considered as a radially-localizedphenomenon. However, the situation differs in the ver-tical direction.In realistic PPDs, dust settles into a layer of thickness H d , which can be related to the midplane dust-to-gasmass density ratio (cid:15) and the metallicity Z ≡ Σ d Σ g (cid:39) (cid:15) H d H g (3)(Johansen et al. 2014), where Σ d,g are the dust and gassurface densities, respectively. The SI operates on dy-namical timescales when (cid:15) (cid:38) (Youdin & Goodman tratified streaming instabilities Z (cid:39) . , dust layers should be thin, H d (cid:39) . H g . For a SI mode with vertical lengthscale ηr to exist, it should fit inside the dust layer, or χ ≡ ηrH d = ˆ η (cid:16) (cid:15) Z (cid:17) (cid:46) , (4)where ˆ η ≡ η/h g and PPDs typically have ˆ η (cid:39) . .For a settled dust layer with (cid:15) (cid:39) in a disk withstandard metallicity Z (cid:39) . we find χ (cid:39) , violatingthe above condition. Moreover, when gas viscosity andparticle diffusion are considered, SI modes have verticallengthscales comparable to H g (Umurhan et al. 2020);while further restricting vertical lengthscales to (cid:46) H d results in negligible growth rates (Chen & Lin 2020).These findings call for stratified analyses.2.2. Energetic considerations
Dusty PPDs possess vertical shear: the settled, dust-rich midplane rotates closer to the Keplerian speed r Ω than the gas-dominated, pressure-supported disk awayfrom the dust layer, which has a sub-Keplerian rotationof (1 − η ) r Ω because η > usually. Such a verticalvariation of the disk’s rotation speed is an importantsource of free energy (as borne out of our calculations).Now, the difference in the azimuthal velocity of thedusty midplane and the overlaying gas is ∆ v φ ∼ ηr Ω .For a dust layer thickness H d we can thus estimate thevertical shear rate within the layer as ηr Ω /H d = χ Ω .For small grains, this vertical shear rate is larger thanthe relative radial drift rate v drift /ηr by a factor of St − (cid:29) .We can also compare vertical shear to the verticalsettling of grains, as the latter can trigger a ‘dust set-tling instability’ (DSI, Squire & Hopkins 2018b; Krappet al. 2020). Taking the typical settling speed of a dustgrain to be | v dz | ∼ St H d Ω (Takeuchi & Lin 2002), wefind | ∆ v φ /v dz | ∼ ˆ η(cid:15) / (St Z ) . For PPDs, with (cid:15) ∼ , ˆ η ∼ Z ∼ O (10 − ) , this ratio is ∼ St − , i.e. large forsmall grains.The above estimates suggest that for small grains, ver-tical shear is a much larger energy source than the rela-tive radial drift or dust settling. Indeed, in the limit of St → the classic SI and DSI are suppressed and verticalshear drives non-axisymmetric Kelvin-Helmholtz insta-bilities (KHI, Chiang 2008; Lee et al. 2010).In this work we consider axisymmetric disturbancesso that KHIs are not applicable. Lin & Youdin (2017)showed that axisymmetric dusty disks are generally sta-ble in the limit St → , however large the vertical shearrate. This is due to stabilization by dust-induced, effec-tive buoyancy forces (see also Lin 2019). However, the above result ceases to be valid for St (cid:54) = 0 ,because in this case dust and gas are no longer perfectlycoupled and they can stream past one another, whichdiminishes the stabilizing effect of dusty buoyancy. Thisis similar to how rapid cooling can enable the ‘verticalshear instability’ (VSI) in gaseous PPDs (Nelson et al.2013; Lin & Youdin 2015) by eliminating gas buoyancy.(For the gaseous VSI, vertical shear originates from thedisk’s radial thermal structure.)We can therefore expect in a stratified, dusty disk thefree energy associated with vertical shear, here a resultof dust settling, to be accessible through an instabilitywith non-vanishing particle sizes. Indeed, we will findsuch instabilities are the dominant modes in stratifieddisks. We refer to them as vertically-shearing streaminginstabilities (VSSIs), since both vertical shear and dust-gas streaming motions are necessary. BASIC EQUATIONSWe consider a non-self-gravitating, unmagnetizedPPD comprised of gas and a single species of dust grainsaround a central star of mass M ∗ . The gas componenthas density, pressure, and velocity fields ( ρ g , P, V g ) . Weassume an isothermal gas so that P = c s ρ g with a con-stant sound-speed c s = H g Ω .We treat the dust population as a pressureless fluidwith density and velocity ( ρ d , V d ) . The dust and gasfluids interact via a drag force parameterized by a stop-ping time τ s (see below). The fluid approximation fordust is then valid for well-coupled, small grains such that τ s (cid:46) Ω − (Jacquet et al. 2011).3.1. Two-fluid, radially local, axisymmetric disk model
We study radially-localized disturbances in the afore-mentioned dusty disk using the shearing box framework(Goldreich & Lynden-Bell 1965). The shearing box iscentered about a fiducial point ( r , φ , in cylindricalco-ordinates on the star, which rotates at the referenceKeplerian frequency Ω( r ) ≡ Ω , i.e. φ = Ω t . Carte-sian co-ordinates ( x, y, z ) in the box correspond to theradial, azimuthal, and vertical directions in the globaldisk. The radial extent of the box is assumed to be muchsmaller than r , so that curvature terms from the cylin-drical geometry can be neglected. Keplerian rotation isthen approximated as the linear shear flow − (3 / x ˆ y .We define v d , g ≡ V d , g − ( r − x/ ˆ y as the local dustand gas velocities in the shearing box relative to thislinear shear flow. We assume axisymmetry throughout,so that ∂ y ≡ . M.-K. Lin
The governing equations for the dust component inthe shearing box are ∂ρ d ∂t + ∇ · ( ρ d v d ) = ∇ · (cid:20) Dρ g ∇ (cid:18) ρ d ρ g (cid:19)(cid:21) , (5) ∂ v d ∂t + v d · ∇ v d = 2Ω v d y ˆ x − Ω v d x ˆ y − Ω z ˆ z − τ s ( v d − v g ) , (6)where D is a constant diffusion coefficient defined be-low. The third term on the right-hand-side (RHS) ofEq. 6 corresponds to the vertical component of the stel-lar gravity in the thin-disk limit. The last term on theRHS corresponds to gas drag, the strength of which ischaracterized by the stopping time τ s .For the gas, we include the effect of a global radialpressure gradient in the shearing box by writing ∇ P → ∇ P − η r Ω ρ g ˆ x , (7)where η = η ( r = r , z = 0) and η is defined by Eq. 2.That is, the global radial pressure gradient is modeledas a constant forcing. We then re-interpret P as pres-sure fluctuations in the shearing box. The governingequations for the gas component are then ∂ρ g ∂t + ∇ · ( ρ g v g ) = 0 , (8) ∂ v g ∂t + v g · ∇ v g = 2Ω v g y ˆ x − Ω v g x ˆ y − ∇ Pρ g + 2 η Ω r ˆ x + 1 ρ g ∇ · T − Ω z ˆ z − (cid:15)τ s ( v g − v d ) . (9)The fifth term on the RHS of Eq. 9 represent viscousforces, where T = ρ g ν (cid:20) ∇ v g + ( ∇ v g ) † − I ∇ · v g (cid:21) (10)is the viscous stress tensor and ν is a kinematic viscosity,prescribed later. The final term on the RHS is the back-reaction of dust drag onto the gas.The basic equations 5–6 and 8–9 extend those used byChen & Lin (2020) with the addition of vertical grav-ity, which themselves are extensions of that in Youdin& Johansen (2007) with the addition of dust diffusionand gas viscosity. We solve Eqs. 5–6 and 8–9 in fullto obtain equilibrium states, then solve their linearizedversions to study the stability of said equilibria. Bothproblems are one-dimensional in z . For numerical solu-tions we consider the half-disk z ∈ [0 , z max ] by imposingsymmetry conditions at the midplane. Details are given in §3.5.1–3.5.2 and §4. Hereafter we drop the subscript‘0’ for clarity. Below, H g refers to the pressure scaleheight at the reference radius.3.2. Dust-gas drag
The stopping time τ s is the timescale for a dust par-ticle to reach velocity equilibrium with its surroundinggas. In this work we take τ s to be a constant parameterfor simplicity. It is convenient to define the dimension-less stopping time or Stokes number, St = τ s Ω . (11)We consider well-coupled, or small dust grains with St (cid:28) .Physically, St depends on the particle and gas prop-erties, such as grain size ( a ), internal density ( ρ • ), andgas density (Weidenschilling 1977). To put our calcula-tions in context, consider grains in the Epstein regimein a Minimum Mass Solar Nebula-like disk (MMSN) de-scribed in Chiang & Youdin (2010). We then find St = 0 . F − (cid:16) r au (cid:17) / (cid:18) ρ • gcm − (cid:19) (cid:16) a mm (cid:17) , (12)where F is a mass scale relative to the standard MMSN( F = 1 ). We are mostly interested in mm or sub-mm-sized grains with internal density gcm − at tens of auin MMSN-like disks.3.3. Dust diffusion
We include dust diffusion to allow a stratified equilib-rium state to be defined, in which dust settling is bal-anced by dust diffusion. Without dust diffusion, parti-cles would continuously settle and no steady state can beestablished for standard stability analyses. Dust diffu-sion is usually attributed to gas turbulence (e.g. Youdin& Lithwick 2007; Laibe et al. 2020), which is often mod-eled as a gas viscosity. We thus parameterize dust dif-fusion in terms of a gas viscosity, although for the mostpart we neglect viscosity in the gas equations.We model dust diffusion via the constant parameter δ such that D = δc s H g , (13)with δ given by δ = 1 + St + 4St (cid:0) (cid:1) α (14)(Youdin & Lithwick 2007; Youdin 2011), where α is aninput constant turbulent viscosity parameter defined be-low. In practice, δ (cid:39) α since we consider small grains. tratified streaming instabilities Turbulent viscosity
When considered, we model gas turbulence via a vis-cous stress tensor (see Eq. 10) and adopt the standardalpha prescription (Shakura & Sunyaev 1973) such thatthe kinematic viscosity is ν = αc s H g ρ g , eqm ρ g , (15)where ρ g , eqm ( z ) denotes the equilibrium gas density, de-rived below. We use this prescription so that the dy-namic viscosity ρ g ν is a fixed function of space, whichavoids viscous overstabilities that could complicate re-sults (Latter & Ogilvie 2006; Lin & Kratter 2016).3.5. Two-fluid equilibria
We seek steady, horizontally uniform equilibria with ∂ t = ∂ x = 0 . The gas continuity equation then imply v g z = 0 . All other are quantities are non-zero and z -dependent, e.g. ρ g = ρ g , eqm ( z ) . For clarity, hereafter wedrop the subscript ‘eqm’ on the equilibrium fields.3.5.1. Vertical equilibrium
The equilibrium dust mass and vertical momentumequations are: d ln (cid:15)dz = v d z D , (16) c s d ln ρ g dz = (cid:15) ΩSt v d z − Ω z, (17) v d z dv d z dz = − Ω z − ΩSt v d z , (18)where we recall (cid:15) = ρ d /ρ g . For constant St these maybe solved exactly to yield (cid:15) ( z ) = (cid:15) exp (cid:18) − β δ z H (cid:19) , (19) ρ g ( z ) = ρ g0 exp (cid:20) δ St ( (cid:15) − (cid:15) ) − z H (cid:21) , (20) v d z ( z ) = − βz Ω , (21)where (cid:15) , ρ g0 is the mid-plane dust-to-gas ratio and gasdensity, respectively; and β ≡ (cid:16) − (cid:112) − (cid:17) . (22)We thus require St ≤ . . Note that β (cid:39) St for thesmall particles considered in this work ( St (cid:28) ). We This can been seen by performing a Taylor expasion of the nu-merator in Eq. 22 also consider weak diffusion such that δ (cid:28) St . The dustlayer thickness can then be approximated by H d = (cid:114) δδ + St H g (23)(Dubrulle et al. 1995; Zhu et al. 2015). The localmetallicity is Z ≡ (cid:82) ∞−∞ ρ d dz (cid:82) ∞−∞ ρ g dz = Σ d Σ g . (24)In practice, we adjust (cid:15) until a specified value of Z isobtained. However, Eq. 3 and Eq. 23 also gives an ade-quate estimate, (cid:15) (cid:39) Z (cid:112) St /δ . The vertical structure isthen completely determined.3.5.2. Horizontal equilibrium
The equilibrium horizontal momentum equations are v d z dv d x dz = 2Ω v d y − ΩSt ( v d x − v g x ) , (25) v d z dv d y dz = − Ω2 v d x − ΩSt ( v d y − v g y ) , (26) v g y + 2 η Ω r − (cid:15) ΩSt ( v g x − v d x ) + νρ g ddz (cid:18) ρ g dv g x dz (cid:19) , (27) − Ω2 v g x − (cid:15) ΩSt ( v g y − v d y ) + νρ g ddz (cid:18) ρ g dv g y dz (cid:19) , (28)with (cid:15) ( z ) , ρ g ( z ) , and v d z ( z ) given by Eqs. 19 – 21. Thehorizontal velocity profiles must, in general, be solvednumerically subject to appropriate boundary conditions.However, for | z | → ∞ and thus (cid:15) → , the dust and gasequations decouple and we obtain lim (cid:15) → v d x = − ηr Ω1 + St , (29) lim (cid:15) → v d y = − ηr Ω1 + St , (30) lim (cid:15) → v g x = 0 , (31) lim (cid:15) → v g y = − ηr Ω , (32)which are constants. These correspond to a sub-Keplerian gas flow that does not feel the dust drag, whilethe dust drifts inwards in response to gas drag. Eqs.29–32 are consistent with the unstratified solutions ofNakagawa et al. (1986).When gas viscosity is neglected, we impose Eqs. 29–30at a finite height z = z max such that (cid:15) (cid:28) . When gasviscosity is included, we impose Eqs. 29–32 at z = z max ,as well as v (cid:48) g x (0) = v (cid:48) g y (0) = 0 , where (cid:48) denotes d/dz . M.-K. Lin
One-fluid models
We also employ the ‘one-fluid’ description of dustygas (Laibe & Price 2014; Price & Laibe 2015; Lin &Youdin 2017) to confirm selected results. In this frame-work, we work with the total mass ρ and the center-of-mass velocity v c of the dust-plus-gas mixture, which istreated as a single, ideal fluid subject to a special cool-ing function. This approximation is valid for small par-ticles with St (cid:28) . Our one-fluid formulation includesdust diffusion, but without gas viscosity (cf. Lovascio &Paardekooper 2019). Details are given in Appendix B. LINEAR PROBLEMWe consider axisymmetric Eulerian perturbations ofthe form δρ g ( z ) exp (i k x x + σt ) , (33)where k x is a (real) radial wavenumber taken to be pos-itive without loss of generality; and σ is the complexfrequency or eigenvalue, σ ≡ s − i ω, (34)where s is the real growth rate and ω is the oscillationfrequency. We also refer to the complex amplitudes suchas δρ g ( z ) and their normalized versions (e.g. δρ g /ρ g ) asthe eigenfunctions. The initial perturbation in real spaceis then obtained by taking Re [ δρ g exp (i k x x )] . Similardefinitions apply to other variables.The linearized equations for the dust fluid read: σ δρ d ρ d + i k x (cid:18) v d x δρ d ρ d + δv d x (cid:19) + ρ (cid:48) d ρ d (cid:18) v d z δρ d ρ d + δv d z (cid:19) v d z (cid:18) δρ d ρ d (cid:19) (cid:48) + v (cid:48) d z δρ d ρ d + δv (cid:48) d z = − Dk x δ(cid:15)(cid:15) + D(cid:15) (cid:34) ρ (cid:48) g ρ g (cid:18) (cid:15) (cid:48) δρ g ρ g + δ(cid:15) (cid:48) (cid:19) + (cid:15) (cid:48)(cid:48) δρ g ρ g + (cid:15) (cid:48) (cid:18) δρ g ρ g (cid:19) (cid:48) + δ(cid:15) (cid:48)(cid:48) (cid:35) , (35) σδv d x + i k x v d x δv d x + v (cid:48) d x δv d z + v d z δv (cid:48) d x = 2Ω δv d y − ΩSt ( δv d x − δv g x ) , (36) σδv d y + i k x v d x δv d y + v (cid:48) d y δv d z + v d z δv (cid:48) d y = − Ω2 δv d x − ΩSt ( δv d y − δv g y ) , (37) σδv d z + i k x v d x δv d z + v (cid:48) d z δv d z + v d z δv (cid:48) d z = − ΩSt ( δv d z − δv g z ) , (38) and that for the gas equations are: σ δρ g ρ g + i k x (cid:18) v g x δρ g ρ g + δv g x (cid:19) + ρ (cid:48) g ρ g δv g z + δv (cid:48) g z = 0 , (39) σδv g x + i k x v g x δv g x + v (cid:48) g x δv g z = 2Ω δv g y − i k x c s δρ g ρ g + δF br x + δF visc x , (40) σδv g y + i k x v g x δv g y + v (cid:48) g y δv g z = − Ω2 δv g x + δF br y + δF visc y , (41) σδv g z + i k x v g x δv g z = − c s (cid:18) δρ g ρ g (cid:19) (cid:48) + δF br z + δF visc z , (42)where the linearized back-reaction force is δ F br ≡ − (cid:15) ΩSt (cid:20) ( v g − v d ) δ(cid:15)(cid:15) + ( δ v g − δ v d ) (cid:21) , (43)and the components of the linearized viscous forces are δF visc x = ν (cid:20) δv (cid:48)(cid:48) g x − k x δv g x + 13 i k x δv (cid:48) g z + ρ (cid:48) g ρ g (cid:0) δv (cid:48) g x + i k x δv g z (cid:1)(cid:21) − ν (cid:18) v (cid:48)(cid:48) g x + ρ (cid:48) g ρ g v (cid:48) g x (cid:19) δρ g ρ g , (44) δF visc y = ν (cid:20) δv (cid:48)(cid:48) g y + ρ (cid:48) g ρ g δv (cid:48) g y − k x δv g y (cid:21) − ν (cid:18) v (cid:48)(cid:48) g y + ρ (cid:48) g ρ g v (cid:48) g y (cid:19) δρ g ρ g , (45) δF visc z = ν (cid:20) δv (cid:48)(cid:48) g z − k x δv g z + 13 i k x δv (cid:48) g x + ρ (cid:48) g ρ g (cid:18) δv (cid:48) g z −
23 i k x δv g x (cid:19)(cid:21) (46)(Lin & Kratter 2016). We remark that one can differen-tiate the gas continuity equation (39) to eliminate δv (cid:48)(cid:48) g z from the expression of δF visc z .In practice, we solve for the perturbation to the dust-to-gas ratio instead of the dust density perturbation,which are related by δρ d ρ d = δ(cid:15)(cid:15) + δρ g ρ g ≡ Q + W. (47)Note that for the strictly isothermal gas we consider W ≡ δρ g /ρ g = δP/P .The linearized equations may be written in the form L g = σ g , (48) tratified streaming instabilities L is an × matrix of linear differential oper-ators and g = [ W, δ v g , Q, δ v d ] † . When supplementedwith appropriate boundary conditions (see below) thisconstitutes an eigenvalue problem for L .4.1. Boundary conditions
We consider modes symmetric about the mid-planesuch that W (cid:48) (0) = Q (cid:48) (0) = δv (cid:48) c ,x (0) = δv (cid:48) c ,y (0) = δv c ,z (0) = 0 , (49)where v c is the center of mass velocity, see Eq. B2. Atthe top boundary we impose W (cid:48) ( z max ) = Q ( z max ) = 0 . (50)When gas viscosity is included we additionally impose δv (cid:48) g x (0) = δv (cid:48) g y (0) = δv (cid:48) g x ( z max ) = δv (cid:48) g y ( z max ) = 0 . (51)We generally find dominant modes have amplitudesthat maximize off the disk-midplane and decay towardsthe domain boundaries, as found by Ishitsu et al. (2009)in direct simulations. As such, boundary conditions areunlikely to modify our main findings.4.2. Numerical method
We use dedalus (Burns et al. 2019), a general-purpose spectral code for solving partial differentialequations, including the linear eigenvalue problem de-scribed above. The eigenfunctions are expanded inChebyshev polynomials T n up to order n = N − , andthe domain is discretized into N points correspondingto the roots of T N . Unless otherwise stated, we take z max (cid:39) H d .We use N = 1024 for computing the background diskstructure and N = 384 for the linearized equations.For the latter, dedalus transforms Eq. 48 into a gen-eralized matrix eigenvalue problem and solves it via theSciPy package (see Section 9D of Burns et al. 2019).This directly yields the eigenvalues σ and the associatedeigenfunctions.For the eigenvalue problem we also use the eigen-tools package to filter out spurious numerical solu-tions due to the discretization. This is done by com-paring eigenvalues obtained from different vertical res-olutions and only keeping those within some tolerance In the one-fluid formulation we use approximate analytic equilib-rium solutions, see Appendix B. https://bitbucket.org/jsoishi/eigentools. (here − ), i.e. only physical solutions that convergewith respect to N z are kept. See Barker & Latter (2015)for a similar treatment. We also filter out unphysical so-lutions with large growth rates compared to Ω (Lin &Youdin 2015).Example source codes for used this work may be ob-tained from the author’s GitHub repository .4.3. Units and notation
We use normalized units such that c s = H g = Ω = 1 ,and quote the dimensionless radial wavenumber K x ≡ k x H g . It turns out that only the reduced pressure-gradient parameter ˆ η = ηr/H g is relevant. Eigenfunc-tions are normalized such that δρ d /ρ d = 1 at its maxi-mum amplitude. For clarity, in plot labels we drop thesubscripts ‘d’, ‘g’, and ‘c’ when collectively referring tothe dust, gas, and center-of-mass velocity fields.We quote δ to distinguish models in which only dustdiffusion is included, from models wherein correspondingviscous terms are also included in the gas momentumequations, in which case we quote α . RESULTSWe present results for a high dust density layer (§5.1),a low dust density layer (§5.2), and a viscous disk (§5.3).For a given set of disk parameters and K x , solving thelinearized equations accounts for all vertical structurespermitted by the boundary conditions and the finite res-olution. This can result in a large number of modes. Weare interested in unstable modes as they will dominateover decaying ones in a real disk. Thus, solutions with s < are discarded and we focus on those with thelargest growth rate s = s max at a given K x . Table 1lists, for each case, approximately the most unstablemode over ≤ K x ≤ (based on two-fluid calcula-tions). 5.1. Case A: High dust density layer
We first present a fiducial case with well-settled dust.Here, we neglect gas viscosity but include particle diffu-sion. This enables a comparison between the two-fluidand one-fluid frameworks, since the latter does not in-clude gas viscosity. We choose Z = 0 . , ˆ η = 0 . , St = 10 − , and δ (cid:39) − . This gives H d (cid:39) . H g . Theequilibrium disk profiles are shown in Fig. 1.Fig. 2 shows the maximum growth rate and corre-sponding frequencies for case A for K x = 100 to .Growth rates increase with K x and maximizes around K x ∼ . We find two classes of unstable modes: https://github.com/minkailin/stratsi. Due to the finite range and sampling in K x -space. M.-K. Lin
Table 1.
Selected unstable modes in stratified dusty disks.Case ˆ η Z St α ( (cid:39) δ ) Viscosity a k x H g / s max / Ω ω/ Ω CommentA 0.05 0.03 − − no . . − . high dust density layer, Fig. 20.01 0.03 − − no . . − . smaller pressure gradient, Fig. 60.1 0.03 − − no . . − . larger pressure gradient, Fig. 60.05 0.03 − − no
10 0 . − . smaller particle, Fig. 70.05 0.03 − − no . . − . larger particle, Fig. 70.05 0.01 − − no . . − . smaller metallicity, Fig. 90.05 0.1 − − no . . − . larger metallicity, Fig. 9B 0.05 0.03 − − no
10 0 . − . low dust density layer, Fig. 11C 0.05 0.01 − − yes . . − . viscous disk, Fig. 15 a Denotes whether or not viscous terms are included the gas momentum equations.
Figure 1.
Two-fluid equilibrium for case A with a highdust density layer. From top to bottom: dust-to-gas ratio,gas density, vertical dust velocity, radial velocities, and az-imuthal velocities. for K x (cid:46) , the oscillation frequency ω ∼ Ω and isconstant; while for K x (cid:38) oscillation frequencies arenegative and increase in magnitude. We find good agree-ment between the one- and two-fluid results, giving con-fidence that these are physical solutions.Example eigenfunctions from the above modes areshown in Fig. 3. We find that with increasing K x , un-stable modes become increasingly localized about z ∼ . H g . Notice there is little perturbation in the gasdensity for either mode, which indicate they are nearlyincompressible.5.1.1. Pseudo-energy decomposition
Figure 2.
Maximum growth rate (top) and correspondingoscillation frequency (bottom) for unstable modes in case A(high dust density layer), as a function of the dimensionlessradial wavenumber K x . In order to identify physical origin of the above in-stabilities, we follow Ishitsu et al. (2009) and examinethe energy-like quantity U tot = (cid:80) i =1 U i associated witheach mode, as described in Appendix C. The contribu-tions U i include: vertical shear in the equilibrium veloc-ity field ( U , which is dominated by the azimuthal com-ponent U y ), vertical dust settling ( U ), pressure forces( U ), dust-gas relative drift ( U ), buoyancy ( U ), andviscosity ( U ). Note that U ≡ for inviscid disks, asconsidered here.Fig. 4 shows the pseudo-energy decomposition of thetwo main type of modes we find. For K x = 100 the modeis driven by a mixture of relative dust-gas drift ( U , red),itself dominated by radial drift, and the vertical shear inthe azimuthal velocity ( U y , crosses). However, the high K x = 3594 mode is entirely driven by vertical shear. tratified streaming instabilities Figure 3.
Normalized eigenfunctions of the most unstable modes found in case A (high dust density layer) with K x = 100 (left)and K x (cid:39) (right, also the most unstable over K x ). Perturbations from top to bottom: relative dust density, relative gasdensity, radial velocity, azimuthal velocity, and vertical velocity. For clarity we plot the amplitudes of the velocity eigenfunctions. Figure 4.
Pseudo-energy decomposition for the modes shown in Fig. 3.
In Fig. 5 we show the vertically-integrated pseudo-energy contributions as a function of K x . We confirmthat the abrupt change in oscillation frequencies around K x = 200 (see Fig. 2) is due to a change in the char-acter of the most unstable mode. For K x (cid:46) , modesare destabilized by a combination of dust-gas drift andvertical shear in the azimuthal velocity; while the latterdominates entirely for K x (cid:38) . Interestingly, dust-gas drift becomes a stabilizing effect (its contributionbecomes negative) for high- K x modes. On the otherhand, vertical shear is always destabilizing.Fig. 5 show that dust settling is always destabiliz-ing, which is consistent with Squire & Hopkins (2018b),who find dust settling alone can lead to instability in vertically-local disk models. However, this effect is sub-dominant in our stratified models because vertical shearis much more significant. We also find that buoyancyforces are always stabilizing, as dust-gas coupling in-creases the mixture’s inertia (Lin & Youdin 2017).We find pressure forces provide increasing stabiliza-tion with increasing K x , which is expected since pressureacts on small scales. The increased restoring force frompressure may explain the increasing magnitude of oscil-lation frequencies (Lubow & Pringle 1993; Balbus 2003).Ishitsu et al. (2009) showed that in the limit of an incom-pressible gas, pressure forces do not contribute to modegrowth or decay. This indicates that gas compressibil-ity becomes non-negligible for the high- K x modes in our0 M.-K. Lin case. These modes have short vertical wavelengths andare localized to regions of largest vertical shear (see Fig.3, right panel).This is reminiscent of ‘surface modes’ ofthe gaseous VSI (Nelson et al. 2013), which are also sta-bilized by gas compressibility (McNally & Pessah 2014).These similarities motivate us to interpret the high- K x ,dust-driven, vertical-shear modes as dusty analogs of thegaseous VSI, see §6.1.5.1.2. Dependence on the global pressure gradient
Fig. 6 shows the effect the global radial pressure gra-dient, as measured by ˆ η ≡ ηr/H g . Note that dust-gasdrift and vertical shear in the azimuthal velocity, whichare the destabilizing effects for the modes we find, bothscale with ˆ η . This is consistent with the drift-drivenclassic SI, which also grows faster with increasing ˆ η ata fixed spatial scale (Jacquet et al. 2011, see their Eq.29). For the vertical shear-driven unstable modes, thediscussion in Appendix 6.1 also indicate that a minimum ˆ η is needed for instability (see Eq. 58). Hence, we findgrowth rates increase with ˆ η .We also find that the transition to modes purely drivenby vertical shear occurs at smaller K x for larger ˆ η : K x > for ˆ η = 0 . and K x > for ˆ η = 0 . . This isagain similar to the gaseous VSI as a weaker verticalshear requires larger radial wavenumbers to destabilize(Latter & Papaloizou 2018).5.1.3. Dependence on particle size
Here we consider Stokes numbers
St = 10 − and St = 0 . . The midplane dust-to-gas ratios are then ∼ and ∼ , respectively (see §3.5.1). Growth rates and fre-quencies are shown in Fig. 7. The trend for St = 10 − isqualitatively similar to our fiducial case with St = 10 − ,but with reduced growth rates. The modes again tran-sition from drift-dominated to vertical-shear dominatedas K x increases, here beyond ∼ , somewhat higherthan the fiducial case. The oscillation frequencies for the St = 10 − vertical-shear modes are also much smaller inmagnitude than that for St = 10 − .For St = 0 . the curves in Fig. 7 are truncated at K x (cid:38) — as we were unable to find convergedtwo-fluid solutions; and one-fluid eigenfunctions werefound to have large, unphysical oscillations near the diskboundary. Notice also the one-fluid model over-predictgrowth rates for K x (cid:38) . This is not surprising as the In the two-fluid calculation with
St = 0 . , numerical artifactsdeveloped near the disk surface, which were remedied by taperingthe equilibrium vertical dust velocity to zero near the upper of the domain. However, this had negligible effects on the growthrates, oscillation frequencies, or the eigenfunctions in the diskbulk. modes have St | σ | (cid:38) Ω , which can invalidate the one-fluid approximation (Lin & Youdin 2017; Paardekooperet al. 2020).Nevertheless, we find St = 0 . growth rates exceedthat for St = 10 − . Moreover, all the St = 0 . modesare driven by vertical shear. However, for K x (cid:38) the two-fluid modes have nearly constant growth ratesand were found to be centered around z = 3 . H d , unlikethe K x (cid:46) modes which are centered around H d ,including the most unstable mode.The two-fluid results show that the most unstablemodes in the St = 10 − and St = 0 . disks occurat K x (cid:38) and K x ∼ , respectively. That is,instability with larger particles occur on larger radialscales. Fig. 8 compares the vertical profiles in the rel-ative dust density perturbations for the most unstablemodes. Here, we re-scale the vertical co-ordinate to ac-count for different dust scale heights in these cases. For St = 10 − the mode is vertically-localized with a char-acteristic lengthscale l z (cid:39) H d / ; whereas for St = 0 . we find l z (cid:39) H d . Thus instability with larger particlesare also more global in the vertical direction.We interpret the above results as looser dust-gas cou-pling provides more rapid ‘cooling’ to mitigate buoyancyforces, which then allows destabilization of disturbanceson longer lengthscales, similar to the gaseous VSI (Lin& Youdin 2015, see also §6.1).5.1.4. Dependence on dust abundance
In Fig. 9 we plot the maximum growth rates and cor-responding frequencies for modes in disks with differ-ent metallicities. We find growth rates are modestly in-creased with increasing solid abundance, but overall theresults are insensitive to Z . In particular, the transitionfrom mixed-modes to vertical shear-dominated modesdoes not depend on Z .5.2. Case B: Low dust density layer
We now consider a low dust density layer with (cid:15) < throughout the disk column by choosing a strongerdiffusion coefficient, δ (cid:39) − , and smaller particles, St = 10 − . Other parameters are the same as the fidu-cial setup in case A. The equilibrium disk profile for caseB is shown in Fig. 10.Growth rates and oscillation frequencies for case B areshown in Fig. 11. We again find two distinct classes ofunstable modes: for K x < growth rates are small( s (cid:46) − Ω ) with oscillation frequencies O (Ω) ; whilefor K x > modes are nearly purely growing with s saturating around . . Case B is much more stablethan case A owing to the smaller dust-to-gas ratio andparticle size. tratified streaming instabilities Figure 5.
Vertically-integrated pseudo-energy contributions for the most unstable modes in case A, as a function of K x . Notethat we plot the cube root for improved visualization. Figure 6.
Maximum growth rate (top) and correspondingoscillation frequency (bottom) for unstable modes in case A(high dust density layer), with different values of the radialpressure gradient ˆ η = 0 . (solid) and ˆ η = 0 . (dashed). Fig. 12 shows the vertically-integrated pseudo-energies. Modes with K x (cid:46) are dominated by dust-gas drift with minor contributions from vertical shearand dust settling. This is unlike for case A where low- K x modes have equal contributions from dust-gas driftand vertical shear. However, for K x (cid:38) modes aredriven by vertical shear, as observed for case A. Figure 7.
Maximum growth rate (top) and correspondingoscillation frequency (bottom) for unstable modes in caseA (high dust density layer), with different Stokes numbers:
St = 10 − (solid) and St = 0 . (dashed). Fig. 13 shows the mode with K x = 100 , primarilydriven by dust-gas drift, involves ultra-short vertical os-cillations of characteristic lengthscale − H g , which ismuch smaller than H d (cid:39) . H g . This should be com-pared to the case A mode in the left panel of Fig. 3,which is driven by a combination of relative dust-gas ra-dial drift and vertical shear, and is more global, i.e. it2 M.-K. Lin
Figure 8.
Magnitude of the relative dust density perturba-tion for the most unstable mode (over K x , based on two-fluidgrowth rates) in case A with St = 10 − (top) and St = 0 . (bottom). Figure 9.
Maximum growth rate (top) and correspondingoscillation frequency (bottom) for unstable modes in case A(high dust density layer), with different metallicities: Z =0 . (solid) and Z = 0 . (dashed). varies on a scale comparable to the dust layer thickness.This suggests that vertical-shear drives a more globaldisk response.For the high- K x modes driven mostly by vertical shearwe find similar behaviors in the eigenfunctions betweencase B and A: modes become increasingly localized withincreasing K x .5.3. Case C: Viscous disk
We briefly examine a viscous disk. To obtain appre-ciable growth rates in the presence of viscosity we set α = 10 − . This value is much smaller than that ex-pected in PPDs, but is sufficient to demonstrate the im-pact of viscosity. We discuss this issue further in §6.3.Here, we use fiducial values of ˆ η = 0 . and St = 10 − ,but set Z = 0 . so that the midplane dust-to-gas ra-tio (cid:15) (cid:39) is similar to case A. We also use a largerdomain with z max = 7 H d as we find viscous modes at Figure 10.
Two-fluid equilibrium for case B with a lowdust density layer. From top to bottom: dust-to-gas ratio,gas density, vertical dust velocity, radial velocities, and az-imuthal velocities.
Figure 11.
Maximum growth rate (top) and correspondingoscillation frequency (bottom) for unstable modes in case B(low dust density layer), as a function of the dimensionlessradial wavenumber K x . the smaller K x values tend to be vertically extended, aresult already hinted by unstratified calculations (Chen& Lin 2020; Umurhan et al. 2020).Fig. 14 show that the equilibrium structure for thisviscous disk is qualitatively similar to case A (Fig .1),except in the radial velocities, which is noticeably non-monotonic away from the mid-plane.Fig. 15 shows the growth rates of the most unstablemodes as a function of K x and corresponding oscillation tratified streaming instabilities Figure 12.
Vertically-integrated pseudo-energy contributions for most unstable modes in case B as a function of K x . Notethat we plot the cube root for improved visualization. Figure 13.
Normalized eigenfunctions of the most unstablemode found in case B (low dust density layer) with K x =100 . Perturbations from top to bottom: relative dust density,relative gas density, radial velocity, azimuthal velocity, andvertical velocity. For clarity, we plot the amplitudes of thevelocity eigenfunctions. Figure 14.
Two-fluid equilibrium for the viscous case C.From top to bottom: dust-to-gas ratio, gas density, verticaldust velocity, radial velocities, and azimuthal velocities. frequencies. For comparison, we also plot results foran inviscid disk. We find viscosity strongly suppressesdust-gas instabilities. In the viscous disk, growth ratesmaximize at K x (cid:39) with s ∼ . and is essentially4 M.-K. Lin
Figure 15.
Growth rates (top) and oscillation frequencies(bottom) for modes found in the viscous case C (black circles)as a function of the dimensionless radial wavenumber K x .Corresponding results for an inviscid disk are also shown(red crosses). quenched for K x (cid:38) ; while growth rates continue toincrease with radial wavenumber in the inviscid disk.In the left panel of Fig. 16 we show a meridional vi-sualization of the most unstable mode found for caseC. The corresponding flow in the inviscid disk is shownin the right panel. As expected, viscosity tends pro-duce vertically-elongated disturbances, here with lengthscales ∼ -2 H d . This mode is again predominantlydriven by vertical shear, as demonstrated by its pseudo-energy decomposition shown in Fig. 17.Fig. 18 shows the vertically-integrated pseudo-energycontributions to the unstable modes. Note that there isnow a viscous contribution (brown curve), see AppendixC. For K x (cid:46) unstable modes are driven by the rela-tive dust-gas radial drift with minor contributions fromdust settling. These modes have relatively small growthrates ( s (cid:46) . ). For (cid:46) K x (cid:46) , modes aredriven by vertical shear with contributions from dust-gas drift. For K x (cid:38) , modes are mostly driven byvertical shear, but their growth rates decline rapidly dueto viscosity, which is more effective at stabilizing smallerlengthscales.As expected, buoyancy and viscous forces always actto stabilize the system. On the other hand, dust set-tling is always destabilizing (Squire & Hopkins 2018b),although here its effect is small. Gas pressure has neg-ligible effects, which reflects the incompressible natureof all the unstable modes presented. Like the inviscidcases A and B, we find dust-gas drift becomes stabiliz-ing at high radial wavenumbers (here (cid:38) ). However,unlike those inviscid cases where vertical shear is always destabilizing, in the viscous disk we find that for low K x modes ( (cid:46) ) vertical shear becomes stabilizing. DISCUSSION6.1.
Vertically-shearing streaming instabilities
Our numerical results show that the most unstablemodes in stratified dust layers occur on radial length-scales (cid:46) − H g and are driven by the vertical gra-dient in the dusty disk’s azimuthal velocity combinedwith partial dust-gas coupling. This is similar to theVSI in gaseous PPDs (Nelson et al. 2013). To interpretthese vertically-shearing streaming instabilities (VSSIs),we invoke the analogy between isothermal dusty gas anda pure gas subject to cooling as developed by Lin &Youdin (2017).For the gaseous VSI, the destabilizing vertical shearresults from the global radial temperature gradient.However, in PPDs vertical gas buoyancy is strongly sta-bilizing. The gaseous VSI thus requires rapid cooling toremove the effect of buoyancy. In terms of these physi-cal quantities, Lin & Youdin (2015) found the instabilitycriterion t cool < | ∂ z v y | N z , (52)where ∂ z v y is the vertical shear rate, N z is the verti-cal buoyancy frequency, and t cool is the thermal coolingtimescale such that the linearized cooling rate is δ Λ = − t cool (cid:18) δP − Pρ g δρ g (cid:19) . (53)We can obtain a criterion analogous to Eq. 52 fordusty disks as follows. We treat the isothermal dustygas as a single fluid subject to a special cooling function,as described in Appendix B. The dusty disk’s azimuthalvelocity profile is given by Eq. B15. The vertical shearrate is thus ∂v y ∂z = ηr Ω (cid:15) (cid:48) (1 + (cid:15) ) . (54)Next, Lin & Youdin (2017) showed that the square ofthe vertical buoyancy frequency in a dusty disk is N z = c s ∂ ln ρ g ∂z ∂f d ∂z , where f d = (cid:15)/ (1 + (cid:15) ) is the dust fraction. Using theequilibrium condition (Eq. B8) we find N z = − z Ω (cid:15) (cid:48) (cid:15) . (55) tratified streaming instabilities Figure 16.
Structure of the most unstable mode found in the viscous case C (left, with K x (cid:39) ) and the correspondingstructure for an unstable mode in an inviscid disk with the same radial wavenumber. Streamlines correspond to the perturbeddust velocity field and colors correspond to the relative dust density perturbation. Figure 17.
Pseudo-energy decomposition of the most un-stable mode found in the viscous disk (case C).
To estimate the appropriate ‘cooling time’ in anisothermal dusty disk, we examine the one-fluid effec-tive energy equation in Appendix B (Eq. B21). We as-sume modes have small lengthscales and write ∂ z → i k z ,where k z is a real vertical wavenumber. Assuming both k x and k z are large in magnitude, the RHS of Eq. B21,which can be interpreted as a cooling rate after multi-plying the equation by P , can be approximated by itsleading term δ Λ d ≡ − c s St (cid:15)k (1 + (cid:15) ) Ω δP, (56) where k = k x + k z . Comparing Eq. 56 with Eq. 53motivates the identification t cool , d ≡ (1 + (cid:15) ) Ω c s St (cid:15)k = (1 + (cid:15) ) (cid:15) St K Ω − (57)as the cooling timescale of an isothermal dusty gas,where K = kH g .Inserting Eq. 57, 55, and 54 into Eq. 52 gives theminimum wavenumber needed to trigger the dust-drivenVSI, K > (1 + (cid:15) ) (cid:15) Stˆ η (cid:18) H d H g (cid:19) zH d . Note that the RHS is a function of height. To obtain amore practical criterion, we evaluate it at z = H d andapproximate (cid:15) ∼ ZH g /H d (see Eq. 3). For settled dustlayers with St (cid:29) δ , as considered throughout this work,we have H d /H g (cid:39) (cid:112) δ/ St . These approximations thengive K (cid:38) (cid:16) Z (cid:112) St /δ (cid:17) δZ ˆ η St . (58)Evaluating Eq. 58 for the fiducial case A ( Z = 0 . , St = 10 − , δ = 10 − , ˆ η = 0 . ) suggest K (cid:38) isneeded for finite dust-gas decoupling to destabilize thedisk through vertical shear. Indeed, for case A we find6 M.-K. Lin
Figure 18.
Vertically-integrated pseudo-energy contributions to the most unstable modes found in the viscous case C, as afunction of K x . Note that we take the cube root for improved visualization. vertical shear contributes to the instability for all K x considered ( ≥ ), see Fig. 5.On the other hand, for case B we have St = 10 − and δ (cid:39) − , giving K (cid:38) . We thus expect a muchlarger K x is needed to tap into the free energy associatedwith vertical shear. Indeed, Fig. 12 show that verticalshear is sub-dominant for modes below the transition at K x ∼ . This is similar to the gaseous VSI (Lin &Youdin 2015): a slower ‘cooling’ rate, here associatedwith stronger dust-gas coupling, means it is can onlyeffectively destabilize smaller length scales.We remark that Eq. 58 can also applied to estimatethe minimum radial pressure gradient needed to triggerVSSIs at a given spatial scale, which may explain theincreasing growth rates with ˆ η seen in §5.1.2.6.2. Comparison to Ishitsu et al. (2009)
Ishitsu et al. (2009) carried out direct simulations toinvestigate the effect of a vertical density gradient on thestability of dust layers. Their disk models were initial-ized with a prescribed, non-uniform vertical dust den-sity distribution and corresponding horizontal velocityprofiles given by Nakagawa et al. (1986). This is equiva-lent to stacking layers of unstratified disk models. Theyneglected vertical gravity, physical diffusion, viscosity,and assumed incompressible gas. We instead consid-ered compressible gas (although this has negligible ef- fects) and solve for the steady vertical disk structureself-consistently, which requires one to at least includedust diffusion.Our results are broadly consistent with Ishitsu et al..Their simulation with
St = 10 − yield disturbances withcharacteristic wavenumber k x ηr ∼ (or K x ∼ as-suming ˆ η = 0 . ) and growth rate ∼ Ω . This is com-parable to our case A with St = 10 − in §5.1.3 (seeFig. 7), for which we find growth rates of . – . for K x (cid:38) . Importantly, Ishitsu et al. also find that thevertical shear in the azimuthal velocity of the dust-gasmixture is the main driver of instability and that dis-turbances are centered off the disk midplane, similar tothat observed in our VSSI eigenfunctions.6.3. Gas viscosity
The example in §5.3 (case C) showed that even a smallamount of gas viscosity of α = 10 − can reduce growthrates significantly. However, in PPDs one may expectup to α ∼ − , for example due to turbulence driven bythe gaseous VSI (Manger et al. 2020). Here, we discusstwo supplementary calculations to further explore therole of gas viscosity.We first repeat the fiducial case A, but enable vis-cosity ( α = 10 − ) and use z max = 7 H d . We find themost unstable mode, shown in Fig. 19, has K x (cid:39) , s (cid:39) × − Ω , ω (cid:39) − × − Ω , and is driven by tratified streaming instabilities Figure 19.
Structure of the most unstable mode in a viscousdisk with α = 10 − and other disk parameters taken fromcase A (see §5.1). Streamlines correspond to the perturbeddust velocity field and colors correspond to the relative dustdensity perturbation. vertical shear. Unsurprisingly, the mode has a muchlarger spatial scale than that for case C (see Fig. 16),with a vertical lengthscale several times larger than H d ( = 0 . H g ). However, the non-negligible mode ampli-tudes near z = z max indicates that boundary conditionsmay be important.Next, we increase viscosity to α = 10 − , which alonewould ‘puff up’ the dust layer such that (cid:15) < . . Wethus compensate by setting Z = 0 . to obtain the samemid-plane dust-to-gas ratio as the above case ( (cid:15) (cid:39) ).We find the most unstable mode has K x (cid:39) . , s (cid:39) × − Ω , and s (cid:39) − Ω . Interestingly, the mode is found tohave comparable contributions from vertical shear anddust settling. However, its radial lengthscale of order r (for h g = 0 . ) calls the shearing box framework itselfinto question.Nevertheless, these calculations are instructive toshow that viscosity is strongly stabilizing and producesdisturbances that exceed the dust layer thickness, simi-lar to the classic SI in unstratified, viscous disk models(Chen & Lin 2020; Umurhan et al. 2020).6.4. Classic streaming instability
In unstratified disks the classic SI of Youdin & Good-man (2005), driven by the dust-gas relative radial drift,is usually the only linearly unstable mode (but see Jau- part & Laibe 2020), which can lead to dust-clumping inthe non-linear regime (Johansen & Youdin 2007).In stratified disks, we find that the relative dust-gasdrift can provide a significant (although not total) con-tribution to the most unstable modes on wavelengthsof O (10 − H g ) or larger. Thus we still refer to them asclassic SI modes. If their nonlinear evolution is similarto their unstratified counterparts, then we can expectdust-clumping on such scales.However, classic SI modes are not the most unstableacross all scales, which are the VSSI modes that occuron radial scales of O (10 − H g ) or smaller, as discussedabove. The nonlinear evolution of classic SI modes willthus be affected by small-scale VSSIs that develop first.If VSSIs saturate in turbulence, as the early work ofIshitsu et al. appear to suggest (albeit based on simu-lations with a limited parameter range and integrationtimes – see §6.7) then we expect a reduced efficiency ofdust-clumping via classic SIs.In this case, we suggest that low-resolution (e.g.global) simulations that are biased towards classic SImodes should include a physical dust diffusion to mimicthe effect of unresolved VSSI turbulence.6.5. Dust settling instability
The dust settling instability (DSI, Squire & Hopkins2018b; Zhuravlev 2019, 2020) is an analog of the classicSI, except the DSI is driven by the vertical drift of dustrelative to the gas (i.e. dust settling), rather than theirrelative radial drift. The DSI has been proposed to seedplanetesimal formation by acting as a dust-clumpingmechanism, although recent simulations show this effectmay be weak in practice (Krapp et al. 2020). Studies ofthe DSI have so far adopted vertically-local disk models,which neglect vertical shear, but permit equilibria to bedefined without dust diffusion. This is not possible inour vertically-global disk models.Our vertically-global models confirm that dust settlingacts to destabilize stratified dusty disks, which suggestthat the DSI is present. However, it is generally sub-dominant to vertical shear, relative radial drift, or both.We can crudely understand this by comparing the dustsettling velocity v d z to the azimuthal velocity differenceacross the dust layer, H d v (cid:48) y . Using Eqs. 19, 21, 54, andconsidering St (cid:28) , we find (cid:12)(cid:12)(cid:12)(cid:12) H d v (cid:48) y v d z (cid:12)(cid:12)(cid:12)(cid:12) (cid:39) (cid:15) (1 + (cid:15) ) ˆ ηδ H d H g ∼ Z ˆ η (1 + (cid:15) ) δ , (59)where we used Z ∼ (cid:15)H d /H g (see Eq. 3). We can fur-ther use δ (cid:39) St Z /(cid:15) to rewrite Eq. 59 as ∼ ˆ η/ (St Z ) ,assuming (cid:15) (cid:38) .For the fiducial case A we find Eq. 59 gives ∼ i.e. dust settling is dwarfed by vertical shear. Hence for8 M.-K. Lin well-settled dust (small δ ), instability is primarily due tovertical shear, provided its free energy can be accessed(see §6.1).Dust settling may dominate over vertical shear if δ (cid:38) Z ˆ η . For PPDs with Z and ˆ η both of O (10 − ) thisrequirement translates to δ (cid:38) − . However, such alarge diffusion parameter is likely to be strongly stabi-lizing (Chen & Lin 2020; Umurhan et al. 2020; Krappet al. 2020).Similarly, we can compare settling to the dust-gas rel-ative radial drift (Eq. 1), say at z = H d , to find (cid:12)(cid:12)(cid:12)(cid:12) v drift v d z (cid:12)(cid:12)(cid:12)(cid:12) ∼ (cid:15) (cid:15) ˆ ηZ , (60)again assuming St (cid:28) . In PPDs the last factor is O (1) .Then for settled dust layers with (cid:15) of O (1) we expectsettling to be at most comparable to radial drift. Forwell-mixed dust layers with (cid:15) ∼ Z ∼ O (10 − ) , the aboveratio is O (ˆ η ) (cid:28) so that dust settling can dominate.However, having such an equilibrium requires a largediffusion coefficient ( δ (cid:29) St ), which may provide com-plete stabilization.6.6. Applicability of RDI theory
Both the classic SI (for (cid:15) (cid:28) ) and the DSI are ‘reso-nant drag instabilities’ (RDI Squire & Hopkins 2018a,b).RDIs arise when the background relative dust-gas mo-tion resonates with a wave in the gas. The local condi-tion for an RDI is k · ( v d − v g ) = ω gas ( k ) (61)(Squire & Hopkins 2018a), where ω gas ( k ) is the fre-quency of a wave mode in the gas (when there is nodust) with local wavenumber k . The small- (cid:15) classic SIand the DSI occur when radial dust drift and verticaldust settling resonates with inertial waves in the gas,respectively. It is then natural to ask whether or notvertically-global modes in our stratified disks can be alsointerpreted as RDIs.The first step of the RDI recipe given by Squire &Hopkins (2018b) is to choose a gas mode in the absenceof dust. Fortunately, analytic dispersion relations canbe obtained for stratified gas disks (Lubow & Pringle1993; Lin & Youdin 2015). Specifically, inertial waves inan isothermal Keplerian disk satisfy ω = LK x + L Ω , (62)where L is an integer and it is assumed that L (cid:29) ω (Barker & Latter 2015; Lin & Youdin 2015).Let us consider inertial waves with L (cid:29) K x . Then ω gas (cid:39) Ω , as observed in the oscillation frequency for modes with K x = 100 in cases A and B (see Figs. 2and 11, respectively). Using the radial drift and dustsettling velocities given by Eq. 1 and 21, respectively,Eq. 61 becomes − ηK x − zH g K z = 1St , where we have assumed St , (cid:15) (cid:28) and K z = k z H g . Wecan use this condition to estimate the vertical wavenum-ber K z required for resonance.Consider the K x = 100 mode in case A with St = 10 − and z = 0 . H g or case B with St = 10 − and z =0 . H g , where radial drift and dust settling contributesto instability. The heights are chosen where mode ampli-tudes maximize, see Figs. 3 (left panel) and 13. We thenfind | K z | ∼ . However, the actual global eigenfunc-tions are better characterized by | K z | ∼ , indicatingsuch modes do not reflect a RDI, at least locally.This discrepancy may be related to the fact that thesemodes are not purely associated with radial drift anddust settling: vertical shear also contributes (see Figs.5 and 12), but this effect does not enter local RDI theory.On the other hand, the dominant VSSI modes areunrelated to the relative dust-gas motion in the back-ground disk. Instead, it is associated with the singleazimuthal velocity of the dust-plus-gas disk. It is thusunclear if VSSI modes can be interpreted as RDIs.It will be necessary to develop a global RDI theory toaddress the above issues.6.7. Implications for planetesimal formation
In non-linear simulations, Ishitsu et al. (2009) showedthat VSSI modes first lead to turbulence. They foundlarge grains with unit St then underwent clumping,possibly due to the classic SI, which is most effective formarginally-coupled solids (Youdin & Goodman 2005).However, small grains ( St = 10 − ) were dispersed bythe turbulence and did not clump, but this may be dueto insufficient metallicities and integration times.Recent simulations carried out by Yang et al. (2017)show that the clumping of small grains ( St = 10 − – − ) require sufficient metallicities ( Z (cid:38) . – . )and integration times ( (cid:38) – orbits). They alsofound dust-clumping occurs after the disk saturates in aturbulent state. It is worth noting that the resolutionsadopted in their simulations, of O (10 − H g ) , should re-solve VSSI modes, which we find to dominate on radialscales of O (10 − H g ) . This regime cannot be probed in our disk models because noequilibrium can be defined, see §3.5.1. tratified streaming instabilities
Caveats and outlooks
Analytical models
We have relied on full numerical solutions to the lin-earized equations. Although our subsequent analyseshint at the physical origin of the various instabilitiesuncovered, a true understanding of the instability mech-anisms require more rigorous mathematical modeling(Jacquet et al. 2011; Squire & Hopkins 2018b; Jaupart& Laibe 2020; Pan & Yu 2020; Pan 2020).To this end, it is desirable to derive an algebraic dis-persion relation for modes in a stratified dusty disk. Thiswill allow us to classify modes and explain their growthas well as oscillation frequencies. This might be pos-sible for the VSSI by exploiting the analogy betweendust-laden flows and pure gas subject to cooling (Lin &Youdin 2017), as in the latter case analytic solutions forthe gaseous VSI can be obtained (Lin & Youdin 2015).6.8.2.
Multiple dust species
We have only considered a single dust species. How-ever, a distribution of particle sizes is expected in reality (Mathis et al. 1977; Birnstiel et al. 2012). Recent gen-eralizations of the classic SI in unstratified disks showthat having multiple dust species can significantly re-duce growth rates when (cid:15) (cid:46) (Krapp et al. 2019; Zhu& Yang 2020; Paardekooper et al. 2020). For the DSI,though, a particle size distribution has a limited effect(Krapp et al. 2020).The VSSI is associated with the vertical shear of thedust-plus-gas system. Considering small, tightly cou-pled grains, all dust species and the gas share the sameazimuthal velocity to O (St) . We may thus naively ex-pect the VSSI to be qualitatively similar for single andmultiple dust species, if the two systems have the samedust-to-gas ratio profile and average Stokes number.This should be checked with explicit calculations.One approach is to add vertical gravity and dust dif-fusion to the one-fluid model of a dusty gas with acontinuous particle size distribution recently developedby Paardekooper et al. (2020). This is equivalent toadding one extra equation for the particle size-densityto those presented in Appendix B, which can then beimplemented in the codes developed for this study.6.8.3. Dust diffusion model
We adopted a simple dust diffusion model so thatstratified equilibria can be defined and standard lin-ear stability analyses can be carried out. Physically,this model assumes there exists an underlying, externalmechanism that stirs up dust grains, such as turbulence.In our implementation this is characterized by a single,constant diffusion coefficient. However, realistic turbu-lence may depend on the disk structure. For example,turbulence driven by the gaseous VSI can be reduced bydust-loading (Lin 2019; Schäfer et al. 2020). In this case,particle diffusion within the dusty midplane should beweaker than the dust-free gas above and below it.Particle stirring may also result from dust-gas insta-bilities itself. Consider, for example, an initially laminardisk. As grains settle, it may (instantaneously) meet theconditions for the classic SI, dust-driven VSI, KHIs, orothers. However, to properly describe how a settlingdust layer becomes unstable, one needs to perform sta-bility analyses with respect to non-steady backgrounds(e.g. Garaud & Lin 2004), which is beyond the scope ofthis work. Nevertheless, such instabilities are thought todrive turbulence that prevents further settling (e.g. Jo-hansen et al. 2009) and maintain a quasi-steady state.In the above contexts, our study should be interpretedas the stability of dust layers whose equilibrium state ismaintained by turbulence driven by pre-existing dust-gas instabilities. To better reflect realistic PPDs, ourmodels should thus be generalized to diffusion (and pos-0
M.-K. Lin sibly gas viscosity) coefficients with strength and spa-tial dependencies based on explicit simulations of quasi-steady, turbulent dust layers (e.g. Bai & Stone 2010b;Yang et al. 2017).A more fundamental issue with the adopted diffu-sion model, though common, is that it can lead tonon-conservation of total angular momentum (Tomi-naga et al. 2019). However, we suspect this will notto qualitatively affect the VSSI as it is unrelated to dustdiffusion: Ishitsu et al. observe the same instabilities,but only included a small diffusion term for numericalstability. Nevertheless, it would be useful to examinethe VSSI in the angular momentum-conserving formal-ism introduced by Tominaga et al..6.8.4.
Non-axisymmetry
Finally, we considered axisymmetric perturbations ex-clusively, which preclude non-axisymmetric KHIs (Chi-ang 2008; Lee et al. 2010). The growth of KHIs and theaxisymmetric instabilities presented in this work shouldbe compared to assess which is more relevant in PPDs.However, in the shearing box framework, linear, non-axisymmetric disturbances may only undergo transientor algebraic growth (e.g. Balbus & Hawley 1992; John-son & Gammie 2005). Describing them requires one tosolve an initial value problem, rather than the eigen-value problem herein. Alternatively, one can forgo theplane wave ansatz and compute the full radial structureof linear disturbances with non-periodic boundary con-ditions (e.g. Adams et al. 1989; Savonije & Heemskerk1990; Lin & Papaloizou 2011a,b). However, this wouldresult in a partial differential equation eigenvalue prob-lem (e.g. Lin 2013), which is significantly more complexthan that considered in this work . SUMMARYIn this paper, we study the axisymmetric linear sta-bility of vertically stratified dust layers in protoplane-tary disks (PPDs). Our disk models extend those usedto study the classic streaming instability (SI, Youdin &Goodman 2005), namely unstratified disks, by account-ing for the vertical structure of dust and gas in PPDs,as solids are expected to settle near the disk midplane. We find the dominant instability in stratified disks isone driven by the vertical gradient in the dusty-gas’ az-imuthal velocity. The large vertical shear within a set-tled dust layer is a significant source of free energy, whichcan be accessed via partial dust-gas coupling. This al-lows unstable modes to grow on orbital timescales. Ourfindings are consistent with earlier non-linear simula-tions carried out by Ishitsu et al. (2009).In PPDs, these vertically-shearing streaming instabil-ities (VSSIs) occur on radial scales (cid:46) − H g , where H g is the local gas scale height. On the other hand, clas-sic SI modes, associated with the relative radial driftbetween dust and gas, occur on radial length scales (cid:38) − H g , but have much smaller growth rates thanVSSIs.However, the non-linear evolution of VSSIs may driveturbulence that mixes up the dust layer (Ishitsu et al.2009), rather than dust clumping like the classic SI (Jo-hansen & Youdin 2007). Given their dynamical growthrates, we suggest VSSI turbulence may have alreadymanifested in some simulations (e.g. Bai & Stone 2010b;Yang et al. 2017), which show that stratified dust lay-ers first settle into a quasi-steady, turbulent state beforeclumping.If VSSIs are inherent to PPDs and its primaryoutcome is turbulence, then planetesimal formationthrough the classic SI may be less efficient than pre-viously thought, as clumping will always be hinderedby small-scale VSSI-turbulence. High-resolution simu-lations that fully resolve VSSI scales will be necessaryto clarify this issue.ACKNOWLEDGMENTSI thank the anonymous referee for a thorough re-port that considerably improved the connection betweenthis work and the literature. I thank Volker Elling,Pin-Gao Gu, Ming-Chih Lai, Yueh-Ning Lee, Te-ShengLin, Jack Ng, Debanjan Sengupta, Ryosuke Tominaga,Orkan Umurhan, and David C.C. Yen, for useful discus-sions, tips, and advice. This work is supported by Tai-wan’s Ministry of Science and Education through grant107-2112-M001-043-MY3.APPENDIX A. LIST OF SYMBOLSTable 2 summarizes frequently used and related symbols in the main text. tratified streaming instabilities Table 2.
Frequently used symbolsNotation Definition Description ρ d , g Dust and gas densities v d , g Dust and gas velocities in the shearing box, relative to Keplerian flow (cid:15), (cid:15) ρ d /ρ g , (cid:15) ( z = 0) Local dust-to-gas ratio, midplane dust-to-gas ratio Σ d , g (cid:82) ∞−∞ ρ d , g dz Dust and gas surface densities Z Σ d / Σ g Metallicity c s Constant gas sound-speed Ω (cid:112) GM ∗ /r Keplerian rotation frequency H g , h g c s / Ω , H g /r Gas disk pressure scale-height, aspect-ratio
P c s ρ g Pressure in the global disk or local pressure fluctuations in the shearing box η , ˆ η − (cid:0) r Ω ρ g (cid:1) − ∂ r P , η/h g Dimensionless global pressure gradient, reduced pressure gradient parameter α Dimensionless gas viscosity δ (cid:0) (cid:1)(cid:14) (cid:0) (cid:1) α Dimensionless dust diffusion coefficient St τ s Ω Stokes number with particle stopping time τ s H d (cid:112) δ/ ( δ + St) H g Dust scale-height δρ g , etc. Complex amplitude of Eulerian perturbations (eigenfunctions) s , ω Re( σ ) , − Im( σ ) Growth rate, oscillation frequency of the complex growth rate σK x k x H g Dimensionless radial wavenumberB.
ONE-FLUID MODEL OF DUSTY GAS IN THE SHEARING BOXIn the ‘one-fluid’ description of dusty-gas we work with the total density ρ = ρ g + ρ d (B1)and center-of-mass velocity v c = ρ g v g + ρ d v d ρ . (B2)Furthermore, by considering small, tightly-coupled dust particles we relate the gas and dust velocities by the ‘terminalvelocity approximation’, v d = v g + t s (cid:18) ∇ Pρ g − ηr Ω ˆ x (cid:19) (B3)(Youdin & Goodman 2005; Laibe & Price 2014). Here, t s = τ s ρ g /ρ is the relative stopping time. Recall that v d , g are dust and gas velocities in the shearing box relative to the Keplerian flow, respectively, and P is the local pressurefluctuation. The term ∝ η represents the radial pressure gradient in the global disk. Thus, in the unperturbed statewith vanishing P , dust drifts radially relative to the gas.The dust-gas mixture is modeled as a single, adiabatic fluid with a special cooling function (Lin & Youdin 2017;Lovascio & Paardekooper 2019). Dropping the subscript ‘c’ for clarity, our one-fluid model equations in the localshearing box to O ( t s ) are ∂ρ∂t + ∇ · ( ρ v ) = ∇ · ( Dρ g ∇ (cid:15) ) , (B4) ∂ v ∂t + v · ∇ v = − ∇ Pρ + 2 ηr Ω ρ g ρ ˆ x + 2Ω v y ˆ x − Ω2 v x ˆ y − Ω z ˆ z , (B5) ∂P∂t + ∇ · ( P v ) = c s ∇ · (cid:2) t s f d (cid:0) ∇ P − ηr Ω ρ g ˆ x (cid:1)(cid:3) . (B6)2 M.-K. Lin (Laibe & Price 2014; Lin & Youdin 2017; Lovascio & Paardekooper 2019; Chen & Lin 2020; Paardekooper et al. 2020),where f d ≡ ρ d /ρ is the dust fraction. Note that ρ g = P/c s for the isothermal gas we consider, so f d = 1 − P/c s ρ . Thesecond term on the RHS of Eq. B5 vanishes in a particle disk where ρ g → , since solids do not feel pressure gradients.Conversely, for a dust-free gas disk ( ρ g /ρ → ) we recover the full pressure support from the global disk. The secondterm in the parenthesis on the RHS of Eq. B6 arises from the contribution to the terminal velocity approximationfrom the large-scale radial pressure gradient in Eq. B3.B.1. Approximate equilibria
As in the two-fluid model we seek steady, horizontally uniform solutions. The mass, energy, and vertical momentumequations are then ρv z = Dρ d (ln (cid:15) ) (cid:48) , (B7) v z v (cid:48) z = − P (cid:48) ρ − Ω z (cid:39) , (B8) P v z = c s t s f d P (cid:48) , (B9)where in Eq. B8 we neglect the O ( v z ) term a posterior for consistency with the small t s approximation used to derivethe one-fluid model. These equations may then be solved for constant τ s = t s /f g = St / Ω to yield (cid:15) = (cid:15) exp (cid:18) − St2 δ z H (cid:19) , (B10) v z = − (cid:15) (cid:15) St z Ω , (B11) P = P exp (cid:20) δ St ( (cid:15) − (cid:15) ) − z H (cid:21) . (B12)From here it is clear that v z = O (St) , so it is self-consistent to neglect the O ( v z ) term. Eq. B10–B12 are in fact thesame solutions as in the full two fluid model in the limit St (cid:28) (see Eq. 19–21, recall β → St for St → ; and note v z = f d v d z ).The horizontal momentum equations are v z v (cid:48) x = 2 ηr Ω ρ g ρ + 2Ω v y (cid:39) , (B13) v z v (cid:48) y = − Ω2 v x , (B14)where we neglect the quadratic term in Eq. B13 to obtain v y = − ηr Ω1 + (cid:15) . (B15)This is expected on physical grounds for tightly-coupled dust ( St → ). In this limit the mixture behaves close to asingle fluid with orbital velocity depending on the level of dust enrichment. For (cid:15) → we have a pressure-supportedgas disk at sub-Keplerian velocity (assuming η > ); while (cid:15) → ∞ corresponds to a particle disk on exactly Keplerianorbits, since then v y → .Next, we use Eq. B15, Eq. B14, and Eq. B11 to obtain v x = − v z v (cid:48) y = 2 ηr Ω (cid:15)(cid:15) (cid:48) (1 + (cid:15) ) St z. (B16)Notice the radial velocity of the dust-gas mixture’s center-of-mass depends on height. It is only zero at the midplaneand for | z | → ∞ where (cid:15) → . For η > the specific angular momentum decreases with increasing | z | : the mixturegains pressure support as it becomes more gas-rich away from the midplane. This means that as a parcel of the mixturesettles, it finds itself having an angular momentum deficit compared to its surrounding; it thus drifts inwards ( v x < ),as indicated by Eq. B16. tratified streaming instabilities Linearized equations
We linearize the one-fluid equations about the above basic state, with non-uniform v x ( z ) , v y ( z ) , and v z ( z ) . As in themain text we assume axisymmetric perturbations in the form of δρ ( z ) exp ( σt + i k x x ) , and similarly for other variables.The linearized equations are σ δρρ + i k x (cid:18) δv x + v x δρρ (cid:19) + ρ (cid:48) ρ (cid:18) v z δρρ + δv z (cid:19) + v z (cid:18) δρρ (cid:19) (cid:48) + v (cid:48) z δρρ + δv (cid:48) z = Df d (cid:34) Q (cid:48)(cid:48) − k x Q + ρ (cid:48) d ρ d (cid:18) (cid:15) (cid:48) (cid:15) δρ d ρ d + Q (cid:48) (cid:19) + (cid:15) (cid:48) (cid:15) (cid:18) δρ d ρ d (cid:19) (cid:48) + (ln (cid:15) ) (cid:48)(cid:48) δρ d ρ d (cid:35) , (B17) σδv x + i k x v x δv x + v (cid:48) x δv z + v z δv (cid:48) x = − i k x Pρ W − ηr Ω (cid:15) (1 + (cid:15) ) Q + 2Ω δv y , (B18) σδv y + i k x v x δv y + v (cid:48) y δv z + v z δv (cid:48) y = − Ω2 δv x , (B19) σδv z + i k x v x δv z + v (cid:48) z δv z + v z δv (cid:48) z = P (cid:48) ρ (cid:15)Q (cid:15) − Pρ W (cid:48) , (B20) σW + i k x ( δv x + v x W ) + P (cid:48) P ( v z W + δv z ) + v (cid:48) z W + v z W (cid:48) + δv (cid:48) z = K P (cid:34) W (cid:48)(cid:48) − k x W + K (cid:48) K (cid:18) P (cid:48) P δ KK + W (cid:48) (cid:19) + P (cid:48) P (cid:18) δ KK (cid:19) (cid:48) + (ln P ) (cid:48)(cid:48) δ KK (cid:35) − k x ηr Ω K c s P (cid:20)(cid:18) − (cid:15) (cid:15) (cid:19) Q + W (cid:21) , (B21)where K ≡ c s P St (cid:15) Ω(1 + (cid:15) ) , (B22)and recall Q = δ(cid:15)/(cid:15) and W = δρ g /ρ g = δP/P . B.3. Mode energetics
Following Ishitsu et al. (2009), we multiply the x , y , and z momentum equations (Eq. B18–B20) by δv ∗ x,y,z ,respectively, combine them appropriately, then take the real part. We also scale the overall result by a factor of (1 + (cid:15) ) for easier comparison with the corresponding two-fluid treatment in Appendix C. The one-fluid result is E tot ≡ (1 + (cid:15) ) (cid:16) | δv x | + 4 | δv y | + | δv z | (cid:17) = (cid:88) i =1 E i , (B23)with sE = − (1 + (cid:15) ) (cid:104) v (cid:48) x Re ( δv z δv ∗ x ) + 4 v (cid:48) y Re (cid:0) δv z δv ∗ y (cid:1) + v (cid:48) z | δv z | (cid:105) (B24) ≡ sE x + sE y + sE z sE = − v z (1 + (cid:15) ) Re (cid:0) δv (cid:48) x δv ∗ x + 4 δv (cid:48) y δv ∗ y + δv (cid:48) z δv ∗ z (cid:1) (B25) sE = c s [ k x Im (
W δv ∗ x ) − Re ( W (cid:48) δv ∗ z )] (B26) sE = − ηr Ω (cid:15) (1 + (cid:15) ) Re ( Qδv ∗ x ) = (cid:15) Ω ( v d x − v g x )St Re ( Qδv ∗ x ) (B27) sE = − (cid:15)z Ω Re (