© The Authors 2026

1 Introduction

The first stage of planet formation involves the coagulation of micron-sized grains into millimeter-sized dust (Dullemond & Dominik 2005). However, the bouncing (Zsom et al. 2010) and fragmentation (Blum & Wurm 2008) barriers tend to prevent the further growth of solids, such that the way in which particles grow beyond the millimeter-size range remains an unanswered question. The streaming instability (hereafter SI, Youdin & Goodman 2005; Youdin & Johansen 2007), which arises from angular momentum exchange between the gas and dust components, has emerged as a leading mechanism to overcome these barriers. It has indeed been shown that gravitationally bound 100–1000 km sized bodies can form directly through the gravitational collapse of the filaments produced during the nonlinear evolution of the SI (Johansen et al. 2009; Simon et al. 2016; Abod et al. 2019). The original linear stability analysis, restricted to inviscid discs, indicates that the growth rate of the SI is highest for solids that are marginally coupled to the gas and for local dust-to-gas ratios ϵ higher than unity.

Recent works have extended this original linear analysis of the instability to examine whether the SI remains robust in the presence of external turbulence in the disk (Chen & Lin 2020; Umurhan et al. 2020). In the particular case when turbulence is driven by the magneto-rotational instability (hereafter MRI, Balbus & Hawley 1991), Johansen et al. (2007) found that the SI can still operate and lead to the formation of structures that concentrate into high-density regions; while Balsara et al. (2009) and Tilley et al. (2010) also reported significant dust clumping through the SI during the dust-settling process. However, except in its innermost regions, where thermal ionization guarantees that MRI operates, the bulk of the disk is expected to remain laminar because the ionization level of the gas is low. In these regions, nonideal MHD effects play an important role in the disk dynamics, as they can drive magnetically driven disk winds, which are thought to be the dominant mechanism for angular momentum transport in these discs (Bai & Stone 2013; Gressel et al. 2015; Béthune et al. 2017).

For these nonideal effects on the SI, Yang et al. (2018) considered a disk model that was sandwiched between a turbulent surface layer and a dead zone dominated by ohmic resistivity. The authors found that strong clumping of solids can arise in the dead zone due to a weak radial diffusion of particles near the midplane. On the other hand, by examining the stability of a dusty magnetized gas, Lin & Hsu (2022) found that in resistive discs, the SI can be stabilized by magnetic perturbations for low values of the dust-to-gas ratio. The consequence of the Hall effect on the SI was investigated by Wu et al. (2024). By employing a combination of linear analysis and numerical simulations, these authors found a new instability resulting from the advection of magnetic perturbations relative to the dusty gas, which can dominate the SI for low dust-to-gas ratios and weak magnetic fields.

We focus on the effect of ambipolar diffusion on the SI. Given the importance of ambipolar diffusion in the context of wind-driven protoplanetary discs, this is an important issue to consider. We find evidence for a strong resonant drag instability (hereafter RDI; Squire & Hopkins 2018a,b; Zhuravlev 2019) resulting from the destabilization of an ambipolar-modified Alfv’en wave by the relative drift between the gas and dust components. Compared to a standard Alfvén wave RDI (Hopkins & Squire 2018; Lin & Hsu 2022), ambipolar diffusion can trigger RDI over a wider range of wavelengths because it can modify the Alfvén frequency. We also find that dust can stabilize the unstable oblique modes of the MRI (Kunz & Balbus 2004; Desch 2004), giving rise to the ambipolar-shear instability.

The paper is organized as follows. In Sect. 2 we present the governing evolution equations of the gas and dust components. We then describe in Sect. 3 the linear theory and present numerical solutions to the linear stability problem in Sect. 4. In Sect. 5, we discuss the effect of ambipolar diffusion on the Alfvén wave RDI. In Sect. 6 we present the results of nonlinear simulations of the instability. We finally summarize our findings in Sect. 7.

2 Physical model

2.1 Shearing sheet approximation

We modeled a local patch in a protoplanetary disk within the shearing-box framework (Goldreich & Lynden-Bell 1965; Latter & Papaloizou 2017). A Cartesian coordinate system with the origin located at an arbitrary distance R₀ from the central star corotates with angular velocity Ω = Ω_k(R₀), with Ω_k the Keplerian angular velocity. The x- and y-axes were oriented radially outward and along the orbital direction, respectively, while the z-axis was directed in the vertical direction. We assumed that the domain is small compared to the orbital distance so that Keplerian rotation appears as a linear shear flow with velocity ${\mathbf{U}}_{\mathbf{k}}=-\frac{3}{2}x\mathrm{\Omega }\hat{\mathbf{y}}$. We assumed that the system is axisymmetric and neglected the vertical component of stellar gravity. In these limits, the continuity and momentum equations for the gas and dust components are given by $$\frac{\mathrm{\partial }{\rho }_g}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_g\mathbf{v}\right)=0$$(1) $$\begin{array}{rl}\frac{\mathrm{\partial }\mathbf{v}}{\mathrm{\partial }t}+\mathbf{v}\cdot \mathrm{\nabla }\mathbf{v}-2v_y\mathrm{\Omega }\hat{\mathbf{x}}+v_x\frac{\mathrm{\Omega }}{2}\hat{\mathbf{y}}=& 2\eta R{\mathrm{\Omega }}^2\hat{\mathbf{x}}-\frac{1}{{\rho }_g}\mathrm{\nabla }(P+\frac{B^2}{2{\mu }_0})\\ & +\frac{1}{{\rho }_g{\mu }_0}(\mathbf{B}\cdot \mathrm{\nabla })\mathbf{B}+\frac{ϵ}{{\tau }_s}(\mathbf{w}-\mathbf{v})\end{array}$$(2) $$\frac{\mathrm{\partial }{\rho }_d}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\rho }_d\mathbf{w}\right)=0$$(3) $$\frac{\mathrm{\partial }\mathbf{w}}{\mathrm{\partial }t}+\mathbf{w}\cdot \mathrm{\nabla }\mathbf{w}-2w_y\mathrm{\Omega }\hat{\mathbf{x}}+w_x\frac{\mathrm{\Omega }}{2}\hat{\mathbf{y}}=-\frac{1}{{\tau }_s}(\mathbf{w}-\mathbf{v}),$$(4)

respectively, where ρ_g and ρ_d are the gas and dust densities, respectively, v and w are their velocities measured relatively to the background Keplerian shear, and μ₀ is the magnetic permeability. We adopted an isothermal equation of state with the pressure given by $P={\rho }_g{c}_s^2$, where c_s = H_gΩ is the (constant) sound speed, with H_g the gas-pressure scale height. In Eq. (2), the term ∝ η corresponds to an outward force acting on the gas through a background radial pressure gradient and is determined by the dimensionless parameter $$\eta =-\frac{1}{2{\rho }_g{\mathrm{\Omega }}^2{R}_0}\frac{\mathrm{\partial }P}{\mathrm{\partial }r}.$$(5)

Finally, ϵ = ρ_d/ρ_g is the dust-to-gas ratio, and τ_s is the stopping time that we characterized in the following in terms of the dimensionless Stokes number: St = τ_sΩ.

The induction equation is given by $$\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=\mathrm{\nabla }\times (\mathbf{v}\times \mathbf{B})-\frac{3}{2}\mathrm{\Omega }{B}_x\hat{\mathbf{y}}+\mathrm{\nabla }\times [\frac{1}{{\gamma }_{in}{\rho }_i{\rho }_g}(\mathbf{J}\times \mathbf{B})\times \mathbf{B}],$$(6)

where $\mathbf{J}\equiv {\mu }_0^{-1}\mathrm{\nabla }\times \mathbf{B}$ is the current density, γ_in is the ion-neutral drag coefficient, and ρ_i is the ion mass density.

2.2 Equilibrium state

