© The Authors 2025

1 Introduction

Visible matter in the universe mostly consists of gas, with 0.99 being the canonical fraction (see e.g. Draine et al. 2007). The remaining 1% is in solid form, and consists of particles of various sizes, from sub-micron size to planetary size, that are often mixed with the gas. Examples of dusty gas flows include AGB star winds (Höfner & Olofsson 2018), the interstellar medium (Draine 2003), and protoplanetary discs (Bae et al. 2023). The dust component can play an important role in the chemistry and opacity of the mixture, for example, but also in its dynamics. In particular, adding dust to a gas in an otherwise stable configuration can make the system go hydrodynamically unstable.

One particular class of instabilities in gas-dust mixtures are the so-called resonant drag instabilities (RDIs, see Squire & Hopkins 2018a,b). These instabilities rely on the fact that gas and dust are usually drifting relative to each other. If a wave in the gas can propagate at the same velocity as the drift velocity of the dust, the result is often a very strong instability that grows fast even if the total amount of dust is very small (Squire & Hopkins 2018b). Of particular interest is the streaming instability (SI, Youdin & Goodman 2005), which can play an important role in planet formation, and was shown to be an inertial wave RDI (Squire & Hopkins 2018a; Zhuravlev 2019; Magnan et al. 2024b). If the system of interest allows for different gas waves, different types of RDI can manifest, such as the acoustic RDI (Hopkins & Squire 2018) or RDIs in magnetised gas (Hopkins et al. 2020).

Dust drift speeds usually depend on particle size. When accelerated for example due to gravity, dust particles eventually reach an equilibrium terminal velocity where the acceleration is balanced by (size-dependent) gas drag. The smallest particles, tightly coupled to the gas, drift slowest, while larger particles drift faster. Since RDIs rely on the drift speed being equal to the gas wave speed, it is natural to ask what happens to an RDI when there is a distribution of dust sizes, and therefore a distribution of drift velocities.

In this paper, we focus on the acoustic resonant drag instability. This instability occurs in systems where drift is fast enough to couple to sound waves, and where magnetic forces are subdominant (see Figure 2 of Hopkins et al. 2022). A recent physical picture of the acoustic RDI was presented in Magnan et al. (2024a). The instability occurs for example in dust-driven outflows, where radiation pressure drives dust particles away from AGB stars (e.g. Höfner & Olofsson 2018), around supernova remnants (Micelotta et al. 2018), or around star clusters in starburst galaxies (Menon et al. 2023). The nonlinear outcome of the acoustic RDI can be virulent turbulence with strong dust clumping (Moseley et al. 2019). This has strong implications for the chemistry and opacity in the system (see also Seligman et al. 2019; Hopkins et al. 2020).

The purpose of this paper is twofold: first of all, we want to study the linear phase of the acoustic RDI with a dust size distribution for its own sake. In this sense, we are continuing the work of Squire et al. (2022), who mostly focused on the case where all sizes have the same drift velocity (but see their appendix a3). We exclusively deal with the case where the drift speed is different for different grain sizes. The second purpose of this paper is to provide a framework for interpreting RDIs in the polydisperse regime. Since the acoustic RDI (hereafter ARDI) is probably the simplest RDI (in terms of the physics required as well as the number of spatial dimensions), it lends itself well for this. For a companion paper, we used the framework to deal with more complex cases of the streaming instability (Youdin & Goodman 2005) and settling instability (Squire & Hopkins 2018b,a), which are likely important in protoplanetary discs during planet formation.

The plan of this paper is as follows. We introduce the basic equations in Section 2, where we also define our equilibrium state and explain how we derived equations for linear perturbations. Results are presented in Section 3, and we conclude in Section 4.

2 Governing equations

