Streaming Instabilities in Accreting and Magnetized Laminar Protoplanetary Disks

The streaming instability (SI) is one of the most promising pathways to the formation of planetesimals from pebbles. Understanding how this instability operates under realistic conditions expected in protoplanetary disks (PPDs) is therefore crucial to assess the efficiency of planet formation. Contemporary models of PPDs show that magnetic fields are key to driving gas accretion through large-scale, laminar magnetic stresses. However, the effect of such magnetic fields on the SI has not been examined in detail. To this end, we study the stability of dusty, magneftized gas in a protoplanetary disk. We find the SI can be enhanced by passive magnetic torques and even persist in the absence of a global radial pressure gradient. In this case, instability is attributed to the azimuthal drift between dust and gas, unlike the classical SI, which is driven by radial drift. This suggests that the SI can remain effective inside dust-trapping pressure bumps in accreting disks. When a live vertical field is considered, we find the magneto-rotational instability can be damped by dust feedback, while the classic SI can be stabilized by magnetic perturbations. We also find that Alfvén waves can be destabilized by dust–gas drift, but this instability requires nearly ideal conditions. We discuss the possible implications of these results for dust dynamics and planetesimal formation in PPDs.


Introduction
The formation of 1-100 km sized planetesimals is a key stage in the core accretion theory of planet formation (Chiang & Youdin 2010;Johansen et al. 2014;Birnstiel et al. 2016) but is still not well understood. Specifically, the growth of solids from micron-sized grains through sticking is limited to millimeter to centimeter sized pebbles, beyond which collisions result in bouncing or fragmentation (Blum & Wurm 2008;Blum 2018). Furthermore, gas drag can lead to a rapid inwards drift of solids (Whipple 1972;Weidenschilling 1977), which introduces a radialdrift barrier (Birnstiel et al. 2010(Birnstiel et al. , 2012. One way to circumvent these growth barriers is the collective, self-gravitational collapse of a particle swarm directly into planetesimals (Goldreich & Ward 1973;Youdin & Shu 2002). To do so, the particle swarm must attain a dust-to-gas mass density ratio, ò, well above unity (Shi & Chiang 2013), which should be compared to the typical value of ò ∼ 0.01 in the interstellar medium and is a reasonable expectation in protoplanetary disks (PPDs; Testi et al. 2014), at least initially. Thus, some other mechanism is needed to enhance the local dust-to-gas ratio.
The SI is a linear instability in rotating flows of dust and gas when the two components interact through mutual drag (Jacquet et al. 2011;Lin & Youdin 2017;Jaupart & Laibe 2020;Pan 2020Pan , 2021Pan & Yu 2020; and is powered by their relative radial drift, which itself is usually driven by a global radial pressure gradient that offsets the gas from Keplerian rotation. More recently, the SI was shown to be a member of a broader class of resonant drag instabilities (RDIs; Squire & Hopkins 2018a, 2018bZhuravlev 2019) in dusty gas, where instability arises from a resonance between neutral waves in the gas and the relative motion between dust and gas.
On the other hand, the influence of a magnetic field on the SI has been less well explored, but PPDs are expected to be magnetized, albeit subject to nonideal magneto-hydrodynamic (MHD) effects (Lesur 2021). Early studies focused on the impact of MHD turbulence sustained by the magneto-rotational instability (MRI; Balbus & Hawley 1991) on dust dynamics, including their vertical settling and radial diffusion and migration (Fromang & Nelson 2005;Johansen & Klahr 2005;Fromang & Papaloizou 2006;Johansen et al. 2006). Subsequent simulations do find dust clumping in MHD-turbulent disks, which has been attributed to the SI Balsara et al. 2009;Tilley et al. 2010), weak radial diffusion in resistive disks , or concentration by zonal flows or pressure bumps induced by the MRI (Johansen et al. 2009a(Johansen et al. , 2011Dittrich et al. 2013;Xu & Bai 2022). The latter two effects may facilitate the SI as it requires ò  1 for dynamical growth, although pressure bumps present a dilemma as the classic SI does not formally operate without a radial pressure gradient.
In addition to driving small-scale (perhaps weak) turbulence, magnetic fields are also expected to regulate the large-scale gas dynamics of PPDs. Modern numerical simulations often find magnetized disk winds that drive laminar accretion (Gressel et al. 2015;Bai 2017;Béthune et al. 2017;Wang et al. 2019;Gressel et al. 2020). These windy PPDs are the new paradigm for planet formation and evolution. In some of these models, namely those including the Hall effect with the field and disk rotation being aligned, large-scale horizontal magnetic stresses can develop in the disk midplane, which leads to radial gas accretion. Such inflows can strongly modify the orbital migration of protoplanets compared to conventional, alphatype viscous disk models (McNally et al. 2017;Kimmig et al. 2020).
A natural question is how do planets form in these laminar, accreting PPDs in the first place? Related to this issue is the formation of rings or pressure bumps commonly seen in such simulations (Béthune et al. 2016;Suriano et al. 2017Suriano et al. , 2018Suriano et al. , 2019Riols & Lesur 2019;Cui & Bai 2021), which can act as dust traps (Krapp et al. 2018;Riols et al. 2020) and may be sites of preferential planetesimal formation, for example through the SI. These pressure bumps may also explain dust rings observed in bright PPDs (e.g., Andrews et al. 2018;Long et al. 2018).
In this work, we examine the basic properties of dusty-gas dynamics in a magnetized PPD. Motivated by the aforementioned studies, we are particularly interested in the effect of a magnetically driven accretion flow on the SI, as well as how the MRI and SI interact, and their implications for planetesimal formation. We find that a laminar gas accretion flow can enhance the SI and lead to instability even without radial pressure gradients, which suggests that pressure bumps remain feasible sites for the SI. We also find dust loading can reduce MRI growth rates, while magnetic perturbations can be effective in stabilizing the classic SI. Finally, we demonstrate a RDI unique to dusty and magnetized gas, in which Alfvén waves are rendered unstable by dust-gas drift, although it may have limited relevance to realistic PPDs. This paper is organized as follows. In Section 2 we list the basic equations for a dusty, magnetized PPD and specify the physical setups under consideration. We describe the linear problem in Section 3 and give an overview of analytic results, some of which are developed in the Appendixes. We present numerical results in Section 4. In Section 5 we verify the main findings from our linear-stability analyses with direct integration of the dusty MHD equations. We discuss the implications of our results to dust dynamics in PPDs in Section 6 before summarizing in Section 7.

Basic Equations
We consider a three-dimensional (3D) PPD comprised of gas and a single species of uncharged dust grains in orbit around a central star of mass M * . Cylindrical coordinates (R, f, z) are centered on the star. The gas component has density, pressure, and velocity fields (ρ g , P, V), respectively, and is threaded by a magnetic field, B. We consider small grains (defined below) and approximate the dust population as a pressureless fluid with density and velocity fields (ρ d , W), respectively. The frictional drag between dust and gas is characterized by the stopping time, τ s .
We make several simplifications for tractable analyses. We consider a strictly isothermal gas so that r = P C s 2 g , where C s = H g Ω K is a constant sound-speed. Here H g is the pressure scale height, W º * GM R K 3 is the Keplerian frequency, and G is the gravitational constant. We neglect gas viscosity and dust diffusion (Dubrulle et al. 1995) except in some calculations for regularization and comparison purposes. As a proxy for nonideal MHD effects, we consider ohmic resistivity with a constant diffusion coefficient, η O . For the dust, we assume a constant Stokes number, St ≡ τ s Ω K . Then the fluid approximation applies to small grains with St = 1 (Jacquet et al. 2011), which we assume throughout. We focus on dynamics close to the disk midplane (z = 0) and thus neglect the vertical component of stellar gravity, so our models are unstratified (∂ z = 0 in equilibrium). We also impose axisymmetry from the outset (∂ f ≡ 0).
Under these approximations, the governing equations for the magnetized gas are as follows: O together with ∇ · B = 0, μ 0 is the magnetic permeability, and ò ≡ ρ d /ρ g is the dust-to-gas ratio. The dust equations are

Steady-state Drift with a Magnetic Field
We consider axisymmetric (∂ f = 0), unstratified (∂ z = 0) basic states with vanishing vertical velocities (V z = W z = 0). The vertical component of the magnetic field B z is taken to be a constant, which does not affect the equilibrium solutions below. The solenoidal condition implies B R ∝ R −1 . Hence, we write where R 0 is a reference radius and B R0 is a constant. For a thin disk we expect f W V R K . Inserting this into the induction equation, we find where Ω 0 = Ω K (R 0 ); see also McNally et al. (2017). We next seek approximate solutions to the horizontal velocity fields given the above field configuration. We write and assume the disk is nearly Keplerian so that ¢ W | |  V R K , and similarly for W. Neglecting terms quadratic in the primed variables, the equilibrium equations become where F R,f are the radial and azimuthal components of the Lorentz force: with B R,f given by Equations (6)- (7), and is a dimensionless measure of the radial gas pressure gradient. Solving Equations (8)- (11) give is a dimensionless measure of the total (gas plus magnetic) radial pressure gradient. Equations (14)-(17) generalize the standard, two-fluid steady-state drift solutions for a dusty gas (Nakagawa et al. 1986) to include the effect of horizontal magnetic fields in a resistive disk, where magnetic torques drive gas accretion. See Umurhan et al. (2020) for a similar set of equations that account for viscous gas accretion. We remark that the above velocity fields do not produce global steady states for arbitrary density fields. This requires the mass fluxes r ¢ RV R g and r ¢ RW R d to be global constants, which constrains the density profiles. We explore this in Appendix A for the case of constant St and ò. However, small-scale instabilities such as the SI are not expected to be sensitive to the global density profile. Indeed, it is possible to model them in a local disk model, in which case a strict steady state can be obtained for constant densities, velocities, pressure gradients, and magnetic torques, as we do below.

Dust Drift Indirectly Induced by Magnetic Fields
In the limit ò → 0, or negligible dust feedback, we find so the gas rotates at a sub-Keplerian speed (assuming η tot > 0) and drifts inward since F f < 0. The gas accretion rate is a constant: In this limit the dust radial drift is Thus the inwards drift of dust is enhanced by the magnetic torque. This is expected as the dust is partially coupled to the inwardly accreting gas through drag. The dust, therefore, feels the magnetic field indirectly.

Shearing Box Approximation
We study the local dynamics of the above system by focusing on a small patch of the disk using the shearing box framework (Goldreich & Lynden-Bell 1965). The box is anchored at a fiducial point (R 0 , f 0 , 0) that rotates around the star at an angular frequency of Ω 0 , so f 0 = Ω 0 t. Cartesian coordinates (x, y, z) in the box correspond to the radial, azimuthal, and vertical directions in the global disk. For a sufficiently small box size (=R 0 ) we can ignore curvature effects and approximate Keplerian rotation as = -WÛ y qx K 0 with q = 3/2. We then define v and w as the dust and gas velocities in the shearing box relative to this linear shear flow. The total gravitational and centrifugal force in the box is Wx x 3 0 2 . In terms of velocity deviations, the axisymmetric shearing box equations read Yang et al. 2018), where subscript 0 here and below denotes evaluation at the reference radius in steady state. In Equation (23), we include the effect of a global gas-plusmagnetic pressure gradient through the term ∝ η tot,0 and that from a large-scale horizontal magnetic torque through F f0 . In the shearing box approximation these are modeled as constant forcing terms that do not respond to the dynamics in the box. Note that P here refers to pressure fluctuations and is zero in equilibrium. An exact equilibrium state consists of constant densities and velocities, with the latter given by Equations (14)

Problem Specification
We study two specific roles of magnetic fields and simplify Equations (22)-(26) accordingly: 1. Case I. The effect of gas accretion flows induced by large-scale, horizontal fields. Here, the magnetic field is only included in the equilibrium disk and is neglected in the perturbed state. We thus set the Lorentz force to zero in Equation (23)  which enables the MRI and other magnetic modes. A vertical field does not modify the equilibrium from the hydrodynamic limit. However, magnetic forces are fully active in the perturbed state, i.e., we include Lorentz forces and the induction equation. Note that a live radial field cannot be included initially as it would be sheared apart by differential rotation so no equilibrium can be constructed. For simplicity, we also neglect background horizontal fields, thus η tot,0 → η 0 and F f0 → 0, and there is no magnetically induced gas accretion.

Physical Parameters
Our magnetized, dusty disks are characterized by several parameters, some of which are specific to Cases I and II. For clarity, we drop the subscript 0 notation with the understanding that all variables below are evaluated at the reference radius in the equilibrium state.
The disk opening angle at the reference radius, is a measure of the disk temperature. In this work we fix h g = 0.05. The global pressure gradient, η tot , is typically of ( ) O h g 2 and is also a measure of gas compressibility when considering the SI (Youdin & Goodman 2005;.
For convenience we also define the reduced pressure gradient parameter: We take h  0  unless otherwise stated. We will consider small values of h  to explore how the SI behaves around pressure extrema.
For the magnetic field, we define the azimuthal and vertical Alfvén speeds: A z z , , 0 g and the plasma beta parameters: to quantify the (inverse) strength of the magnetic field in Cases I and II, respectively. We use the Elsässer numbers , to quantify the effect of ohmic resistivity. A given disk model is parameterized by h , St, ò, β f,z , and Λ f,z . However, for Case I the magnetic parameters only appear in the azimuthal force, F f , which physically represents Maxwell stresses from the background disk. It is therefore convenient to define the dimensionless stress Our results for Case I do not depend on values of β f and Λ f separately, but only on the magnetic stress, α M . Note that α M should not be interpreted as the viscous α SS parameter of Shakura & Sunyaev (1973) that attempts to mimic small-scale turbulence. Here, α M characterizes radial angular momentum transported by large-scale, horizontal magnetic torques, although our results are also applicable to gas accretion driven by other means (see Section 6.1).

Radial and Azimuthal Drifts
We define a dimensionless measure of the relative radial dust-gas drift as The first term in the square brackets of Equation ( is attributed to the magnetic torque. For the same parameters as above, the critical h  is O(10 −3 ), implying that as one approaches a pressure extremum the magnetically induced azimuthal drift dominates over that due to pressure gradients, before the same occurs for the radial drift. We remark that while the radial pressure gradient causes an O(St) radial drift and an O(St 2 ) azimuthal drift, the magnetic torque has the opposite effect in driving an O(St) azimuthal drift and an O(St 2 ) radial drift. When h  = 0, the azimuthal drift is larger than the radial drift by a factor of O(St −1 ).

Linear Theory
We consider axisymmetric Eulerian perturbations for any variable A such that where δA is a complex amplitude; k x and k z are real radial and vertical wavenumbers, respectively; and the complex growth rate σ = s − iω, where s is the real growth rate and ω is the oscillation frequency. For the linear problem, we take k x,z > 0 without loss of generality. Explicit expressions of the linearized equations are given in Appendix B, where we also reproduce the standard MRI and SI.
The linearized system constitutes the eigenvalue problem where M is the matrix representation of the right-hand side of Equations (B1)-(B11) and q is the eigenvector of the complex amplitudes. Note that for Case I, M is a 8 × 8 matrix and since the induction equation is dropped. For Case II, M is an 11 × 11 matrix and q additionally includes magnetic field perturbations. We solve Equation (36) using standard numerical methods in MATLAB. In addition to the physical parameters described in Section 2.5, each calculation also depends on the wavenumbers k x,z .

Streaming Instabilities Driven by Azimuthal Drift
The classic SI of Youdin & Goodman (2005) is driven by the radial drift between dust and gas, ζ x . There also exists an azimuthal drift, ζ y , between dust and gas (Equation (34)), but this is subdominant to the radial drift provided that magnetic torques are weak compared to pressure gradients, which is usually the case.
However, as discussed in Section 2.6, azimuthal drifts can dominate over radial drifts near pressure extrema. Indeed, for Case I and sufficiently small h , we find the SI can still operate, but is now related to the azimuthal drift, and thereby to the magnetically induced accretion flow. In Appendix C we present a simplified model of this new form of SI. We find they have growth rates 3 7 x y for K z = 0, where K x,z = k x,z H g are normalized wavenumbers. This is distinct from the classic SI that requires K z ≠ 0 (Youdin & Goodman 2005). This azimuthal-drift-driven SI also differs from the RDI described below, which would necessarily require nonaxisymmetric disturbances.

Resonant Drag Instabilities
A powerful description of instabilities driven by mutual dust-gas drag is the RDI theory developed by Squire & Hopkins (2018a, 2018b. This family of instabilities arise from a resonance between waves in the gas and the relative drift between dust and gas when ω gas (k) = k · (w − v), where ω gas is the neutral frequency of a wave mode in the gas when there is no dust, and k is the wavevector.
In our unstratified, axisymmetric shearing box the resonance condition is .
3 8 x z x x x gas Equation (38) gives the relation between k x and k z that maximizes growth rates. Since a variety of waves can be supported in a gas, RDIs are generic so that dusty gas is generally unstable. Note that the azimuthal drift (w y − v y ) cannot cause an RDI in axisymmetric disks.
Although RDI theory formally applies to the limit ò = 1, we find Equation (40) is still a useful guide in identifying classic SI modes even when ò > 1. Note, also, that RDIs do not distinguish between different origins of dust-gas drift. Thus, for example, at a pressure extremum (h  = 0) the radial drift caused by the magnetic torque alone can, in principle, also drive RDIs. However, in practice, we find the azimuthal drift dominates in this limit and leads to the non-RDI instability described in Section 3.1.

Magnetic Fields and the Classic SI
When a live magnetic field is present, for example in Case II, we find that magnetic perturbations can stabilize the classic SI, most notably in dust-poor disks, even with large resistivities. In Appendix D we present a toy model based on a one-fluid description of a dusty, magnetized gas.
to leading order in Λ z , where s s = Ŵ . One can then show that at where the resonant K z → ∞ , SI growth rates decrease as Λ z increases from zero, presumably due to stabilization by magnetic tension.
A magnetized gas can support additional MHD modes: the fast and slow magneto-sonic and Alfvén waves. In a dusty media these MHD waves can resonate with the dust-gas drift and produce a variety of instabilities (Hopkins & Squire 2018;Seligman et al. 2019;Hopkins et al. 2020). However, their relevance to PPDs may be limited by nonideal effects (Hopkins & Squire 2018), a conclusion consistent with our numerical calculations below.
We demonstrate with Case II that Alfvén waves can drive streaming-type instabilities. When rotation is neglected, the dispersion relation for an Alfvén wave with a purely vertical background field is w = ( ) C k 42 Az z gas 2 2 2 (e.g., Ogilvie 2016). Using Equations (42), (33), and (38), the resonance condition between Alfvén waves and dust-gas drift becomes where the second expression applies to small grains in dustpoor disks without a horizontal field. We can therefore expect instabilities in a magnetized, dusty disk with K z ∝ K x . For St = 0.1, h =  0.05, β z = 10 4 , we find K z ; K x , while smaller β z reduces the resonant K z for a given K x . That is, a stronger vertical field produces more vertically elongated disturbances.
For Case II we only consider an initially vertical field. If a radial field is also present, Equation (42) where C AR is the Alfvén speed associated with the radial field. The RDI condition then generalizes to 2 , implying the resonant K z is reduced in magnitude. However, a live radial field will be wound up by the background shear so that β R →∞, and this reduction eventually vanishes. The azimuthal field generated by the shear would take no part in axisymmetric RDIs.

Dynamical Effect of Dust on the MRI
The presence of tightly coupled dust is expected to reduce the Alfvén speed since the total density increases with dust loading, ρ g → ρ g (1 + ò). One can then define as the effective vertical Alfvén speed of a dusty gas, and similarly for the azimuthal speed. In ideal MHD this reduction increases the most unstable k z ( ∝ 1/C Az,eff ), while in resistive disks dust loading reduces growth rates (µC Az,eff 2 ) and the most unstable k z ( ∝ C Az,eff ); see Sano & Miyama (1999) for expressions in the pure gas limit. We demonstrate these effects in Sections 4.3 and 4.4.

Numerical Results
In this section, we solve Equation (36) numerically to obtain the dispersion relation s s , . To connect our results to previous studies, we base our setups on that for the LinA (St = 0.1, ò = 3) and LinB (St = 0.1, ò = 0.2) SI eigenmodes described in . We normalize growth rates by Ω and wavenumbers by H g . Note that this differs from the conventional wavenumber normalization by η tot R, since we will consider disks with η tot = 0. Otherwise, the fiducial value of the reduced pressure gradient . For LinA and LinB the fiducial wavenumbers are, then, K x,z = 600 and K x,z = 120, respectively. We apply a lower limit to the growth rates at s = 10 −5 Ω.
In Table 1 we list the selected modes in our dusty, magnetized disks. These include SI modes modified by a background accretion flow, SI modes at pressure extrema that are driven by azimuthal drift, and the SI of Alfvén waves.
To help identify the dominant source of the instabilities uncovered, we also examine the pseudo-kinetic energy of the modes, Ishitsu et al. 2009; Lin 2021), constructed from the eigenvectors, where U 1 corresponds to thermal or pressure contributions, U 2 to that from dust-gas drift, and U 3 from the magnetic field. In all of the cases examined, U 1 is negligible as the modes are nearly incompressible. We thus focus on U 2 and U 3 (the latter only applicable to Case II). We further decompose U 2 into components associated with the background radial and azimuthal drifts, and the drift in the perturbed velocities. Details are given in Appendix E. In plots we normalize U i by U tot .  Note. All modes employ St = 0.1. The classic SI modes LinA and LinB of  are reproduced for comparison. The limit of zero initial toroidal and vertical fields are denoted by β f,z → ∞, while the ideal MHD limit is attained for Λ f,z → ∞ . Disks initialized with a purely horizontal field (β z = ∞ ) do not evolve the induction equation, and magnetic perturbations are assumed to be zero. These models possess a nonzero magnetic stress, α M (Equation (32)). Disks initialized with a purely vertical field (β f = ∞ ) include the full induction equation and Lorentz forces in the perturbed state.

Case I: Classic SI in Magnetically Accreting Disks
We begin by adding a horizontal magnetic field with β f = 10 to LinA and LinB. We vary the drift speeds induced by the field, which remains passive, through Λ f . Similar values of β f have been considered (e.g., Krapp et al. 2018), although these are rather strong fields compared to those typically found in numerical simulations (e.g., Cui & Bai 2021, who find a gas-tomagnetic pressure ratio of order 100). However, as emphasized in Section 2.5, it is the magnetic stress, α M , that is relevant to Case I. For nominal values of h g = 0.05, β f = 10, and Λ f = 10 −3 , we obtain α M = 0.0125, which is consistent to nonideal MHD disk simulations (e.g., Bai 2017; Béthune et al. 2017;Riols et al. 2020). Figure 1 shows that SI growth rates increase with decreasing Λ f , while oscillation frequencies become more negative and can exceed Ω in magnitude. For LinB, growth rates can increase by 50% for Λ = 10 −4 (α M ∼ 0.1) compared to the unmagnetized limit. Growth rates for other field strengths can be inferred from the equivalent α M , as shown in the Figure. For example, if β f = 100 (i.e., a factor 10 increase), then α M would decrease by a factor of 100 at fixed Λ z (see Equation (32)). Thus for Λ z  10 −4 the new growth rates would be equal to those with β f = 10 and Λ z  10 −2 , or α M  10 −3 , in which case the magnetically induced accretion has no effect.
Recall that the magnitude of the dust-gas radial drift decreases with increasing strength of the magnetic torque, while that of the azimuthal drift increases (Section 2.2). This suggests that the increasing growth rates are attributed to the azimuthal drift. We confirm this in the bottom panels of Figure 1 by plotting the contributions to the modes' pseudokinetic energy. For LinA, modes are driven by the radial drift, but its contribution drops slightly at Λ f = 10 −4 (α M ∼ 0.1), whereas that from the azimuthal drift increases.
On the other hand, for LinB, azimuthal drift has twice the contribution as radial drift even for negligible magnetic stresses, while the radial-drift contribution drops noticeably (even becoming slightly negative) for Λ f  10 −3 (α M  0.01). It is also the case for LinB that the increase in growth rates is more pronounced. This indicates that an accretion flow primarily affects SI modes in which the azimuthal drift plays a dominant role.
We conclude that in dust-poor disks SI will be moderately boosted by an underlying, magnetically induced accretion flow if the corresponding α M  10 −2 , while in dust-rich disks the SI only becomes slightly more unstable. The enhancement is limited due to the fact that with a typical h  of O(h g ), the drifts induced by the magnetic field remain small compared to those directly induced by pressure gradients. However, this picture changes when we consider regions with vanishing h , as explored next.

Case I: SI with Vanishing Pressure Gradients
Here we examine the SI with vanishing pressure gradients (h   0). These models can be considered as representing regions near to or at a pressure bump. In this section, we include a small gas viscosity and dust diffusion characterized by α SS (see Appendix B) to regularize the problem at large wavenumbers. Otherwise, we find unstable modes can develop on arbitrarily small scales with frequencies |σ| ? Ω, which are likely unphysical. Unless otherwise stated, we set α SS = 10 −9 .
In Figure 2 we plot unstable modes as a function of wavenumbers and pressure gradients for the magnetized LinA setup. We fix β f = 10 and Λ f = 10 −3 , or α M = 0.0125. For h  0.005  , we find classic SI modes with growth rates up to Ω for K x ∼ 10 3 -10 4 and K z ∼ 10 4 . The most unstable wavenumbers translate to k x η tot R ∼ 50 and k z η tot R ∼ 500 for h =  0.05 and to about 50 for h =  0.005. As shown in the pseudo-energy plots, these modes are attributed to the radial drift.
We find unstable modes change character as h  drops to 5 ×10 −4 . As shown in the right two columns of Figure 2, the instability maintains a dynamical growth rate of O(0.1Ω) with a characteristic K x ∼ 10 4 , independent of K z , except at large K z where the viscous cut-off applies. In fact, instability persists even for K z ≡ 0, indicating these modes are one-dimensional (1D) with negligible perturbations in the gas and dust vertical velocities, as well as that of the gas radial velocities (see Appendix C).
We remark that for h =´- 5 10 4 , K x = 10 4 corresponds to k x η tot R = 5. In this sense, these modes have longer radial lengthscales than the standard SI observed for h  0.005  . Most notably, however, is that as h   0, unstable modes are driven by the azimuthal drift with negligible contributions from the Similarly, in Figure 3 we show how modes in the dust-poor LinB setup evolve as the pressure gradient decreases to zero. Again we find the classic SI modes dominate for h  0.005  , here corroborated by the coincidence of the most unstable wavenumbers with the RDI resonance condition for the SI as given by Equation (40). For h =´- 5 10 4 , however, the RDI condition no longer predicts the most unstable modes, although a cluster of subdominant, classic SI modes are still found around the RDI condition. For h´- 5 10 4  , we again find the most unstable modes are completely dominated by azimuthal drift and can reach growth rates of O(0.1Ω). However, for large K x (here ∼ 10 4 ) we find the modes have large oscillation frequencies, |ω|  10 Ω, implying that t w | | 1 s  , which may violate the fluid approximation of dust (Jacquet et al. 2011).
In Figure 4 we examine how growth rates vary with β f at fixed Λ f = 10 −3 , or more precisely as a function of α M . Here, we set K z = 0 and maximize growth rates over K x . Per the simplified model in Appendix C, growth rates decrease with α M and consequently with azimuthal drift. Note that the cut-off at finite α M is due to dust diffusion, without which growth would persist for decreasing α M by going to arbitrarily small scales. For α SS = 10 −9 , the instability is quenched when α M  10 −3 for ò = 3 and when α M  10 −4 for ò = 0.2. Although the corresponding critical magnetic fields of β f ; 40 and β f ; 100 are strong, these β f values scale as (32)), and so would increase with increasing h g , decreasing Λ f , or both. Thus the instability is more easily realized in highly resistive disks. Larger α SS requires stronger azimuthal drifts for instability.
The results in this section suggests that in accreting disks the SI can still operate in regions of weak or even zero pressure gradients, provided there is sufficient magnetic stress to drive gas accretion with a corresponding α M  10 −4 -10 −3 . However, its nature differs from the classic SI: here the SI is driven by the relative azimuthal drift between dust and gas, which is largely induced by the torque acting on gas that is responsible for the underlying accretion flow.

Case II: Dust Feedback on the MRI and the SI of Alfvén Waves
We now consider disks in the limit of ideal MHD with an initially vertical field and account for Lorentz forces in the perturbed state, i.e., a live field. Recall that a purely vertical field does not modify the background drift velocities from the hydrodynamic limit, since β f → ∞ and hence F f → 0 (see Equations (14)- (17)): there is no magnetically induced accretion flow as in Case I. Here we fix h =  0.05. We begin by examining how dust loading affects the MRI in disks with ò = 0.2 and ò = 3. The ensuing MRI turbulence would stir up dust grains, so considering dust-to-gas ratios of order unity may not be self-consistent with a settled dust layer; e.g., Yang et al. (2018) find ò ∼ 0.3 under ideal MRI turbulence  (40)). Bottom: contributions to the pseudo-energy from the background radial drift (dashed), background azimuthal drift (solid), and perturbed drift velocities (dotted). These are shown for a fixed K z as indicated by the dashed-dotted line in the top panel and limited to modes with growth rates > 0.01 Ω. Left to right: decreasing total radial pressure gradients, h , with the right-most column showing results for vanishing pressure gradient. A small gas viscosity and dust diffusion with α SS = 10 −9 is included to regularize the problem at large wavenumbers. and solar metallicities (or without feedback). However, orderunity values of ò can be reached at supersolar metallicities because of feedback  or from radial concentrations by zonal flows (Dittrich et al. 2013) and so are still relevant to explore. Figure 5 shows MRI growth rates as a function of K z at fixed K x = 0 and β z = 10 4 . The maximum growth rate of 0.75 Ω is unaffected by the dust. However, with increasing ò the MRI shifts to smaller vertical scales. In a pure gas disk the most unstable MRI mode has = k C C 15 4 z s A z (Sano & Miyama 1999). However, in a dusty gas we expect  + C C 1 Az Az (see Section 3.3) and thus the most unstable µ +  k 1 z increases with dust loading. This prediction is confirmed in the figure.
Next, Figure 6 show growth rates for the LinA (ò = 3) and LinB (ò = 0.2) setups in wavenumber space and for different levels of magnetization. Starting with the leftmost column with β z = 100, the rectangular region with K x,z  20-30 corresponds to the MRI. However, the most prominent feature here is the cluster of new modes around the straight solid line for K x  10 2 (ò = 3) and K x  20 (ò = 0.2), which corresponds to the resonance condition between Alfvén waves and the dust-gas radial drift (Equation (44)), indicating they are a new class of RDI modes unique to dusty, magnetized gas.
These Alfvén wave streaming instabilities (AwSI) extend the parameter space of unstable modes to much smaller lengthscales compared to a pure gas disk, wherein only the MRI occurs (see Figure 16). However, the maximum growth rates of the AwSI modes are typically of O(10 −1 Ω), which is slower than the MRI. 4 Going rightwards to weaker magnetizations, the AwSI modes shift to smaller vertical scales, which is consistent with the resonant condition Equation (44).
In the top panels of Figure 6 we also overplot the RDI condition for the classic SI (Equation (40)) as the dotted curves. However, for β z 10 4 we do not find modes to cluster around this resonance. For example, with β z = 100, at the LinA (LinB) wavenumbers K x,z = 600 (K x,z = 120) we find growth rates 10 −4 Ω (2 × 10 −4 Ω), which is much smaller than the classic, unmagnetized SI growth rate of 0.42 Ω (0.015 Ω). This  . Growth rates of the azimuthal-drift-driven streaming instability in disks without pressure gradients (h  = 0) as a function of the plasma beta parameter, β f , at fixed Λ f = 10 −3 , which corresponds to a range of magnetic stresses, α M . Modes have no vertical structure (K z = 0) and growth rates are maximized over the radial wavenumbers, K x , at each β f (or α M ). The black (red) curve show results for a dust-rich (dust-poor) disk with ò = 3 (ò = 0.2). The solid (dashed-dotted) curves are obtained with a viscosity and diffusion coefficient of α SS = 10 −9 (10 −5 ).
indicates that the classic SI is suppressed by magnetic perturbations. It is only with an extremely weak field (here for β z 10 6 ) do we observe classic SI modes, but even then the most unstable modes occur at high K z and are attributed to the MRI and the AwSI.

Case II: MRI and SI in Resistive Disks
We now add a constant resistivity, η O , to Case II to mimic a dead zone in the disk midplane (Gammie 1996a). As before, we first examine the effect of dust loading on the MRI. Here, the resulting MRI turbulence would be dampened by ohmic resistivity. In this case, while grains do not settle much more than in ideal MHD (Fromang & Papaloizou 2006;Yang et al. 2018), weakened horizontal diffusion coupled with dust feedback can drive strong radial dust concentrations that lead to order-unity (or larger) values of ò in the dead zone . 5 This motivates our consideration of high dust-togas ratios. Figure 7 shows that increasing ò decreases both the MRI growth rates and the corresponding K z , so modes develop on larger vertical scales in dust-rich disks. This can again be explained by the fact that the effective Alfvén speed decreases with ò as the total density of the system increases. In resistive, gaseous disks with Λ z = 1, however, the most unstable MRI Figure 5. Ideal magneto-rotational instability growth rates in dusty disks threaded by a vertical magnetic field with strength β z = 10 4 and dust-to-gas ratio ò = 3 (blue) and ò = 0.2 (brown). The grain size is fixed to St = 0.1. The most unstable vertical wavenumber in the pure gas limit is marked by the solid vertical line. Assuming dust loading reduces the Alfvén speed by increasing the total density of the system, the resulting most unstable K z are marked by the dashed and dashed-dotted vertical lines. Figure 6. Growth rates of unstable modes in a dusty, actively magnetized disk initialized with a vertical field in the limit of ideal magneto-hydrodynamics. The grain size is fixed to St = 0.1. From left to right: decreasing magnetic field strength (increasing plasma beta parameter β z ). Upper row: dust-rich disks with ò = 3. Bottom row: dust-poor disks with ò = 0.2. Solid lines correspond to the resonant drag instability (RDI) condition between Alfvén waves and the dust-gas radial drift (Equation (44)). Dotted lines correspond to the RDI condition for the classic streaming instability (Equation (40)). (Sano & Miyama 1999). Again, in a dusty gas C Az → C Az,eff (see Equation (45)) so now k z ∝ (1 + ò) −1/2 while s ∝ (1 + ò) −1 . These reductions are indeed observed for ò = 3 compared to ò = 0.2 in the figure. Figure 8 shows growth rates with fixed β z = 10 4 as a function of K x,z and Λ z . We find that AwSI modes are easily damped by resistivity: even with Λ z = 10 there is only a hint of their presence around K x,z ∼ 10 2 as the "flared" region to the top-right of the MRI modes (this is more pronounced for ò = 0.2). Note that the RDI condition for the AwSI modes (solid lines) appears to have no relevance to the features in the figure; although coincidentally for ò = 0.2 and Λ z = 0.1 they mark the maximum K z beyond which classic SI modes are strongly stabilized.

mode has
In the ò = 3 disk the classic SI appears for Λ z 1 and its properties are similar in this regime. However, in the ò = 0.2 disk the SI only appears for Λ z 0.1, indicating that slowly growing SI modes are more easily stabilized by the magnetic field. For Λ z = 0.1, notice also the SI modes are damped for K z  10 2 , as mentioned above.
Thus in resistive, dust-poor disks, there is a lower bound to the scales at which the SI operates, here being  0.01 H g in both radial and vertical directions. This is unlike the Λ z = 10 −2 disk, which is effectively unmagnetized, where the SI can operate on much smaller vertical scales as growth rates increase and plateau at large K z .
To confirm the magnetic stabilization of the SI, we show in Figure 9 the pseudo-energy decomposition for the SI modes satisfying the RDI condition for ò = 0.2 and Λ z = 0.01, 0.1. Note that the resonant K z → ∞ as K x → 1/|ζ x | (∼100 here) according to Equation (40), as expected drag forces (U 2 ) are destabilizing, but magnetic forces are stabilizing (for ease of comparison we show − U 3 ). The latter effect intensifies with Figure 7. Resistive magneto-rotational instability growth rates in dusty disks threaded by a vertical magnetic field with strength β z = 10 4 and dust-to-gas ratio ò = 3 (blue) and ò = 0.2 (brown). The grain size is fixed to St = 0.1. The most unstable vertical wavenumber in the pure gas limit is marked by the solid vertical line. Assuming dust loading reduces the Alfvén speed by increasing the total density of the system, the dashed and dashed-dotted vertical lines mark the resulting most unstable K z , and the red open circles mark the corresponding maximum growth rates. increasing wavenumber, but for Λ z = 0.01 it is always subdominant, so the SI persists. However, for Λ z = 0.1, both magnetic and drag contributions increase rapidly with K x,z , here resulting in a near cancellation and the SI being effectively quenched. The toy model presented in Appendix D also indicates magnetic stabilization at low dust-to-gas ratios for resonant modes with K z ? K x .

Direct Simulations
To verify the new instabilities uncovered above-namely, the azimuthal-drift-driven SI without pressure gradients (Section 4.2) and the Alfvén wave SI (Section 4.3)-we also solve the full shearing box equations directly. To this end, we have developed a finite-difference code to evolve Equations (22)-(26) in their conservative form. We approximate spatial derivatives with a sixth-order central finitedifference scheme and integrate in time with a fourth-order Runge-Kutta method. We follow the ATHENA code to treat the source terms related to tidal and Coriolis forces, the large-scale radial pressure gradient, and torques due to a horizontal magnetic field if applicable (Stone & Gardiner 2010). We integrate the dust-gas drag term explicitly.
In Appendix F we test the code by reproducing the standard LinA and LinB SI growth rates, as well as MRI growth rates with and without resistivity. We have found this simple code to be adequate for our primary goal of confirming linear theory.

Numerical Setup
We perform axisymmetric simulations in the (x, z) plane but include all three components of the velocity and magnetic fields (the latter for Case II). The domain x ä [ − L x /2, L x /2] and z ä [ − L z /2, L z /2] is discretized into N x = 100 and N z = 100 uniform cells, respectively. To test an eigenmode with wavenumbers k x and k z , we set the domain size to be one wavelength in each direction, i.e., L x,z = 2π/k x,z . We adopt L x as a unit of length and Ω −1 as the unit of time. The velocity normalization is therefore L x Ω ≡ u norm . The mass scale is arbitrary for a non-self-gravitating disk, but, for convenience, we take ρ g = 1 ≡ ρ norm in the equilibrium.
We initialize each simulation in a steady state consisting of constant densities and drift velocities, as discussed in Section 2.3. Following Youdin & Johansen (2007), we add a pair eigenmodes with the same k x but oppositely signed k z as a perturbation. The eigenmode amplitude is scaled such that δρ d = 0.01ρ norm . We apply periodic boundary conditions in x and z.

Azimuthal-drift-driven SI without Pressure Gradients
Here, we test for the SI driven by the azimuthal drift induced by magnetic torques from a horizontal field. We adopt the setup of Case I with β f = 10 and Λ f = 10 −3 , or equivalently α M = 0.0125. We set h  = 0 so the equilibrium drift, which is dominated by the azimuthal drift, is entirely due to magnetic torques. For the perturbation we choose K x = 5000 and K z = 100, which results in a "tall" shearing box as our domain. We consider ò = 3 and 0.2, for which the linearly unstable eigenmodes are labeled "LinAIeta0" and "LinBIeta0" in Table 1. For this test we do not evolve the magnetic field, assuming the resistivity is sufficiently large such that the field remains passive. Figures 10 and 11 show the evolution in the perturbed dust density. The upper panel shows the maximum value measured in the domain, while the lower panel shows w z at a fixed grid point. The measured growth rates 0.1360 Ω for ò = 3 and 0.1380 Ω for ò = 0.2 are close to that calculated from linear theory (s = 0.1358 Ω and s = 0.1377 Ω, respectively). Furthermore, for ò = 3 we measure an oscillation period of T ; 10 Ω −1 , which is close to that expected from linear theory as 2π/|ω| = 10.233 Ω −1 . Similarly, for ò = 0.2 we measure T ; 2.6 Ω −1 , compared to the theoretical period of 2.56 Ω −1 . This shows that in dust-poor disks the azimuthal-drift-driven SI is highly oscillatory. We find growth rates begin to deviate from linear theory after a few periods. Nevertheless, the close agreement in the early phase of the simulations confirms that the SI can operate without pressure gradients if the gas is subject to passive magnetic torques.
We note that for large K x the azimuthal-drift-driven SI is insensitive to K z (see Figures 2-3) and in fact persists with Figure 9. Pseudo-energy decomposition of the classic streaming instability (SI) in a magnetized disk with β z = 10 4 , ò = 0.2, and St = 0.1. Left: Λ z = 0.01. Right: Λ z = 0.1. Blue: the destabilizing contribution from dust-gas drag. Red: stabilizing contribution from magnetic perturbations (the negative of the magnetic pseudo-energy is shown). Here, we consider modes with K x and K z satisfying the resonant drag instability condition for the SI, Equation (40), i.e., along the dotted curve in the corresponding panels of Figure 8. K z = 0 (Appendix C). Indeed, we have also performed 1D simulations in x only and recovered similar growth rates as above.

Alfvén Wave SI
Here we adopt the setup of Case II with a live magnetic field to test for the AwSI, which results from a resonance between Alfvén waves and the radial dust drift. The induction equation and Lorentz forces are evolved. We consider ideal MHD with an initially vertical field with β z = 100 and revert to the nominal pressure gradient h =  0.05. The eigenmodes are labeled as "LinAII" and "LinBII" in Table 1 for ò = 3 and 0.2, respectively. We fix K x = 1000 but choose K z = 25 for LinAII and K z = 83 for LinBII in accordance with the RDI condition given by Equation (44), i.e., these are resonant wavenumbers. Note that our domain size of one mode wavelength is sufficiently small to exclude the MRI, which operates on larger scales. However, the MRI has the larger growth rate (see Figure 6), which in reality would dominate and drive rapid disk evolution. The examples here are chosen to confirm the Alfvén wave SI. Figures 12 and 13 show the evolution in the perturbed dust density for the above cases. For LinAII the measured growth rate and oscillation period are s = 0.1250 Ω and T = 12.8 Ω −1 , which compare well with the theoretical values of 0.1248 Ω and 12.79 Ω −1 . Similarly, for LinBII we measure s = 0.228 Ω and T = 0.94 Ω −1 , compared with 0.2215 Ω and 0.9353 Ω −1 from linear theory.

Planetesimal Formation in Accreting Pressure Bumps
We have shown that the SI persists in regions of weak or vanishing pressure gradients, provided that there is sufficient azimuthal drift between dust and gas. This suggests that the SI can remain effective inside pressure bumps. Although we considered the case of an azimuthal drift resulting from largescale, horizontal magnetic torques acting on the gas, our results are not limited to this context. The constant azimuthal force, F f , which represents the background magnetic torque in our models, may also represent other causes of radial gas accretion. For example, as argued by McNally et al. (2017), if a vertical angular momentum loss (e.g., due to disk winds) resulted in a gas accretion equivalent to that induced by a corresponding F f , then our results still apply.
Using Equations (37) and (34), we can express the growth rate of the azimuthal-drift-driven SI at a pressure extremum in terms of the laminar, horizontal stress, α M , where we assumed St = 1. We can use this equation to assess the relevance of the SI around pressure bumps, for example in those found by Riols et al. (2020) in global nonideal MHD simulations of dusty PPDs. Recall from Equation (32) that 1 when other parameters are fixed. In their nonideal MHD simulations, Riols et al. find that gas accretion is largely due to angular momentum removed vertically by a magnetized wind (and therefore occurs near the disk surface), which is a few times larger than that transported radially. The latter is dominated by horizontal magnetic fields with dimensionless stresses of order 10 −2 (see also Béthune et al. 2017). Hence, we consider α M ∼ 10 −2 below. Riols et al. (2020) also find the spontaneous formation of pressure bumps due to a wind instability (Riols & Lesur 2019). In one of their disk models, for a pressure bump at ∼20 au they find that 3 mm sized dust grains (with St ∼ 0.04) accumulates into a ring of width ∼ 0.6 H g and settles to a thin layer of thickness 0.25 H g . In the ring, the vertically integrated dust-togas ratio (or metallicity, Z) of all four dust species they considered is about 0.3. We therefore take ò ∼ 0.3. We then find~Ws K 10 x 3 for h g = 0.05, as used in their disk models. The ring width limits K x  10, which implies a growth timescale 50 orbits.
Of course, the relative importance between the azimuthaldrift-driven SI and the classic, radial-drift-driven SI depends on the local pressure gradient, h . It is true that the classic SI formally ceases only where h  = 0 exactly; elsewhere in the pressure bump it can still operate, though on vanishing scales as one approaches the bump center. However, as discussed in Section 2.6, the torque-induced azimuthal drift is expected to become important when h a  h St M g  , which is consistent with the numerical results in Section 4.2. We therefore expect a finite region in which the azimuthal-drift-driven SI operates.
To make a crude estimate, we consider a Gaussian pressure bump in the disk midplane, µ -- 2 2 , where a is the ring width. Neglecting magnetic pressure, the magnitude of the reduced pressure gradient near the bump center is where we took a = H g in the second equality. For the above parameters we find azimuthal drift is significant within |R − R 0 |  0.025 H g . Hence the longest permissible radial wavelength for the azimuthal-drift-driven SI is 0.05 H g , or K x  10 2 , for which the growth timescale is 16 orbits. We can therefore expect rapid instability even close to the bump center. We emphasize that the SI discussed above requires azimuthal drift. However, it is possible to have a pressure bump without such drift, e.g., those in geostrophic balance, which have been considered in planetesimal formation simulations (e.g., Carrera et al. 2021aCarrera et al. , 2021b. In this case, both gas and dust orbit at the Keplerian velocity at the bump center, for which we expect neither the classic nor the azimuthal-drift-driven SI.

Comparison with Secular Gravitational Instabilities
The secular gravitational instability (SGI; Youdin 2011; Takahashi & Inutsuka 2014;Latter & Rosca 2017) is another dust-clumping mechanism (Tominaga et al. 2018;Pierens 2021) that may facilitate planetesimal formation (Abod et al. 2019). Here, dust-gas friction provides an effective "cooling" (Lin & Youdin 2017) that removes pressure support against gravitational instability (Lin & Kratter 2016). 6 The SGI does not require large-scale pressure gradients, which also makes it a candidate for planetesimal formation at pressure bumps. It is thus of interest to compare the SGI with the azimuthal-drift-driven SI discussed above.
To this end, in Figure 14 we calculate SGI growth rates using Equation (13) of Takahashi & Inutsuka (2014) and compare it with that of the azimuthal-drift-driven SI. We consider two disk masses corresponding to Toomre parameters, Q T = 10 and 100. A dimensionless dust diffusion coefficient of α SS = 10 −9 is used without a corresponding gas viscosity as it is neglected in Takahashi & Inutsuka's Equation (13).
In the Q T = 10 disk, the SGI dominates on most scales except the largest ( ∼ H g ), which are stabilized by rotation. For Q T = 100, the SGI still has the largest overall growth rate, but only dominates for K x  60. These results suggest that the azimuthal-drift-driven SI is likely more relevant in low-mass disks and at large scales compared to the SGI. When a larger dust diffusion coefficient of α SS = 10 −6 is adopted (not shown), we find the instabilities have comparable maximum growth rates ( ∼ 10 −2 Ω), but the SI dominates over the SGI for K x  100, while both are cut-off for K x  200-300.
Note that the SI model does not include disk self-gravity and the SGI model does not include magnetically driven azimuthal drifts. A proper comparison between these instabilities in a unified framework is beyond the scope of this paper but should be conducted in the future.

SI in Off-midplane Accreting Layers
Another region of vanishing pressure gradients is off of the disk midplane. To see this, consider a standard Gaussian atmosphere; then, the 3D pressure profile is The corresponding dimensionless pressure gradient is   Figure 14. Comparison between the secular gravitational instability (SGI, red) and the azimuthal-drift-driven streaming instability (black) as a function of radial wavenumbers, K x . SGI growth rates are computed from Equation (13) of Takahashi & Inutsuka (2014) with Toomre parameters Q T = 10 (solid) and Q T = 100 (dashed-dotted). The dust-to-gas ratio and Stokes number are fixed to ò = 0.2 and St = 0.1, respectively. Neither disk models have radial pressure gradients (h  = 0). Modes have no vertical structure, K z = 0.
If h >  0, i.e., a negative midplane pressure gradient, then at some height η 3D will vanish since H g increases with R. For H g ∝ R and h =  0.05, we find η 3D → 0 as  | | z H 2 g . Of course, the dust-to-gas ratio will be small at this altitude due to dust settling, especially for large grains. We may thus only expect small grains there. However, Equation (46) shows that growth rates are not strongly dependent on these parameters. Indeed, even for ò = 0.01 and St = 0.01 we find W s K 10 x 4 (assuming α M = 0.01 and h g = 0.05), so for K x ∼ 100 the growth timescale ∼160 orbits is still short compared to disk lifetimes. However, a small dust-to-gas ratio may only lead to order-unity turbulent fluctuations rather than dust clumping .
Note that for z ≠ 0 complications can arise from the dustsettling instability (DSI): an RDI associated with a resonance between inertial waves and the dust's vertical settling motion (Squire & Hopkins 2018b;Krapp et al. 2020). However, the DSI requires K z ≠ 0, while the azimuthal-drift-driven SI can operate for K z = 0. We thus expect vertically extended modes to be unaffected by settling.
We remark that the SGI is probably unimportant in the disk layers away from the midplane due to low values of ò and high values of Q T locally (Takahashi & Inutsuka 2014).

Dust Feedback and the Ideal MRI
In a realistically stratified disk and considering ideal MHD, MRI modes can only develop if they fit into the gas disk thickness, 2H g . Now, the characteristic MRI vertical wavenumber is k z ∼ Ω/C Az,eff , where the effective Alfvén speed C Az,eff = (1 + ò) −1/2 C Az decreases with dust loading since the total density increases. Then for a given field strength, β z , in the dust-free limit, we require k z > π/H g , or b + ( )  10 1 50 z  upon dust loading for the ideal MRI to operate. Thus if the MRI operates without dust, it will continue to operate in the presence of dust, although it will shift to smaller vertical scales. While the maximum growth rate remains at 0.75 Ω, growth rates for modes with k z  Ω/C Az are somewhat reduced by dust loading, as shown in Figure 5.
Recently, Yang et al. (2018) carried out stratified MHD simulations including dust and found that the particles sediment to a thinner layer when the metallicity, Z, is increased. They attributed this to the reduction in the effective sound-speed of the dusty gas (Lin & Youdin 2017). We point out another possible contributing factor from dust feedback onto the MRI.
For ideal MHD simulations without particle feedback, Yang et al. (2018) find a turbulence strength of α SS ∼ 10 −2 . Using the diffusion model of Dubrulle et al. (1995) . Then, for β z = 10 4 we find λ z,opt ; 0.06 H g = 2 H d . This means that, upon the introduction of dust feedback, MRI modes with vertical lengthscales comparable to the dust layer thickness (i.e., those having k z  Ω/C Az ; see Figure 5) will be damped, which may promote further settling.
In practice, however, feedback is only important if the local ò is O(1) or larger, so this effect should only be appreciable at high Z. Indeed, Yang et al. (2018) find H d is only reduced by a factor of 1.5 for an eight-fold increase in Z to 0.08 from its canonical value of 0.01.

Dust Feedback and the Resistive MRI
In resistive disks the characteristic MRI vertical wavelength is k z ∼ C Az,eff /η O . Thus for fixed η O , defined through Λ z in the dust-free limit, requiring the MRI to fit inside the gas disk yields This condition becomes more difficult to fulfill as ò increases, which reflects the increasing vertical length scale. Furthermore, growth rates are damped on all length scales by dust feedback (Figure 7). The result above suggests the possibility of self-sustained dust settling as follows. Consider a dusty disk undergoing resistive MRI turbulence but without dust feedback. For reference, in such a simulation Yang et al. (2018) measured a H d = 0.2 H g for St = 0.1 and β z ∼ 10 4 in quasi-steady state as settling is balanced by turbulent stirring. Switching on feedback weakens the MRI, so the dust grains begin to settle, which further increases ò and weakens the MRI. This positive feedback loop cannot continue indefinitely, however, as the gaseous layers above and below the dusty midplane can still undergo full (resistive) MRI turbulence and stir the particle layer (Fromang & Papaloizou 2006). Nevertheless, we can expect thinner dust layers in higher-metallicity disks because the MRI is progressively suppressed. This scenario is similar to that outlined for the vertical shear instability in the presence of dust (Lin 2019).
A shortcoming of the above discussion (including Section 6.4) is that the estimates given in Equations (50) and (51) assume a uniform dust-to-gas ratio. In the stratified context, ò could be taken as the average value over the vertical column, in which case ò ∼ Z = 1 and thus feedback should always be negligible in an averaged sense.
On the other hand, we can apply a similar argument for fitting MRI modes within the dust layer, wherein ò may reach order unity or larger. Doing so, we find the critical implying it is difficult to have MRI modes confined to a highly settled dust layer. This does not prevent MRI modes from operating in the gas outside the dust layer, however. Stratified analyses are required to study MRI modes with vertical wavelengths that exceed the dust layer thickness yet still fit within the gas disk.

Magnetic Stabilization of the Classic SI
In Sections 4.3-4.4 we found that a live vertical magnetic field can suppress the classic SI, which is a stronger statement than MRI modes outgrowing the SI. For ideal MHD, there is no trace of the classic SI for β z 10 4 . While the classic SI does appear at β z = 10 6 , the most unstable modes are mixed with the MRI. It is only with β z = 10 8 , i.e., an extremely weak field, and with ò = 3, do we find distinct classic SI modes that dominate the system. This observation, together with the fact that we considered rather large grains with St = 0.1, which should favor the SI, suggests that the classic SI is sensitive to magnetic stabilization.
As a consequence, for typical values of β z = 10 4 adopted for PPDs, the classic SI is only present with sufficient resistivity. For ò = 0.2 we only partially recover the SI with Λ z 0.1 (high vertical wavenumbers are still damped), and Λ z = 0.01 is needed for SI modes to outgrow the MRI. At high dust-to-gas ratios, magnetic stabilization is less effective as we find the SI dominates over the MRI with Λ z 1 when ò = 3. However, should such a high ò be attained through dust settling, then (MRI) turbulence cannot be too strong in the first place, which sets an upper bound on Λ z . Therefore, whether or not the SI can occur in magnetized disks is still determined by the value of Λ z below which either the MRI wavelengths exceed the disk scale height (and therefore does not operate), or when SI growth rates exceed the resistive MRI.

Dust Clumping in MHD-turbulent Disks
At present, the connection between the linear SI and nonlinear clumping-and hence planetesimal formation-is still unclear . Nevertheless, our results suggest that dust clumping directly via the SI may be more difficult in magnetized disks due to magnetic stabilization. This is in addition to stirring by MRI turbulence, which may prevent dust sedimentation.
On the other hand, planetesimal formation has been successfully simulated in MRI-turbulent disks (Johansen et al. , 2011, but this is not necessarily in contradiction with our results. This is because of the formation of zonal flows or pressure bumps in MRI-turbulent disks (Johansen et al. 2009a) that can trap dust (Dittrich et al. 2013;Xu & Bai 2022). Similar radial dust concentrations can develop in models with dead zones owing to the weak radial diffusion of particles . In either case, the enhanced dust-to-gas ratio in these dust traps should, according to our results, reduce the effect of the MRI and favor further concentration and eventually allow the classic SI to operate, though perhaps not at the exact center of pressure bumps.

Relevance of the Alfvén Wave SI in PPDs
Our numerical results indicate that the Alfvén wave SI requires nearly ideal MHD conditions to operate. This makes their relevance to PPDs questionable, since if ideal MHD conditions are met, for example in the inner disk (Flock et al. 2016), the MRI would generate vigorous turbulence and overwhelm the Alfvén wave SI.
One possible exception is in the disk atmosphere where β z increases relative to the midplane as the gas density drops. If sufficiently magnetized, i.e., if β z  10 (Equation (50)), then the MRI is stabilized, leaving AwSI modes to survive. As an example, suppose β z ∼ 10 4 at z = 0, then at |z| = 4 H g we find β z ∼ 3 for a Gaussian atmosphere. Using Equation (49), at this height the reduced pressure gradient parameter is , which is much larger in magnitude than the midplane value. For this estimate we assumed H g = h g R with h g = 0.05.
In Figure 15 we plot growth rates for the above example under ideal MHD conditions with ò = 0.01 and St = 0.1. As expected, the MRI is excluded from geometric considerations, but the atmosphere is unstable to the Alfvén wave SI with maximum growth rates ;0.03 Ω, which may lead to turbulent activity. However, the DSI also operates off of the midplane (Squire & Hopkins 2018b; Krapp et al. 2020), but is neglected here. How the DSI interacts with the Alfvén wave SI, or how it is affected by magnetic fields, is beyond the scope of this work.

Caveats and Outlook
In the shearing box formalism, the effect of the radial pressure gradient from the global disk is treated through the constant parameter, h . We therefore cannot model pressure bumps as a whole. The persistence of the SI with h   0 only indicates instability localized to regions of vanishing pressure gradients, for example near the center of a pressure bump. A proper stability analysis of pressure bumps requires a variable h h = ( )   x . For consistency with global disks, one may also need a variable magnetic torque, F f = F f (x). These generalizations inevitably result in a global eigenvalue problem, i.e., ordinary differential equations in x. This should be tackled in a future study.
Our results demonstrate the potential importance of the azimuthal drift between dust and gas, especially in disks torqued by a laminar magnetic field. According to RDI theory, nonaxisymmetric disturbances may grow from a resonance between the azimuthal drift and neutral waves in the gas, for example, Rossby waves (Pan & Yu 2020). This hypothetical RDI is excluded in our axisymmetric models but should be investigated and compared with the axisymmetric, azimuthaldrift-driven SI discussed in this work (which is not an RDI). However, nonaxisymmetric modes may not grow exponentially in the shearing box due to differential rotation (Johnson & Gammie 2005), so a radially global or cylindrical treatment would be more appropriate.
We have only considered ohmic resistivity as the sole nonideal MHD effect. For passive fields, this is not a limitation since only the effective body force, F R,f , acting on the gas is relevant, regardless of how it arises. For live fields, however, one should also consider ambipolar diffusion and the Hall effect, which may, in fact, dominate the midplane of PPDs (see Lesur 2021, and references therein). The Hall effect leads to qualitatively new phenomena such as whistler waves and the Hall-shear instability (Kunz 2008); while ambipolar diffusion can make unstable modes with k x ≠ 0 dominant (Kunz & Balbus 2004). The present work should be generalized to study how dust interacts with these modes and instabilities.
Finally, our results are based on the fluid approximation of dust grains. However, in some calculations, we find modes with high frequency or growth rates such that |στ s |  O(1), which may put the fluid treatment into question (Jacquet et al. 2011). Ultimately, our results need to be checked against particle-gas simulations (e.g., .

Summary and Conclusions
In this paper, we study the stability of dusty, magnetized PPDs. We are mainly motivated by recent models of PPDs that highlight the dominant role of large-scale magnetic fields in driving disk accretion (e.g., Gressel et al. 2015;Bai 2017;Béthune et al. 2017) and how this would affect planetesimal formation through the SI (Youdin & Goodman 2005;. We apply standard linear-stability analyses and verify key results with direct MHD simulations including dust. We extend previous local shearing box models of the SI to account for large-scale magnetic stresses from the global disk, which is treated as a constant, passive torque that modifies the equilibrium dust and gas drift velocities. Here, we find the magnetic torque primarily induces an azimuthal drift between dust and gas, which becomes important in regions where the global radial pressure gradient, h , is small. This results in instability even when h   0 (in which case the classic SI of Youdin & Goodman ceases to operate). The existence of azimuthal-drift-driven SIs suggests that pressure bumps in accreting disks remain viable sites for accelerated planetesimal formation.
We also consider how a live magnetic field interacts dynamically with dust. Under ideal MHD we find the radial drift between dust and gas can destabilize Alfvén waves at sufficiently high wavenumbers. However, these modes are likely overwhelmed by the MRI and therefore may be of limited relevance to realistic PPDs, unless the MRI becomes subdominant, for example in strongly magnetized disks. In resistive disks, we find dust feedback can stabilize the MRI by reducing the effective Alfvén speed of a dusty gas compared to pure gas. Conversely, we find the classic SI can be stabilized by magnetic perturbations, especially at low dust-to-gas ratios. Whether clumping observed in simulations of dusty, magnetized gas can be directly attributed to the SI is, therefore, an open question.
In a follow-up work, we will use numerical simulations to explore the nonlinear evolution of the azimuthal-drift-driven SI, as well as the classic SI and MRI in dusty, magnetized disks.
We are grateful to Xuening Bai and Pinghui Huang for helpful advice on nonlinear simulations. We thank the anonymous referee for an insightful report that prompted us to compare our results with the SGI. This research is supported by the Ministry of Science and Technology of Taiwan

Appendix A Global Steady State
We derive the global pressure profile in a steady-state disk with velocity fields given by Equations (14)-(17) at each R. We assume St and ò are both constants. Then the condition for a constant mass flux in dust and gas is equivalent. From Equation (14), the gas mass flux is From Equations (6), (7), and (12), we see that F f ∝ R −3/2 /Rρ g , so the second term on the right-hand side is a constant. That is, the magnetic torque induces a constant mass flux (see also Equation (20)). We therefore require the first term to be constant. This implies, using Equation (18) with (12), that are dimensionless constants and subscript "0" denotes evaluation at R 0 . Using the definition of η (Equation (13)), Equation (A2) is an ordinary differential equation for P, This can be integrated with the boundary condition P → 0 as R → ∞ to give The corresponding gas density profile can be obtained through r = P C s g 2 for an isothermal disk. Note that we require h >  0 to ensure P > 0 at all radii. This is satisfied in a thin, weakly magnetized disk wherein η 0 is ( ) O h g0 2 and β f ? 1.

Appendix B Linearized Equations
We linearize Equations (22)-(26) assuming a uniform, vertical background field, B z (Case II). We define . The result is Az z x x z x 2 s s 2 Figure 15. Unstable modes in a strongly magnetized disk atmosphere with β z = 3 and grains with St = 0.1 and ò = 0.01. Modes with large wavenumbers along the solid line correspond to Alfvén wave streaming instabilities. The horizontal dashed-dotted line (K z = π) is the minimum wavenumber for modes to fit within a pressure scale height, which excludes the block of magnetorotational instability modes (in red) in the lower left.
hydrodynamic model with the magnetic field only modifying the background drift speeds. First, we note that this new instability can operate even when K z = 0 (contrary to the classic SI, which requires K z ≠ 0; Youdin & Goodman 2005). Henceforth we consider K z = 0. Second, as compressibility is negligible in all of the modes examined, we replace the gas continuity equation with the incompressible condition, ∇ · v = 0. Then the linear perturbations satisfy k x δv x + k z δv z = 0. Hence for k z = 0 the gas has no radial motion, δv x = 0. Incompressibility also means letting the gas density perturbations W → 0 but C s → ∞ , such that the enthalpy perturbation d º h C W s 2 remains finite and is determined from the gas' radial momentum equation. Furthermore, the linearized vertical momentum equations for gas and dust can be satisfied with δv z = δw z = 0. We include dust diffusion but ignore gas viscosity.
The linear eigenvalue problem now only involves the gas' azimuthal momentum equation and the dust's continuity and horizontal momentum equations. In the limit considered, these can be written as  Figure 16. Left: the standard magneto-rotational instability in a pure gas disk with β z = 10 4 in the ideal magneto-hydrodynamic limit. Middle: classic streaming instability (SI) in an unmagnetized disk with ò = 3 and St = 0.1. Right: similar to the middle panel but with ò = 0.2. In the SI panels the dotted curves correspond to the resonant drag instability condition given by Equation (40).
and recall ζ y is the dimensionless azimuthal drift (Equation (34)). The full dispersion relation is a quartic in σ and its direct solutions are unwieldy. This quadratic equation for στ s + μ d , and hence σ, can be solved readily.

C.1. Explicit Solutions without Dust Diffusion
When dust diffusion is negligible (D = 0), a further simplification can be made in the limit St 2 = |μ d |, which means K x cannot be too small. The term ∝ St 2 in Equation (C7) can then be neglected. Seeking growing solutions, we find The growth rate and oscillation frequency are then (Recall that we defined σ ≡ s − iω.) We see from Equation (C9) that instability is directly powered by azimuthal drift. Furthermore, for azimuthal drifts driven entirely by the magnetic torque, we have ζ y ∝ St (Equation (34) with h  = 0), thus µ s St . Note that both s and ω are unbounded with increasing K x , but this limit eventually violates the assumption of small |μ d | and |μ g | used to derive Equations (C9)-(C10). In any case, the fluid approximation probably breaks down when |στ s |  1 (Jacquet et al. 2011).
In Figure 17 we test the above model by comparing it with the solution to the full eigenvalue problem. We compute unstable modes with K z = 0 in disks with h  = 0, β f = 10, and Λ f = 10 −3 (or α M = 0.0125), for a range of stopping times. The analytic model (Equation (C7)) accurately reproduces the growth rates, while the fully analytic solutions (Equation (C9)) work best for small St, large K x , or both.

C.2. Physical Interpretation
The origin of the instability can be examined through the linearized Equations (C1)-(C4). Under the same approximations as that used to derive Equation ( The instability mechanism can now be read from Equations (C11), (C12), and (C13). Suppose the dust experiences an azimuthal acceleration, moves outward and is slowed down by gas drag (Equation (C11)). Then the local dust density increases (Equation (C12)), which can be seen by inserting Equation (C8) to give Q ∝ (1 − i)δw x with a positive proportionality constant. For ζ y > 0 this increases the feedback onto the gas that accelerates it in the azimuthal direction as Figure 17. Streaming instability in disks without pressure gradients (h  = 0) but including a magnetically induced accretion flow with a dimensionless Maxwell stress of α M = 0.0125. No viscosity or diffusion is used (ν = D = 0). Solid curves are computed from the full eigenvalue problem, the circles from the analytic model (Equation (C7)), and the dashed curve is the fully analytic solution (Equation (C9)). The dust-to-gas ratio is ò = 3. Three Stokes numbers are shown: St = 0.1 (black), St = 10 −2 (red), and St = 10 −3 (blue). Modes have no vertical structure, K z = 0. δw y ∝ (1 + i)Q, but for tightly coupled dust the latter is dragged along (Equation (C13)). This leads to a positive feedback and hence instability.

C.3. Effect of Dust Diffusion
We briefly examine the influence of dust diffusion. Figure 18 shows growth rates for α SS = 10 −9 , 10 −8 , and 10 −7 . As expected, diffusion reduces growth rates and sets a minimum scale of instability, which increases with α SS . The analytic model (circles) reproduces results from the full eigenvalue problem with slight deviations at large K x . This is not surprising since the model assumes sufficiently small m µ | | K x d .

Appendix D Single-fluid Model of a Magnetized, Dusty Gas
Our numerical results show the disappearance of classic SI modes with large k z as one increases Λ z , even when Λ z = 1 (Section 4.4, e.g., the bottom left two panels in Figure 8). This is curious because, for large k z , ohmic diffusion should render magnetic fields ineffective. To check this result, we analyze a one-fluid model of a dusty, magnetized gas following the framework developed by Laibe & Price (2014). We assume the gas is incompressible (∇ · v = 0) and ρ g is constant. We neglect gas viscosity and dust diffusion.
In the single-fluid approach, one works with the total density where f g = 1/(1 + ò) is the gas fraction and f d = ò/(1 + ò) is the dust fraction. We define the dust-gas drift as Δu ≡ w − v. Then v = u − f d Δu and w = u + f g Δu.
For sufficiently small grains, dust follows the gas with a correction from the differential force between dust and gas. In our case, this arises from pressure gradients and Lorentz forces that only act on the gas. ). While the TVA simplifies the modeling of dusty gas considerably, its shortcomings were recently analyzed by Pan (2021), where it was shown that the TVA underestimates linear growth rates when ò = 1, k z ? k x , or both, and is particularly significant at resonant wavenumbers (i.e., those satisfying Equation (40))-regimes that we will consider below. The remedy to the TVA's deficiency involves additional contributions to Δu; see Pan (2021) for details. We leave this to future work and proceed with the caution that the following discussion is aimed at capturing qualitative effects, rather than producing a quantitative replacement of the twofluid treatment.
The gas incompressibility condition, dust continuity equation, and the center-of-mass momentum equation for the dust-gas mixture are See also Lin (2021, Appendix B) for the case of a compressible, unmagnetized gas. Strictly speaking, one should express the induction equation in terms of u and Δu. We neglect this complication and set v → u in Equation (24), so the induction equation retains the form as in the two-fluid model; see Fromang & Papaloizou (2006) for a similar treatment. This approximation is applicable when f d = 1, but to connect to

( ) Code Tests
To test our finite-difference code used in Section 5, we first reproduce the classic SI eigenmodes LinA and LinB, as described in , and summarized in Table 1. For this test the magnetic field is switched off. Figure 20 shows the evolution of the maximum dust density perturbation in the domain, which shows a good agreement between the growth measured in the simulation and the expected growth rate calculated from linear theory. We remark that for LinB the measured growth rate in δρ d is slightly higher than the theoretical value, but for δw z we do measure the same growth rate of s = 0.0154 Ω as in linear theory.
Similarly, we reproduce the standard magneto-rotational instability with a purely vertical field in Figure 21. For this test we disable the dust component and set h  = 0. We fix β z = 100 , K x = 0, and choose the most unstable = K 5 15 2 z ( = K 3 2 z ) in the ideal (resistive cases) according to linear theory. For ideal MHD we obtain a linear growth rate of 0.75 Ω as expected. For the nonideal case, we set Λ z = 0.1 and obtain a growth rate of 0.074 Ω, which is close to the analytic value of 0.75Λ z Ω in the limit Λ z = 1 (Sano & Miyama 1999).