In the following, we assumed a constant magnetic field and null radial component. In practice, any radial component is sheared by differential rotation such that the radial component rapidly becomes much smaller that the azimuthal component. Under these conditions, there is no deviation from Keplerian rotation because the magnefic field and the equilibrium conditions correspond to the traditional steady-state drift solutions for a dusty disk obtained by Nakagawa et al. (1986). For constant values of ρ_g, ρ_d, and St, equilibrium solutions to Eqs. (1)–(4) can be obtained, leading to velocity deviations from Keplerian rotation given by $$v_x=\frac{2ϵ\mathrm{S}\mathrm{t}}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(7) $$v_y=-\frac{(1+ϵ+{\mathrm{S}\mathrm{t}}^2)}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(8) $$w_x=\frac{-2\mathrm{S}\mathrm{t}}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(9) $$w_y=-\frac{(1+ϵ)}{{\mathrm{\Delta }}^2}\eta R_0\mathrm{\Omega },$$(10)

where Δ² = St² + (1 + ϵ)².

3 Linear theory

3.1 Linearized equations

In this section, we perturb Eqs. (1)–(6) assuming a background field $\mathbf{B}=B_y\hat{\mathbf{y}}+B_z\hat{\mathbf{z}}$, with B_y (B_z) the azimuthal (vertical) component of the magnetic field. For any variable A, we also assumed Eulerian perturbations such that $$A\to A+\delta A\exp [i(k_{x}x+k_{z}z)+\sigma t)],$$(11)

where linearized quantities are indicated by the δ notation, and where k_x (k_z) is the radial (vertical) wavenumber. Under these conditions, the linearized equations for our dusty magnetized disk are given by $$\sigma W=-i(k_x\delta v_x+k_z\delta v_z)-i{k}_x{v}_{x}W$$(12) $$\begin{array}{rl}\sigma \delta v_x=& 2\mathrm{\Omega }\delta v_y-i{k}_x{v}_x\delta v_x-\frac{ϵ}{{\tau }_s}(v_x-w_x)(Q-W)-\frac{ϵ}{{\tau }_s}(\delta v_x-\delta w_x)\\ & -i{k}_x{c}_s^{2}W+i{c}_{A,z}^2(k_z\delta b_x-k_x\delta b_z-k_x{B}_y^{\mathrm{\prime }}\delta b_y)\end{array}$$(13) $$\begin{array}{rl}\sigma \delta v_y=& -\frac{\mathrm{\Omega }}{2}\delta v_x-i{k}_x{v}_x\delta v_y-\frac{ϵ}{{\tau }_s}(v_y-w_y)(Q-W)\\ & -\frac{ϵ}{{\tau }_s}(\delta v_y-\delta w_y)+i{c}_{A,z}^2{k}_z\delta b_y\end{array}$$(14) $$\sigma \delta v_z=-i{k}_x{v}_x\delta v_z-\frac{ϵ}{{\tau }_s}(\delta v_z-\delta w_z)-i{k}_z{c}_s^{2}W-i{c}_{A,z}^2{k}_z{B}_y^{\mathrm{\prime }}\delta b_y$$(15) $$\sigma Q=-i(k_x\delta w_x+k_z\delta w_z)-i{k}_x{w}_{x}Q$$(16) $$\sigma \delta w_x=2\mathrm{\Omega }\delta w_y-i{k}_x{w}_x\delta w_x-\frac{1}{{\tau }_s}(\delta w_x-\delta v_x)$$(17) $$\sigma \delta w_y=-\frac{\mathrm{\Omega }}{2}\delta w_x-i{k}_x{w}_x\delta w_y-\frac{1}{{\tau }_s}(\delta w_y-\delta v_y)$$(18) $$\sigma \delta w_z=-i{k}_x{w}_x\delta w_z-\frac{1}{{\tau }_s}(\delta w_z-\delta v_z)$$(19) $$\sigma \delta b_x=i{k}_z\delta v_x-i{k}_x{v}_x\delta b_x+{\eta }_{\mathrm{A}\mathrm{D}}(-k^2\delta b_x+k_x{k}_z{B}_y^{\mathrm{\prime }}\delta b_y)$$(20) $$\begin{array}{rl}\sigma \delta b_y=& i{k}_z\delta v_y-i(k_x\delta v_x+k_z\delta v_z)B_y^{\mathrm{\prime }}-i{k}_x{v}_x\delta b_y-\frac{3}{2}\mathrm{\Omega }\delta b_x\\ & +{\eta }_{\mathrm{A}\mathrm{D}}[k^2\frac{k_x}{k_z}{B}_y^{\mathrm{\prime }}\delta b_x-(k^2{B}_y^{\mathrm{\prime }2}+k_z^2)\delta b_y]\end{array}$$(21) $$\begin{array}{rl}\sigma \delta b_z=& -i{k}_x\delta v_x-i{k}_x{v}_x\delta b_z\\ & +{\eta }_{\mathrm{A}\mathrm{D}}(k_x{k}_z\delta b_x-k_x^2{B}_y^{\mathrm{\prime }}\delta b_y-k_x^2\delta b_z),\end{array}$$(22)

where $k=\sqrt{k_x^2+k_z^2},W=\delta {\rho }_g/{\rho }_g,Q=\delta {\rho }_d/{\rho }_d,\delta b_{x,y,z}=\delta B_{x,y,z}/B_z,B_y^{\mathrm{\prime }}=B_y/B_z$, and where $c_{A,z}=B_z/\sqrt{{\mu }_0{\rho }_g}$ is the Alfvén speed and ${\eta }_{\mathrm{A}\mathrm{D}}=\frac{B_z^2}{{\gamma }_{in}{\rho }_i{\rho }_g}$ is the ambipolar resistivity. The linearized system consisting of the previous equations corresponds to an eigenvalue problem, $$\mathbf{M}\mathbf{q}=\sigma \mathbf{q},$$(23)

where M is the matrix representation of the right-hand side of Eqs. (12)–(22), and q is the eigenvector of complex amplitudes. Here, M is an 11 × 11 matrix and q = [W, Q, δv, δw, δb]^T. We solved the stability problem numerically using the function eig of NumPy to find the dimensionless growth rate s = σ/Ω as a function of the following parameters:

• St: The Stokes number of the particle size

• ϵ: The dust-to-gas ratio

• ${\beta }_z=\frac{c_s^2}{c_{A,z}^2}$: the plasma β parameter

• ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}=\frac{c_{\mathrm{A},z}^2}{{\eta }_{\mathrm{A}\mathrm{D}}\mathrm{\Omega }}$: the ambipolar Elsasser number

• $B_y^{\mathrm{\prime }}$: the azimuthal field B_y relative to the vertical field B_z

• ${\tilde{K}}_{x,z}=k_{x,z}\eta R$: the dimensionless radial and vertical wavenumbers.

3.2 Resonant drag instabilities

It is generally expected that the relative drift between dust and gas can render a gas-dust mixture unstable to RDIs (Squire & Hopkins 2018a,b). An RDI can be triggered whenever the gas can support waves that propagate at the same speed as the dust drift velocity. The RDI condition is then obtained by equating the frequency ω_gas of the wave involved in the RDI to k · (w − v), with k the wavevector. In the particular case of an unstratified and axisymmetric shearing box, the RDI condition then becomes $${\omega }_{\mathrm{g}\mathrm{a}\mathrm{s}}(k_x,k_z)=k_x(w_x-v_x).$$(24)

In this context, the classical streaming instability is an RDI driven by the destabilization of inertial waves with frequencies given by $${\omega }_{\mathrm{g}\mathrm{a}\mathrm{s}}^2=\frac{k_z^2}{k_x^2+k_z^2}{\mathrm{\Omega }}^2,$$(25)

such that the RDI condition given by Eq. (24) reads $${\tilde{K}}_z^2=\frac{{\tilde{K}}_x^2{\zeta }_x^2}{1-{\tilde{K}}_x^2{\zeta }_x^2},$$(26)

with $${\zeta }_x^2=\frac{4{\mathrm{St}}^2(1+ϵ{)}^2}{{\mathrm{\Delta }}^4}\eta .$$(27)

In a magnetized disc, the gas also supports Alfvén waves that can be eventually destabilized by the dust-gas relative velocity when condition (24) is fulfilled. When we neglect the effect of rotation for the time being and assume a purely vertical magnetic field, the Alfvén waves have frequencies given by $${\omega }_{\mathrm{g}\mathrm{a}\mathrm{s}}^2=k_z^2{c}_{A,z}^2,$$(28)