The equations governing a mixture of gas and dust, where the dust component consist of a continuum of sizes, were presented in Paardekooper et al. (2021): ${\mathrm{\partial }}_t\sigma +\mathrm{\nabla }\cdot (\sigma \mathbf{u})=0,$(1) ${\mathrm{\partial }}_t\mathbf{u}+(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{u}={\alpha }_{\mathrm{d}}-\frac{\mathbf{u}-{\mathbf{v}}_{\mathrm{g}}}{{\tau }_{\mathrm{s}}}.$(2)

Here, σ the size density (Paardekooper et al. 2020), u the size-dependent dust velocity, and v_g the gas velocity. We specialise drag law of the form (u – v_g)/τ_s to couple gas and dust, with stopping time τ_s (see e.g. Woitke & Helling 2003). While the stopping time can in principle depend on gas density and |u – v_g| (Squire & Hopkins 2018b), we consider only the simple case where the stopping time does not depend on any dynamical quantity. Any additional accelerations acting on the dust are contained in α_d, which we assume to depend on u only. We note that under this assumption, the dust momentum equation is completely decoupled from the evolution of the size density.

The gas component, similarly, has its continuity and momentum equation, where the backreaction of drag on the gas enters the latter: ${\mathrm{\partial }}_t{\rho }_{\mathrm{g}}+\mathrm{\nabla }\cdot \left({\rho }_{\mathrm{g}}{\mathbf{v}}_{\mathrm{g}}\right)=0,$(3) ${\mathrm{\partial }}_t{\mathbf{v}}_{\mathrm{g}}+({\mathbf{v}}_{\mathrm{g}}\cdot \mathrm{\nabla }){\mathbf{v}}_{\mathrm{g}}=-\frac{\mathrm{\nabla }p}{{\rho }_{\mathrm{g}}}+{\mathit{\alpha }}_{\mathrm{g}}+\frac{1}{{\rho }_{\mathrm{g}}}\int \sigma \frac{\mathbf{u}-{\mathbf{v}}_{\mathrm{g}}}{{\tau }_{\mathrm{s}}}\mathrm{d}{\tau }_{\mathrm{s}},$(4)

where ρ_g is the gas density, p the pressure, and all additional accelerations on the gas are contained in α_g. At this point, we can remain agnostic about any other equations governing the gas; there might be an energy equation, equations coupling the magnetic field, etc. The only assumption here is that the only way gas and dust interact is through the drag force. Nevertheless, in this paper, we are concerned with only very simple gas physics: no magnetic fields, and a barotropic equation of state, so that (3) and (4) completely describe the evolution of the gas component.

It is worth noting that, as long as we do not consider any interaction between the dust particles themselves, we have complete freedom in letting σ depend on multiple dust properties. Apart from size, one could consider shape, porosity, chemical composition, etc. Given dust parameters p₁, p₂, …, p_n, we write $\sigma =\sigma (p_1,p_2,\dots ,p_n),$(5) ${\tau }_{\mathrm{s}}={\tau }_{\mathrm{s}}(p_1,p_2,\dots ,p_n),$(6)

and demand that ${\rho }_{\mathrm{d}}=\int \sigma \mathrm{d}{p}_1\mathrm{d}{p}_2\dots \mathrm{d}{p}_n,$(7)

where ρ_d is the mass volume density of dust. The drag integral appearing in Eq. (4) will also be over all dust parameters. This is useful if, for example, α_d depends on dust parameters, or if one considers a more complex drag law. In this paper, we do not consider these cases, but since the only relevant dust parameter is the stopping time, we can simplify notation by choosing p₁ = τ_s, so that the drag integral is performed over stopping time rather than size¹.

The size distribution is often taken to be a power law, at least initially. Most often, the MRN distribution is considered (Mathis et al. 1977), which has the number density of dust grains ∝a^−7/2, where a is the dust size, and therefore $\sigma \propto {\tau }_{\mathrm{s}}^{-1/2}$ . Different power laws were explored, for example, in Krapp et al. (2019); Zhu & Yang (2021), while in McNally et al. (2021) power laws were augmented with a bump at larger sizes to account for grain growth and in particular cratering (Birnstiel et al. 2011). In this paper, to make the various integrals more tractable, we concern ourselves with very simple size distributions only. In particular, we focus on a constant size distribution between a minimum and a maximum stopping time: $\sigma \left({\tau }_{\mathrm{s}}\right)=\{\begin{array}{ll}\frac{{\rho }_{\mathrm{d}}}{{\tau }_{\mathrm{s},\text{max}}-{\tau }_{\mathrm{s},\text{min}}}& {\tau }_{\mathrm{s},\text{min}}<{\tau }_{\mathrm{s}}<{\tau }_{\mathrm{s},\text{max}}\\ 0& \text{otherwise}\end{array}.$(8)

It is straightforwardly verified that the integral over stopping time yields the total dust density, ρ_d.

2.1 Equilibrium state

The ARDI can be set up in 1D Cartesian geometry. We fix the gas equation of state to be isothermal with constant sound speed c_g, and consider a drag law where the stopping time τ_s is independent of gas density. The dust size distribution is parametrised using the size density σ(τ_s), defined so that the dust density is ${\rho }_{\mathrm{d}}=\int \sigma \left({\tau }_{\mathrm{s}}\right)\mathrm{d}{\tau }_{\mathrm{s}}.$(9)

The equations governing the dynamics of the coupled dust-gas system are given by conservation of mass and momentum of gas and dust: ${\mathrm{\partial }}_t{\rho }_{\mathrm{g}}+{\mathrm{\partial }}_x\left({\rho }_{\mathrm{g}}{v}_{\mathrm{g}}\right)=0,$(10) ${\mathrm{\partial }}_t{v}_{\mathrm{g}}+v_{\mathrm{g}}{\mathrm{\partial }}_x{v}_{\mathrm{g}}=-\mu \alpha -\frac{c_{\mathrm{g}}^2{\mathrm{\partial }}_x{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{g}}}+\frac{1}{{\rho }_{\mathrm{g}}}\int \frac{\sigma }{{\tau }_{\mathrm{s}}}(u-v_{\mathrm{g}})\mathrm{d}{\tau }_{\mathrm{s}},$(11) ${\mathrm{\partial }}_t\sigma +{\mathrm{\partial }}_x(\sigma u)=0,$(12) ${\mathrm{\partial }}_{t}u+u{\mathrm{\partial }}_{x}u=\alpha -\frac{u-v_{\mathrm{g}}}{{\tau }_{\mathrm{s}}},$(13)

where the dust velocity u = u(τ_s). A constant acceleration α was chosen to generate drift, similar to Hopkins & Squire (2018) and Magnan et al. (2024a). In reality, this may for example come from radiation pressure, see, for example, Steinwandel et al. (2022). An equilibrium state can be constructed where all quantities are constant in space and time, with drift velocity $u^{(0)}-v_{\mathrm{g}}^{(0)}=\alpha {\tau }_{\mathrm{s}}$ . For simplicity, we choose $v_{\mathrm{g}}^{(0)}=0$ , or, equivalently, we work in a frame that moves with the background gas velocity. We denote the equilibrium gas density ${\rho }_{\mathrm{g}}^{(0)}$ , equilibrium dust size density σ⁽⁰⁾(τ_s), and equilibrium dust-to-gas ratio $\mu =\int {\sigma }^{(0)}\mathrm{d}{\tau }_{\mathrm{s}}/{\rho }_{\mathrm{g}}^{(0)}$ .

2.2 Linear perturbations

We consider small perturbations around the equilibrium state, denoting perturbed quantities with superscript ‘1’: ${\mathrm{\partial }}_t{\rho }_{\mathrm{g}}^{(1)}=-{\rho }_{\mathrm{g}}^{(0)}{\mathrm{\partial }}_x{v}_{\mathrm{g}}^{(1)},$(14) $\begin{array}{rl}{\mathrm{\partial }}_t{v}_{\mathrm{g}}^{(1)}=& -\frac{c_{\mathrm{g}}^2{\mathrm{\partial }}_x{\rho }_{\mathrm{g}}^{(1)}}{{\rho }_{\mathrm{g}}^{(0)}}+\frac{1}{{\rho }_{\mathrm{g}}^{(0)}}\int \frac{{\sigma }^{(0)}}{{\tau }_{\mathrm{s}}}(u^{(1)}-v_{\mathrm{g}}^{(1)})\mathrm{d}{\tau }_{\mathrm{s}}\\ & +\int {\left(\frac{\sigma }{{\rho }_{\mathrm{g}}}\right)}^{(1)}\frac{u^{(0)}-v_{\mathrm{g}}^{(0)}}{{\tau }_{\mathrm{s}}}\mathrm{d}{\tau }_{\mathrm{s}},\end{array}$(15) ${\mathrm{\partial }}_t{\sigma }^{(1)}=-{\sigma }^{(0)}{\mathrm{\partial }}_x{u}^{(1)}-u^{(0)}{\mathrm{\partial }}_x{\sigma }^{(1)},$(16) ${\mathrm{\partial }}_t{u}^{(1)}=-u^{(0)}{\mathrm{\partial }}_x{u}^{(1)}-\frac{u^{(1)}-v_{\mathrm{g}}^{(1)}}{{\tau }_{\mathrm{s}}}.$(17)

We take perturbations ${\rho }_{\mathrm{g}}^{(1)}={\hat{\rho }}_{\mathrm{g}}\exp (\mathrm{i}kx-\mathrm{i}\omega t)$, and use a similar form for other perturbed quantities, $\omega {\hat{\rho }}_{\mathrm{g}}=k{\rho }_{\mathrm{g}}^{(0)}{\hat{v}}_{\mathrm{g}},$(18) $\begin{array}{rl}\omega {\hat{v}}_{\mathrm{g}}=& \frac{k{c}_{\mathrm{g}}^2{\hat{\rho }}_{\mathrm{g}}}{{\rho }_{\mathrm{g}}^{(0)}}+\frac{\mathrm{i}}{{\rho }_{\mathrm{g}}^{(0)}}\int \frac{{\sigma }^{(0)}}{{\tau }_{\mathrm{s}}}(\hat{u}-{\hat{v}}_{\mathrm{g}})\mathrm{d}{\tau }_{\mathrm{s}}\\ & +\mathrm{i}\int \left(\frac{\hat{\sigma }}{{\rho }_{\mathrm{g}}}\right)\frac{u^{(0)}-v_{\mathrm{g}}^{(0)}}{{\tau }_{\mathrm{s}}}\mathrm{d}{\tau }_{\mathrm{s}},\end{array}$(19) $\omega \hat{\sigma }=k{\sigma }^{(0)}\hat{u}+k{u}^{(0)}\hat{\sigma },$(20) $\omega \hat{u}=k{u}^{(0)}\hat{u}-\mathrm{i}\frac{\hat{u}-{\hat{v}}_{\mathrm{g}}}{{\tau }_{\mathrm{s}}},$(21)

where dust quantities $\hat{\sigma }$ and û depend on the stopping time. We note that the perturbation in dust-to-gas ratio $\left(\frac{\hat{\sigma }}{{\rho }_{\mathrm{g}}}\right)=\frac{\hat{\sigma }}{{\rho }_{\mathrm{g}}^{(0)}}-\frac{{\sigma }^{(0)}}{{\rho }_{\mathrm{g}}^{(0)}}\frac{{\hat{\rho }}_{\mathrm{g}}}{{\rho }_{\mathrm{g}}^{(0)}},$(22)

so that $\mathrm{i}\int \left(\frac{\hat{\sigma }}{{\rho }_{\mathrm{g}}}\right)\frac{u^{(0)}-v_{\mathrm{g}}^{(0)}}{{\tau }_{\mathrm{s}}}\mathrm{d}{\tau }_{\mathrm{s}}=\mathrm{i}\alpha \int \frac{\hat{\sigma }}{{\rho }_{\mathrm{g}}^{(0)}}\mathrm{d}{\tau }_{\mathrm{s}}-\mathrm{i}\alpha \mu \frac{{\hat{\rho }}_{\mathrm{g}}}{{\rho }_{\mathrm{g}}^{(0)}}.$(23)

Equations (18)–(21) form a continuous eigenvalue problem in τ_s, with eigenvalue ω, with a linear operator involving integral terms. These integrals can cause significant difficulties when determining eigenvalues (Paardekooper et al. 2021). These difficulties are related to the presence of the resonance that in the monodisperse regime gives the RDIs their name. In the case of the ARDI, we see from Eq. (20) that if ω = ku⁽⁰⁾ the dust density perturbation diverges. If in the limit of μ → 0 we are dealing with a gas sound wave, that is ω = kc_g, the resonant condition requires dust to drift at the sound speed. This is where fast growth in the monodisperse regime is expected (e.g. Squire & Hopkins 2018b,a). In the polydisperse case, we often have to integrate over this resonance, the presence of which makes the integral nearly singular (for examples see Paardekooper et al. 2021, see also the discussion below Equation (49)). This is where standard integration techniques fail (Paardekooper et al. 2021).

2.3 Dispersion relation

We can combine (18)–(21) into a single equation for the temporal frequency ω: $\begin{array}{rl}& {\omega }^2-c_{\mathrm{g}}^2{k}^2+\int \frac{{\sigma }^{(0)}}{{\rho }_{\mathrm{g}}^{(0)}}\frac{\omega (\omega -u^{(0)}k)}{1-\mathrm{i}{\tau }_{\mathrm{s}}(\omega -u^{(0)}k)}\mathrm{d}{\tau }_{\mathrm{s}}\\ & -\mathrm{i}\alpha \int \frac{{\sigma }^{(0)}}{{\rho }_{\mathrm{g}}^{(0)}}\frac{k\omega }{\omega -k{u}^{(0)}}\frac{1}{1-\mathrm{i}{\tau }_{\mathrm{s}}(\omega -u^{(0)}k)}\mathrm{d}{\tau }_{\mathrm{s}}+\mathrm{i}\alpha k\mu =0.\end{array}$(24)

The monodisperse limit is obtained by considering σ⁽⁰⁾ to be a Dirac delta function at stopping time τ_s, $\begin{array}{rl}& ({\omega }^2-k^2{c}_{\mathrm{g}}^2)(\omega -u^{(0)}k)(\frac{\mathrm{i}}{{\tau }_{\mathrm{s}}}+\omega -u^{(0)}k)\\ & +\frac{\mathrm{i}\mu }{{\tau }_{\mathrm{s}}}[(\omega +\alpha k{\tau }_{\mathrm{s}}){(\omega -u^{(0)}k)}^2-\mathrm{i}\alpha k^2{u}^{(0)}]=0,\end{array}$(25)

which agrees with Magnan et al. (2024a). A non-dimensional version, with timescale c_g/α and length scale $c_{\mathrm{g}}^2/\alpha$ , reads as: ${\tilde{\omega }}^2=K^2+\mu \int \mathrm{\Sigma }(s)\frac{\mathrm{i}K\tilde{\omega }-\tilde{\omega }(\tilde{\omega }-Ks{)}^2}{(\tilde{\omega }-Ks)(1-\mathrm{i}s(\tilde{\omega }-Ks))}\mathrm{d}s-\mathrm{i}K\mu ,$(26)

with $\tilde{\omega }=c_{\mathrm{g}}\omega /\alpha ,K=k{c}_{\mathrm{g}}^2/\alpha ,s={\tau }_{\mathrm{s}}\alpha /c_{\mathrm{g}}$, and a nondimensional size distribution $\mathrm{\Sigma }=c_{\mathrm{g}}{\sigma }^{(0)}/\left(\alpha {\rho }_{\mathrm{d}}^{(0)}\right)$ , so that ∫ Σ(s) ds = 1. It is worth noting the various contributions on the right hand side: pressure, dust density (first term in the numerator under the integral), relative velocity (second term in the numerator under the integral), and gas density (last term, outside the integral). The stability of the system is then determined by μ, K, and the size distribution. The non-dimensional version of the size density Eq. (8) that we use for the ARDI reads as: $\mathrm{\Sigma }(s)=\{\begin{array}{llll}\frac{1}{s_{\text{max}}-s_{\text{min}}}& s_{\text{min}}& <s& <s_{\text{max}}\\ 0& \text{otherwise}.\end{array}$(27)

Numerical results were obtained by solving (26) with the techniques detailed in Paardekooper et al. (2021); McNally & Paardekooper (2021).

3 Linear ARDI modes

In the analysis of RDIs, the dust-to-gas ratio μ is assumed to be small, and a series solution is developed. From Eq. (25), we see that in the monodisperse case, there are two possible (non-damped) solutions in the limit μ → 0: a gas sound wave with ${\omega }^2=k^2{c}_{\mathrm{g}}^2$ , and a dust advection mode with ω = u⁽⁰⁾k. Since the gas does not feel the dust for μ = 0, the sound wave propagates independent of any dust perturbation, while for the dust advection mode, a dust density perturbation is advected with the unperturbed dust velocity u⁽⁰⁾, without any dust velocity perturbation.

3.1 Non-resonant growth

We first focus on the non-resonant case, when dust drift is different from the gas sound speed (see also Hopkins & Squire 2018). If $\tilde{\omega }=\pm K\ne Ks$ , we can have a non-resonant growing sound wave in the monodisperse regime. Trying a solution $\tilde{\omega }=K+\mu {\tilde{\omega }}_1$ in Eq. (25) yields ${\tilde{\omega }}_1=\frac{1}{2}\frac{K(1+s)(s-1{)}^2-\mathrm{i}s}{(s-1)(1+\mathrm{i}Ks(s-1))},$(28)

from which we can deduce that growing modes exist when s < 1 and $K^2<\frac{1}{(s+1)|s-1{|}^3}.$(29)

We note that, for our ordering in μ ≪ 1 to be valid, we can not have K smaller than μ. Interestingly, these growing modes can survive in the polydisperse case. In the limit of K → 0, with a size distribution (27) such that s_max < 1, we find that $\mathfrak{J}\left({\tilde{\omega }}_1\right)=-\frac{1}{2}(\log \left(\frac{1-s_{max}}{1-s_{min}}\right)+1).$(30)

Physically, we get growth if the destabilising effect of the dust density perturbation, represented by the logarithm, overcomes the stabilising effect of the gas density perturbation (the constant term, originating from the last term in Eq. (26)). This can happen if s_max is close enough to the resonant size s = 1. For finite K, we get an additional stabilising effect from the dust velocity perturbation, which means we may need to approach the resonant size even closer. From Eq. (28), it is clear that the approximation breaks down at the resonant size where s = 1, but as long as we stay away from this resonance, we can have growing modes ῶ = K + μ ῶ₁. The condition for being far enough away from the resonance is simply |ῶ₁| ≪ K/μ.

If we choose a backward propagating sound wave we can avoid the resonance altogether (see also Squire et al. 2022). Trying a solution $\tilde{\omega }=-K+\mu {\tilde{\omega }}_1$ in Eq. (26), we find that ${\tilde{\omega }}_1=\frac{\mathrm{i}}{2}+\frac{1}{2}\int \mathrm{\Sigma }(s)\frac{K(1+s{)}^2-\mathrm{i}}{(1+s)(1+\mathrm{i}Ks(1+s))}\mathrm{d}s.$(31)

The first term is due to the gas density perturbation, while the first term in the numerator under the integral sign is due to the relative velocity perturbation, and the second term is due to the dust density perturbation. The integrand is well-behaved for s > 0, and while the closed form solution does not provide much insight, a numerical solution is straightforward to calculate. In the limits of small and large K, ω₁ becomes independent of K: ${\tilde{\omega }}_1=\frac{\mathrm{i}}{2}\{\begin{array}{cl}1-\frac{1}{s_{\text{max}}-s_{\text{min}}}\log \left(\frac{1+s_{\text{max}}}{1+s_{\text{min}}}\right),& K\ll 1\\ 1-\frac{1}{s_{\text{max}}-s_{\text{min}}}\log \left(s_{\text{max}}/s_{\text{min}}\right)& K\gg 1.\end{array}$(32)

For K ≪ 1, we always get growth for our constant size distribution, while for K ≫ 1 we usually have $\mathfrak{J}\left({\tilde{\omega }}_1\right)<0$ , as long as s_min is sufficiently smaller than unity. The transition from growth to damping at K ∼ 1 was also found by Squire et al. (2022). As a reminder, note that we need K > μ for our ordering scheme to remain valid.

If instead of a sound wave we choose a non-resonant dust advection mode with ῶ = Ks ≠ K, an expansion ῶ = Ks + μ ῶ₁ yields in the monodisperse case ${\tilde{\omega }}_1=\mathrm{i}\frac{s}{s^2-1},$(33)

with growth for s > 1. In this case, however, the generalisation to the polydisperse case is less straightforward. We consider an expansion ῶ = Ks_* + μ ῶ₁ for some dimensionless, non-resonant stopping time s_*. Equation (26) gives, up to first order in μ: $\begin{array}{rl}& 2\mu K{s}_{\ast }{\tilde{\omega }}_1=K^2(1-s_{\ast }^2)\\ & +\mu \int \mathrm{\Sigma }(s)\frac{\mathrm{i}K\tilde{\omega }-\tilde{\omega }(\tilde{\omega }-Ks{)}^2}{(K{s}_{\ast }+\mu {\tilde{\omega }}_1-Ks)(1-\mathrm{i}s(\tilde{\omega }-Ks))}\mathrm{d}s-\mathrm{i}K\mu ,\end{array}$(34)

where for clarity of notation we have only expanded the relevant terms of ῶ. At zeroth order in μ, in order to have a non-resonant mode with s_* ≠ 1, the first term on the right hand side needs to be balanced by a term in the integral. This happens if K|s−s_*| ≪ μ ῶ₁ for all s. In this case, we simply find at zeroth order in μ that $K^2(1-s_{\ast }^2)+\int \mathrm{\Sigma }(s)\frac{\mathrm{i}{K}^2{s}_{\ast }}{{\tilde{\omega }}_1}\mathrm{d}s=0,$(35)

which immediately leads to the equivalent of the monodisperse result ${\tilde{\omega }}_1=\frac{\mathrm{i}{s}_{\ast }}{s_{\ast }^2-1},$(36)

with growth for s_* > 1, which corresponds to supersonic drift. The condition that K|s − s_*| ≪ μῶ₁ is a condition on the width of the size distribution: away from the resonance, ῶ₁ = O(1), so that we need |s − s_*| ≪ μ/K. Only for such a very narrow size distribution does the dust advection mode exist (remember μ ≪ 1): if the size distribution is too wide, there can be no coherent gas reaction to sustain the mode. This was also noted in Paardekooper et al. (2020) for the streaming instability, which, away from the resonance, also relies on a dust advection mode to go unstable (Youdin & Goodman 2005).

3.2 Resonant modes

We now turn to the resonant case. As was shown in Squire & Hopkins (2018b), if a gas wave propagates with the same velocity as the dust is drifting, a very strong reaction results even for μ ≪ 1. Mathematically, this is manifested as a growth rate proportional to $\sqrt{\mu }$ instead of μ, which for μ ≪ 1 usually leads to stronger growth. In other words, at resonance, the perturbation of the system in the small parameter μ is singular, and in the limit of μ → 0 multiple distinct solution collapse onto each other. The associated defective nature of the matrix eigenvalue problem lies at the heart of the resonant drag instability formalism (Squire & Hopkins 2018b).

Following this reasoning, we try an expansion ῶ = K(1+ ϵ ῶ₁ + ϵ² ῶ₂, with $\mu ={ϵ}^2\overline{\mu }$ with $\overline{\mu }=O(1)$ : $\begin{array}{rl}& 2ϵ{\tilde{\omega }}_1+{ϵ}^2({\tilde{\omega }}_1^2+2{\tilde{\omega }}_2)\\ & -\frac{\overline{\mu }{ϵ}^2}{K}\int \frac{\mathrm{\Sigma }(s)(\mathrm{i}-K{(1+ϵ{\tilde{\omega }}_1+{ϵ}^2{\tilde{\omega }}_2-s)}^2)}{(1+ϵ{\tilde{\omega }}_1+{ϵ}^2{\tilde{\omega }}_2-s)(1-\mathrm{i}sK(1+ϵ{\tilde{\omega }}_1+{ϵ}^2{\tilde{\omega }}_2-s))}\mathrm{d}s\\ & +\frac{\mathrm{i}\overline{\mu }{ϵ}^2}{K}=0,\end{array}$(37)

where we formally assume K = O(1). We note that in the absence of dust (μ = 0), we have a sound wave with ῶ = K, or, equivalently, ω = kc_g. The resonant condition for an RDI is ῶ = Ks, which means at the resonance we need s = 1, which corresponds to exactly sonic drift. We first make contact with the monodisperse results by considering a narrow size distribution.

3.2.1 Narrow size distribution

We consider a narrow size distribution (to be defined below) around the resonance at s = 1. Under the integral, we only need to keep terms up to order ϵ. Let $\begin{array}{rl}I& \equiv {ϵ}^2\int \frac{\mathrm{\Sigma }(s)(\mathrm{i}-K{(1+ϵ{\tilde{\omega }}_1-s)}^2)}{(1+ϵ{\tilde{\omega }}_1-s)(1-\mathrm{i}sK(1+ϵ{\tilde{\omega }}_1-s))}\mathrm{d}s\\ & =ϵI^{(1)}+{ϵ}^2{I}^{(2)},\end{array}$(38)

allowing for the singularity at s = 1 to lead to a term I⁽¹⁾. It is convenient to split the integral: $\begin{array}{rl}I\equiv & \mathrm{i}{ϵ}^2\int \frac{\mathrm{\Sigma }(s)}{(1+ϵ{\tilde{\omega }}_1-s)(1-\mathrm{i}sK(1+ϵ{\tilde{\omega }}_1-s))}\mathrm{d}s\\ & -{ϵ}^{2}K\int \frac{\mathrm{\Sigma }(s)(1-s+\mathrm{i}sK(1-s{)}^2)}{1+s^2{K}^2(1-s{)}^2}\mathrm{d}s.\end{array}$(39)

The first term is due to the dust density perturbation, the second to the perturbation in relative velocity. The last term is O(ϵ²) for all values of K. The first integral can be analysed by noting the resonance at s = 1, and developing a series around s = 1: $\begin{array}{rl}& \mathrm{i}{ϵ}^2\int \frac{\mathrm{\Sigma }(s)}{(1+ϵ{\tilde{\omega }}_1-s)(1-\mathrm{i}sK(1+ϵ{\tilde{\omega }}_1-s))}\mathrm{d}s\\ & =\mathrm{i}{ϵ}^2\int \frac{\sum _{n=0}^{\mathrm{\infty }}{A}_n(1-s{)}^n}{1+ϵ{\tilde{\omega }}_1-s}\mathrm{d}s,\end{array}$(40)

noting that apart from the resonant factor, the integrand is well-behaved around s = 1 for all values of ῶ₁ as long as ϵ ≪ 1. The A_n are simply the coefficients of the Taylor expansion of the function f(s) ≡ Σ(s)/(1−isK (1+ϵῶ₁ − s)) at s = 1: $A_n=\frac{f^{(n)}(1)}{n!},$(41)

and for our constant size distribution (27) the first few terms read: $A_0=\frac{\mathrm{\Sigma }(1)}{1-\mathrm{i}Kϵ{\tilde{\omega }}_1},$(42) $A_1=\frac{iK(ϵ{\tilde{\omega }}_1-1)\mathrm{\Sigma }(1)}{{(\mathrm{i}+Kϵ{\tilde{\omega }}_1)}^2},$(43) $A_2=\frac{\mathrm{i}K(K({ϵ}^2{\tilde{\omega }}_1^2-ϵ{\tilde{\omega }}_1+1)+i)}{{(Kϵ{\tilde{\omega }}_1+\mathrm{i})}^3}.$(44)

If the size distribution is narrow enough so that |1 − s| ≪ ϵ|ῶ₁| for all s, we can approximate the integrand by a constant: $\begin{array}{rl}\mathrm{i}{ϵ}^2\int \frac{\mathrm{\Sigma }(s)}{(1+ϵ{\tilde{\omega }}_1-s)(1-\mathrm{i}sK(1+ϵ{\tilde{\omega }}_1-s))}\mathrm{d}s& \approx \\ \mathrm{i}{ϵ}^2(s_{max}-s_{min})\frac{\mathrm{\Sigma }(1)}{ϵ{\tilde{\omega }}_1(1-\mathrm{i}Kϵ{\tilde{\omega }}_1)}& =\\ \frac{\mathrm{i}ϵ}{{\tilde{\omega }}_1(1-\mathrm{i}Kϵ{\tilde{\omega }}_1)}& \approx \frac{\mathrm{i}ϵ}{{\tilde{\omega }}_1}.\end{array}$(45)

Taking s_max = 1 + δ/2 and s_min = 1 − δ/2, the condition that the integration domain is narrow enough, translates into δ ≪ ϵ|ῶ₁|. This is going to be our definition of a narrow size distribution.

We therefore find that, for a narrow size distribution as defined above, that I⁽¹⁾ = i/ω₁, and therefore ${\tilde{\omega }}_1={\overline{\mu }}^{1/2}\frac{1+\mathrm{i}}{2\sqrt{K}},$(46)

which is equivalent to the monodisperse result. This limit holds for δ ≪ ϵ|ῶ₁|, or $K\ll \frac{\mu }{4{\delta }^2}.$(47)

We can recast this as $\alpha ({\tau }_{max}-{\tau }_{min})\ll \frac{\omega }{k}-c_{\mathrm{g}}.$(48)

Therefore, the system acts as monodisperse if the perturbation in wave velocity is much larger than the width of the background dust velocity distribution. We call this regime quasi-monodisperse. The resonance can tolerate a wider size distribution for larger dust-to-gas ratios and smaller wave numbers, both of which increase the perturbed wave velocity.

The quasi-monodisperse regime is illustrated schematically in Figure 1. We start with a gas sound wave (the solid red arrow), propagating in the absence of dust at a speed indicated by its horizontal position. The resonant drift speed is indicated by the blue arrow, where the dust drifts at the same velocity as the propagation speed of the gas wave. The backreaction of the dust on the gas leads to a perturbed wave speed for the gas, a perturbation that is proportional to μ^1/2 (see Eq. (46)). The perturbed gas wave is indicated by the dashed red arrow. A dust size distribution gives rise to a range of drift speeds, indicated by the blue rectangle. As long as this blue rectangle stays far away from the dashed red arrow, we are in the quasi-monodisperse regime, see Eq. (48).

Under the conditions of a narrow size distribution, the integral in the full dispersion relation (26) can be approximated by a constant in general, resulting in the exact same dispersion relation as in the monodisperse case. This immediately means that, in the quasi-monodisperse regime, we also find the ‘long wavelength pressure-free mode’ (Hopkins & Squire 2018) for K ≪ μ, with a growth rate ∝K^2/3 μ^1/3, and a ‘short wavelength resonant quasi-drift mode’, with a growth rate ∝K^1/3 μ^1/3 (Hopkins & Squire 2018; Magnan et al. 2024a). We note, however, that the latter requires K ≫ 1/μ and K ≪ μ/δ², which means that this short wavelength mode can only operate for extremely narrow size distributions that have δ ≪ μ.

These results are illustrated in Figure 2, where we consider the case with μ = 0.01 and the size distribution given by Eq. (27). The monodisperse limit, indicated by the red dotted line, shows growth at all wavenumbers, with the highest wave numbers growing fastest. In the polydisperse case, there is a cutoff at high wave numbers, which is where the system leaves the quasimonodisperse regime. We will see below that growth is more difficult if the size distribution is too wide. For wave numbers far below the cutoff Eq. (47) and K = O(1), growth rates follow the monodisperse result ∝K^1/2. For K≪ μ, the growth rates follow that of the long wavelength pressure-free mode ∝K^2/3.

Fig. 1 Schematic representation of the quasi-monodisperse regime. The top half of the figure deals with the gas, while the bottom half deals with the dust component. The horizontal axis indicates (wave) speed. The red solid arrow denotes the propagation speed of a gas sound wave in the absence of dust, and the blue arrow indicates the resonant dust drift speed. The backreaction of the dust modifies the propagation speed of the sound wave, leading to the dashed red arrow. The blue rectangle indicates a range of drift speeds if there is a distribution of dust sizes. As long as this rectangle does not come close to the dashed red arrow, we are in the quasi-monodisperse regime.

Fig. 2 Growth rate of the acoustic resonant drag instability as a function of the wave number, which was obtained by numerically solving Eq. (26). The dust-to-gas ratio μ = 0.01, and the size density is constant, centred on the resonant size s = 1, with width δ. The red dotted line indicates the monodisperse result, while the vertical dotted lines indicate the cutoff of Eq. (47).

3.2.2 Wide size distribution

We now return to Eq. (40) but do not make the assumption of a narrow size distribution. The series expansion around s = 1 does not necessarily converge far away from s = 1, depending on the value of K: for K = 1, for example, we need |s − 1| < 0.5. For larger values of K, the size distribution needs to be more narrow in order for the series expansion to converge for all values of s in the domain, while smaller values of K means we can tolerate a wider size distribution. Regardless of the convergence of the series, we can always split off the resonant term involving n = 0.

If we do not make the assumption of a narrow size distribution, for n = 0 in Eq. (40) we get from the full integral: $\begin{array}{rl}& \mathrm{i}{ϵ}^2{A}_0\int \frac{1}{1+ϵ{\tilde{\omega }}_1-s}\mathrm{d}s\\ & =-\mathrm{i}{ϵ}^2{A}_0\int \frac{x-\mathrm{i}}{x^2+1}\mathrm{d}x\\ & =\mathrm{i}{ϵ}^2{A}_0{[-\frac{1}{2}\log (1+x^2)+\mathrm{i}{\tan }^{-1}x]}_{x_{min}}^{x_{max}},\end{array}$(49)

where $x\equiv \frac{1-s+ϵ\mathfrak{R}\left({\tilde{\omega }}_1\right)}{\mathfrak{J}\left({\tilde{\omega }}_1\right)}.$(50)

As an aside, note that the integral diverges in the limit where the imaginary part of ῶ₁ vanishes. This is why standard numerical integration techniques have difficulty evaluating integrals of this type, in particular for growth rates that are close to zero (Paardekooper et al. 2021).

Define a wide size distribution by taking s_max =1+Δ_R and s_min = 1 − Δ_L, with Δ_L,R ∼ 1: $\begin{array}{rl}\int \frac{1}{1+ϵ{\tilde{\omega }}_1-s}\mathrm{d}s=& -\frac{1}{2}\log \left(\frac{{ϵ}^2\mathfrak{J}{\left({\tilde{\omega }}_1\right)}^2+{(ϵ\mathfrak{R}\left({\tilde{\omega }}_1\right)-{\mathrm{\Delta }}_R)}^2}{{ϵ}^2\mathfrak{J}{\left({\tilde{\omega }}_1\right)}^2+{(ϵ\mathfrak{R}\left({\tilde{\omega }}_1\right)+{\mathrm{\Delta }}_L)}^2}\right)\\ & +\mathrm{i}{\text{\hspace{0.17em}tan}}^{-1}\left(\frac{ϵ\mathfrak{R}\left({\tilde{\omega }}_1\right)-{\mathrm{\Delta }}_R}{ϵ\mathfrak{J}\left({\tilde{\omega }}_1\right)}\right)\\ & -\mathrm{i}{\text{\hspace{0.17em}tan}}^{-1}\left(\frac{ϵ\mathfrak{R}\left({\tilde{\omega }}_1\right)+{\mathrm{\Delta }}_L}{ϵ\mathfrak{J}\left({\tilde{\omega }}_1\right)}\right)\\ =& \log \left(\frac{{\mathrm{\Delta }}_L}{{\mathrm{\Delta }}_R}\right)-\mathrm{i}\pi +O(ϵ).\end{array}$(51)

For a wide size distribution, therefore, Eq. (40) is O(ϵ²), and we must have ῶ₁ = 0, and we can subsequently neglect all ϵ terms under the integral: ${\tilde{\omega }}_2=\frac{\overline{\mu }}{2K}\int \frac{\mathrm{\Sigma }(s)(\mathrm{i}-K(1-s{)}^2)}{(1-s)(1-\mathrm{i}sK(1-s))}\mathrm{d}s-\frac{\mathrm{i}\overline{\mu }}{2K}.$(52)

The term proportional to K under the integral is due to the perturbation in relative velocity, and it is straightforward to show that the contribution of this term to the imaginary part of ω₂ is always negative. Together with the last term, representing the gas density perturbation, the relative velocity perturbation is therefore a stabilising effect. The dust density perturbation has to overcome both in order for the instability to grow.

For simplicity and avoidance of cluttered notation, we have only considered n = 0 in Eq. (40). It is straightforward to show, by defining x = 1+ϵ ῶ₁ − s: $\int \frac{(1-s{)}^n}{1+ϵ{\tilde{\omega }}_1-s}\mathrm{d}s=-\int \frac{{(x-ϵ{\tilde{\omega }}_1)}^n}{x}\mathrm{d}x,$(53)

and expanding the term in the numerator, that all terms are O(ϵ²) and that the conclusion that we need ῶ₁ = 0 holds in general.

If we split off the resonant term in Eq. (52), we get that ${\tilde{\omega }}_2=\frac{\overline{\mu }}{2K}\int \frac{\mathrm{i}\mathrm{\Sigma }(s)}{1-s}\mathrm{d}s-\frac{\overline{\mu }}{2}\int \frac{\mathrm{\Sigma }(s)}{1-\mathrm{i}sK(1-s)}\mathrm{d}s-\frac{\mathrm{i}\overline{\mu }}{2K},$(54)

or, integrating the first term as above, and including the full frequency: $\begin{array}{rl}\omega & =k{c}_{\mathrm{g}}+\frac{\mu {\sigma }^{(0)}\left({\tau }_{\text{res}}\right)}{2{\rho }_{\mathrm{d}}^{(0)}}[\mathrm{i}\text{\hspace{0.17em}log\hspace{0.17em}}\left(\frac{{\mathrm{\Delta }}_L}{{\mathrm{\Delta }}_R}\right)+\pi ]\\ & -\frac{\mu k{c}_{\mathrm{g}}}{2{\rho }_{\mathrm{d}}^{(0)}}\int \frac{{\sigma }^{(0)}\left({\tau }_{\mathrm{s}}\right)}{1-\mathrm{i}{\tau }_{\mathrm{s}}k{c}_{\mathrm{g}}(1-{\tau }_{\mathrm{s}}/{\tau }_{\text{res}})}\mathrm{d}{\tau }_{\mathrm{s}}-\frac{\mathrm{i}\mu }{2{\tau }_{\text{res}}},\end{array}$(55)

with resonant stopping time τ_res = c_g/α. The contribution of the resonance (the second term on the right hand side), is proportional to the size density at the resonant stopping time, therefore increasing the importance of the resonance if the size distribution is dominated by the resonant size. In order for the resonance to promote instability, we need log (Δ_L/Δ_R) > 0, so Δ_L/Δ_R > 1: the size distribution should not extend too far towards larger sizes. This particular effect is due to our assumption that Σ is constant: for more general power law size distributions there can be growing modes for Δ_R ≫ 1. These two issues make the ARDI sensitive to the exact size distribution: the importance of the resonant size, as well as the maximum and minimum size.

It is important to realise that the fact that ῶ₁ = 0 means that growth is now proportional to μ rather than $\sqrt{\mu }$ which is the hallmark of an RDI (Squire & Hopkins 2018b). Mathematically speaking, integrating over the singularity that is the resonance gives a finite result, essentially removing the resonant nature of the instability. Physically speaking, the dust density contribution to the instability is smaller in the polydisperse case, to the extent that it now has to compete with stabilising effects of velocity perturbations and gas density perturbations.

The remaining integral in Eq. (55) contains contributions from both the (non-resonant) density and relative velocity perturbation. While the integral can be solved in terms of a complex inverse tangent function, in practice it is easier to solve the integral numerically. The non-resonant contribution to the density perturbation is usually much smaller than the resonant contribution, unless Δ_R ≈ Δ_L. We will assume this is the case in what follows.

As in the case of a narrow size distribution, there is again an important wave number dependence: the stabilising effect of the velocity perturbations becomes stronger with increasing K. If we approximate the contribution of the relative velocity perturbation as $-\frac{\mathrm{i}K\mathrm{\Sigma }(1)}{2}\int \frac{s(1-s{)}^2}{1+s^2{K}^2(1-s{)}^2}\mathrm{d}s\approx -\frac{\mathrm{i}K\mathrm{\Sigma }(1)}{2}\int s(1-s{)}^2\mathrm{d}s,$(56)

then the imaginary part of ω₂ is, ignoring non-resonant contributions from the density perturbation, $\mathfrak{J}\left({\tilde{\omega }}_2\right)=\frac{\overline{\mu }\mathrm{\Sigma }(1)}{2K}\log \left(\frac{{\mathrm{\Delta }}_L}{{\mathrm{\Delta }}_R}\right)-\frac{\overline{\mu }K\mathrm{\Sigma }(1)}{2}\int s(1-s{)}^2\mathrm{d}s-\frac{\overline{\mu }}{2K},$(57)

which means there is a cut-off wave number $K^2<\frac{\mathrm{\Sigma }(1)\log \left(\frac{{\mathrm{\Delta }}_L}{{\mathrm{\Delta }}_R}\right)-1}{\mathrm{\Sigma }(1)\int s(1-s{)}^2\mathrm{d}s}=O(1),$(58)

which means that for K ≫ 1 we expect no growth. This is purely because the relative velocity perturbations provide stability at small scales.

In the limit of K ≪ μ we can find a polydisperse version of the long wavelength pressure-free mode, as identified in Squire et al. (2022). Rescaling the wavenumber to $K=\overline{K}{ϵ}^4$ , with $\overline{K}=O(1)$ , and at the same time setting $\tilde{\omega }=\overline{\omega }{ϵ}^3$ , we find from Eq. (26) that $ϵ{\overline{\omega }}^2={ϵ}^3{K}^2+\overline{\mu }\int \mathrm{\Sigma }(s)\frac{\mathrm{i}\overline{K}\overline{\omega }-{ϵ}^2\overline{\omega }(\overline{\omega }-ϵ\overline{K}s{)}^2}{(\overline{\omega }-ϵKs)(1-\mathrm{i}s{ϵ}^3(\overline{\omega }-ϵKs))}\mathrm{d}s-\mathrm{i}\overline{K}\overline{\mu }.$(59)

At lowest order in ϵ, this holds straightforwardly. At order ϵ, we need ${\overline{\omega }}^3=\mathrm{i}\overline{\mu }{\overline{K}}^2\int \mathrm{\Sigma }(s)s\mathrm{d}s=\mathrm{i}\overline{\mu }{\overline{K}}^2\frac{s_{max}+s_{min}}{2}.$(60)

A similar result was obtained by Squire et al. (2022) (their Appendix A).

Fig. 3 Growth rate for the backward propagating sound wave, with ω0 = −kcg, for μ = 0.001, smin = 0.5 and smax = 1.2. The solid blue curve represents the exact solution to Eq. (26), while the orange curve is the approximation (31) for μ ≪ 1, which overlaps with the blue curve everywhere except for K ≲ μ. The dotted line is the approximation (32) for small K.

3.3 Interpretation

We summarise the results of the previous subsections using Figures 2, 3, 4 and 5. We start with the non-resonant case in Figure 3. This is the case of a backward propagating soundwave with ῶ₀ = −K going unstable. We focus on μ = 0.001, and a constant size distribution with 0.5 < s < 1.2. The exact result, which was obtained by solving Eq. (26), indicated by the blue curve, shows growth for K < 1 and stability for K > 1 (see also Squire et al. 2022). The approximation (31) for μ ≪ 1, indicated by the orange curve, is indistinguishable from the exact result except for K ≲ μ, which is where our ordering scheme becomes invalid. The limit (32) of K ≪ 1 does a good job of predicting the plateau of maximum growth. The maximum growth rate is comparable to the monodisperse result for s = 1, in which case ῶ₁ = i/4 for K ≪ 1. However, in the monodisperse case, this growth rate of O(μ) is usually dwarfed by the resonant case which has growth of $O(\sqrt{\mu })$ . In the polydisperse case, however, this mode tends to dominate because of the reduced growth rates of the resonant modes (Squire et al. 2022), which we explore below.

In the previously discussed Figure 2, we consider narrow size distributions for a dust-to-gas ratio μ = 0.01. For small K, the results follow the monodisperse result of Eq. (46). We note that for very small wave numbers, the growth rate starts to differ from Eq. (46), which represents the regime of the ‘long wavelength pressure-free mode’ (Hopkins & Squire 2018; Magnan et al. 2024a). This regime occurs for K ≪ μ, and occurs most readily in the quasi-monodisperse case given by Eq. (47). For a wide size distribution, it is hard to enter the K ≪ μ regime without entering the quasi-monodisperse regime at the same time (see also below).

Towards larger K, growth disappears as we leave the quasimonodisperse regime. In this simplified system, this means that the width of the size distribution actually sets the maximum growth rate. In reality, there is probably competition from diffusive processes killing off growth at high wave numbers. We briefly look at the effects of diffusion in Section 3.4. As noted in Hopkins & Squire (2018), subsonic drift also leads to a cut-off at high wave numbers, see Eq. (29). For μ = δ = 0.01, the cut-off occurs at a similar wave number as for monodisperse subsonic drift with s ≈ 0.9.

In Figure 4, we show the growth rate at fixed K = 1 and a size distribution that is constant between s_min = 0.875 and s_max = 1.05. These limits were chosen to get the transition between quasi-monodisperse and the polydisperse regime at reasonable wave numbers. When μ is large enough, the quasi-monodisperse condition (47) is satisfied, and we recover a growth rate $\propto \sqrt{\mu }$ , characteristic of resonant drag instabilities. Although the size distribution is asymmetric around s = 1, if we take δ = 0.0875 as an average width, from Eq. (47) we expect the quasi-monodisperse regime to hold for μ ≫ 0.03. At lower values of μ, we leave the quasi-monodisperse regime and thereby the dominance of the resonance, which leads to a growth rate ∝μ.

In Figure 5, we show results for a wide size distribution, which is constant between s_min = 0.5 and s_max = 1.2. We note that we need the size distribution to be asymmetric around the resonant size s = 1, otherwise the second order contribution from the dust density perturbation vanishes, see Eq. (55). This is also where the shape of the size distribution can make a difference. At large wave numbers K ≳ 1, growth vanishes, see Eq. (58). The plateau in growth towards the right side of the curves, most notably for μ = 10⁻⁵, is the polydisperse mode, and it can be readily verified that growth in this regime is indeed ∝μ. The size distributions considered in Figure 2 are too narrow for polydisperse growth to show up. The dotted curves show the relation Eq. (55), which does a good job in predicting the plateau and the cut-off at high wave numbers. The bump towards the left of the curves is the tail of the quasi-monodisperse regime, which is not part of the derivation of Eq. (55). If we take δ = 0.35 for this size distribution, from Eq. (47) we expect the quasimonodisperse regime to hold for K ≪ 2μ, consistent with what is seen in Figure 5.

In summary, the polydisperse version of the ARDI is different from its monodisperse counterpart in two important ways. First, a distribution of dust sizes lead to a cutoff of growth at high wave numbers, even for very narrow size distributions. Second, for wider size distributions, there exists an intermediate range of wave numbers, where the resonance plays only a minor role, as the integration over size dilutes its contribution. Here, the growth rates scale with μ, rather than $\sqrt{\mu }$ , which is the case for monodisperse RDIs. In other words, a wide size distribution converts resonant growth into ordinary growth, with usually reduced growth rates for μ ≪ 1. At smaller wave numbers, we find the quasi-monodisperse regime, which behaves similar to the monodisperse ARDI.

Fig. 4 Growth rate of the acoustic resonant drag instability as a function of the dust-to-gas ratio μ, calculated numerically from Eq. (26). The dimensionless wavenumber K = 1, and the size density is constant between smin = 0.875 and smax = 1.05. The green dotted line indicates the monodisperse result (46), while the orange dotted line indicates the polydisperse limit (55).

Fig. 5 Growth rate of the acoustic resonant drag instability as a function of the wave number for different dust-to-gas ratios μ, calculated numerically from Eq. (26). The size distribution is a constant between smin = 0.5 and smax = 1.2. The dotted curves indicate the approximation (55).

3.4 Acoustic drag instability with diffusion

Before we conclude, we briefly look at diffusive effects by including gas viscosity and dust diffusion. We only consider the monodisperse version of the ARDI. It is expected that diffusion stabilises the instability at large wave numbers, and it is interesting to compare the diffusive cut-off to the cut-off due to having a size distribution. Realistic levels of viscosity and associated dust diffusion will strongly depend on the background state, in particular whether the underlying gas flow is turbulent (see e.g. Woitke 2006, for the case of AGB star winds). Here, we do not aim at modelling a particular physical situation in terms of the input level of viscosity, but rather ask the question at what level of viscosity do dissipative effects associated with gas viscosity and dust diffusion compete with effects due to the finite with of the size distribution discussed in the previous subsections.

The monodisperse equations, including gas viscosity and dust diffusion, read, adopting the formulation of Lin (2021): ${\mathrm{\partial }}_t{\rho }_{\mathrm{g}}+{\mathrm{\partial }}_x\left({\rho }_{\mathrm{g}}{v}_{\mathrm{g}}\right)=0,$(61) $\begin{array}{rl}{\mathrm{\partial }}_t{v}_{\mathrm{g}}+v_{\mathrm{g}}{\mathrm{\partial }}_x{v}_{\mathrm{g}}=& -\mu \alpha -\frac{c_{\mathrm{g}}^2{\mathrm{\partial }}_x{\rho }_{\mathrm{g}}}{{\rho }_{\mathrm{g}}}+\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{u-v_{\mathrm{g}}}{{\tau }_{\mathrm{s}}}\\ & +\frac{4}{3{\rho }_{\mathrm{g}}}v{\mathrm{\partial }}_x\left({\rho }_{\mathrm{g}}{\mathrm{\partial }}_x{v}_{\mathrm{g}}\right),\end{array}$(62) ${\mathrm{\partial }}_t{\rho }_{\mathrm{d}}+{\mathrm{\partial }}_x\left({\rho }_{\mathrm{d}}u\right)=D{\mathrm{\partial }}_x\left(({\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}}){\mathrm{\partial }}_x\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{d}}+{\rho }_{\mathrm{g}}}\right)\right),$(63) ${\mathrm{\partial }}_{t}u+u{\mathrm{\partial }}_{x}u=\alpha -\frac{u-v_{\mathrm{g}}}{{\tau }_{\mathrm{s}}},$(64)

where v is the gas kinematic viscosity and D the dust diffusion coefficient. The background flow is constant in space and time, and exactly the same as the inviscid case. Linear perturbations are governed by: ${\mathrm{\partial }}_t\mathcal{G}+{\mathrm{\partial }}_x{v}_{\mathrm{g}}^{(1)}=0,$(65) $\begin{array}{rl}{\mathrm{\partial }}_t{v}_{\mathrm{g}}^{(1)}=& -c_{\mathrm{g}}^2{\mathrm{\partial }}_x\mathcal{G}+\mu \frac{u^{(1)}-v_{\mathrm{g}}^{(1)}}{{\tau }_{\mathrm{s}}}\\ & +\mu (\mathcal{D}-\mathcal{G})\frac{u^{(0)}}{{\tau }_{\mathrm{s}}}+\frac{4v}{3}{\mathrm{\partial }}_x^2{v}_{\mathrm{g}}^{(1)},\end{array}$(66) ${\mathrm{\partial }}_t\mathcal{D}+{\mathrm{\partial }}_x{u}^{(1)}+u^{(0)}{\mathrm{\partial }}_x\mathcal{D}=\frac{D}{\mu +1}{\mathrm{\partial }}_x^2(\mathcal{D}-\mathcal{G}),$(67) ${\mathrm{\partial }}_t{u}^{(1)}+u^{(0)}{\mathrm{\partial }}_x{u}^{(1)}=-\frac{u^{(1)}-v_{\mathrm{g}}^{(1)}}{{\tau }_{\mathrm{s}}},$(68)

with $\mathcal{G}={\rho }_{\mathrm{g}}^{(1)}/{\rho }_{\mathrm{g}}^{(0)}$ and $\mathcal{D}={\rho }_{\mathrm{d}}^{(1)}/{\rho }_{\mathrm{d}}^{(0)}$ . Taking perturbations ∝ exp (ikx − iωt) as before: $k\hat{v}=\omega \hat{\mathcal{G}},$(69) $c_{\mathrm{g}}^{2}k\hat{\mathcal{G}}+\mathrm{i}\mu \frac{\hat{u}-\hat{v}}{{\tau }_{\mathrm{s}}}+\mathrm{i}\mu \alpha (\hat{\mathcal{D}}-\hat{\mathcal{G}})-\frac{4\mathrm{i}v}{3}{k}^2\hat{v}=\omega \hat{v},$(70) $k\hat{u}+\alpha {\tau }_{\mathrm{s}}k\hat{\mathcal{D}}-\frac{\mathrm{i}D}{\mu +1}{k}^2(\hat{\mathcal{D}}-\hat{\mathcal{G}})=\omega \hat{\mathcal{D}},$(71) $\alpha {\tau }_{\mathrm{s}}k\hat{u}-\mathrm{i}\frac{\hat{u}-\hat{v}}{{\tau }_{\mathrm{s}}}=\omega \hat{u},$(72)

This is a straightforward eigenvalue problem with matrix: $\mathrm{A}=\left(\begin{array}{cccc}0& k& 0& 0\\ c_{\mathrm{g}}^{2}k-\mathrm{i}\mu \alpha & -\frac{\mathrm{i}\mu }{{\tau }_{\mathrm{s}}}-\frac{4\mathrm{i}{k}^2}{3}& \mathrm{i}\mu \alpha & \mathrm{i}\mu /{\tau }_{\mathrm{s}}\\ \frac{\mathrm{i}D{k}^2}{1+\mu }& 0& \alpha {\tau }_{\mathrm{s}}k-\frac{\mathrm{i}D{k}^2}{1+\mu }& k\\ 0& \mathrm{i}/{\tau }_{\mathrm{s}}& 0& \alpha {\tau }_{\mathrm{s}}k-\frac{\mathrm{i}}{{\tau }_{\mathrm{s}}}\end{array}\right).$(73)

We can construct a non-dimensional version, using a length scale $c_{\mathrm{g}}^2/\alpha$ and timescale c_g/α: $K\hat{V}=\mathrm{\Omega }\hat{\mathcal{G}},$(74) $K\hat{\mathcal{G}}+\mathrm{i}\mu \frac{\hat{U}-\hat{V}}{s}+\mathrm{i}\mu (\hat{\mathcal{D}}-\hat{\mathcal{G}})-\frac{4\mathrm{i}{\alpha }_{\mathrm{g}}}{3}{K}^2\hat{V}=\mathrm{\Omega }\hat{V},$(75) $K\hat{U}+sK\hat{\mathcal{D}}-\frac{\mathrm{i}{\alpha }_{\mathrm{d}}}{\mu +1}{K}^2(\hat{\mathcal{D}}-\hat{\mathcal{G}})=\mathrm{\Omega }\hat{\mathcal{D}},$(76) $sK\hat{U}-\mathrm{i}\frac{\hat{U}-\hat{V}}{s}=\mathrm{\Omega }\hat{U},$(77)

with nondimensional diffusion coefficients² ${\alpha }_{\mathrm{g}}=v\alpha /c_{\mathrm{g}}^3$ and ${\alpha }_{\mathrm{d}}=D\alpha /c_{\mathrm{g}}^3$ .

The resulting growth rates for μ = 0.01 are shown in Figure 6, for three different levels of diffusion. In all cases, we focus on the resonant case with s = 1. For simplicity, we keep α_g = α_d, noting that dust diffusion is far more important than gas viscosity for damping the unstable modes. This choice corresponds to a Schmidt number of unity. The levels of viscosity were chosen so that the cutoff wave numbers roughly match those of Figure 2. We see that a size distribution of width δ = 0.1 leads to a similar cutoff wavenumber as a dimensionless diffusion coefficient of α_g = α_d = 0.005. In other words, for levels of diffusion lower than α_g = α_d = 0.005, it is the width of the size distribution that determines the maximum growth rate, rather than diffusion. Narrower size distributions need lower levels of diffusion to be dominant.

Fig. 6 Growth rate of the monodisperse acoustic resonant drag instability as a function of the wave number, including gas viscosity and dust diffusion, calculated from the eigenvalues of Eq. (73). The dust-to-gas ratio μ = 0.01.

4 Discussion and conclusion

We have presented linear calculations of the acoustic resonant drag instability in the polydisperse regime. We have shown that a quasi-monodisperse regime exists, with growth rates comparable to the resonant monodisperse case, if the relative width of the size distribution around the resonant size is small: $\delta =\frac{{\tau }_{\mathrm{s},\text{max}}-{\tau }_{\mathrm{s},\text{min}}}{{\tau }_{\text{res}}}\ll \frac{1}{2}\sqrt{\frac{\mu }{k{c}_{\mathrm{g}}{\tau }_{\text{res}}}}.$(78)

There is a straightforward correspondence between the analysis presented here, and the discussion on detuning of the resonance in Magnan et al. (2024a), their appendix A. They studied how the growth rate varies just away from the resonance, possibly because the resonant modes do not fit exactly in the numerical domain. Our relative width of the size distribution, δ, plays the role of the detuning parameter. A narrow size distribution corresponds to small detuning, while a wide size distribution corresponds to large detuning. The effect of this ‘detuning’ is strongest for the highest wave numbers, which is also where the monodisperse ARDI grows fastest. The fastest growing monodisperse modes require extremely narrow size distributions to survive in the polydisperse case.

For wider size distributions, genuinely polydisperse modes exist. Realistic widths of the dust size distribution can, for example, be based on the ISM: Squire et al. (2022) take a minimum size of 5 nm and a maximum size of 0.25 μm, based on the work of Mathis et al. (1977). This gives a spread in sizes, and, for our simple drag law, stopping times of two orders of magnitude, or δ ≈ 1, definitely outside the quasi-monodisperse regime. Moreover, the true ISM size distribution may well include larger sizes (e.g. Weingartner & Draine 2001). In dynamical environments such as AGB star winds it is more difficult to establish a size distribution, but hydrodynamical simulations including dust nucleation find that dust particles can grow rapidly over at least an order of magnitude (Simis et al. 2001), suggestive of a wide size distribution in these environments as well.

Even though growing polydisperse modes exist, their growth rates are ∝μ rather than $\propto \sqrt{\mu }$ for the corresponding monodisperse case, which means they usually grow slower for μ ≪ 1. Integrating over the resonance removes the resonant nature of the instability, but given the right size distribution, the contribution of the resonance can still promote growth. For the simple, constant size distribution considered here, we found that the size distribution should be asymmetric on either side of the resonant size.

We have focused on the case of a continuous size distribution. That is, we have never explicitly discretised the integral appearing in the backreaction term in the gas momentum equation. As a result, the dispersion relation (26) contains an integral, which makes numerically solving for the frequencies not easy. A different approach is to discretise the integral from the start, and solve the resulting matrix eigenvalue problem with standard techniques (McNally et al. 2021; Squire et al. 2022). In the limit where the number of dust species N goes to infinity, one expects the frequencies found to converge to the results obtained from Eq. (26). However, as noted in McNally et al. (2021) for the streaming instability, if the integral is nearly singular, this approach of discretising the integral leads to extremely slow convergence. In Figure 7, we illustrate the problem for the acoustic drag instability. The orange curves indicate the growth rates for the non-resonant backward propagating sound wave, and here we find that for N < 100 the results are indistinguishable from the continuum case. In the case where the integration domain includes the resonance, which happens for the blue curves that have ῶ₀ = K, convergence is much harder to obtain. The top panel shows that we reach convergence for μ = 0.0004 at N < 1024, but if we lower the dust-to-gas ratio to 0.0001 the result does not converge for N up to 2048. Over the whole range of N considered, the fastest growing mode is not the backward propagating sound wave, but rather a spuriously growing resonant wave. We note that, for the relatively narrow size distribution chosen here, we need the dust-to-gas ratio to remain small or we enter the quasi-monodisperse regime (see Figure 5). These spuriously growing modes could pose a problem for numerical simulations, for which it is usually too expensive to include 1000s of dust species. Even more, when we increase K so that there should be no growing modes in the polydisperse regime (see Figures 3 and 5), the direct matrix method finds spurious growing modes with growth rates $\mathfrak{J}(\tilde{\omega })\sim \mu$ up to N = 2048. These can easily come to dominate a numerical simulation.

Several important simplifications were made to make the problem tractable. First of all, we have only considered very simple (i.e. constant) size distributions. While most results are likely to be insensitive to the exact form of the size distribution (i.e. the existence of the quasi-monodisperse regime and polydisperse modes), in some cases details will matter. In particular, the contribution of the resonance to the total growth rate depends on the limits of the size distribution, see Eq. (55). More complicated size distributions, either inherited from the ISM or self-consistently generated through nucleation and collisional growth, are likely to affect the growth rates of the acoustic resonant drag instability.

We have considered only a very simple drag law with constant stopping time. A more general form of Epstein drag was given in Hopkins & Squire (2018), where the stopping time depends both on gas density and relative velocity between gas and dust. This form was also used in Squire et al. (2022). For the fast drift necessary to trigger the ARDI, in particular the dependence of the stopping time on relative velocity will be important. A comparison in the monodisperse case was made in Magnan et al. (2024a), but for the polydisperse case there is still work to be done.

Despite the differences in setup, our results are in agreement with Appendix A3 of Squire et al. (2022). In particular the existence of the unstable backward propagating sound wave for small K, and its dominance over the resonant mode (see Figure 7). With our continuum approach, we were able to study the resonant modes in more detail. Apart from the constant stopping time mentioned above, another difference with the work of Squire et al. (2022) is that we have worked in one spatial dimension only, greatly simplifying the analysis. As was mentioned in Magnan et al. (2024a), the physics of the instability is well-captured by a one-dimensional approach.

In conclusion, for realistically wide size distributions, the acoustic resonant drag instability mostly shows non-resonant growth in the limit μ ≪ 1. The highest wave number modes are rendered stable by having a size distribution, which takes the maximum growth rate down substantially, much like diffusion. Therefore, the acoustic resonant drag instability survives in the polydisperse regime, but with growth rates that are typically 𝔍(ω) ∼ μα / c_g, which is a fraction $\sim \sqrt{\mu /K}$ of the monodisperse growth rate.

Fig. 7 Convergence of the direct method with number of dust species, for a case with smin = 0.5 and smax = 1.2, all at K = 0.1. Top panel: μ = 0.0004, bottom panel: μ = 0.0001. The blue curves indicate a mode with unperturbed frequency ῶ0 = K, which contains the resonance, and the orange curves a non-resonant mode with ῶ0 = −K. The dotted lines indicate the result of Eq. (26).

Data availability

All data used to create figures in this work are available at https://doi.org/10.4121/df27e782-0a37-43a9-8fdd-bd4081b96dfb

Acknowledgements

We thank the anonymous referee for a thorough and insightful report. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (Grant Agreement No. 101054502). This work made use of several open-source software packages. We acknowledge numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020), and psitools (McNally et al. 2021). We acknowledge 4TU.ResearchData for supporting open access to research data.

