text:
\section*{Acknowledgments}\end{small}{\noindent #1}\vspace{10pt}}
\newcommand\altaffilmark[1]{$^{#1}$}
\newcommand\altaffiltext[1]{$^{#1}$}
\newcommand\wsj{w_{s,j}}
\newcommand\ws{w_{s}}
\newcommand\ad{\epsilon_{\rm grain}}
\newcommand\ab[1]{\epsilon_{{\rm gr},{#1}}}
\newcommand\abc[1]{\epsilon_{{\rm gr},{#1}}}
\newcommand\amin{\epsilon_{\rm grain}^{\rm min}}
\newcommand\amax{\epsilon_{\rm grain}^{\rm max}}
\newcommand\adj{\epsilon_{{\rm grain},j}}
\newcommand\tsj{t_{s,j}}
\newcommand\ddt[1]{\frac{\partial #1}{\partial t}}
\newcommand\aext{{\rm a}^{\rm ext}}
\newcommand\baext{{\bf a}^{\rm ext}}

\newcommand\hs{HS18}
\newcommand\msh{MSH19}

\newcommand\tfigh{4.7cm}
\voffset=-0.6in

\title[The RDI with a spectrum of grains]{The Acoustic Resonant Drag Instability with a Spectrum of Grain Sizes\vspace{-0.5cm}}

\author[Squire et al.]{
\parbox[t]{\textwidth}{
	Jonathan Squire\altaffilmark{1}, Stefania Moroianu\altaffilmark{2}, \&\ Philip~F.~Hopkins\altaffilmark{3}
}
\vspace*{6pt} \\
\altaffiltext{1}{Physics Department, University of Otago, Dunedin 9010, New Zealand} \\
\altaffiltext{2}{Department of Applied Physics, Stanford University, Stanford, CA 94305, USA}\\
\altaffiltext{3}{TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA\vspace{-0.3cm}}
}

\date{Submitted to MNRAS, 2021\vspace{-0.6cm}}
\begin{document}
\maketitle

\begin{abstract}
We study the linear growth and nonlinear saturation of the ``acoustic Resonant Drag Instability'' (RDI) when the dust grains, which drive the instability, have a wide, continuous spectrum of different  sizes. This physics is generally applicable to dusty winds driven by radiation pressure, such as occurs around red-giant stars, star-forming regions, or active galactic nuclei. Depending on the physical size of the grains compared to the wavelength of the radiation field that drives the wind, two qualitatively different  regimes emerge. In the case of  grains that are larger than the radiation's wavelength -- termed the \emph{constant-drift} regime -- the grain's equilibrium drift velocity through the gas is approximately independent of grain size, leading to strong correlations between differently sized grains that persist well into the saturated nonlinear turbulence.  For grains that are smaller than the radiation's wavelength  -- termed the \emph{non-constant-drift} regime -- the linear instability grows more slowly than the single-grain-size RDI and only the larger grains exhibit RDI-like behavior in the saturated state. A detailed study of grain clumping and grain-grain collisions shows that outflows in the constant-drift regime may be effective sites for grain growth through collisions, with large collision rates but low collision velocities.    \end{abstract}

\begin{keywords}
instabilities -- turbulence -- ISM: kinematics and dynamics -- galaxies: formation -- stars: winds, outflows -- dust, extinction
\vspace{-1.0cm}
\end{keywords}

\vspace{-1.1cm}

\section{Introduction}

Cosmic dust is ubiquitous across the universe and vital to a wide range of astrophysical processes. By mass, it makes
up around $\sim\!1\%$ of the  interstellar medium (ISM) of galaxies, but its strong coupling to radiation
fields implies it can nonetheless strongly influence gas dynamics and cooling in many situations   \citep{Draine2010}.
More generally, because around half of the metal content of our galaxy is locked up in dust, it plays
crucial roles in any process that requires metals or solids \revchng{\citep{Whittet1992,Draine2003}}. Notably, dust is almost certainly the key ingredient
for planet formation and life, supplying the necessary reservoir of solids that provide the seeds to make
planetesimals in protostellar disks \revchng{\citep{Chiang2010}}.

This paper deals with the physics of dust moving through gas, with the interaction between the species mediated by drag forces. Such conditions occur, for example, in dust-radiation-pressure driven winds, where an outflow of dusty gas
 is driven by an anisotropic radiation field that couples strongly to the dust. Such outflows are thought
 to be important in  the evolution of asymptotic giant-branch (AGB) stars (which also produce large quantities
 of dust; e.g., \citealp{Habing1996,Norris2012}), in feedback processes that help regulate star  formation
 and/or active-galactic nuclei \citep{Scoville2003,Thompson2005,Ishibashi2015}, around supernovae \citep[e.g.,][]{Micelotta2018},
 and in the bulk ISM \citep{Weingartner2001a} and circum-galactic medium (CGM; \citealp{Menard2010}).
As shown by \cite{Squire2018}, this situation -- specifically, when the radiation pressure on the dust is stronger than that on the gas, such that
 the gas outflow is driven indirectly through the drag force from the dust -- is unstable to the ``{Resonant Drag Instability}'' (RDI): small perturbations in the gas or dust will grow in time until they become large, driving turbulence  in the gas
 and strong dust clumping.  \cite{Hopkins2018} (hereafter\defcitealias{Hopkins2018}{\hs}\hs) studied the linear features and growth rates of the RDI
 for the case of nearly neutral grains and neutral gas in outflows, while \cite{Hopkins2018a} generalized these
 results to charged grains in magnetized gas.
 These results were then extended to the nonlinear regime by  \cite{Moseley2019} (hereafter\defcitealias{Moseley2019}{\msh}\msh) and  \citealt{Seligman2019,Hopkins2020} (in the magnetized regime), who
 studied the turbulence induced by the RDI, constructing simple estimates for
 its saturation amplitude and other properties.

 However, each of these studies  has allowed for the dust grains to have only
 a single size.
\revchng{ As shown by \cite{Krapp2019,Paardekooper2020,Zhu2021} for} the streaming instability in protoplanetary disks (which is
 part of the RDI family; \citealp{Squire2018a}), a range of grain sizes can
 strongly influence the instabilities in some regimes.   Thus, for more realistic application to astrophysical fluids and
 outflows, where the
 range in grain sizes easily spans two orders of magnitude or more  \citep{Draine2010},
 we must relax the single-grain-size assumption and better understand the influence
 of a spectrum of grain sizes on the growth rate and saturation of the RDI. This is the first purpose of the present paper.
 We study both the linear growth rates and nonlinear saturation of the ``acoustic RDI'' (\hs) that applies to uncharged
 grains and neutral gas, and involves the driving of compressive shocks and sound waves
 by the drifting dust. We find that depending on whether grains are smaller or larger
 than the wavelength of the accelerating radiation field, the presence of a spectrum of grains either has
 little effect on the acoustic RDI   or reduces the growth rate and saturation of smaller-scale motions.
 In both cases,  key features of the linear-instability structure persist well into the highly turbulent saturated state.

 The second purpose of this paper is to better understand  dust clumping and collisions in
 RDI turbulence (i.e., the saturated state of the \revchng{acoustic RDI}). This is important because outflows are highly
 dynamic and  often thought to be key sites for grain condensation, coagulation, and fragmentation, and the
  latter two of these processes are strongly influenced by turbulence. While well
 developed theories exist to describe the how standard gas turbulent motions influence dust collisions
 \cite[e.g.,][]{Ormel2007,Pan2013,Pumir2016}, the structure of the turbulence and clumping driven by the
acoustic RDI is quite different to standard turbulence in many ways, because the instability operates
across all scales of the system simultaneously (\msh). In this context, it is important to consider the RDI with a
wide spectrum of grain sizes (as opposed to the single-grain-size RDI) because the nature of the instability suggests that there
could exist interesting
correlations between differently sized grains in some regimes, and grain clumping and collision  statistics
 depend strongly on grain sizes \citep{Pan2014a,Blum2018,Mattsson2019}.
 With this in mind, our nonlinear study is designed to compare RDI saturation to forced turbulence with passively advected dust. We  do this by designing ``equivalent'' forced turbulence simulations  (i.e., simulations with parameters
 chosen to match the saturated RDI as closely as possible), allowing an explicit comparison of  the statistics of dust in
 RDI turbulence with those of passive dust in forced turbulence.
 Given the rather detailed nature of these comparisons, we focus on just two RDI case
 studies at the high numerical resolution, but with parameters that can be applicable
 to a range of astrophysical scenarios. Depending on the regime, we find that the RDI involves
 a  significantly faster rate of  lower-velocity collisions than forced turbulence, particularly  between
 grains of different sizes. It also exhibits far stronger clumping of the smallest grains, even when
 simple estimates suggest these small grains should be very well coupled to the gas.

The paper is split in two: first the main exposition; second an extended appendix that
studies analytically the linear behavior of the acoustic RDI with a spectrum of grain sizes. This split was chosen because
the linear calculations, in which we derive simple analytic expressions for the RDI growth rate
in most of the important regimes,  are necessarily rather technical. They do provide useful
understanding of the nonlinear behavior, however, so we will refer to those results throughout.
The main paper starts in  \cref{sec: rdi description} with a detailed description
of the problem, model, and numerical setup, particularly focusing  on
important differences that arise in the dust (quasi-)equilibrium depending on the wavelength
of the radiation compared to the dust size (\cref{sub: acceleration} and \cref{tab: regime}). The numerical results are presented in \cref{sec: results},
starting with a discussion of the broad morphological features of the turbulence and how this differs between
regimes and driving (\cref{sub: results: time evolution}), then followed by a more detailed analysis of
the grain clumping and collisions (\cref{sub: results: statistics,sub: results: collisions}).

\section{Numerical model and physical setup}\label{sec: rdi description}

\subsection{Dust and gas model}\label{sub: basic model}

We model gas dynamics using the standard neutral fluid equations,
ignoring magnetohydrodynamic effects and charged grains for simplicity in this study, although such effects are important
for many astrophysical regimes (\hs). Dust is modelled numerically
by treating it as a population of individual particles (the super-particle approach), which interact with the
gas through drag forces that depend on the grain size.
We use $f_{d}(\ad;\bm{x},\bm{v},t)$ to denote the phase-space density of grains with radius $\ad$, velocity $\bm{v}$, at position
$\bm{x}$, such that the equations of motion are
\begin{gather}
\ddt{\rho_{g}} + \nabla\cdot(\rho_{g}\bm{u}_{g}) = 0,\label{eq: dt cont}\\
\ddt{\bm{u}_{g}} + \bm{u}_{g}\cdot\nabla\bm{u}_{g} = -c_{s}^{2}\frac{\nabla\rho_{g} }{\rho_{g}}+ \frac{1}{\rho_{g}}\int d\ad d\bm{v}\, f_{d}(\bm{x},\bm{v},a) \frac{\bm{v}-\bm{u}_{g}}{t_{s}(\ad,\bm{v})},\label{eq: dt mom}\\
\ddt{f_{d}} + \bm{v}\cdot\nabla f_{d} + \frac{\partial }{\partial\bm{v}} \cdot \left[ \left( -\frac{\bm{v}-\bm{u}_{g}(x)}{t_{s}(\ad,\bm{v})}+ \baext (\ad)\right) f_{d}\right]=0.\label{eq: f dt}
\end{gather}
Given the super-particle approach, the final equation is equivalent to
modeling individual grains of size $\adj$, with  velocity $\bm{v}_{j}$ and position $\bm{x}_{j}$, with
\begin{equation}
\ddt{\bm{x}_{j}}=\bm{v}_{j},\quad\ddt{\bm{v}_{j}} = -\frac{\bm{v}_{j} - \bm{u}_{g}(\bm{x}_{j})}{t_{s}(\adj,\bm{v}_{j})} + \baext (\adj)\label{eq: dt v indiv}
\end{equation}
as they move through the gas velocity field, then constructing $f_{d}$ by counting dust particles at each gas position.
Here $\rho_{g}$ is the gas density, $\bm{u}_{g}$ is the gas velocity, $c_{s}$ is the sound speed (the equation of state is taken as isothermal), $t_{s}(\ad,\bm{v})$ is the stopping time (see \cref{eq: stopping time epstein} below), and $\baext$ is an external force from radiation pressure on the grains, which we will take (arbitrarily) to be in the $\hat{\bm{z}}$ direction $\baext=\aext\hat{\bm{z}}$ (see \cref{sub: acceleration}). The
final term in \cref{eq: dt mom} is the dust ``backreaction'' force -- it the force on the gas from the dust -- which
is neglected in most studies of dust dynamics. We also use $\langle\cdot\rangle$ to denote the volume average, and
the subscripts $\perp$ and $\|$  to denote velocities in the directions perpendicular and parallel to the mean drift, respectively (e.g., $\langle u_{g,\|}^{2}\rangle=\langle u_{g,z}^{2}\rangle$, $\langle u_{g,\perp}^{2}\rangle=\langle u_{g,x}^{2}+u_{g,y}^{2}\rangle$). A $\delta$ indicates that the mean is subtracted from a quantity before averaging (e.g., $ \delta v_{d,\|}=v_{d,\|}-\langle v_{d,\|}\rangle$). All quantities will be measured in the frame where the gas is stationary $\langle \bm{u}_{g}\rangle=0$.

To understand dust dynamics, analyze simulations, and compute linear growth rates, it is  helpful to take the zeroth and first velocity moments
of $f_{d}$, defining
\begin{equation}
\rho_{d}=\int_{\abc{l}}^{\abc{h}} d\ad d\bm{v}\, f_{d},\qquad \bm{v}_{d} = \frac{1}{\rho_{d}}\int_{\abc{l}}^{\abc{h}} d\ad d\bm{v}\,\bm{v} f_{d},\label{eq: vd rhod defs}
\end{equation}
where the integration limits $\abc{l}$ and $\abc{h}$ can be taken over just a subset of grains (if specified as such; see \cref{sub: dust binning}), or the full distribution (if unspecified).
We use $\mu_{0}$ to denote the
total average dust-to-gas-mass ratio $\mu_{0}=\langle \rho_{d}\rangle/\langle \rho_{g}\rangle$,
and
\begin{equation}
w_{s}(\ad) = \frac{\langle |\bm{v}_{d}(\ad)-\bm{u}_{g}|\rangle}{c_{s}},\label{eq: ws def}
\end{equation}
to denote the mean drift between dust of size $\ad$ and the gas (here $\bm{v}_{d}(\ad)$ is computed from $\abc{l}=\ad$,
$\abc{h}=\ad+d\ad$ in \cref{eq: vd rhod defs}).

Throughout, we assume Epstein drag, using the simple approximation of  \citep{Draine1979a},
\begin{equation}
t_{s}(\ad,\bm{v}) = \sqrt{\frac{\pi}{8}} \frac{\bar{\rho}_{d,{\rm int}}\ad}{\rho_{g}c_{s}}\left(1+\frac{9\pi}{128}\frac{|\bm{v}-\bm{u}_{g}|^{2}}{c_{s}^{2}} \right)^{-1/2},\label{eq: stopping time epstein}
\end{equation}
where $\bar{\rho}_{d,{\rm int}}$ is the internal grain density. Epstein drag is generally appropriate for astrophysical conditions
in which MHD and charging effects can be neglected as done here (generally, for cooler, denser gas; see, e.g., \hs).

\subsubsection{Grain mass distribution}
 We will assume a simple power law distribution of grain sizes between $\ad=\amin$ and $\ad=\amax$. In
 order to reduce the number of
free parameters, we use the standard MRN distribution \citep{Mathis1977}, which postulates that the
 mass of grains $d\mu$ within logarithmic range of sizes $d\ln\ad$ is
$d\mu/d\ln\ad\propto \ad^{0.5}$, such that most of the mass is in the largest grains, along with the total dust-to-gas-mass ratio $\mu_{0}\approx0.01$. The original distribution postulates that $\amin\approx5{\rm nm}$ and
$\amax \approx 0.25 {\rm \mu m}$, although subsequent works have suggested that
this underestimates significantly the population of small grains and misses a population of larger grains, even in the diffuse ISM \citep{Weingartner2001,Zubko2004,Draine2009}. Less is known about the grain distribution in more dynamic environments (e.g., around AGB stars or the AGN dusty torus\revchng{, see, e.g., \citealp{Hoefner2018,Murray2005}}), where our simulations can apply by virtue of their dimensionless nature (see \cref{{sub: mapping to astro}} below). Cursory lower-resolution tests of different
grain-mass distributions (e.g., small-grain dominated, $d\mu/d\ln\ad\propto \ad^{-0.5}$) have not revealed significant differences, so we shall not explore this in detail. \revchng{It is, however, worth noting that for the gas of the protoplanetary-disk streaming instability, the grain
distribution can affect important details of the linear instability \citep{Paardekooper2020,McNally2021,Zhu2021}, so this issue may be worth revisiting in more detail in future work.}

\begin{table*}
\centering
 \begin{tabular}{|| m{0.48\textwidth} m{0.48\textwidth} ||}
 \hline
 \emph{Constant drift regime} & \emph{Non-constant drift regime}\\
 \hline
 $ \lambda_{\rm rad} < \ad$, $Q_{\rm abs}\sim 1$ (radiative force on the dust) & $ \lambda_{\rm rad} > \ad $, $Q_{\rm abs}\sim \ad/\lambda_{\rm rad}$  (force on the dust) \emph{or} acceleration of the gas \\
  $\aext \propto \ad^{-1}$ & $\aext \sim {\rm const.}$\\
   $\ws\sim {\rm const.}$, $t_{s}\propto \ad$ & $\ws\propto \ad^{1/2}$, $t_{s}\propto \ad^{1/2}$ (supersonic drift);  $\ws\propto \ad$, $t_{s}\propto \ad$ (subsonic drift)\\
   Linear instability similar to single-grain-size case & Range of resonant angles changes character of linear instability\\
   Strong  correlations between grains of different sizes in saturated state& Grain correlations broadly similar to externally forced turbulence\\
 \hline
 \end{tabular}
 \caption{Summary of some basic   properties of the grains and the RDI in  the constant- and non-constant-drift regimes.}\label{tab: regime}
\end{table*}

\subsection{Grain acceleration regimes}\label{sub: acceleration}

The acoustic RDI  requires a net drift between grains and the gas as their energy source.  As discussed extensively
in \hs\ and \citet{Hopkins2018a},
such a net drift is expected to occur generically in the presence of
 radiation fields, which couple to the grains more strongly than to the gas under most conditions (e.g., \citealt{Weingartner2001a,Murray2005}), sourcing the $\baext$ term in \cref{eq: f dt}.
In this radiatively driven situation, two acceleration regimes naturally emerge, applying  respectively to grains smaller or larger than the
wavelength of the accelerating radiation $\lambda_{\rm rad}$. The different scaling of $\ws(\ad)$
has important implications for the development of instabilities. Some basic properties of the different regimes are summarized in \cref{tab: regime} for quick reference.

A coherent radiation field of energy density $e_{\rm rad}$ causes a grain acceleration $\aext \sim Q_{\rm abs} e_{\rm rad} \ad^{2}/m_{d} c^{2}$, where $m_{d}=4/3\, \pi \bar{\rho}_{d}\ad^{3}$ is the grain's mass. The factor $Q_{\rm abs}$ is the absorption efficiency; $Q_{\rm abs}\sim 1$ if $\ad\gg \lambda_{\rm rad}$, while $Q_{\rm abs}\sim \ad/\lambda_{\rm rad}$ if $\ad\ll \lambda_{\rm rad}$ \citep{Weingartner2001a}. Thus,
large grains in shorter wavelength radiation fields feel an acceleration $\aext \propto \ad^{-1}$, while small
grains in longer wavelength radiation fields feel an acceleration that is independent of $\ad$.  In such a driven situation, the (quasi-)equilibrium  occurs when the acceleration and drag forces on grains balance,\footnote{Note
that this is not actually a true equilibrium because it occurs with the
system as a whole (gas and dust) accelerating linearly at a constant rate, unless $\aext$ can be balanced
by another external force on both gas and dust (e.g., gravity). However, this subtlety does not make
a difference to the arguments here or to our simulations. The issue is discussed in detail in  \hs. } \emph{viz.,} when \begin{equation}
\frac{\ws(\ad,v_{d})}{t_{s}(\ad,v_{d})}=\aext(\ad).\label{eq: ws to solve}
\end{equation}
 Because $t_{s}(\ad,v_{d})\propto \ad$,
this implies that $\ws$ is independent of grain size if  $\ad\gg \lambda_{\rm rad}$, which we term the \emph{constant drift} regime;
the opposite regime, where $\ws(\ad)$ is a function of grain size (for $\ad\ll \lambda_{\rm rad}$) is termed the \emph{non-constant drift} regime. In the latter case, we see from \cref{eq: stopping time epstein,eq: ws to solve} that
for subsonic drift $\ws\ll1$, $t_{s}$ in the saturated state is independent of $\ws$ so that $t_{s}\propto\ad$ and $\ws\propto\ad$, while
for supersonic drift $\ws\gg1$ (the case of more relevance to this article), the decrease in $t_{s}$ with $\ws$ implies that $t_{s}\propto \ad^{1/2}$ and $\ws\propto \ad^{1/2}$.
Of course, in many physical situations, the distribution of grain sizes could fall around $\lambda_{\rm rad}$ (i.e., $\amin<\lambda_{\rm rad}<\amax$), in which case the larger grains will lie in the
constant-drift regime and the smaller grains in the non-constant-drift regime. However, given that
our goal here is to explore the basic physics of the multi-grain RDI, we will not consider such situations in detail in this work.

A net drift of grains can also be set up through a force on the gas (and not the dust), due to radiation pressure absorbed by gas (e.g., line pressure)  or gravity (e.g., in a stratified medium). In this case, the situation is identical to the non-constant-drift regime. \revchng{This physics is applicable, albeit with a
different linear instability, to the  polydisperse streaming
instability in protoplanetary disks \citep{Krapp2019}. }

\revchng{Note that \cite{Hopkins2021} also explore the RDI with a spectrum of grain sizes, explicitly
including vertical stratification and more complex radiation-MHD effects. For some cases, they also consider simulation variants in the $\ad\gg \lambda_{\rm rad}$  and $\ad\ll \lambda_{\rm rad}$
 regimes. Rather than using the labels ``constant-drift'' and ``non-constant-drift,'' they label the $\ad\ll \lambda_{\rm rad}$ (non-constant-drift) regime simulations with ``\textbf{-Q},'' to signify that $Q_{\rm abs}$ depends on $\ad$ in this regime.
 This different nomenclature is used their because  most of their simulations use
 explicit radiation-MHD effects, which leads to more complex relations between $\ws$ and $\ad$ that depend on the dynamics of
 the radiation field.   }

\subsection{Behavior of the acoustic Resonant Drag Instability}\label{sub: RDI regimes}

As discussed in detail in \citet{Squire2018}, \hs, and \msh, the acoustic RDI
exhibits different behaviors based on a dimensionless scale parameter $k \,c_{s} t_{s}$, where
$k$ is the wavenumber of the mode.
In \cref{app: linear}, we cover in detail the linear
behavior of the acoustic RDI with a spectrum of grain sizes, showing, as expected, that
the same parameter (which now covers a range of values because of the range in $\ad$) has a similar influence on the
instability.
Specifically, three regimes emerge: the instability is in the low-$k$ regime if  $k\,c_{s} t_{s}\lesssim \mu$; the mid-$k$ regime if  $\mu\lesssim k\,c_{s} t_{s}\lesssim \mu^{-1}$; and the high-$k$ regime if $\mu^{-1}\lesssim k\,c_{s} t_{s}$.
With a spectrum of grains, there is ambiguity surrounding these delineations ($t_{s}$ and $\mu$ depend on $\ad$),
but they are nonetheless useful for general understanding (see \cref{app: linear} for more precise estimates).
In the low-$k$ regime, the fastest growing modes are generally non-resonant and grow in
the direction aligned with the drift, involving strong perturbations to both the gas and dust densities. We show
in \cref{sub: low k linear} that such modes are generally agnostic to the presence of a spectrum of grain sizes or the drift regime (constant versus non-constant).
In contrast, modes in the mid- and high-$k$ regimes are  fastest growing at a specific ``resonant'' mode direction,
where the projection of the dust drift speed  onto that direction is equal to the sound speed. They have a much stronger dust-density response than gas-density  response. Such
modes behave similarly in the constant-drift regime with a spectrum of grain sizes, but are significantly modified in the non-constant-drift regime
because each grain size resonates with a different mode angle (see \cref{app: grain indep accel}).

In addition to the linear behavior, the three regimes  control the \revchng{acoustic RDI's} nonlinear evolution (\msh).
While gas motions and dust clumping driven by the \revchng{acoustic RDI} in the low-$k$ regime broadly resemble standard
driven supersonic  turbulence (although there are distinct differences),  the mid- and high-$k$ regimes  are very  different, with the
resonant mode structure remaining clear well into the saturated state and  across
all scales. In larger-$\ws$ cases, this manifests itself through thin dust filaments, which ``draft'' on the
nearby dust and never reach a saturated turbulent steady state.
For subsonic drift, the linear and nonlinear behavior is most similar to the
low-$k$ regime (indeed, the subsonic instability at mid- to high-$k$ is non-resonant
and depends on details of the equation of state and drag law; \hs).

Finally, it is worth reiterating from previous works that the acoustic RDI generally has no preferred scale
in any of the three regimes. Rather, modes at smaller scales grow faster, with the growth rate $\Im(\omega)$
scaling as $\Im(\omega)\sim k^{2/3}$, $\Im(\omega)\sim k^{1/2}$, and $\Im(\omega)\sim k^{1/3}$, in the
low-, mid-, and high-$k$ regimes, respectively. Thus, simulations cannot be converged in the conventional sense,
in that a higher-resolution simulation will resolve faster-growing modes (in the absence of a small-scale dissipative effect such as viscosity). However, as shown by \msh\ (appendix B3), the bulk properties of the saturated state are effectively resolution independent  once box-scale  modes saturate nonlinearly.

\subsection{Simulation design}\label{sub: simulation set up}

We use the code GIZMO,\footnote{A public version of the code, including all methods used in this paper, is available at \href{http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html}{http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html} }
which solves the fluid equations using the second-order Lagrangian ``Meshless Finite Volume'' (MFV) method \citep{Hopkins2015}. Dust is included via the super-particle approach \citep{YoudinA2007,Hopkins2016a},
using a random sampling of grain sizes $\ad$ across the full continuous distribution (i.e., we do not use a
set number of grain-size bins). The backreaction force of the dust on the gas is computed using a standard momentum-conserving scheme \citep{YoudinA2007}, with details of the scheme and a variety of numerical tests described in appendix B of \msh\ (although \msh\ considers
only a single grain size, there are no significant  numerical complications that arise
from the use of a spectrum of grains).

The range of available parameter space for simulations using a grain-size spectrum is
extensive, even without magnetization and grain charge: it includes the drift-velocity regimes (constant versus non-constant, supersonic versus
subsonic, and mixtures of each), stopping-time distributions (which could straddle the different $k$ regimes),  the distribution of grain masses, and the total dust-to-gas-mass ratio.
Because  of this, and motivated by the goal of better understanding RDI-turbulence physics rather than
detailed matching of specific astrophysical situations, we choose
in this article to focus on the detailed understanding of just two sets of RDI parameters.
We supplement this by comparing these directly to simulations without
dust backreaction ($\mu_{0}=0$), where turbulence
is driven by large-scale external forcing to have a similar  velocity dispersion.
\revchng{The purpose of this comparison is two fold: firstly, and most importantly,
it enables us to probe the physics of dust clumping in RDI-generated turbulence by direct comparison to the better-understood
case of standard (Kolmogorov) turbulence, revealing clearly their most important differences.
Secondly, in application to AGB-star winds or AGN outflows, the forced-turbulence simulation
could provide a reasonable model for dust clumping   if the RDI were not operating. Specifically, one might expect larger-scale global instabilities (e.g.,
 Rayleigh-Taylor-like instabilities; \citealp{Krumholz2012}) to drive turbulence, which would then, through a turbulent cascade, drive
 fluctuations on the small scales considered here. Thus, the simulations act as a benchmark for how such a dust-driven wind might
 clump dust in the absence of RDIs (e.g., at extremely small dust-to-gas ratios), although we caution that the
 magnitude of the turbulent driving is not at all tuned to explore this in detail (rather it is tuned to address the first point and probe the physics). }

The overall approach complements and builds on that of \msh\ and \citet{Hopkins2020}, which
surveyed a wide range of parameters to understand how the RDI behaves in different regimes. Specifically,
 the results of \msh\ tell us that the most interesting computationally accessible RDI behaviour -- i.e.,
 that which exhibits the most interesting differences compared to standard turbulence -- occurs in the ``mid-$k$''
 range.
 As discussed above (\cref{sub: RDI regimes}) this regime is also expected to show
 more interesting differences  between the grain-spectrum and single-grain RDIs, so is, overall, the
 most obvious candidate for further study.\footnote{Unfortunately, reaching the true high-$k$ regime
has proved to be computationally challenging, because the width of the resonant wavelengths become increasingly narrow at increasing $k$, necessitating excessively high resolution. \msh\ considered some cases
around the boundary between the mid- and high-$k$ regimes, which were  similar to the mid-$k$ cases (see their appendix A for more information).}
 Thus we will ignore the low-$k$ regime and/or subsonic drift in this study. Although they are
 potentially astrophysically relevant in many situations \citep{Hopkins2021}, these regimes
 can likely be mostly adequately understood using a combination
 of RDI-related understanding from \msh\  and \hs\, and theories of collisions/clustering in turbulence
with a spectrum of sizes without dust backreaction \citep[e.g.,][]{Pan2014a,Pan2014b,Hopkins2016a,Mattsson2019,Li2020}.

Based on the discussion of the previous paragraph, we focus on four simulations with a grain-size spectrum covering a factor of 100 ($\amax=100\amin$), each in a cubic box
of size $L^{3}$.  These are:
\begin{description}
\item[\textbf{Constant-Drift RDI}] This simulation sets $\aext\propto1/\ad$ such that $\ws$ is independent of $\ad$ and  $t_{s}\propto \ad$. The acceleration and grain size range (see below) is chosen such that $\ws\approx1.5$, which
is convenient because the resonant angle at which the modes grow fastest ($\cos\theta_{k}=1/\ws$) is neither
too oblique nor too parallel, as well as being astrophysically reasonable. The outer-scale normalized wavenumbers of grains range from $(2\pi/L)c_{s}t_{s}\approx0.005$ to $(2\pi/L)c_{s}t_{s}\approx0.5$.
\item[\textbf{No-Backreaction Constant Drift}] This simulation has parameters (including $\aext$) that match exactly Constant Drift RDI, but with
no dust backreaction, and thus no RDI. Instead, the turbulence is externally driven, with the amplitude of the forcing chosen so that the resulting turbulence has a velocity dispersion that matches (as closely as possible) the Constant Drift RDI simulation.
\item[\textbf{Non-Constant-Drift RDI}] This simulation sets $\aext$ constant such that $\ws$ and $t_{s}$ both depend
on $\ad$ ($Q_{\rm abs}\sim\ad/\lambda_{\rm rad}$). The acceleration and grain size range  is chosen such that $\ws\sim \ad^{1/2}$ ranges between $\ws\approx 0.73$ (for the smallest grains) and $\ws\approx 12.7$ (for the largest grains). The outer-scale normalized wavenumbers of grains range from $(2\pi/L)c_{s}t_{s}\approx0.003$ to $(2\pi/L)c_{s}t_{s}\approx0.05$ (note that $t_{s}\sim \ad^{1/2}$ for the larger/faster grains).
\item[\textbf{No-Backreaction Non-Constant Drift}] This simulation has parameters (including $\aext$)  that match exactly Non-Constant Drift RDI, but with no dust backreaction and using external driving to match the turbulence amplitude of Non-Constant Drift RDI.
\end{description}
Note that in each case, $(2\pi/L)c_{s}t_{s}<1$ is chosen so that smaller-scale modes in the box do not move too
far into the high-$k$ regime, where they may be adversely affected by numerical resolution. In the non-constant drift cases, the range in $\ws$ is chosen so
that most grains are supersonic, but the largest grains are not highly supersonic ($\ws$ not too much larger than $\simeq 10$), which was shown
by \msh\ to cause particularly extreme behavior (drafting) that does not necessarily converge in time (owing in part due to the finite periodic box assumption that we adopt). The equilibrium $\ws$ and $(2\pi/L)c_{s}t_{s}$ are shown
below in the insets of \cref{fig:time.evolution.sda,{fig:time.evolution.ca}}.

\revchng{\subsubsection{Numerical parameters}}

\revchng{In each case, it is convenient to use code units defined by $\langle \rho_{g}\rangle=1$, $c_{s}=1$, $\bar{\rho}_{d,{\rm int}}=1$, and $L=1$. Then,
$\ws(\ad)$ and $t_{s}(\ad)$ are numerically determined by the combination of $\ad$ and $\aext$ (in code units); in order to obtain
the simulation parameters discussed above, we need
$\amax=0.15=100\amin$ with $\aext = 3/\ad$ for the constant-drift simulations, and  $\amax=0.0878=100\amin$ with $\aext =1400$ for the non-constant-draft simulations.}
Grains are initialized with their equilibrium drift velocity ($\ws(\ad)/t_{s}(\ad)=\aext$), which is also
present in the driven cases without backreaction (unlike previous studies of turbulent grain dynamics). \revchng{We use periodic boundary conditions with a resolution of
$256^{3}$ gas particles and $4\times256^{3}$ dust particles in all simulations}. This was chosen based on the scaling tests in \msh\ (appendix B), which showed  reasonable convergence in the mid-$k$ regime between $128^{3}$ and $256^{3}$
 (although RDI simulations of this type can never be truly converged because the instability generally grows fastest at the smallest scales of the box).\footnote{\revchng{The numerical time step in all simulations is strongly limited by the
integration of the smallest grains, which have a very short stopping time. This increases their computational cost by a factor of ${\simeq}10$ compared
to an equivalent system without dust, which, combined
with the relatively slow saturation of the RDI (up to $t\simeq 30 L/c_{s}$),  makes the simulations relatively computationally expensive despite their modest resolution.}} \msh\ also showed  convergence in RDI dust probability-density functions (PDFs) once the dust resolution approached that of the gas (a factor $4$ lower gives good results, but not a factor 16 lower) in line with previous works \citep[e.g.,][]{Bai2010b}; our choice of $4\times256^{3}$ dust particles is made due to the wide range of grain sizes in our simulations and is explored further in  \cref{app: linear gizmo comparison}.
 Although we use an MRN  dust-mass distribution in all cases $d\mu/d\ln\ad\propto \ad^{1/2}$, we
use an equal number of numerical super-particles randomly sampled across each logarithmic size range (i.e., the super-particles change in their ``total'' super-particle mass proportionally to $\ad^{1/2}$)
in order to not under-resolve the small particles.
\revchng{We note that, by spreading the grain-size distribution across the full range of particles, this method appears to avoid some of the
convergence problems that are well known in the numerical computation of polydisperse dust instabilities \citep{Paardekooper2021}.
We study this issue further in \cref{app: linear gizmo comparison} by comparing the early phases of GIZMO simulations to linear results.  }

In the no-backreaction, driven cases, we use time-correlated incompressible (solenoidal) forcing at the  largest scales (Fourier modes with $k\leq4\pi/L$; \citealp{Hopkins2016a}). The  correlation time is $0.5L/c_{s}$ in the constant-drift simulation and $0.2L/c_{s}$ in the non-constant-drift simulation, because the turbulence in the non-constant-drift case
has higher Mach number, so a shorter box-crossing time.
We use the total dust-to-gas-mass ratio $\mu_{0}=0.01$ for both of the RDI simulations ($\mu_{0}=0$ for the no-backreaction simulations).

\subsubsection{Grain-size bins for diagnostics}\label{sub: dust binning}

Although our simulations involve a continuous distribution of grain sizes, for most diagnostics -- for instance,
any quantity involving dust density or velocity as per \cref{eq: vd rhod defs} -- it is necessary
to bin the dust size distribution. We choose to do this across $8$ logarithmically spaced bins, which we label
with $\ab{1}$ for the smallest grains, through to $\ab{8}$ for the largest grains. More precisely, the label $\ab{i}$
refers to grains with $\ad$ between $10^{b_{i}}\amax$ and $10^{b_{i+1}}\amax$, where $b_{i}=(-2,-1.75,-1.5,-1.25,-1,-0.75,-0.5,-0.25,0)$. So, for example, grains of size $\ab{1}$  have a density and bulk velocity  given
by \cref{eq: vd rhod defs}, with $\abc{l}=10^{-2}\amax=\amin$ and $\abc{h}=10^{-1.75}\amax\approx 1.78 \amin$.

\subsection{Mapping to astrophysical applications}\label{sub: mapping to astro}

Our simulations are not intended to map to a specific astrophysical object, but
rather study the generic behavior of RDI-generated turbulence. Here, we
outline how the simulation units -- with $\langle \rho_{g}\rangle=1$, $c_{s}=1$, $\bar{\rho}_{d,{\rm int}}=1$, and $L=1$ (the box size) -- translate into various astrophysical situations and processes of interest. The
two important properties are  the grains' drift velocity -- which depends on the radiation field, grain size, and other
gas properties -- and the physical box scale $L_{\rm phys}$, which depends on the gas density and grain sizes.
Following \hs, we  briefly consider as examples asymptotic giant-branch (AGB) stars, Active Galactic Nucleii (AGN), and star-forming regions (giant molecular clouds; GMCs),
which are the contexts most relevant to uncharged dust; more
detailed estimates for a wider range of situations are given in \citet{Hopkins2018a} (see their figure 6).

A simulation can be mapped to a given physical situation by comparing the box size $L=1$ to the grain-drift length at $\ws=0$, $L_{\rm drift}=c_{s} t_{s}|_{\ws=0} = \sqrt{\pi/8}\ad \bar{\rho}_{d,{\rm int}}/\rho_{g}$. In order
to obtain the desired $\ws$ and $(2\pi/L)c_{s}t_{s}$, as discussed above, the
constant-drift simulations use $\amax=0.15=100\amin$, while the non-constant-drift simulations use $\amax=0.0878=100\amin$. Thus, the box scale translates into $L_{\rm phys}\approx6.7\, \amax \bar{\rho}_{d,{\rm int}}/\rho_{g}$ and $L_{\rm phys}=11.4\, \amax  \bar{\rho}_{d,{\rm int}}/\rho_{g}$ for the constant- and non-constant-drift
cases, respectively.
\begin{description}
\item[\textbf{Star-forming regions}]
\hs\ estimates GMC conditions  under reasonable assumptions (e.g., a cloud of $\sim\!10{\rm pc}$ that has converted $\sim\!0.1$ of its mass into stars), finding
$\ws\sim10 $ for larger grains with $Q_{\rm abs}\sim 1$ ($\ws$ is larger closer to a more luminous source and at lower gas density). Depending on the wavelength of the radiation and the dominant grain sizes, this is  in reasonable agreement to both the constant- and non-constant-drift simulations. The above estimates yield $L_{\rm phys}\simeq 20 {\rm pc}$, showing that the smaller-scale, shorter-time dynamics in our simulations will
apply to larger scales around a GMC.
\item[\textbf{AGB-star winds}]
The envelopes/winds of AGB stars are dust laden and driven by radiation pressure. Using similar estimates to
\hs\ (a stellar luminosity of $\sim\!10^{5}M_{\sun}$ with  mass-loss rate $\sim\!10^{-4}M_{\sun}{\rm yr}^{-1}$) yields
$\ws \sim 0.1\rightarrow 10$, which is in  the range probed by either simulation. Similarly, estimating the gas density at $\sim\!1{\rm AU}$ from a star with a wind velocity $\sim\!10{\rm kms}^{-1}$ yields $L_{\rm phys}\sim 5\times10^{-5}{\rm AU}$ showing that our boxes are probing relatively small scales inside the outflowing wind, where the RDI would seed small-scale clumping.
\item[\textbf{AGN-driven outflows}]
\hs\ and \citet{Hopkins2018a} estimate very fast drift velocities, from $\ws\sim 1$ for small grains to $\ws\gtrsim 100 $ for larger grains, in the ``dusty torus'' region around an AGN with luminosity $\sim\!10^{46}{\rm erg}{\rm s}^{-1}$. Using $n\sim 10^{6}{\rm cm}^{-3}$, as appropriate to a   denser closer-in  region, yields $L_{\rm phys}\sim 2\times 10^{-4}{\rm pc}$, showing that (as for the AGB wind), the numerical box represents a small patch of the larger system, representing fast dust dynamics on small scales.
\end{description}
Overall, we see that our  simulation parameters are  likely most  applicable to  AGB winds; GMC-related  applications would be better served by a somewhat smaller-scale box, while the drift velocities in the AGN are somewhat more extreme than simulated. \revchng{Also, grain charging and MHD effects are likely
more important for many GMC-like conditions, which would change our results substantially \citep{Hopkins2021}.} However, parameters in all three cases vary enormously, and the physical ideas we explore are generally relevant to different
regions of all three cases. Our chosen dust-to-gas-mass ratio of $\mu_{0}=0.01$
broadly applies to most situations, although this can vary significantly both above and below this value (see, e.g., \citealt{Knapp1985,Dharmawardena2018,Wallstroem2019} for some examples of AGB winds).

\begin{figure*}
\begin{center}
\includegraphics[width=1\textwidth]{{DustGasPanels_s}.pdf}
\caption{Dust and gas visualizations for the constant-drift simulations. Top images show the gas density on the surface of the cube with medium-sized ($\ab{3}$) grains shown as black dots, while bottom images show volume renderings of the dust density for small ($\ab{1}$) and large ($\ab{8}$) grains (the dust binning is described in   \cref{sub: dust binning}).  Panel (a) shows the near-linear phase of the RDI at early times and  panel (b) shows show its saturated turbulent state at late times (see also \cref{fig:time.evolution.sda}). Panel (c)  illustrates the saturated turbulence in the no-backreaction constant-drift simulation, where the turbulence is externally forced (only large grains are shown in the bottom image).
}
\label{fig:panels.sda}
\end{center}
\end{figure*}
\begin{figure*}
\begin{center}
\includegraphics[height=\tfigh]{{sda-dust-evolution}.pdf}~\includegraphics[height=\tfigh]{{sda-gas-evolution}.pdf}\includegraphics[height=\tfigh]{{sturb-gas-evolution}.pdf}
\caption{Time evolution of bulk quantities in the constant-drift simulations. The left-hand panel shows the evolution
of dust-related quantities in the RDI simulation for different-sized grains from small ($\ab{1}$; blue) to large ($\ab{8}$; red), showing the measured perpendicular dust velocity dispersion $\langle \delta v_{d,\perp}\rangle= \langle v_{d,\perp}^{2}\rangle$ with solid lines, and the drift velocity $\langle \ws\rangle=\langle v_{d,\|}\rangle$ with dashed lines.  The inset shows the
equilibrium $\ws(\ad)$ (black, left axis) and $t_{s}(\ad)$ (red, right axis). Similar quantities for the gas and dust averaged over all bins
are shown in the middle panel for the RDI simulation, and in the right-hand panel the no-backreaction case with externally driven turbulence. Their comparison some clear differences in the structure of RDI-driven turbulence (larger parallel velocity fluctuations and stronger relative density fluctuations) that cannot be matched by isotropic incompressible driving. In taking averages, gas quantities are weighted by the gas mass and dust  quantities by the dust mass.
The black dotted lines in the middle panel show the linear growth rates of resonant modes with $k\approx 32\pi/L$ (approximately $1/16$ of the box scale) and $k\approx 2\pi/L$ (approximately the box scale); see   \cref{eq: size dep final ans} and \cref{fig:linear adep}.
}
\label{fig:time.evolution.sda}
\end{center}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[width=1\textwidth]{{DustGasPanels_c}.pdf}
\caption{As for \cref{fig:panels.sda}, but for the non-constant-drift simulations.
}
\label{fig:panels.ca}
\end{center}
\end{figure*}

\begin{figure*}
\begin{center}
\includegraphics[height=\tfigh]{{ca-dust-evolution}.pdf}~\includegraphics[height=\tfigh]{{ca-gas-evolution}.pdf}\includegraphics[height=\tfigh]{{cturb-gas-evolution}.pdf}
\caption{As for \cref{fig:time.evolution.sda}, but for the non-constant drift simulation. Because the larger grains, with most of the mass,
move faster  in this case, the resulting turbulence is trans-sonic, with $\langle u_{g}^{2}\rangle\sim c_{s}^{2}$
and larger density fluctuations than the constant-drift simulation. The black dotted line in the middle panel indicates the linear growth
rate of \revchng{the box-scale mode seen in the left panels of \cref{fig:panels.ca}; the details and convergence of this linear growth phase are assessed
in detail in \cref{app: linear gizmo comparison}}.
}
\label{fig:time.evolution.ca}
\end{center}
\end{figure*}

\section{Results}\label{sec: results}

In this section we present results of the four GIZMO simulations outlined above (\cref{sub: simulation set up}). We
start with a general exploration of the time evolution and turbulence  structure (\cref{sub: results: time evolution,sub: results: statistics}) then consider more detailed statistics related to dust clumping and collisions in \cref{sub: results: collisions}.

\subsection{General morphology and time evolution}\label{sub: results: time evolution}

In \cref{fig:panels.sda,fig:time.evolution.sda} we show  three-dimensional visualizations of the
turbulence structure and the time evolution of important quantities for the constant-drift simulations. The
same is shown for the non-constant-drift simulations in \cref{fig:panels.ca,fig:time.evolution.ca}.
In each case, the turbulence structure is illustrated  during  the instability-growth phase (left-hand panels)  and once it saturates nonlinearly (middle panels)  for the RDI simulations,
 and compared to the gas and dust structure of the externally forced runs without dust backreaction (right-hand panels).
 The bottom panels visualize the smallest and largest grains in the RDI, to compare differences in their structure.
 The time-evolution panels (\cref{fig:time.evolution.sda,fig:time.evolution.ca}) show how
 velocity dispersions and drift velocities vary with grain size (left-hand panels; insets show the equilibrium dust parameters) and gas and dust
 velocity dispersions integrated over all grain sizes (middle and right-hand panels), comparing
 the RDI cases to the driven-turbulence ones (right-hand panels).
 These plots illustrate the basic time evolution of the RDI  through the linear phase and  saturation, and allow
 simple comparison to the saturated state of driven turbulence.

The early-time growth in the constant- and non-constant-drift cases is rather different and well explained
by the linear mode structure in each case. In the constant-drift case, the growth
is initially rapid but slows in time, broadly matching the predicted linear growth rates  (dotted lines in the middle panel of \cref{fig:time.evolution.sda}), which increase monotonically with scale (see  \cref{fig:linear adep}).
This behavior is expected and discussed in previous works on RDI evolution (\citealp{Seligman2019}; \msh); it results from
smaller-scale modes growing and saturating nonlinearly more rapidly than large-scale modes, so that the growth of a bulk quantity
(such as $\langle u_{g}^{2}\rangle$) is first dominated by the faster smaller scales, then, at later times (once the small-scales saturate), by the slower
larger scales. Such evolution is also clearly seen in the morphology in \cref{fig:panels.sda}: there is strong clumping of small grains at
small scales by $t\approx5 L/c_{s}$ even though the turbulence is far from full saturation at this point. The
structures in the gas and dust, which clearly show the RDI resonant angle ($\theta_{k}=\cos^{-1}(\ws^{-1})\approx 48^{\circ}$)
are broadly similar to those that develop in saturation at the box scale (middle panel). We see that the smallest grains
exhibit the strongest clumping (note higher-density patches in the bottom panels)  and have a modestly higher velocity dispersion in the
saturated state. However, it is also clear that large and small grains are undergoing similar dynamics:
 $w_{s}$ remains remarkably similar for all grains even as it fluctuates in time (left-hand panel of \cref{fig:time.evolution.sda}), and
high-density  regions of large and small grains are clearly correlated spatially (bottom middle panels of \cref{fig:panels.sda}). Thus, as suggested by linear calculations (\cref{app: grain dep accel}),
 the constant-drift RDI  involves different grain sizes interacting with the gas in similar ways,
 driving resonant modes that cause strong clumping for all sizes concurrently.

The non-constant-drift RDI is more complex, with significant differences between the dynamics
of small and large grains. As discussed in detail in \cref{app: grain indep accel},
a number of  different linear-instability mechanisms can operate and/or dominate in the non-constant-drift RDI, and this
is also true for the chosen scale and parameters of the simulation. In particular, there exists a large-scale parallel
mode that resembles a backward-propagating sound wave, which has a similar growth rate to a smaller-scale, more oblique resonant mode
that predominantly effects the largest grains (see \cref{fig:linear aindep}).
\revchng{In \cref{app: linear gizmo comparison}, we test the detailed linear growth of these modes across different scales,
showing generally good agreement with linear predictions (see \cref{fig: nl lin comparison large,fig: nl lin comparison k theta}).}
The \revchng{effect} of this linear mode  structure can be seen in the turbulence morphology that develops, as well as
in the time evolution of the RDI. We see the formation of a box-scale strong shock at early times (left-hand panel of \cref{fig:panels.ca}), which
creates a strong  density contrast in the gas and the smaller dust grains (because they are better
coupled to the gas).  It is also clear in the gas and dust  time evolution (middle panel of \cref{fig:time.evolution.ca}),
manifesting  as a large density dispersion that develops until $t\approx 5L/c_{s}$, at which point in breaks up and grows  larger
turbulent velocity fluctuations into the true steady state for $t\gtrsim 8L/c_{s}$.
The visual turbulent morphology in the saturated state at late times (middle panel of \cref{fig:panels.ca}) is quite different from the driven
turbulence  (right-hand panel), with
smaller-scale structures in the gas caused by highly elongated dust velocity filaments. We interpret this behavior as
being due to the quasi-resonant mode, which affects only the largest grains and is highly oblique because of their fast drift velocity  (\cref{fig:linear aindep} inset \revchng{and \cref{fig: nl lin comparison k theta}}); indeed, there is strong, oblique clumping of large grains but not small grains (middle-lower panel).
Further evidence for the different
dynamics of small and large grains is seen in the evolution of $\ws$ (left-hand panel of \cref{fig:time.evolution.ca}),
which exhibits much larger fluctuations for small grains than large grains in the saturated turbulence, even
though their perpendicular dust velocity dispersions (solid lines) are similar.

\begin{figure}
\begin{center}
\includegraphics[width=0.9\columnwidth]{{spectra}.pdf}
\caption{Top panel: gas velocity spectra from all four simulations, with blue lines showing the constant-drift simulations, red lines showing the non-constant-drift simulations, solid lines showing RDI simulations, and dashed lines showing no-backreaction (externally forced) cases. We see similar behavior in all cases, with kinetic energy spectra $\sim\!k^{-5/3}$ that broadly match
previous results for subsonic and trans-sonic turbulence. Bottom panel:  gas density spectra, which also broadly match expectations, with a significantly flatter spectrum in the higher-Mach-number non-constant-drift case.
}
\label{fig: spectra}
\end{center}
\end{figure}

\subsubsection{The level of turbulent driving}\label{subsub: turb driving level}

It is worth briefly commenting on some differences between RDI turbulence and the externally forced runs without dust backreaction.
The level of forcing in the no-backreaction runs was chosen to match the RDI saturation as best as possible; however,
 the choice of isotropic forcing in the driven runs means that
$u_{g,\perp}^{\rm rms}\approx 2 u_{g,\|}^{\rm rms}$ (where $u_{g,\cdot}^{\rm rms}= \langle u_{g,\cdot}^{2}\rangle^{1/2}$ is the root-mean-squared gas velocity), while the saturated state of the RDI runs are distinctly non-isotropic, with much larger parallel velocity dispersions
$u_{g,\|}^{\rm rms}\gg u_{g,\perp}^{\rm rms}/2$. Thus, although both RDI cases have a modestly larger total velocity dispersion than
their equivalent no-backreaction runs ($u_{g}^{\rm rms}\approx 0.9c_{s}$ versus $u_{g}^{\rm rms}\approx 0.6c_{s}$ for
the constant-drift case; $u_{g}^{\rm rms}\approx 3.6c_{s}$ versus $u_{g}^{\rm rms}\approx 1.7c_{s}$ for
the non-constant-drift case), their perpendicular velocity dispersions are smaller ($u_{g,\perp}^{\rm rms}\approx 0.3c_{s}$ versus $u_{g,\perp}^{\rm rms}\approx 0.5c_{s}$ for
the constant-drift case; $u_{g,\perp}^{\rm rms}\approx 0.8c_{s}$ versus $u_{g,\perp}^{\rm rms}\approx 1.4c_{s}$ for
the non-constant-drift case). In addition, we see that externally forced turbulence has a somewhat lower density dispersion
than RDI turbulence in the constant-drift case, while the opposite is true for the non-constant-drift simulations.
While it might be possible to rectify some of these discrepancies using non-isotropic driving, such an
exercise would not necessarily be helpful: the differences arise from fundamental differences in the structure and evolution
of RDI turbulence and driven turbulence, which is exactly the physics that we wish to study. However, it is
important to keep them in mind as we explore some of the differences in more detail.

\subsubsection{Turbulent Stokes numbers of grains}\label{subsub: stokes numbers}
A parameter of particular importance for turbulent grain dynamics is the ratio of stopping time
to the  turnover time of the turbulence, known as the Stokes number, ${\rm St}$. For comparison to previous
results  on dust dynamics in turbulence without backreaction  \citep[e.g.,][]{Ormel2007,Pan2013,Pumir2016,Mattsson2019} it is helpful to estimate ${\rm St}$ at the largest and smallest resolvable scales in our
simulations. The large-scale turnover time is estimated as $\tau_{L}\sim L/u_g^{\rm rms}$, which (using
the externally forced turbulence values from above to avoid dealing with anisotropic turbulence) is $\tau_{L}\simeq 1.7 L/c_{s}$ for the
constant-drift case, and $\tau_{L}\simeq 0.6 L/c_{s}$  for the non-constant-drift case. As expected, the saturation of the turbulent runs, seen in the right-rand panels of \cref{fig:time.evolution.sda,fig:time.evolution.ca}, occurs over a timescale $\sim\!\tau_{L}$.
The fastest time scale in the turbulence, which occurs at small scales and is termed $\tau_{\nu}$,
 can be estimated from basic Kolmogorov arguments. For a $\sim\!k^{-5/3}$ velocity spectrum, which is at least roughly
 valid  for the trans-sonic turbulence here (see \cref{fig: spectra}), the turnover time of structures of lengthscale $l$ scales as $l^{2/3}$. Using twice the average point spacing at our fiducial resolution, $l_{\nu}\simeq L/128$,
 as an estimate for the smallest scales available before numerical viscosity
 starts damping fluid motions, we estimate $\tau_{\nu}\simeq 0.07L/c_{s}$ and $\tau_{\nu}\simeq 0.02L/c_{s}$ for the constant- and non-constant-drift cases, respectively.
Denoting the outer scale and smallest-scale Stokes numbers as ${\rm St}_{L}=t_{s}/\tau_{L}$ and  ${\rm St}_{\nu}=t_{s}/\tau_{\nu}$, respectively, we thus find
\begin{gather}
{\rm St}_{L}\simeq 0.01\left(\frac{\ad}{\amax}\right)^{1/2},\: {\rm St}_{\nu}\simeq 0.4\left(\frac{\ad}{\amax}\right)^{1/2},\quad \text{Non-constant drift},\nonumber\\
{\rm St}_{L}\simeq 0.05\frac{\ad}{\amax},\: {\rm St}_{\nu}\simeq 1.1\frac{\ad}{\amax},\quad \text{Constant drift}.
\label{eq: stokes in simulations}
\end{gather}
We see that these simple estimates would naively suggest that nearly all grains in our simulations are  ``well coupled'' (${\rm St}\lesssim1$) to the gas turbulence across all scales, meaning that they should passively trace gas motions. We will show below that small grains in the constant-drift RDI are actually not ``well coupled'' to the gas turbulence in this sense, despite the fact that
${\rm St}_{\nu}\simeq 0.01$ for this population. This is likely because the compressive and dust-drift motions associated with the instability are both faster (thus increasing ${\rm St}$)
and more effective at driving dust clumping.
Finally, it is worth noting that the RDI quasi-linear saturation estimate of  \msh, $u_{g}^{\rm rms}\sim \mu^{1/2} (c_{s}t_{s}/L)^{-1/2}$, which compared well
against single-grain-size simulations, yields similar estimates for ${\rm St}$; ${\rm St}_{\nu}\sim \mu^{1/2}(c_{s}t_{s}/L)^{1/2}(l_{\nu}/L)^{-2/3}$, for the Stokes number on scale $l_{\nu}$, which is generally $\ll1$ for outer scales in the low- or mid-$k$ regime with $\mu\ll1$, except at very small scales.
Thus, it seems that simple (Kolmogorov) estimates
of turbulent Stokes numbers in RDI turbulence generically yield  rather small values (${\rm St}\ll1$) compared to what
is usually  necessary to strongly clump grains in turbulent flows (i.e., ${\rm St}\gtrsim 1$; see, e.g., \citet{Hopkins2016a}, who see only very weak clumping
for ${\rm St}\ll 1$ in compressible turbulence with the same numerical methods).

\revchng{The discussion above also suggests that another clumping regime could emerge at extremely small scales: at scales where ${\rm St}_{\nu}\gg 1$,
the turbulent cascade may start to dominate the clumping for grains that  are not strongly affected directly by the RDI.
Given the efficiency of the constant-drift RDI at clumping all grain sizes, such an effect would likely be important only for smaller grains in
the non-constant-drift regime, and would manifest as an enhanced clumping at larger resolution. Using
the estimate \eqref{eq: stokes in simulations}, we see that accessing this regime --  ${\rm St}_{\nu}\gtrsim 1$ for the smallest grains -- requires
${\rm St}_{\nu}$ to be ${\simeq100}$ times larger than the simulations presented here, thus requiring ${\simeq}1000$ times more resolution
elements in each direction (assuming $\tau_{\nu}\propto l^{2/3}$). This is clearly unattainable with present resources, although it is possible
that a similar regime could be accessed with a different set up.}

\subsection{The turbulent steady state}\label{sub: results: statistics}

In this and the following section, we explore more detailed statistics of the saturated states of RDI turbulence, comparing
to the turbulence runs as a reference.

\paragraph*{Spectra}
The simplest measure of the scale-dependent turbulent structure is the spectrum, which
is illustrated for all four cases in \cref{fig: spectra}. The kinetic energy spectra (panel a) of RDI turbulence and externally forced (no-backreaction) turbulence are
seen to be quite similar.  Spectral slopes in the inertial range ($20/L\lesssim k\lesssim 200/L$)
are approximately $\sim\!k^{-5/3}$, as expected for a standard Kolmogorov cascade, and are slightly steeper in the
non-constant-drift runs, as expected because of their higher Mach numbers (velocity spectra steepen to $\sim\!k^{-2}$ for highly supersonic flows; e.g., \citealt{Federrath2013}).
The consistency of the RDI and forced velocity spectra is interesting, given the clear differences in their structures observed in \cref{fig:panels.sda,fig:panels.ca}. For example, we see that although the non-constant-drift RDI turbulence  looks quite
different to its  forced counterpart, with smaller-scale features, the scaling of the spectra at smaller wavenumbers $k\gtrsim 20/L$
 is very similar. The density spectra (panel b) are also consistent with previous works \citep{Konstandin2016},
 although we see more interesting differences between the RDI and externally forced turbulence runs.
The constant-drift runs exhibit similar, relatively steep ($\sim \!k^{-2}$) spectral scaling, but
with stronger compressive motions in the RDI (see also \cref{fig:time.evolution.sda}); however, the non-constant-drift density spectra
are quite different, perhaps because of the higher-mach-number  parallel flows in the RDI case.

\begin{figure}
\begin{center}
\includegraphics[width=0.9\columnwidth]{{sda-2D-hist}.pdf}
\caption{Two-dimensional volume-weighted  PDF of the local (resolution-scale) dust density $\rho_{d}(\ad)$ and gas density $\rho_{g}$ for the constant-drift simulations. Solid lines show
contours at $0.1$ of the maximum, while dashed lines show the contour of $0.01$ of the maximum, with
each color showing the PDF for different grain-size bins from the smallest grains $\ab{1}$ in blue, to the the
largest $\ab{8}$ in dark red. The PDFs are averaged over 3 snapshots during the saturated phase.  The two cases are markedly different: in  externally forced turbulence case without backreaction (bottom panel)
the smallest grains mostly trace the gas and larger grains exhibit more clumping; in the RDI turbulence, the smallest
grains instead exhibit the strongest clumping, with a long tail extending to very small dust densities.
}
\label{fig:2d hist sda}
\end{center}
\end{figure}

\paragraph*{Dust and gas distribution}
A useful measure of the level and structure of
dust clumping is the joint Probability Density Functions (PDFs) of dust and gas density, the contours of which
 are shown in \cref{fig:2d hist sda,fig:2d hist ca} for the
constant- and non-constant-drift cases, respectively.
These illustrate how regions of high gas density correlate with those of high dust density, and likewise for low-density regions. The black dotted lines illustrate the one-to-one correlation
that would be observed if dust were perfectly coupled to the gas.

Let us first consider the constant-drift case (\cref{fig:2d hist sda}), which shows a significant difference between the RDI-generated turbulence (top panel)
and the  forced turbulence without dust backreaction (bottom panel). This is surprising, given
the similarity of their spectra. Most clearly, we see that in RDI turbulence the smallest grains (blue contours)
exhibit larger fluctuations than the larger grains (red contours), particularly at low-densities, in stark contrast to the no-backreaction runs. The characteristic shape -- with a high-probability of low dust density in lower-gas-density regions -- was
also seen in \msh\ and seems to be a typical  feature of the saturated state of the mid-$k$ (or high-$k$)  RDI.
The much wider dust density distribution of larger grains in externally forced turbulence is well explained by their relative turbulent Stokes
numbers (which are quite small across all scales for the smallest grains; see \cref{subsub: stokes numbers}, \citealp{Pan2013,Hopkins2016a,Mattsson2019}), so the
fact that this is not true for the RDI turbulence is an important illustration of its different grain-clumping mechanisms.

In contrast, the non-constant-drift RDI PDFs appear more similar to externally forced turbulence (\cref{fig:2d hist ca}), although
the dust-density distribution is wider for the larger grains in the RDI case (i.e., the difference between the large and small grains is larger). This is likely because  the highly oblique resonant instability,
which is driven  only by the large grains (see middle panel of \cref{fig:panels.ca}), is particularly efficient at grain
clumping; the backward-propagating sound-wave mode, in contrast, clumps grains in a similar
way to standard  turbulence, creating a similar (small) density dispersion in the smaller grains.
This explanation is commensurate with the  enhanced clumping of all grains (compared to driven turbulence)
in the constant-drift case, where the resonant instability operates with all grain sizes.

\begin{figure}
\begin{center}
\includegraphics[width=0.9\columnwidth]{{ca-2D-hist}.pdf}
\caption{Same as \cref{fig:2d hist sda}, but for the non-constant-drift simulation. The RDI PDF is  more
similar to externally forced turbulence in this case, although we still see significantly wider dust-density distributions (despite the
somewhat narrower gas-density distribution), particularly for the largest grains that are more strongly affected
by the mid-$k$ RDI (see text).
}
\label{fig:2d hist ca}
\end{center}
\end{figure}

\subsection{Dust clumping and collisions}\label{sub: results: collisions}

The observations above, along with previous results (\msh; \citealp{Seligman2019,Hopkins2020}), show that the grain-clumping
and its dependence on grain size can be very different in the RDI compared to externally forced turbulence without backreaction.
The most obvious question that arises is whether
these differences significantly change grain-collision statistics in the RDI compared to previous theories \citep{Voelk1980,Zaichik2006,Pan2013,Mattsson2019,Li2020}.
Two key properties influence grain collisions and are needed to compute the collision kernel: the first is the relative
clumping of grains in space \citep{Maxey1987},  the second is the relative  velocity of grains that do collide \citep[e.g.,][]{Pumir2016}. Put together,
these can be used to construct estimates for the collision rate and the sticking-bouncing-fragmentation probabilities \citep[e.g.][]{Garaud2013,Pan2014a}, an understanding of which is key for estimating how turbulence influences grain growth
 in astrophysical scenarios (e.g., in an AGB outflow).
A full, careful estimate of these probabilities requires grappling
 with a number of subtle convergence issues; for instance, relative grain velocities depend on particle separation,
 and estimating the true collision velocity (at near-zero particle separation) requires a careful consideration   of
 how different physical effects contribute to the relative velocity \citep{Falkovich2002}
 and how these are affected by numerics \citep{Pan2014b,Haugen2021}. Our goal here
is less detailed -- to compare and contrast the relative clumping and grain-collision velocities
between the constant- and non-constant-drift RDIs and externally forced turbulence without dust backreaction.
\revchng{We find much stronger clumping and reduced collision velocities
of grains in the constant-drift RDI, with qualitatively different trends for small grains and grains of different sizes.
These trends are sufficiently strong to reveal clear differences between grain-grain collisions in the RDI and externally forced turbulence. We thus
forgo a
a careful study of the convergence to the zero-particle-separation limit, which would be a formidable task for RDI turbulence
given the wide range of different regimes \citep{Hopkins2018}.
Overall, our results suggest}
that the RDI could strongly enhance grain growth rates  in outflows, especially in the constant-drift regime ($\ad\gg \lambda_{\rm rad}$).

\begin{figure}
\begin{center}
\includegraphics[width=0.88\columnwidth]{{gfuncs_sda}.pdf}
\caption{Grain-grain density correlation functions $g(r; \abc{l}, \abc{h})$  for the constant-drift
simulations as a function of the ``primary'' (larger) grain size $\abc{h}$, using the $8$ grain-size bins of \cref{sub: dust binning}. This measures the probability of finding another grain
of size $\abc{l}$ a distance $r$ from the primary grain of size $\abc{h}$; i.e., it quantifies the mean local density enhancement between grain pairs. Different line colors show different grain-size ratios $\abc{h}/\abc{l}$, as labeled. Different line styles show the correlation at different particle separations $r$; $r=1/256$ (solid lines), $r=1/512$ (dot-dashed lines), and   $r=1/128$ (dashed lines; see text for further information).  We see much
stronger correlations between grains of all sizes in the RDI simulation, with $g$ rising, rather than falling, with decreasing
grain size for similarly sized grains. The correlations also rise more rapidly
as $r$ decreases in the RDI case, suggesting that the clumping will be even stronger in the $r\rightarrow 0$ limit relevant for grain-grain collisions.
}
\label{fig: gfunc sda}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=0.9\columnwidth]{{gfuncs_ca}.pdf}
\caption{As for \cref{fig: gfunc sda}, but for the non-constant-drift simulations.
The two cases are much more similar than the constant-drift cases, although there
remains a modestly higher correlation between larger grains of different sizes in the RDI turbulence.
}
\label{fig: gfunc ca}
\end{center}
\end{figure}

\subsubsection{Grain clumping}
The key measure of relative grain clumping is the Radial Distribution Function (RDF), $g(r;\abc{l},\abc{h})$,
which measures the relative probability of finding a grain of size $\abc{l}$ a distance $r$ from a grain
of size $\abc{h}$. It is normalized such that a spatially homogenous distribution satisfies $g(r;\abc{l},\abc{h})=1$.
If $g(r;\abc{l},\abc{h})>1$ for small $r$, the collision rate of $\abc{l}$ grains with $\abc{h}$ grains will
be enhanced compared if their distribution were uncorrelated \citep{Maxey1987,Squires1990}. As well as grain collisions, a highly
clumped grain distribution could have interesting implications  for the  opacity, which  may
be significantly reduced  compared to a homogenous grain distribution if photons are primarily scattered around the low-density
regions between dust clumps \citep[see][]{Uli2021RDI}.

We illustrate the RDF for the constant- and non-constant-drift simulations in \cref{fig: gfunc sda,fig: gfunc ca},
comparing the RDI and no-backreaction forced simulations in the top and bottom panels, respectively. We represent $g$
using the same method as \citet{Pan2014b}, setting $\abc{l}<\abc{h}$ and plotting $g(r;\abc{l},\abc{h})$ as a function
of  $\abc{h}$ for a variety of grain-size  ratios  $\abc{h}/\abc{l}$. Grains are binned before
computing $g$ according to the method of \cref{sub: dust binning}. The different line styles show
different $r$, which are computed by including only those grains that lie a distance $r\pm\Delta r$ from
each other\footnote{The turbulence itself is anisotropic with respect to the drift direction, meaning that $g$ can differ
depending on whether $r$ is in the perpendicular or parallel direction. However, we saw only minor
differences, when taking this into account, so show only the isotropic version here.}; the solid lines show $r=L/256$  (the equilibrium gas-particle spacing) and $\Delta r =L/1024$, the dashed lines show $r=L/128$ and $\Delta r =L/4096$,
and the dot-dashed lines show $r=L/512$ and $\Delta r =L/1024$
(note that a relatively wider bin in $r$ is needed when $r$ itself is smaller in order to obtain sufficient particle statistics).
Let us first consider the externally forced turbulence constant-drift case, since this is most directly comparable
to previous work, broadly following the expected behavior (c.f.~figure 1 of \citealt{Pan2014b}). The strongest
turbulent clumping is seen for the largest grains, which is  consistent with the well-documented observation that
clumping is strongest for particles with ${\rm St}_{\nu}\simeq 1$ \citep{Squires1990,Sundaram1997} and the estimate in \cref{eq: stokes in simulations} that
all particles have ${\rm St}_{\nu}\lesssim 1$. The maximum of $g\approx 2$ is less than some previous works, and it
is also clear from the trend that larger grains would clump more strongly. There are a number
of possible causes for this discrepancy: in our driven-turbulence simulations the grains are streaming through the turbulence (unlike previous studies, which used $\aext=0$), which could interfere with the ``sling'' effect that causes turbulent clustering; or the effective viscous scales could be underestimated in the estimates of \cref{subsub: stokes numbers}, thus overestimating ${\rm St}_{\nu}$\footnote{This is supported by examination of the spectra in \cref{fig: spectra}: the velocity spectrum steepens at $k\approx 150/L\approx 24(2\pi/L)$, which is well above the scale of twice the gas particle
spacing as used in the estimates of \cref{eq: stokes in simulations}.}; or, by not including explicit viscosity, we may not be accurately simulating the sub-viscous-scale flows that determine the clumping of the ${\rm St}_{\nu}\simeq 1$ grains. The decrease in clumping between different-sized particles (compared to particles of the same size) is similar to that shown in \citet{Pan2014b}.

 The  constant-drift RDI  contrasts significantly to the externally forced turbulence. Most obviously, there is much higher clumping for all
 grain sizes, but particularly for the smallest grains and in the relative correlations between grains of different sizes. This
 strong clumping of small grains is notable given that our simple estimates suggested they have very small Stokes numbers
  (${\rm St}\lesssim 0.01$; c.f. driven case in \cref{fig: gfunc sda}). We also see that our measurements are far from converged, meaning that grains are increasingly clumped at smaller scales, even at scales well below the gas-particle spacing, strongly suggesting that higher-resolution
  simulations (or reality) would increase $g$ further. Finally, it is worth noting that the general shape of $g(r;\abc{l},\abc{h})$ with $\abc{h}$
  is quite different to those seen in standard hydro turbulence; at least with $\abc{l}\approx\abc{h}$, $g$ becomes independent of $\ad$ for larger $\ad$, rather than decreasing towards $g=1$ at either small or large $\ad$.

  The comparison of the non-constant-drift simulations (\cref{fig: gfunc ca}) tells a less interesting story,
  showing broadly similar grain distributions between the RDI and externally forced turbulence, and less clumping than the
  constant-drift RDI. This is consistent with only the largest grains driving interesting dynamics in the non-constant-drift RDI,
  while the smaller grains are primarily passively advected by the flow. The rather low values of $g$ in these cases are likely
  related to the fact that all grains have a slightly
  different equilibrium streaming velocity $\ws$, meaning that small-scale clumps can be quickly destroyed by grains streaming
  away from each other (this is discussed further below). The convergence to $g\approx1.1$, rather than $g=1$,  is simply because the gas
  is  compressible in these simulations, so that there are always some spatial correlations between grains of different sizes that
  arise due to their mutual  correlation with the gas density (see \cref{fig:2d hist ca}).

\begin{figure}
\begin{center}
\includegraphics[width=1\columnwidth]{{wmag_sda}.pdf}
\caption{Average collision velocity between grain pairs, normalized by the gas velocity dispersion,
in the constant-drift simulations. Solid lines show RDI turbulence and dashed lines show the externally forced turbulence simulation without backreaction (colors show grain-size ratios, as in \cref{fig: gfunc sda}). We see significantly lower collision velocities in RDI turbulence, particularly between grains of different sizes. Coupled with the larger grain clumping factors in RDI-generated turbulence (\cref{fig:2d hist sda,fig: gfunc sda}), this would significantly enhance
grain coagulation compared to standard estimates.
}
\label{fig: wmag sda}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=1\columnwidth]{{wmag_ca}.pdf}
\caption{Average collision velocity between grain pairs, normalized by the gas velocity dispersion,
 in the non-constant-drift RDI simulation (solid lines). The dotted lines show the simple prediction \cref{eq: grain collision velocity cda}. In this case, the collision
velocities are entirely dominated by the different drift velocities of differently sized grains, and the driven and RDI
cases look very similar.
}
\label{fig: wmag ca}
\end{center}
\end{figure}

\subsubsection{Relative velocities}

The second component that is required to estimate the rate and outcome of grain-grain collisions is the relative collision velocity.
Depending on subtle choices related to the definition
of the collision kernel, there are a number of relevant measures of collision velocities that can give slightly different results \citep{Wang2000,Pan2014b}. We measure the mean of the value of the
radial velocity with which particle pairs approach each other, which is needed to compute the ``spherical formulation'' of the collision kernel (unlike, for example, the root-mean-squared relative velocity;  \citealp{Wang2000}). More precisely, for two particles
with velocities $\bm{v}_{1}$ and $\bm{v}_{2}$ and separation vector $\bm{r}$, we define $\bm{w}_{p} = \bm{v}_{2}-\bm{v}_{1}$ and $w_{p,r} = (\bm{v}_{2}-\bm{v}_{1})\cdot\bm{r}/|\bm{r}|$ and take an average $\langle |w_{p,r}|\rangle_{<0}$ over only those pairs that are approaching each other.
As for the particle RDFs, we plot this separately for different ratios of grain sizes, viz., measure the collision velocities
of large and small particles, as well as those of similarly sized particles. Collision velocity
statistics also depend on, and are unconverged in, the particle separation  $r$ (i.e., $\langle |w_{p,r}|\rangle_{<0}=\langle |w_{p,r}(r;\abc{l},\abc{h})|\rangle_{<0}$; \citealp{Pan2014b}). By default we use $r=L/256$
(solid lines in \cref{fig: gfunc sda,fig: gfunc ca}) and comment on this where appropriate. As for the RDFs,
the collision velocity statistics are not isotropic in either $\bm{w}_{p}$ or $\bm{r}$ (in either the RDI or the no-backreaction cases); however, examination of anisotropy of collisions has not yielded any interesting insights, so we show only the isotropic versions here.

The constant-drift simulations are shown in \cref{fig: wmag sda}. In this case all grains drift at the same average speed, so there is no
direct contribution to the collision velocity from the grain's equilibrium drift velocity. We
normalize to gas root-mean-squared velocity in order to make the comparison of the RDI to driven turbulence as apt as
possible (see \cref{subsub: turb driving level}). Overall, we see a modestly lower collision velocity in the RDI for similar
sized grains, but the difference becomes more significant for collisions between grains of different sizes.
In other words, collisions between large and small grains are significantly slower on average in RDI turbulence
compared to standard externally forced turbulence. While this might have been anticipated based on our intuitive understanding
that the constant-drift RDI involves gas motions driven by a wide range of dust grains at once, it could
have interesting implications for the outcome of grain collisions, for instance, by increasing the size to which grains
can grow by sticking  \citep{Blum2018}. Finally, it is worth noting that the collision velocities depend relatively
strongly on the particle separation $r$, but this dependence is similar for the RDI and externally forced turbulence so we do not consider
it in detail here (see \citealt{Pan2014} for extensive discussion).

The non-constant-drift simulations, which are less interesting, are shown in  \cref{fig: wmag ca}. In these cases, $\ws\sim \ad^{1/2}$ is different for each grain size and dominates over the
turbulent dust velocity dispersion (see \cref{fig:time.evolution.ca}). Grains of different sizes thus collide primarily
due to their differing drift velocities, with
\begin{equation}
\bm{w}_{p}=\bm{v}_{2}-\bm{v}_{1}\approx [\ws(\ab{2})-\ws(\ab{1})]\hat{\bm{z}}.\label{eq: grain collision velocity cda}
\end{equation}
For this reason, unrelated to properties of the turbulence, the RDI and externally forced turbulence simulations produce nearly
identical  $\langle |w_{p,r}|\rangle_{<0}$ and we plot only the RDI case in \cref{fig: wmag ca}. The results  match the estimate \cref{eq: grain collision velocity cda} nearly perfectly (dotted lines; there is an extra
geometric factor of $1/2$ because only the radial velocity component is computed). While this
is not surprising, it is nonetheless a potentially relevant physical effect that will significantly enhance the collision
rate and velocities of grains in outflowing winds in the non-constant-drift regime.

In addition
to the mean collision velocities, the PDF of $w_{p}$ is of interest: rare events can have an important impact on the
growth or destruction of grains, for example, by allowing a population to grow beyond particularly important size scales
\citep{Garaud2013,Pan2014a}. We show the PDF of $|\bm{w}_{p}|$, $\mathcal{P}(|\bm{w}_{p}|;\abc{l},\abc{h};r)$ for a variety of grain-size pairs in \cref{fig: wpdf sda},
illustrating their   similar shapes in  the RDI and externally forced turbulence for the constant-drift regime.\footnote{The non-constant-drift
PDFs are dominated by the mean drift velocities, so we do not show them here.} Given their clear differences in small-scale
structure and clumping mechanisms (e.g., \cref{fig: gfunc sda}) this is surprising. One minor difference is a slightly
steeper high-$w_{p}$ tail (and a slightly flatter low-$w_{p}$ tail) in the RDI (see light-blue and yellow curves), indicating that high-velocity collisions between differently
sized particles are even less likely than suggested by the collision-velocity average in \cref{fig: wmag sda}. However,
this seems to be a relatively minor effect.

\begin{figure}
\begin{center}
\includegraphics[width=1\columnwidth]{{wpdfs_sda}.pdf}
\caption{PDF of grain collision velocities $|\bm{w}_{p}|$, at particle separation $r=L/256$, for a selection of different grain size pairs
in the constant-drift simulations. Solid lines show the RDI, dashed lines show the forced turbulence simulation (see
also \cref{fig: wmag sda}).
}
\label{fig: wpdf sda}
\end{center}
\end{figure}

\subsubsection{The collision rate}

The collision rate between grains of size $\ab{i}$ and $\ab{j}$ is given by  $\nu_{c,ij}=\langle n_{d,i}\rangle\langle n_{d,j}\rangle\Gamma_{ij}$,
where $\langle n_{d,i}\rangle$ is the number density of species $i$ and $\Gamma_{ij}$ is the collision
kernel between $i$ and $j$. Using the ``spherical formulation'',\footnote{The alternative ``cylindrical'' formulation gives
very similar results in turbulence simulations, but we have not considered it here \citep{Wang2000,Pan2014b}. } $\Gamma_{ij}=2\pi d^{2} g(d;\ab{i},\ab{j})\langle |w_{p,r}(d;\ab{i},\ab{j})|\rangle_{<0}$, where $d$ is sum of the radii of the grains. This shows
that $\nu_{c}$ and/or $\Gamma_{ij}$  may be inferred (aside from the $d$ dependence) from \crefrange{fig: gfunc sda}{fig: wmag sda}. We see
that in the non-constant-drift RDI, $\nu_{c}$ is large and strongly dominated by the effect of the
mean drift; this situation will involve a large number of very high velocity collisions between grains of different sizes. In contrast,
the constant-drift RDI collision rate is  larger than that of externally  forced turbulence  for similarly sized grains (by a factor $\gtrsim 4$ for small grains), and similar
for differently sized grains (since the larger $g$ cancels with the smaller $\langle w_{p,r}\rangle_{<0}$).
However, the RDI collision velocity is significantly smaller, increasing the probability of slow collisions that lead to grain
growth as opposed to bouncing, cratering, or fragmentation.

\subsection{Discussion: extensions, limitations, and future work}

\revchng{The parameter space of possible  RDIs is very large  \citep{Hopkins2020} and a key limitation
of our study has been the focus on just two parameter sets of the acoustic RDI for simulation. That said, for the RDI  with
neutral gas and grains (acoustic RDI turbulence), most of the likely dependencies on other parameters can be inferred from the results here and those of \msh. For larger effective scales (smaller $(2\pi/L)c_{s}t_{s}$), the linear results in \cref{sub: low k linear} show that
the size distribution of grains becomes unimportant to the instability, suggesting that the non-constant-drift and constant-drift instabilities
will both behave similarly to the larger scales of the non-constant-drift simulation presented here, or similar cases in \msh.
For smaller effective scales (larger $(2\pi/L)c_{s}t_{s}$), the constant-drift cases will likely behave similarly but with stronger relative
clumping for less virulent gas turbulence (\msh), while the non-constant-drift cases will be limited to less virulent quasi-resonant modes involving
only the larger grains (\cref{sub: quasi  modes}). At smaller dust-to-gas-mass ratio. the results of \msh\ suggest that the gas turbulence will become less virulent, but likely cause more {relative} dust clumping (the non-constant-drift RDI will also be less virulent at smaller $\mu$, with the
resonant instability limited to a smaller range of the largest grains; \cref{sub: quasi  modes}). At larger dust-to-gas-mass ratio,
there is no qualitative change to main features the acoustic RDI (as occurs for the streaming instability; \citealp{Youdin2005,Squire2020a});
rather, the larger mass of the dust simply drives stronger gas turbulence (\msh).
Finally, it is also worth mentioning that while GIZMO seems to be able to capture the
linear growth rates of the polydisperse acoustic RDI relatively accurately (see, e.g., \cref{fig: nl lin comparison k theta}),
exploring the detailed convergence to linear predictions in different regimes with different grain-size distributions is a complex task
beyond the scope of this work (see, e.g., \citealt{Paardekooper2021,Zhu2021}). While unlikely to affect our results here, given the
dominance of the large-scale modes in the non-constant-drift simulation, there may be important effects at smaller scales and/or smaller $\mu$,
and a more detailed study of numerical convergence and/or comparison to other codes would be important for exploring such cases.  }

\revchng{We have also not studied in detail the (likely common)
situation where the spectrum of dust grains covers both the constant- and non-constant-drift regimes. Based on linear calculations (see \cref{app: linear}), it is reasonable to
surmise  that this will behave like the constant-drift RDI for the relevant range of grains (albeit with
a somewhat reduced dust-to-gas-mass ratio).  }

\revchng{A rather technical source of uncertainty in our results relating to dust-dust collisions concerns the approach to the
zero-particle-separation limit. As discussed extensively in \citet{Pan2013,Pan2014b,Haugen2021}, a quantitative measurement
of the collision kernel requires a careful convergence study in particle separation, which is not achieved here (e.g., compare the different line styles in
\cref{fig: gfunc sda}). This must take into
account both  explicit dissipative effects, which damp gas motions at small scales, and the different physical contributions
to relative grain-grain velocities (the effect of ``caustics'' on dust motion). This is well beyond the scope of this work and
likely a  difficult task for realistic  application to astrophysical objects with RDI-generated turbulence. In particular,
as well as a better theoretical framework with which to understand the strong clumping even at small Stokes
numbers, the acoustic RDI effectively forces the turbulence down to the viscous scales, changing in character with scale. These issues
make it difficult to apply Stokes-number-based self-similarity  arguments, as usually applied to understand passive grains in turbulence. Thus, the result that  constant-drift RDI turbulence should
be particularly a effective nursery for grain growth is qualitative at this stage, and we refrain from making quantitative
estimates for the collision kernel. }

\revchng{That said, the most  relevant and important uncertainty of our study is our
neglect of grain charge and gas magnetization. The neutral cases explored here provide
a reasonable approximation for cooler, denser regions, for instance around  AGB stars, or in some parts of the cool ISM;
but, in most astrophysical scenarios where RDIs apply, grain charge and gas magnetization are expected to play a key role.
Further, magetized-RDI turbulence  causes significantly stronger clumping than acoustic-RDI turbulence, even
in the regime where drag forces are stronger than Lorentz forces and/or in the non-constant-drift
regime.\footnote{It remains unclear exactly why this is, although
there are many more resonant instabilities available in the magnetic case due to the wider variety of waves
and the dust's gyromotion.} A selection of cases is presented in \cite{Hopkins2021}, mostly focusing
on  the more complex
situation of a stratified wind driven from the base: extremely strong clumping is seen even in non-constant-drift
cases (their ``\textbf{-Q}'' simulations; see e.g., their figure 9) in stark contrast to our results here,  a difference that is related to the magnetization and not the
stratification of the system (see e.g., their figure 13; this is also seen in unstratified magnetized simulations).
From these results it is clear that magnetization will be a key parameter for RDI-turbulence induced grain clumping;
however, given the complexity of these cases and the yet wider parameter space to explore, we leave such studies to future work. }

\revchng{Finally, it is worth mentioning that by neglecting explicit stratification of the system, we are also
neglecting other possible gas and/or radiative instabilities that can occur in dust-driven winds for some systems \citep[see, e.g.,][]{Woitke2006,Krumholz2012}.
Generically, such instabilities operate on much larger scales than the acoustic RDIs considered here, although it is plausible
that they could create turbulence that influences RDI development in some circumstances.
Again, such issues are addressed more explicitly in the stratified simulations of \cite{Hopkins2021}.}

\section{Conclusions}

This paper has presented an in-depth study of the   ``acoustic Resonant Drag Instability''  (\citealp{Squire2018}; \hs), which
is driven by the interaction of an outflowing population of dust grains with compressible gas motions.
The acoustic RDI is expected to operate in a variety of astrophysical scenarios, for example,
in the presence of a radiation source that couples more strongly
to dust grains than to the gas. This accelerates grains outwards, often to supersonic velocities,
which (in addition to driving a gas outflow) destabilizes the RDI  (\hs).
In cooler, denser gas, such as that in molecular clouds, AGB-star winds, or around AGN, the grain charge and MHD effects
are not necessarily dominant, and the system may be well approximated by considering a neutral gas and neutral grains (the ``acoustic'' RDI).
The novel feature of this work has been the inclusion of a wide spectrum of grain sizes -- a factor of $100$ in grain
radius $\ad$ in our numerical simulations -- which has not been included in previous
studies but is clearly an important feature of realistic  systems \citep{Draine2010}.
In our numerical study, rather than surveying a wide parameter space of different simulations, we
have focused in detail on two representative cases that can apply adequately well to a variety of astrophysical processes. By
comparing directly to simulations of externally forced turbulence, this allowed
us to consider in more detail aspects of dust and gas structure in the RDI, and how these might
influence important processes such as grain growth.

With a spectrum of grains, there are two qualitatively different regimes of the RDI depending on
how the grain's acceleration (imparted, for example, by an external radiation field) scales with
the grain radius. In the \emph{constant-drift} regime, applicable for grains larger than the wavelength
of a radiation field ($\ad\gtrsim \lambda_{\rm rad}$), all grains drift through the gas with the
same velocity in the (quasi-)equilibrium. In the  \emph{non-constant-drift} regime, applicable for small
grains in long-wavelength radiation fields ($\ad\lesssim \lambda_{\rm rad}$), or when
the acceleration difference between the dust and gas arises from an external force on the gas, the grain drift velocity increases with grain size, stretching across  a wide range of values in realistic scenarios.

We show in \cref{app: linear}, which presents analytic and numerical calculations of the {linear} RDI growth rate
with a spectrum of grains, that these different regimes strongly
influence the behavior of the acoustic RDI. Generally, the RDI is more virulent, with faster growth rates and behavior
that is very similar to the single-grain-size case, in the constant-drift regime. Our simulations
show that this linear behavior carries over into the nonlinearly turbulent regime also:
 the non-constant-drift RDI, although strongly unstable with a saturated state that shares key features
 of the linear instability, develops into turbulence without strong
 correlations between grains of different size (like externally forced turbulence);
the constant-drift RDI is very different, with much stronger spatial correlations between
small grains and those of different sizes, along with slower grain-grain collision velocities.   These
differences, which imply a high rate of low-velocity collisions in the constant-drift
RDI, suggest that  constant-dust-drift  outflows could be highly effective sites
for dust  growth  through collisions, while the opposite is likely true in the non-constant-drift regime because
grain-grain collisions are dominated  by the (fast) mean drift between grains of different sizes.
\revchng{Another interesting conclusion in the constant-drift regime
is the  strong clumping of small grains, even though their turbulent Stokes numbers remain well below one (see \cref{subsub: stokes numbers}).
This highlights the fact
 that the clumping mechanism in RDI-generated turbulence is quite different -- and much more efficient for similar velocity dispersions --  to standard (Kolmogorov) turbulence \citep{Pumir2016},
 even though the turbulent velocity spectra are relatively similar (\cref{fig: spectra}).  }

\acknowledgements{We thank Eric Moseley and Darryl Seligman for helpful discussion.
Support for JS  was
provided by Rutherford Discovery Fellowship RDF-U001804 and Marsden Fund grant UOO1727, which are managed through the Royal Society Te Ap\=arangi. Support for JS, PFH, and SM was provided by NSF Collaborative Research Grants 1715847 \& 1911233, NSF CAREER grant 1455342, and NASA grants 80NSSC18K0562 and JPL 1589742. Numerical simulations were
run on the Caltech compute cluster ``Wheeler,''  and with  allocation TG-AST130039 from the
 Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562.}

\paragraph*{Data availability}\revchng{The simulation data presented in  this article is available on  request to JS. A public version of the GIZMO code is available at \href{http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html}{http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html}.}

redpajama_set_name: RedPajamaArXiv