resulting in an RDI condition given by $${\tilde{K}}_z^2={\tilde{K}}_x^2{\zeta }_x^2{\beta }_z\frac{\eta }{h^2},$$(29)

where h = H_g/R₀ is the disk aspect ratio. Lin & Hsu (2022) indeed found Alfvén wave RDIs in dusty magnetized discs, but they are easily suppressed by ohmic dissipation. However, we show below that ambipolar diffusion can revive them in the presence of an azimuthal field.

4 Numerical results

4.1 Effect of ambipolar diffusion on the SI

In this section, we describe the numerical solutions to the aforementioned eigenvalue problem. To highlight the effect of ambipolar diffusion on the SI, we first varied ϵ and studied the effect of changing the sign of B_y while keeping the other parameters fixed. We then investigated the effect of changing St and the amplitude of the azimuthal field B_y.

Figure 1 shows the growth rates of unstable modes as a function of wavenumber for St = 0.1, β_z = 10⁴, and for different values of ϵ that increase from the upper to the lower panels. The left column corresponds to an almost ideal MHD configuration, with β_z = 10⁴, whereas the middle and right columns compare the cases $B_y^{\mathrm{\prime }}$ = 1 and $B_y^{\mathrm{\prime }}$ = −1. In this figure, the solid and dashed black lines represent the RDI conditions for the classical and Alfvén wave SI (Squire & Hopkins 2018a,b), given respectively by Eqs. (26) and (29). Finally, the blue line shows the resonance condition for the AD-modified Alfvén wave RDI (see Sect. 5.2).

The upper row corresponds to the dust-free case. For Λ_AD = 10⁴ (left panel), the unstable modes in the lower left corner are the signature of the MRI, which should satisfy $k{c}_{A,z}<\sqrt{3}\mathrm{\Omega }$ (Balbus & Hawley 1991). For Λ_AD = 100 and $B_y^{\mathrm{\prime }}$ > 0, the effect of ambipolar diffusion is expected to be similar to that of ohmic resistivity, with a dissipation rate that decreases faster toward small wavenumbers (Sano & Miyama 1999). Compared to the ideal MHD case, adding ambipolar diffusion in that case is expected to stabilize modes with the highest wavenumbers, which is indeed what is observed. For $B_y^{\mathrm{\prime }}$ < 0, oblique modes with k_x ≠ 0 can be destabilized (Kunz & Balbus 2004) by ambipolar diffusion, which in that case couples with differential rotation to amplify the magnetic field (Desch 2004). Since a value Λ_AD = 100 is considered here, the oblique modes have lower growth rates than axial wavenumbers, although we note in passing that they can be more unstable than the axial modes for Λ_AD < 1 (Lesur 2021).

We now turn to the case of dusty magnetized discs. The middle row of Fig. 1 shows the growth rates for a dust-poor disk with ϵ = 0.2 and St = 0.1, corresponding to the Lin B setup in the terminology of Youdin & Johansen (2007). In the idealMHD case (left panel), the unstable modes clustering along the straight line for ${\tilde{K}}_x,{\tilde{K}}_z\gtrsim 10$ can be clearly identified as Alfvén wave streaming instabilities (AwSI) arising from the resonance between Alfvén waves and the dust-gas radial drift, since they follow the resonance condition given by Eq. (28). For an AD-dominated disk with Λ_AD = 100 and $B_y^{\mathrm{\prime }}$ = 1 (middle panel), however, we find that these RDI unstable modes can be stabilized by AD. Again, as AD-dominated and resistive discs are expected to behave similarly in that case, this is consistent with the results of Lin & Hsu (2022), who found that the AwSI modes are easily damped by ohmic resistivity. For $B_y^{\mathrm{\prime }}$ = −1, a primary effect of dust loading is to stabilize the MRI oblique modes of Kunz & Balbus (2004). Moreover, compared to the $B_y^{\mathrm{\prime }}$ > 0 case, a similar effect of AwSI mode damping is obtained for ${\tilde{K}}_x,{\tilde{K}}_z\gtrsim {10}^3$. For smaller wavenumbers with $10\lesssim {\tilde{K}}_x,{\tilde{K}}_z\lesssim {10}^3$, however, we find that ambipolar diffusion in contrast tends to increase the growth rates of the unstable modes located slightly away from the straight line corresponding to the resonance condition (Eq. (28)). Although ambipolar diffusion does not significantly increase the maximum growth rate compared to the ideal-MHD case, its main effect here is to extend the parameter space of unstable AwSI modes toward larger ${\tilde{K}}_z$.

This effect becomes more pronounced as ϵ increases. The bottom row of Fig. 1 shows the growth rates for a dust-rich disk with ϵ = 3 and St = 0.1 (corresponding to the Lin A setup in the terminology of Youdin & Johansen 2007). For $B_y^{\mathrm{\prime }}$ < 0 (bottom right panel), the flared region of unstable modes with ${10}^2\lesssim {\tilde{K}}_x,{\tilde{K}}_z\lesssim {10}^4$ arises from the effect of ambipolar diffusion. Again, AD does not increase the maximum growth rate of AwSI modes, but rather increases the parameter space of unstable modes for which σ ~ Ω. Compared to the ideal-MHD case, ambipolar diffusion indeed tends to destabilize modes with vertical wavelengths shorter by two orders of magnitude. In the following, we refer to the modes that are destabilized by AD as AmSI modes. We show in Sect. 5.2 that these unstables modes are due to an RDI associated with a resonance between AD-modified Alfvén waves and the dust-gas radial drift.

Fig. 1 Growth rates of unstable modes for different amplitudes of the azimuthal magnetic field for a dust-free disk with ϵ = 0 (top panel), a dust-poor disk with ϵ = 0.2 (middle panel), and a dust-rich disk with ϵ = 3. The Stokes number was fixed to St = 0.1 and the ambipolar Elsasser number to ΛAD = 100. The solid line corresponds to the RDI condition for the streaming instability (Eq. (26)), while the dashed black line corresponds to the RDI condition between Alfvén waves and the dust-gas radial drift (Eq. (29)). The dashed blue line corresponds to the condition for AD-modified Alfvén wave RDI (Eq. (39)).

Fig. 2 Growth rates of unstable modes for different amplitudes of the azimuthal magnetic field. The dust-to-gas ratio was fixed to ϵ = 3, the Stokes number to St = 0.1, and the ambipolar Elsasser number to ΛAD = 10.

Fig. 3 Growth rates of unstable modes for different values of the Stokes number for ϵ = 0.01 (top) and ϵ = 0.2 (bottom). The ambipolar Elsasser number was set to ΛAD = 10 and the azimuthal field to By/Bz = −1.

4.2 Parameter study

In this section, we discuss the dependence of the AD-driven instability described above on parameters, with a particular focus on the effect of varying the azimuthal field and the Stokes number. In Fig. 2 we display growth rates for the Lin A setup (ϵ = 3, St = 0.1) with Λ_AD = 10 and for different values of B_y/B_z. For $B_y^{\mathrm{\prime }}$ = −1 (second panel), the AmSI modes appear as the band of unstable modes enclosed within the region delimited by ${10}^2\lesssim {\tilde{K}}_x\lesssim {10}^3$ and $10\lesssim {\tilde{K}}_z\lesssim {10}^2$. Compared to the case with Λ_AD = 10² (bottom right panel in Fig. 1), the range of unstable modes is therefore truncated at larger wavenumbers because the ambipolar diffusivity is stronger. The smaller Elsasser number also causes the classic SI modes to appear in the blue-green rectangular region corresponding to $1\lesssim {\tilde{K}}_x\lesssim {10}^2$ and ${10}^2\lesssim {\tilde{K}}_z\lesssim {10}^5$, with growth rates that increase with the amplitude of the azimuthal field because the ambipolar resistivity scales with B², such that increasing $\left|B_y^{\mathrm{\prime }}\right|$ pushes the system toward the hydrodynamic limit.