References

All Figures

Fig. 1 Schematic representation of the quasi-monodisperse regime. The top half of the figure deals with the gas, while the bottom half deals with the dust component. The horizontal axis indicates (wave) speed. The red solid arrow denotes the propagation speed of a gas sound wave in the absence of dust, and the blue arrow indicates the resonant dust drift speed. The backreaction of the dust modifies the propagation speed of the sound wave, leading to the dashed red arrow. The blue rectangle indicates a range of drift speeds if there is a distribution of dust sizes. As long as this rectangle does not come close to the dashed red arrow, we are in the quasi-monodisperse regime.

In the text

	Fig. 2 Growth rate of the acoustic resonant drag instability as a function of the wave number, which was obtained by numerically solving Eq. (26). The dust-to-gas ratio μ = 0.01, and the size density is constant, centred on the resonant size s = 1, with width δ. The red dotted line indicates the monodisperse result, while the vertical dotted lines indicate the cutoff of Eq. (47).
In the text

	Fig. 3 Growth rate for the backward propagating sound wave, with ω0 = −kcg, for μ = 0.001, smin = 0.5 and smax = 1.2. The solid blue curve represents the exact solution to Eq. (26), while the orange curve is the approximation (31) for μ ≪ 1, which overlaps with the blue curve everywhere except for K ≲ μ. The dotted line is the approximation (32) for small K.