Figure 3 shows the growth rates for different values of the Stokes numbers and for dust-poor disks with relatively low values ϵ = 0.01, 0.2. For ϵ = 0.01 (upper panel), the growth rates remain relatively modest, typically, with σ ≈ 10⁻³Ω. This case, however, (i) confirms that dust can stabilize the original oblique modes of the MRI, and it shows (ii) that the AmSI modes seem to emerge from the stable branch of these oblique modes. Because the value for ϵ is low, it seems unlikely that this is a consequence of the reduced Alfvén speed caused by dust loading. Although not shown here, we confirmed this by considering a purely magnetized gas disc, but employing an effective Alfvén speed given by $$c_{A,eff}=\frac{c_{A,z}}{\sqrt{1+ϵ}}.$$(30)

For ϵ = 0.2 (lower panel), it is interesting to note that the AmSI modes can have growth rates of σ ~ 0.1Ω, even for Stokes numbers as small as St = 10⁻⁴. This suggests that in comparison with the classic SI modes, the conditions of ϵ and St for the unstable AmSI modes to develop with significant growth rates are less stringent.

5 Unstable AmSI modes

We previously identified a new branch of unstable modes that we referred to as AmSI modes. The aim of this section is to understand the origin of the different properties of these AmSI modes that emerged from our numerical calculations. In the following, we discuss in particular how dust can stabilize the oblique modes of the ambipolar-shear instability. We also demonstrate that the AmSI modes can be interpreted as arising from an RDI driven by an AD-modified Alfvén wave and the relative drift between the gas and dust components.

5.1 Effect of dust on MRI oblique modes

To investigate the effect of dust loading on the stability of the MRI oblique modes, we employed a single-fluid model of a magnetized dusty gas disk that we present in detail in Appendix A. We then used the dispersion relation that we obtained (see Eq. (A.14)) and followed Kunz & Balbus (2004) by assuming that the transition from stability to instability occurs through σ ≈ 0. In the limit σ → 0, Eq. (A.14) becomes $$\sigma =-\frac{C_0}{C_1}(1+\frac{\mathrm{\Re }({\mathcal{D}}_1){\tau }_s}{C_0}),$$(31)

where the coefficients 𝒞₀ and 𝒞₁ > 0 are given in Kunz & Balbus (2004, see their Eqs. (34)–(35)), and where ℜ(𝒟₁) was obtained from Eq. (A.19) and reads $$\mathfrak{R}({\mathcal{D}}_1)=-2f_d\frac{B_y}{B_z}\frac{k_x}{k_z}{\left(c_{A,z}{k}_z\right)}^4\mathrm{\Omega }+\frac{B_z^2}{\gamma {\rho }_i{\rho }_g}{f}_d\frac{B_y^2}{B_z^2}{c}_{A,z}^2{k}_x^2{k}_z^2{f}_d{\mathrm{\Omega }}^2,$$(32)

where f_d = ϵ/(1 + ϵ). The latter expression can also be written as $$\mathrm{\Re }({\mathcal{D}}_1)=-2f_d\frac{B_y}{B_z}\frac{k_x}{k_z}{\left(c_{A,z}{k}_z\right)}^4\mathrm{\Omega }[1-\frac{{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^{-1}}{2}\frac{B_y}{B_z}\frac{k_x}{k_z}].$$(33)

We note that in deriving Eq. (31), 𝒟₀ was not taken into account as it is purely imaginary. For a purely magnetized gas disc, MRI oblique modes are unstable for 𝒞₀ < 0 (Kunz & Balbus 2004), which requires the condition B_yk_x < 0 to be satisfied (Lesur 2021). In this case, it is clear that ℜ(𝒟₁) > 0, which leads from Eq. (31) to a lower growth rate or even to a stabilization of the MRI oblique modes when ℜ(𝒟₁) is large enough. This is fully consistent with our numerical results, which show that adding dust indeed tends to stabilize the oblique modes with high wavenumbers (see, e.g., the middle right panel of Fig. 1). Physically, this can be interpreted as resulting from the dust backreaction on the gas. It is indeed expected that this leads to a decrease in the azimuthal velocity of the gas, such that the effect of the feedback loop between ambipolar diffusion and differential rotation is weakened.

Similarly, for 𝒞₀ > 0 (but still assuming B_yk_x < 0), Eq. (31) confirms that dust has a stabilizing effect because the term in parentheses is positive. This suggests that for small Stokes numbers and low dust-to-gas ratios, the AmSI modes that seem to originate from the stable branch adjacent to the unstable oblique modes (e.g., upper left panel of Fig. 3) are not caused by the effect of dust on the ambipolar shear instability.

5.2 AD-modified Alfvén wave RDI: No rotation

The results presented above also show a clear trend for the instability region to lie along the RDI condition for the streaming instability of Alfvén waves. To examine to which extent ambipolar diffusion causes Alfvén waves to resonate with the dust-gas drift, we first simplified the analysis presented in Sect. 3 by neglecting rotation and also considered the incompressible limit. We discuss the effect of rotation in the next section in more detail. Under these approximations, it is straightforward to show that the perturbed gas momentum and induction equations, projected along the direction of $\hat{\mathbf{n}}=\hat{\mathbf{k}}\times \hat{\mathbf{b}}$, where $\hat{\mathbf{k}}$ and $\hat{\mathbf{b}}$ are unit vectors along the directions of the k and B, respectively, leads to the following dispersion relation for the magnetic perturbations oriented along $\hat{\mathbf{n}}$: $${\omega }^2+i\omega \gamma B^2{k}^2{\cos }^2\theta -k_z^2{c}_{A,z}^2=0,$$(34)

where γ = 1/γ_inρ_iρ_g, θ is the angle between $\hat{\mathbf{k}}$ and $\hat{\mathbf{b}}$ and with σ = −iω, while magnetic perturbations oriented along $\hat{\mathbf{b}}$ obey the following relation dispersion: $${\omega }^2+i\omega \gamma B^2{k}^2-k_z^2{c}_{A,z}^2=0.$$(35)

These expressions are similar to those derived by Desch (2004) and Pandey & Wardle (2008).

In the following, we restrict the study to the case of magnetic perturbations directed along $\hat{\mathbf{n}}$, which correspond to classical Alfvén waves, as we show below that the instability described above can be interpreted as an RDI that is triggered when an AD-modified Alfvén wave resonates with dust. In particular, this is supported by Fig. 4, which shows for the case ϵ = 3, St = 0.1, and B_y/B_z = −1 the instability growth rate as a function of ${\tilde{K}}_z$ for different values of Λ_AD (bottom right panel in Fig. 1), and for ${\tilde{K}}_x={10}^2,{10}^3$. As Λ_AD increases, that is, as the ideal MHD regime is approached, the range of vertical modes that are unstable to pure Alfvén wave streaming instabilities is clearly recovered (bottom left panel in Fig. 1).

To clearly demonstrate that magnetic perturbations oriented along $\hat{\mathbf{n}}$ may be involved in the instability, we plot in Fig. 5 the velocity and magnetic field components along this direction as a function of vertical wavenumber ${\tilde{K}}_z$ for ϵ = 3, St = 0.1, and Λ_AD = 100 (corresponding to the second panel in Fig. 2) and ϵ = 0.2, St = 0.001, and Λ_AD = 10 (corresponding to the lower panel of Fig. 3). For the range of unstable vertical modes with the highest growth rates, the velocity and magnetic perturbations are essentially directed along $\hat{\mathbf{n}}$, which supports the transverse wave interpretation, whose dispersion relation is given by Eq. (34).

We note that this relation dispersion can also be written as $${\omega }^2+i\omega {\eta }_{\mathrm{A}\mathrm{D}}{k}_z^2-k_z^2{c}_{A,z}^2=0.$$(36)

The roots of Eq. (36) have real and imaginary parts that are given by $$\mathrm{\Re }(\omega )=c_{A,z}{k}_z{(1-\frac{{\eta }_{\mathrm{A}\mathrm{D}}^2{k}_z^2}{4c_{A,z}^2})}^{1/2}$$(37)

and $$\mathfrak{J}(\omega )=-i\frac{{\eta }_{\mathrm{A}\mathrm{D}}{k}_z^2}{2}.$$(38)

Solving for the resonance condition ℜ(ω) = k_x(w_x − v_x), where ℜ(ω) is given by Eq. (37) and with v_x, w_x given by Eqs. (7) and (9), respectively, leads to the following relation between the vertical and radial dimensionless wavenumbers: $${\tilde{K}}_z^2=\frac{2{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^2{\eta }^2{\beta }_z}{h^2}(1\pm \sqrt{1-\frac{{\tilde{K}}_x^2{\zeta }_x^2}{{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^2\eta }}).$$(39)

In the limit $D\equiv \frac{{\tilde{K}}_x^2{\zeta }_x^2}{{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^2\eta }\ll 1$, we can expand the square root to O(D²) to obtain by taking the negative solution $${\tilde{K}}_z^2\approx \frac{\eta {\zeta }_x^2{\beta }_z{\tilde{K}}_x^2}{h^2}(1+\frac{{\zeta }_x^2{\tilde{K}}_x^2}{4\eta {\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^2}),$$(40)

and Eq. (29) is then readily recovered in the limit of a vanishing AD resistivity.

The resonance condition given by Eq. (39) is marked by the blue line in Figs. 1, 2 and 3. For St = 0.1 (Fig. 1), the model reproduces the branch of unstable modes that appear in the B_y < 0 case reasonably well. For small Stokes numbers (Fig. 3), however, the agreement is less clear because the branch of unstable modes emerging from the original MRI oblique modes does not closely follow the resonance condition. To determine whether an RDI really operates in this case, we compared for each (k_x, k_z) the oscillation frequencies found numerically by solving the stability problem represented by Eq. (23), that is, corresponding to ℑ(σ), to the frequency expected from the resonance condition k_x(w_x − v_x). The results of this procedure for ϵ = 0.01 and St = 0.001 are displayed in Fig. 6. It is clear that for each unstable (k_x, k_z) pair, the two oscillation frequencies are comparable, which seems to support the RDI interpretation.

Fig. 4 Instability growth rate as a function of vertical wavenumber ${\tilde{K}}_z$ for different values of the ambipolar Elsasser number ΛAD, for ${\tilde{K}}_x={10}^3$ (upper panel) and ${\tilde{K}}_x={10}^2$ (lower panel). The dust-to-gas ratio was fixed to ϵ = 3 and the Stokes number to St = 0.1.

Fig. 5 Velocity and magnetic field components along $\hat{\mathbf{n}}=\hat{\mathbf{k}}\times \hat{\mathbf{b}}$ as a function of vertical wavenumber ${\tilde{K}}_z$ for ϵ = 3, St = 0.1, ΛAD = 100 (upper panel), and ϵ = 0.2, St = 0.001, and ΛAD = 10 (lower panel). The azimuthal field was fixed to By/Bz = −1 and the radial wavenumber to ${\tilde{K}}_x={10}^3$.

Fig. 6 For ϵ = 0.01 and St = 0.001, the ratio of oscillation frequencies ℑ(σ) found by solving the stability problem (Eq. (23)) and frequencies expected from the resonance condition kx(wx − vx).

5.3 Near-resonant growth

Compared with the ideal MHD case, the results presented above suggest that the primary effect of ambipolar diffusion is to widen the range of unstable modes for which Alfvén waves can be involved in an RDI. To explain the significant width of the resonance that is obtained, we note that in presence of ambipolar diffusion, the system is almost identical to a forced oscillator, with the characteristic equation given by $${\omega }^2+i\omega {\eta }_{\mathrm{A}\mathrm{D}}{k}_z^2-k_z^2{c}_{A,z}^2={\omega }_{RDI}^2,$$(41)

where we set ω_RDI = k_x(w_x − v_x), and whose oscillation amplitude 𝒜_RDI would be given by $${\mathcal{A}}_{RDI}\propto \frac{1}{\sqrt{{({\omega }_{RDI}^2-k_z^2{c}_{A,z}^2)}^2+{\left({\eta }_{\mathrm{A}\mathrm{D}}{k}_z^2{\omega }_{RDI}\right)}^2}}.$$(42)

The equivalent quality factor associated with this forced oscillator, which is related to the resonance width Δω through 𝒬 = k_zc_A,z/Δω, reads $$\mathcal{Q}=\frac{c_{A,z}}{{\eta }_{\mathrm{A}\mathrm{D}}{k}_z}=\frac{{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}}{{\tilde{K}}_z}\sqrt{\beta }h,$$(43)

where h is the disk aspect ratio. This corresponds to a resonant width given by $$\mathrm{\Delta }\omega =\frac{{\tilde{K}}_z^2}{{\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}\beta \eta }\mathrm{\Omega }.$$(44)

Setting Λ_AD = 100 and h = 0.05, we obtain Q = 0.5 for ${\tilde{K}}_z$ = 10³, such that we indeed expect the resonance to have a significant width in the presence of ambipolar diffusion. In particular, for the parameters adopted here, k_zc_A,z ≈ 200Ω and the resonance width is estimated to Δω ≈ 400Ω.

5.4 Alfvén-modified RDI: Effect of rotation

Although the instability we found may be interpreted by the destabilization of an AD-modified Alfvén wave by the radial drift between dust and gas, it remains unclear why it disappears for B_yk_x > 0. We suggest that, similar to the unstable oblique modes in nonideal MHD, the instability also originates from the combination of shear and ambipolar diffusion, and thus, it can only operate under certain field geometries. In order to examine how the resonance condition given by Eq. (39) is affected by the effect of shear, we started with the dispersion relation of pure MHD waves in the presence of rotation (Kunz 2008), $$\begin{array}{rl}& ({\omega }^2+i\omega \gamma B_z^2{k}_z^2-k_z^2{c}_{A,z}^2)({\omega }^2+i\omega \gamma B^2{k}^2-k_z^2{c}_{A,z}^2)\\ & =-2A\gamma \left(k_x{B}_y\right)\left(k_z{B}_z\right)[{\omega }^2-k_z^2{c}_{A,z}^2],\end{array}$$(45)

where A = ∂Ω/∂ log R. To extract the effect of rotation, we searched for a correction of the AD-modified Alfvén wave frequency by writing ω = ω₀ + ω₁, where ω₀ is given by Eq. (37) and with ω₁ ∝ A. By evaluating the first derivative of the left term together with the right term of (45) with respect to ω and evaluating it at ω = ω₀, and then equating terms of order A, we obtained $${\omega }_1=2Ai\frac{\left(k_x{B}_y^{\mathrm{\prime }}\right)\gamma B_z^2}{c_{A,z}}\frac{k_z^2}{k_x^2+k^2{B}_y^{\mathrm{\prime }2}}.$$(46)

Because the growth rate is σ = −iω and A < 0, an azimuthal field with orientation such that k_xB_y < 0 results in σ > 0, as expected. Conversely, k_xB_y > 0 leads to a damping of Alfvén waves through ambipolar diffusion. For k_xB_y > 0, we estimated under which conditions Alfvén waves are strongly damped by simply assuming ω₁ = ω₀ in that case. This occurs for $${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}<\frac{2|A|{\mathrm{\Omega }}^{-1}{k}_x{k}_z{B}_y^{\mathrm{\prime }}}{k_x^2+k^2{B}_y^{\mathrm{\prime }2}}.$$(47)

For k_x = k_z and $B_y^{\mathrm{\prime }}$ = 1, this leads to Λ_AD < 1, which is consistent with what is observed in Fig. 1.

6 Direct simulations

We demonstrate the AmSI with selected axisymmetric simulations using the performance-portable Godunov code idefix (Lesur et al. 2023). To solve Eqs. (1)–(4) and (6), we enabled the MHD, dust, and shearing-box modules. We chose piecewise linear reconstruction, the HLLD Riemann solver, and a second-order Runge-Kutta time integration. Ambipolar diffusion was integrated using a Runge–Kutta–Legendre super-time-stepping scheme. We modeled a single dust species with a constant stopping time¹ and included feedback onto the gas.

Our prescription for the ambipolar diffusivity in IDEFIX follows the standard form (e.g., Cui & Lin, 2021), $${\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }t}|}_{\mathrm{A}\mathrm{D},\mathrm{I}\mathrm{D}\mathrm{E}\mathrm{F}\mathrm{I}\mathrm{X}}=\mathrm{\nabla }\times [\frac{{\eta }_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}}{|\mathit{B}{|}^2}(\mathit{J}\times \mathit{B})\times \mathit{B}],$$(48)

with a constant ${\eta }_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}$. This is related to the ambipolar diffusivity used above by ${\eta }_{\mathrm{A}\mathrm{D}}={\eta }_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}{B}_z^2/|\mathit{B}{|}^2$. At the linear level, this is only a constant factor and therefore does not affect the stability properties of the dusty gas. In Appendix B, we present code tests against the AmSI linear modes.