In the text

	Fig. 4 Growth rate of the acoustic resonant drag instability as a function of the dust-to-gas ratio μ, calculated numerically from Eq. (26). The dimensionless wavenumber K = 1, and the size density is constant between smin = 0.875 and smax = 1.05. The green dotted line indicates the monodisperse result (46), while the orange dotted line indicates the polydisperse limit (55).
In the text

	Fig. 5 Growth rate of the acoustic resonant drag instability as a function of the wave number for different dust-to-gas ratios μ, calculated numerically from Eq. (26). The size distribution is a constant between smin = 0.5 and smax = 1.2. The dotted curves indicate the approximation (55).
In the text

	Fig. 6 Growth rate of the monodisperse acoustic resonant drag instability as a function of the wave number, including gas viscosity and dust diffusion, calculated from the eigenvalues of Eq. (73). The dust-to-gas ratio μ = 0.01.
In the text

	Fig. 7 Convergence of the direct method with number of dust species, for a case with smin = 0.5 and smax = 1.2, all at K = 0.1. Top panel: μ = 0.0004, bottom panel: μ = 0.0001. The blue curves indicate a mode with unperturbed frequency ῶ0 = K, which contains the resonance, and the orange curves a non-resonant mode with ῶ0 = −K. The dotted lines indicate the result of Eq. (26).
In the text