For the illustrative cases here, we used η = 0.05h for the radial pressure gradient and initialized the magnetic field with β_z = 10⁴, B_y = −B_z, and an Elsassar number ${\mathrm{\Lambda }}_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}\equiv c_{A,z}^2/{\eta }_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}\mathrm{\Omega }$ = 10. For the dust, we considered an initial ϵ = 0.2 and St = 0.01 or St = 0.1. The initial dust and gas densities are uniform, with velocities given by Eqs. (7)–(10).

To isolate the AmSIs, which have short characteristic length scales, we adopted a small computational domain with x, z ∈ [−0.005, 0.005]H_g to exclude MRI-related modes. We applied periodic boundaries in z and shear-periodic boundaries in x, and we used a resolution of N_x,z = 512. The simulations were made in 3D but were axisymmetric, which was realized by setting the azimuthal domain to one cell. To trigger instability, we added random perturbations to the gas radial velocities with an amplitude of 10⁻³c_s. Computational units were such that H_g = Ω = μ₀ = 1 and the initial ρ_g = 1.

Figure 7 shows the evolution of the maximum gas velocity perturbations and dust-to-gas ratios for the two simulations, and Fig. 8 show snapshots of the dust-to-gas ratio. As expected, the systems are linearly unstable, with growth rates ≃0.05Ω, 0.2Ω and saturation amplitudes max |δv| ~ 3 × 10⁻⁴c_s, 3 × 10⁻³c_s for St = 0.01 and St = 0.1, respectively. For St = 0.01, the dust concentrations are weak, with max(ϵ) ≲ 0.3, or about an increase of 50% from the equilibrium value.

On the other hand, for St = 0.1, the dust densities can reach max (ϵ) ~ 10. The system undergoes strong oscillations, with rapid formation and destruction of transient clumps with ϵ ≳ 1 (see the right panel of Fig. 8 for an example clump). The average is max (ϵ) ≃ 2, that is, ten times the initial value. This run can be compared with the classical SI case of Johansen & Youdin (2007), who also used ϵ = 0.2 initially and St = 0.1 (their run AA), which found much lower dust overdensities of about 25%. This suggests that AmSI is more effective at concentrating dust in magnetized disks than in unmagnetized disks. However, it should be noted that MRI modes are absent by construction. The role of the AmSI in the presence of the MRI needs to be clarified in the future.

Fig. 7 Nonlinear evolution of the AmSI with St = 0.01 (blue) and St = 0.1 (orange). Left: maximum gas velocity perturbation. Right: maximum dust-to-gas ratio.

Fig. 8 Snapshots of the dust-to-gas ratios at the end of the AmSI simulations shown in Fig. 7. Left: St = 0.01. Right: St = 0.1.

7 Discussion and conclusion

We examined the effect of AD on the stability of a dusty magnetized disc. The ionization fraction is expected to be relatively low in the main body of the disc, such that ambipolar diffusion is thought to play an important role in the dynamics and evolution of protoplanetary discs. In particular, it has been shown that the combined effect of ambipolar diffusion and ohmic resistivity can launch magnetically driven disk winds (Bai & Stone 2013; Gressel et al. 2015; Béthune et al. 2017) that drive accretion through the disk surfaces.

Dust feedback tends to stabilize the oblique modes of the MRI identified by Kunz & Balbus (2004). This occurs as a result of a decrease in the gas azimuthal velocity, such that the feedback loop between ambipolar diffusion and differential rotation is weakened. Depending on the orientation of the azimuthal field, our linear stability analysis also showed that ambipolar diffusion can give rise to a vigorous RDI, resulting from the destabilization of an Alfvén wave by dust-gas drift. The main effect of ambipolar diffusion is to strongly modify the Alfvén wave frequency, with the consequence that the range of wavenumbers exhibiting significant growth rates σ ≳ 0.1Ω⁻¹ is significantly extended in comparison with a standard Alfvén wave RDI.

Our analysis showed that ambipolar diffusion has a stronger effect when the azimuthal field amplitude B_y is comparable to the vertical field B_z and for ambipolar Elsasser numbers Λ_AD ~ 10, typically. In the outer regions of protoplanetary discs, where AD is the dominant nonideal MHD effect, Λ_AD ~ 1 in the disk midplane (Bai 2011), while amplitudes of azimuthal fields such that B_y/B_z = 𝒪(10) were found in calculations incorporating AD (Cui & Bai 2021). Under these conditions, the requirements for efficient instability is generally not satisfied. At smaller radii R ~ 1–20 AU, however, global simulations of magnetized disk winds show that the toroidal and vertical fields can have comparable amplitudes (Bai 2017), particularly at heights Z ≲ 2H, with H the gas scale height (Gressel et al. 2015). In these inner disk regions, Λ_AD is still expected to be on the order of unity, but can also be significantly higher in the presence of large grains. For example, at distances R ~ 10 AU, Simon et al. (2015) reported a value Λ_AD ~ 10 at the midplane of a disk that contained 10 μm dust grains. These results therefore suggest that the instability is most likely to operate in the intermediate regions of wind-launching discs, particularly when a substantial population of large grains is present.

Interestingly, we obtained significant growth rates in dust-poor discs with ϵ ~ 0.2 and for Stokes numbers as low as St = 10⁻⁴, although numerical simulations suggest that dust enhancements resulting from unstable AmSI modes might be relatively modest in that case. For higher Stokes numbers with St = 0.1, however, we obtained more significant dust concentrations compared to an unmagnetized disc, which might facilitate planetesimal formation in magnetized discs subject to ambipolar diffusion.

A main limitation of this work is that the effect of vertical stratification was not considered. In the hydrodynamical limit and for St ~0.1, Lim et al. (2025) recently obtained good agreement between stratified and unstratified simulations of the streaming instability, such that the midplane dynamics can be well described using an unstratified setup. In a magnetized disk launching a magnetocentrifugal wind, however, AD tends to generate a flip of the toroidal field close to the disk midplane (Béthune et al. 2017). In this context, numerical simulations investigating the effect of AD on the streaming instability should ideally take into account the vertical structure of the disc. Numerical simulations of vertically global radially unstratified wind-emitting dusty discs will be presented in a future study.

Acknowledgements

Computer time for this study was provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Universite de Bordeaux and by HPC resources of Cines under the allocation A0170406957 made by GENCI (Grand Equipement National de Calcul Intensif). Simulations were also carried out on the Kawas cluster at ASIAA. The authors acknowledge the access to high-performance computing facilities (Theory cluster and storage) provided by ASIAA. MKL is supported by the National Science and Technology Council (grants 113-2112-M-001-036-, 114-2112-M-001-018-, 113-2124-M-002-003-, 114-2124-M-002-003-, 115-2124-M-002-014-), an Academia Sinica Career Development Award (AS-CDA-110-M06), and an Academia Sinica Grand Challenge Seed Grant (AS-GCS-115-M02).

References

Appendix A Single-fluid model

To get insights into the effect of dust on the on the oblique modes of Kunz & Balbus (2004), we employed a more simple model of a dusty, magnetized gas. In this approach, we defined the total dust density as $$\rho ={\rho }_g+{\rho }_d$$(A.1)

and the center of mass velocity as $$\mathbf{u}=f_g\mathbf{v}+f_d\mathbf{w},$$(A.2)

where f_g = 1/(1 + ϵ) is the gas fraction and f_d = ϵ/(1 + ϵ) the dust fraction. We work within the framework of the magnetized terminal velocity approximation, which corresponds to a first order approximation in terms of the Stokes number, and which is therefore better suited to small particles. In this regime, the velocity difference between the dust and gas phases Δu = w − v is given by $$\mathrm{\Delta }\mathbf{u}={\tau }_s[\frac{\mathrm{\nabla }P}{\rho }-2\eta R{\mathrm{\Omega }}^2{f}_g\hat{\mathbf{x}}-\frac{1}{{\mu }_0\rho }(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}]$$(A.3)

In this limit, the gas incompressibility condition, dust continuity equation, and the center-of-mass momentum equation for the dust-gas mixture are $$\mathrm{\nabla }\cdot \mathbf{u}=\mathrm{\nabla }\cdot \left(f_d\mathrm{\Delta }\mathbf{u}\right)$$(A.4) $$\frac{\mathrm{\partial }ϵ}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot (ϵ\mathbf{u})=-\mathrm{\nabla }\cdot \mathbf{u}$$(A.5) $$\frac{\mathrm{\partial }\mathbf{u}}{\mathrm{\partial }t}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}=-\frac{1}{\rho }\mathrm{\nabla }P+2\eta r{\mathrm{\Omega }}^2{f}_g\hat{\mathbf{x}}+2\mathrm{\Omega }{u}_y\hat{\mathbf{x}}-\frac{\mathrm{\Omega }}{2}{u}_x\hat{\mathbf{y}}+\frac{1}{{\mu }_0\rho }(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}.$$(A.6)

In principle, the original induction equation given in Eq. 6 should be expressed in terms of u and Δu. Following Lin & Hsu (2022), we do not take this into account and simply set v → u in Eq. 6, such that the induction equation is leaved unchanged in the single fluid model. We revisit this approximation in Appendix A.4.

A.1 Linearized equations

$$\mathrm{\nabla }\cdot \delta \mathbf{u}=-\frac{\sigma }{1+ϵ}\delta ϵ$$(A.7) $$\begin{array}{rl}\sigma \delta ϵ=& 2i{k}_x{\tau }_s\eta R{\mathrm{\Omega }}^2{f}_g\frac{1-ϵ}{1+ϵ}\delta ϵ+f_d{\tau }_s{k}^2(1+ϵ)\frac{\delta P}{\rho }\\ & +f_d{\tau }_s\frac{B_z}{\rho {\mu }_0}{k}^2(-\frac{k_x}{k_z}\delta B_x+\frac{B_y}{B_z}\delta B_y)\end{array}$$(A.8) $$\sigma \delta u_x=2\mathrm{\Omega }\delta u_y-i{k}_x\frac{\delta P}{\rho }+i\frac{k^2}{k_z}\frac{B_z}{\rho {\mu }_0}\delta B_x-i{k}_x\frac{B_y}{\rho {\mu }_0}\delta B_y-\frac{2\eta R{\mathrm{\Omega }}^2{f}_g}{(1+ϵ{)}^2}\delta ϵ$$(A.9) $$\sigma \delta u_y=-\frac{\mathrm{\Omega }}{2}\delta u_x+i{k}_z\frac{B_z}{\rho {\mu }_0}\delta B_y$$(A.10) $$\sigma \delta u_z=-i{k}_z\frac{\delta P}{\rho }-i{k}_z\frac{B_y}{\rho {\mu }_0}\delta B_y$$(A.11) $$\sigma \delta B_x=i{k}_z{B}_z\delta u_x+\frac{1}{\gamma {\rho }_i{\rho }_g}(-k^2{B}_z^2\delta B_x+k_x{k}_z{B}_z{B}_y\delta B_y)$$(A.12) $$\begin{array}{rl}\sigma \delta B_y=& i{k}_z{B}_z\delta u_y-i(k_x\delta u_x+k_z\delta u_z)B_y-\frac{3}{2}\mathrm{\Omega }\delta B_x\\ & +\frac{1}{\gamma {\rho }_i{\rho }_g}(k^2\frac{k_x}{k_z}{B}_z{B}_y\delta B_x-(k^2{B}_y^2+k_z^2{B}_z^2)\delta B_y)\end{array}.$$(A.13)

A.2 Dispersion relation

We find that the dispersion relation can be written in the form $$\begin{array}{rl}(f_d{\sigma }^6+{\mathcal{D}}_5{\sigma }^5+& {\mathcal{D}}_4{\sigma }^4+{\mathcal{D}}_3{\sigma }^3+{\mathcal{D}}_2{\sigma }^2+{\mathcal{D}}_1\sigma +{\mathcal{D}}_0){\mathit{\tau }}_s\\ & +({\sigma }^4+C_3{\sigma }^3+C_2{\sigma }^2+C_1\sigma +C_0)\sigma =0,\end{array}$$(A.14)

where the coefficients 𝒞_i are those given in Eqs. 32–35 of Kunz & Balbus (2004) and the coefficients 𝒟_i are defined by $${\mathcal{D}}_5=\frac{f_d}{\gamma {\rho }_i{\rho }_g}(k^2{B}^2+k_z^2{B}_z^2)$$(A.15) $$\begin{array}{rl}{\mathcal{D}}_4=& 2f_d{c}_{A,z}^2{k}_z^2+{\mathrm{\Omega }}^2(f_d-2i\eta R{k}_x)+\frac{3}{2}{f}_d\frac{B_z{B}_y}{\gamma {\rho }_i{\rho }_g}{k}_x{k}_z\mathrm{\Omega }\\ & +f_d{\left(\frac{B{B}_z}{\gamma {\rho }_i{\rho }_g}\right)}^2{k}^2{k}_z^2\end{array}$$(A.16) $$\begin{array}{rl}{\mathcal{D}}_3=& -2f_d\frac{B_y}{B_z}{\left(c_{A,z}{k}_z\right)}^2\frac{k_x}{k_z}\mathrm{\Omega }\\ & +\frac{B^2{k}^2+B_z^2{k}_z^2}{\gamma {\rho }_i{\rho }_g}(f_d{\left(c_{A,z}{k}_z\right)}^2+{\mathrm{\Omega }}^2(f_d-2i\eta R{k}_x))\end{array}$$(A.17) $$\begin{array}{rl}{\mathcal{D}}_2=& f_d{\left(c_{A,z}{k}_z\right)}^4-{\left(c_{A,z}{k}_z\right)}^2(3f_d+4i\eta R{k}_x){\mathrm{\Omega }}^2-2i\eta \frac{k_x{k}_z^2}{k^2}{\mathrm{\Omega }}^4\\ & +\frac{B_z{B}_y}{2\gamma {\rho }_i{\rho }_g}{k}_x{k}_z(3(f_d-2i\eta R{k}_x){\mathrm{\Omega }}^2-f_d{\left(c_{A,z}{k}_z\right)}^2)\\ & +{\left(\frac{B{B}_z}{\gamma {\rho }_i{\rho }_g}\right)}^2{k}^2{k}_z^2(f_d-2i\eta R{k}_x){\mathrm{\Omega }}^2\end{array}$$(A.18) $$\begin{array}{rl}{\mathcal{D}}_1=& -2f_d\frac{B_y}{B_z}\frac{k_x}{k_z}{\left(c_{A,z}{k}_z\right)}^4\mathrm{\Omega }\\ & +\frac{B_z^2}{\gamma {\rho }_i{\rho }_g}(\frac{B_y^2}{B_z^2}[{\left(c_{A,z}{k}_x\right)}^2(f_d-2i\eta R{k}_x)-2i\eta k_x{\left(c_{A,z}{k}_z\right)}^2]k_z^2{\mathrm{\Omega }}^2\\ & -2i\eta R{k}_x(1+\frac{k_z^2}{k^2})(c_{A,z}^2{k}^2+{\mathrm{\Omega }}^2))\end{array}$$(A.19)

$$\begin{array}{rl}{\mathcal{D}}_0=& -2i{\left(c_{A,z}{k}_z\right)}^2\eta R{k}_x({\left(c_{A,z}{k}_z\right)}^2-3{\mathrm{\Omega }}^2\frac{k_z^2}{k^2}){\mathrm{\Omega }}^2\\ & -3i\frac{B_y{B}_z}{\gamma {\rho }_i{\rho }_g}\eta R{k}_x^2{k}_z({\left(c_{A,z}{k}_z\right)}^2+{\mathrm{\Omega }}^2\frac{k_z^2}{k^2})\\ & -2i{\left(\frac{B{B}_z}{\gamma {\rho }_i{\rho }_g}\right)}^2\eta R{k}_x{k}_z^4{\mathrm{\Omega }}^4\end{array}$$(A.20)

A.3 Limit of strong nonideal effects

In the limit where ambipolar diffusion is very strong, namely for Elsasser numbers Λ_A ≪ 1, an approximation to the dispersion relation can be obtained by retaining only terms $\propto {\mathrm{\Lambda }}_A^{-2}$. In this regime, we find that the dispersion relation becomes $$f_g{\tau }_s{\sigma }^4+{\sigma }^3+{\tau }_s(f_d{\mathrm{\Omega }}^2-2i\eta k_{x}R{\mathrm{\Omega }}^2)+\frac{k_z^2}{k^2}{\mathrm{\Omega }}^2\sigma -2i{k}_x{\tau }_s\eta R{\mathrm{\Omega }}^4\frac{k_z^2}{k^2}=0$$(A.21)

which is equivalent to classical SI relation dispersion (Jacquet et al. 2011; Laibe & Price 2014; Lin & Youdin 2017). For St ≪ 1 and σ = 𝒪(St), the first term in the previous expression can be neglected, and we recover the cubic dispersion relation of Youdin & Goodman (2005).

A.4 Effective ambipolar diffusion from dust-gas drag

Recall that we neglected Δu in the induction equation for the above analysis. To evaluate the consequence of this approximation, we note that $$\mathbf{v}=\mathbf{u}-f_d\mathrm{\Delta }\mathbf{u}.$$(A.22)

It follows that the induction equation, in the ideal MHD limit, can be rewritten as $$\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=\mathrm{\nabla }\times (\mathbf{u}\times \mathbf{B})-\frac{3}{2}\mathrm{\Omega }{B}_x\hat{\mathbf{y}}-\mathrm{\nabla }\times (f_d\mathrm{\Delta }\mathbf{u}\times \mathbf{B}).$$(A.23)

Using the magnetized terminal velocity approximation given by Eq. A.3, the previous expression becomes $$\begin{array}{rl}\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=& \mathrm{\nabla }\times (\mathbf{u}\times \mathbf{B})-\frac{3}{2}\mathrm{\Omega }{B}_x\hat{\mathbf{y}}-\mathrm{\nabla }\times [f_d{\tau }_s(\frac{\mathrm{\nabla }P}{\rho }-2\eta R{\mathrm{\Omega }}^2{f}_g\hat{\mathbf{x}})\times \mathbf{B}]\\ & +\mathrm{\nabla }\times \{{\tau }_s{f}_d{f}_g{c}_A^2[(\mathrm{\nabla }\times \mathbf{B})\times \hat{\mathbf{b}}]\times \hat{\mathbf{b}}\},\end{array}$$(A.24)

with $c_A^2=B^2/{\mu }_0{\rho }_g$ and $\hat{\mathbf{b}}$ = B/B. Therefore, the relative drift between dust and gas phases has the same effect of the drift between neutral and charged species in MHD, as it also leads to ambipolar-like diffusion in the final term on the right hand side.

From the previous equation, the dust effective ambipolar resistivity would be given by ${\eta }_{A,\text{eff}}={\tau }_s{f}_d{f}_g{c}_A^2$, corresponding to an effective ambipolar Elsasser number: $${\mathrm{\Lambda }}_{A,\mathrm{e}\mathrm{f}\mathrm{f}}=\frac{1}{\mathrm{S}\mathrm{t}{f}_d{f}_g}.$$(A.25)

Fig. B.1 Evolution of the maximum gas velocity perturbations for IDEFIX runs initialized with clean eigenmodes (blue and orange curves). The stars are the expected evolution based on linear theory growth rates.

For small grains with St ≪ 1, we have Λ_A,eff ≫ 1 and thus this effect is likely unimportant compared to the true ambipolar diffusion in protoplanetary discs.

Appendix B Code test

We tested our IDEFIX setup by simulating a clean eigenmode and comparing its growth rate against linear theory. The equilibrium disk parameters are the same as that used for Sect. 6. We perturbed the equilibrium by adding axisymmetric perturbations calculated from linear theory (i.e., Eqs. 12–22 with η_AD replaced by ${\eta }_{\mathrm{A}\mathrm{D}}^{\mathrm{\prime }}$), which are characterized by wavenumbers k_x,z. We normalized the perturbations so that max |δv_x| = 10⁻⁶c_s. We considered two cases: 1) St = 0.01 with k_xH_g = 4000, k_zH_g = 1500; and 2) St = 0.1 with k_x,zH_g = 1000. For these tests, the computational domain is one wavelength, that is, x, z ∈ [−π/k_x,z, π/k_x,z], and we use N_x,z = 128 cells in each direction.

Figure B.1 shows the evolution of the maximum gas velocity perturbations, from which we measure growth rates over the first 8 orbits and also compare them in Table B.1. These show a good agreement between simulation and theoretical growth rates, with a relative error ≲ 1%. Not surprisingly, the St = 0.1 case, with a higher growth rate, is more accurately simulated.

The test here shows that IDEFIX is capable of capturing the AmSI. We have also reproduced growth rates with N_x,z = 64 cells at similar accuracy by using a less diffusive configuration: piecewise parabolic reconstruction, a Roe solver, and third order Runge-Kutta time integration. However, we found this setup is prone to crashing (with vanishing timesteps) when evolving the AmSI into the nonlinear regime. Hence we do not use it in Sect. 6.

All Tables

All Figures

Fig. 1 Growth rates of unstable modes for different amplitudes of the azimuthal magnetic field for a dust-free disk with ϵ = 0 (top panel), a dust-poor disk with ϵ = 0.2 (middle panel), and a dust-rich disk with ϵ = 3. The Stokes number was fixed to St = 0.1 and the ambipolar Elsasser number to ΛAD = 100. The solid line corresponds to the RDI condition for the streaming instability (Eq. (26)), while the dashed black line corresponds to the RDI condition between Alfvén waves and the dust-gas radial drift (Eq. (29)). The dashed blue line corresponds to the condition for AD-modified Alfvén wave RDI (Eq. (39)).

In the text

	Fig. 2 Growth rates of unstable modes for different amplitudes of the azimuthal magnetic field. The dust-to-gas ratio was fixed to ϵ = 3, the Stokes number to St = 0.1, and the ambipolar Elsasser number to ΛAD = 10.
In the text

	Fig. 3 Growth rates of unstable modes for different values of the Stokes number for ϵ = 0.01 (top) and ϵ = 0.2 (bottom). The ambipolar Elsasser number was set to ΛAD = 10 and the azimuthal field to By/Bz = −1.
In the text

	Fig. 4 Instability growth rate as a function of vertical wavenumber ${\tilde{K}}_z$ for different values of the ambipolar Elsasser number ΛAD, for ${\tilde{K}}_x={10}^3$ (upper panel) and ${\tilde{K}}_x={10}^2$ (lower panel). The dust-to-gas ratio was fixed to ϵ = 3 and the Stokes number to St = 0.1.
In the text

	Fig. 5 Velocity and magnetic field components along $\hat{\mathbf{n}}=\hat{\mathbf{k}}\times \hat{\mathbf{b}}$ as a function of vertical wavenumber ${\tilde{K}}_z$ for ϵ = 3, St = 0.1, ΛAD = 100 (upper panel), and ϵ = 0.2, St = 0.001, and ΛAD = 10 (lower panel). The azimuthal field was fixed to By/Bz = −1 and the radial wavenumber to ${\tilde{K}}_x={10}^3$.
In the text

	Fig. 6 For ϵ = 0.01 and St = 0.001, the ratio of oscillation frequencies ℑ(σ) found by solving the stability problem (Eq. (23)) and frequencies expected from the resonance condition kx(wx − vx).
In the text

	Fig. 7 Nonlinear evolution of the AmSI with St = 0.01 (blue) and St = 0.1 (orange). Left: maximum gas velocity perturbation. Right: maximum dust-to-gas ratio.
In the text

	Fig. 8 Snapshots of the dust-to-gas ratios at the end of the AmSI simulations shown in Fig. 7. Left: St = 0.01. Right: St = 0.1.
In the text

	Fig. B.1 Evolution of the maximum gas velocity perturbations for IDEFIX runs initialized with clean eigenmodes (blue and orange curves). The stars are the expected evolution based on linear theory growth rates.
In the text