Simulating the Diverse Instabilities of Dust in Magnetized Gas

TL;DR

This paper uses simulations to explore the complex instabilities of charged dust in magnetized gas, revealing diverse turbulence, structure formation, and magnetic field amplification across various parameter regimes.

Contribution

It provides a comprehensive simulation study of dust-gas instabilities, characterizing their turbulence, structure, and magnetic effects under different conditions, extending prior theoretical work.

Findings

Dust drives turbulence and magnetic field amplification.

Dust forms extreme over-densities and complex structures.

Different regimes produce diverse dust and gas behaviors.

Abstract

Recently Squire & Hopkins showed that charged dust grains moving through magnetized gas under the influence of any external force (e.g. radiation pressure, gravity) are subject to a spectrum of instabilities. Qualitatively distinct instability families are associated with different Alfvenic or magnetosonic waves and drift or gyro motion. We present a suite of simulations exploring these instabilities, for grains in a homogeneous medium subject to an external acceleration. We vary parameters such as the ratio of Lorentz-to-drag forces on dust, plasma $\beta$, size scale, and acceleration. All regimes studied drive turbulent motions and dust-to-gas fluctuations in the saturated state, can rapidly amplify magnetic fields into equipartition with velocity fluctuations, and produce instabilities that persist indefinitely (despite random grain motions). Different parameters produce diverse morphologies and qualitatively different features in dust, but the saturated gas state can be broadly characterized as anisotropic magnetosonic or Alfvenic turbulence. Quasi-linear theory can qualitatively predict the gas turbulent properties. Turbulence grows from small to large scales, and larger-scale modes usually drive more vigorous gas turbulence, but dust velocity and density fluctuations are more complicated. In many regimes, dust forms structures (clumps, filaments, sheets) that reach extreme over-densities (up to $\gg 10^{9}$ times mean), and exhibit substantial sub-structure even in nearly-incompressible gas. These can be even more prominent at lower dust-to-gas ratios. In other regimes, dust self-excites scattering via magnetic fluctuations that isotropize and amplify dust velocities, producing fast, diffusive dust motions.

Figures

Tables

Table 1. Table 1: The default initial conditions for the simulations studied in this paper. Each simulation (by default) adopts an isothermal gas equation-of-state, Epstein drag, constant grain charge, and follows a single population of grains, with resolution 2 × 128 3 2 superscript 128 3 2\times 128^{3} elements (equal number gas and dust), and total dust-to-gas mass ratio μ ≡ ρ d 0 / ρ g 0 = 0.01 𝜇 superscript subscript 𝜌 𝑑 0 superscript subscript 𝜌 𝑔 0 0.01 \mu\equiv\rho_{d}^{0}/\rho_{g}^{0}=0.01 . These choices are varied below. Columns show: (1) Simulation name (used throughout). (2) | 𝐰 s 0 | / c s 0 superscript subscript 𝐰 𝑠 0 superscript subscript 𝑐 𝑠 0 {|{\bf w}_{s}^{0}|}/{c_{s}^{0}} : The initial equilibrium drift velocity, in units of the sound speed. (2) τ ≡ t s 0 / t L 0 𝜏 superscript subscript 𝑡 𝑠 0 superscript subscript 𝑡 𝐿 0 \tau\equiv t_{s}^{0}/t_{L}^{0} : Ratio of Lorentz to drag forces (stopping time to gyro time). (3) | 𝐁 ^ 0 ⋅ 𝐚 ^ | ⋅ subscript ^ 𝐁 0 ^ 𝐚 |\hat{\bf B}_{0}\cdot\hat{{\bf a}}| : Angle between initial magnetic field and direction of differential acceleration/force between dust and gas, 𝐚 ≡ 𝐚 ext , dust − 𝐚 ext , gas 𝐚 subscript 𝐚 ext dust subscript 𝐚 ext gas {\bf a}\equiv{\bf a}_{\rm ext,\,dust}-{\bf a}_{\rm ext,\,gas} . (4) L box / c s 0 ​ t s 0 subscript 𝐿 box superscript subscript 𝑐 𝑠 0 superscript subscript 𝑡 𝑠 0 {L_{\rm box}}/{c_{s}^{0}\,t_{s}^{0}} : Box size in dimensionless units. (5) β ≡ P gas 0 / P B 0 = 2 ​ ( c s 0 / v A 0 ) 2 𝛽 superscript subscript 𝑃 gas 0 superscript subscript 𝑃 B 0 2 superscript superscript subscript 𝑐 𝑠 0 superscript subscript 𝑣 𝐴 0 2 \beta\equiv P_{\rm gas}^{0}/P_{\rm B}^{0}=2\,(c_{s}^{0}/v_{A}^{0})^{2} : Initial ratio of gas thermal to magnetic pressure. (6) Notes: we provide an example physical regime where these parameters are plausible for typical interstellar grains (from Paper I ). We also quote an equivalent set of dimensionless parameters: the “acceleration parameter” a ~ ≡ | 𝐚 | ​ ρ ¯ grain i ​ ϵ grain / ( ( c s 0 ) 2 ​ ρ g 0 ) ~ 𝑎 𝐚 superscript subscript ¯ 𝜌 grain 𝑖 subscript italic-ϵ grain superscript superscript subscript 𝑐 𝑠 0 2 superscript subscript 𝜌 𝑔 0 \tilde{a}\equiv|{\bf a}|\,\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}/((c_{s}^{0})^{2}\,\rho_{g}^{0}) , “charge parameter” ϕ ~ ≡ 3 ​ Z grain 0 ​ e / ( 4 ​ π ​ c ​ ϵ grain 2 ​ ( ρ g 0 ) 1 / 2 ) ~ italic-ϕ 3 superscript subscript 𝑍 grain 0 𝑒 4 𝜋 𝑐 superscript subscript italic-ϵ grain 2 superscript superscript subscript 𝜌 𝑔 0 1 2 \tilde{\phi}\equiv 3\,Z_{\rm grain}^{0}\,e/(4\pi\,c\,\epsilon_{\rm grain}^{2}\,(\rho_{g}^{0})^{1/2}) , and “size parameter” α ~ ≡ ρ ¯ grain i ​ ϵ grain / ρ g 0 ​ L box ~ 𝛼 superscript subscript ¯ 𝜌 grain 𝑖 subscript italic-ϵ grain superscript subscript 𝜌 𝑔 0 subscript 𝐿 box \tilde{\alpha}\equiv\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}/\rho_{g}^{0}\,L_{\rm box} .

Name	$\frac{\left|{𝐰}_s^0\right|}{c_s^0}$ ($\tilde{a}$)	$\tau$ ($\tilde{\varphi }$)	$\left|{\widehat{𝐁}}_0\cdot \widehat{𝐚}\right|$	$\frac{L_{\mathrm{box}}}{c_s^0{t}_s^0}$ ($\tilde{\alpha }$)	$\beta$	Notes
Example	0.89 (25)	29 (50)	$0.05$	0.34 (5)	2	Case study in Seligman et al. (2018)
AGB	3 (8.4)	1e-3 (2.8e-3)	$1/\sqrt{2}$	0.01, 3.1, 920, 3e5	2	AGB wind, $r\sim 100R_{\odot }$, $\dot{M}\sim {10}^{-6}{M}_{\odot }{\mathrm{yr}}^{-1}$
(S/M/L/XL)				(270, 0.9, 0.003, 1e-5)		$L_{\mathrm{box}}\sim$ (5 km, 2000 km, 1 $R_{\odot }$, 300 $R_{\odot }$)
HII-near	4 (20)	2.3 (24)	$1/\sqrt{2}$	0.0088, 2.6, 1000	20	HII region, $r=0.1\mathrm{pc}$, $n={10}^4{\mathrm{cm}}^{-3}$
(S/M/L)				(390, 1.3, 0.0034)		$L_{\mathrm{box}}\sim$ (0.2 au, 50 au, 0.1 pc)
HII-far	0.15 (0.45)	3.5 (16)	$1/\sqrt{2}$	0.0032, 1.0, 330	20	HII region, $r=1\mathrm{pc}$, $n={10}^2{\mathrm{cm}}^{-3}$
(S/M/L)				(500, 1.7, 0.005)		$L_{\mathrm{box}}\sim$ (0.7 au, 200 au, 0.3 pc)
WIM	0.05 (0.12)	100 (160)	$1/\sqrt{2}$	4.8e-4, 0.21, 100	2	WIM, $n\sim 1{\mathrm{cm}}^{-3}$, $T\sim {10}^4$K
(S/M/L)				(3400, 7.6, 0.017)		$L_{\mathrm{box}}\sim$ (1 au, 500 au, 1 pc)
Corona	20 (480)	3200 (1700)	$1/\sqrt{2}$	3.3e-4, 0.067, 20	0.002	Solar corona, $r\sim R_{\odot }$
(S/M/L)				(4.8e4, 240, 0.8)		$L_{\mathrm{box}}\sim$ (10 km, 2000 km, 1 $R_{\odot }$)
CGM	9.5 (100)	1.4e5 (3e7)	$1/\sqrt{2}$	0.29 (25)	2000	CGM at $100$kpc from QSO; $L_{\mathrm{box}}\sim 1$kpc

Table 2. Table 2: Dispersion in various gas and dust quantities in the simulations, during the saturated state. We measure the rms dispersion δ ​ X 𝛿 𝑋 \delta X in quantity X 𝑋 X , averaged over the last several snapshots in time for each run. We show gas velocity δ ​ 𝐮 g x , y , z 𝛿 superscript subscript 𝐮 𝑔 𝑥 𝑦 𝑧 \delta{\bf u}_{g}^{x,\,y,\,z} in each direction (see Fig.  3 for axis convention), dust velocity δ ​ 𝐯 d 𝛿 subscript 𝐯 𝑑 \delta{\bf v}_{d} , magnetic field δ ​ 𝐁 𝛿 𝐁 \delta{\bf B} , gas density δ ​ ln ⁡ ρ g 𝛿 subscript 𝜌 𝑔 \delta\ln{\rho_{g}} and dust density δ ​ ln ⁡ ρ d 𝛿 subscript 𝜌 𝑑 \delta\ln{\rho_{d}} . The dispersions shown are mass-weighted for 𝐮 g subscript 𝐮 𝑔 {\bf u}_{g} and 𝐯 d subscript 𝐯 𝑑 {\bf v}_{d} (so that the kinetic energy of gas is just ( 1 / 2 ) ​ M gas ​ | δ ​ 𝐮 g | 2 1 2 subscript 𝑀 gas superscript 𝛿 subscript 𝐮 𝑔 2 (1/2)\,M_{\rm gas}\,|\delta{\bf u}_{g}|^{2} , and likewise for dust), while dispersions for 𝐁 𝐁 {\bf B} and ρ g subscript 𝜌 𝑔 \rho_{g} , ρ d subscript 𝜌 𝑑 \rho_{d} are volume weighted. Usually the mass and volume weights give similar values; where they differ substantially, the value in parenthesis gives the other (for more details, see Figs.  21 - 29 ). For each parameter set, we list the “default” boxes S/M/L/XL as defined in Table  1 , as well as the variants discussed in 4.2 . For Example , we show variations in β 𝛽 \beta , γ 𝛾 \gamma , and μ 𝜇 \mu from Figs.  11 and 16 . For HII-near L, we compare variations in the dust charge (“L: τ 𝜏 \tau =10”; q grain subscript 𝑞 grain q_{\rm grain} increased by a factor ∼ 4 similar-to absent 4 \sim 4 ), as well as different charging models, including un-saturated collisional charging (“L:CC”, q grain ∝ T proportional-to subscript 𝑞 grain 𝑇 q_{\rm grain}\propto T with γ = 5 / 3 𝛾 5 3 \gamma=5/3 ) and photo-electric charging (“L:PE”; q grain ∝ T / ρ 1 / 2 proportional-to subscript 𝑞 grain 𝑇 superscript 𝜌 1 2 q_{\rm grain}\propto T/\rho^{1/2} , γ = 1 𝛾 1 \gamma=1 ), from Fig.  12 . We also compare variations in the angle, 𝐁 0 ⋅ 𝐚 = cos ⁡ θ 𝐁𝐚 ⋅ subscript 𝐁 0 𝐚 subscript 𝜃 𝐁𝐚 {\bf B}_{0}\cdot{\bf a}=\cos\theta_{{\bf B}{\bf a}} with θ 𝐁𝐚 = ( 70 ∘ ,  20 ∘ ) subscript 𝜃 𝐁𝐚 superscript 70 superscript 20 \theta_{{\bf B}{\bf a}}=(70^{\circ},\,20^{\circ}) (“L:70 ∘ ”, “L:20 ∘ ”, respectively) from Fig.  15 , and different dust-to-gas ratios μ 𝜇 \mu from Fig.  18 . For WIM , we compare variations of S/M/L with 5 × \times lower 𝐚 𝐚 {\bf a} and 𝐰 s subscript 𝐰 𝑠 {\bf w}_{s} (“LoV”) from Fig.  14 . For CGM , we show runs with lower τ 𝜏 \tau ( τ ≈ 4500 𝜏 4500 \tau\approx 4500 implying q grain subscript 𝑞 grain q_{\rm grain} is ∼ 30 × \sim 30\times lower; labeled “ τ low subscript 𝜏 low \tau_{\rm low} ”), and runs with un-saturated photo-electric (PE) charging and γ = 5 / 3 𝛾 5 3 \gamma=5/3 (for both the default and low τ 𝜏 \tau cases), from Fig.  13 .

Name	$\frac{\delta {𝐮}_{g,x}}{c_s^0}$	$\frac{\delta {𝐮}_{g,y}}{c_s^0}$	$\frac{\delta {𝐮}_{g,z}}{c_s^0}$	$\frac{\delta {𝐯}_{d,x}}{c_s^0}$	$\frac{\delta {𝐯}_{d,y}}{c_s^0}$	$\frac{\delta {𝐯}_{d,z}}{c_s^0}$	$\frac{\delta {𝐁}_x}{{\left(4\pi P_0\right)}^{1/2}}$	$\frac{\delta {𝐁}_y}{{\left(4\pi P_0\right)}^{1/2}}$	$\frac{\delta {𝐁}_z}{{\left(4\pi P_0\right)}^{1/2}}$	$\delta \ln \left[\frac{{\rho }_g}{{\rho }_g^0}\right]$	$\delta \ln \left[\frac{{\rho }_d}{{\rho }_d^0}\right]$
Example:
Default	0.075	0.025	0.017	0.086	0.063	0.32	0.075	0.025	6.4e-3	0.013	1.2
$\beta$=1	0.067	0.015	6.6e-3	0.095	0.086	0.20	0.067	0.013	3.1e-3	3.7e-3	1.7 (1.4)
$\gamma$=5/3	0.088	0.039	0.020	0.12	0.11	0.33	0.089	0.037	0.013	8.2e-3	1.2 (1.3)
$\mu$=1e-3	4.3e-3	1.8e-3	6.7e-3	2.6e-3	3.2e-3	0.035	5.1e-3	1.8e-3	4.1e-3	4.9e-3	2.7
$\mu$=0.1	0.18	0.15	0.11	0.32	0.30	0.41	0.17	0.12	0.056	0.055	0.89
AGB:
S	0.052	0.050	0.014	6.9e-4	9.4e-4	7.6e-4	8.6e-3	9.4e-3	0.049	0.052	0.75
M	0.031	0.010	0.038	0.10	6.8e-3	0.11	0.066	4.2e-3	0.071	0.074	0.90
L	1.0	0.50	1.2	1.7	0.50	1.9	1.2	0.51	1.2	0.85	0.88
XL	19	4.3	18	20	3.2	20	11	2.6	12	1.1	1.1
HII-near:
S	0.060	0.057	0.013	3.5e-3	3.6e-3	3.7e-3	5.1e-3	4.7e-3	0.023	0.080	1.2
M	0.15	0.13	0.20	0.76	0.80	0.65	0.15	0.11	0.14	0.058	1.4
L	3.3	1.8	3.6	13	15 (31)	12	3.7	1.9	3.5	1.2	1.5 (2)
L:$\tau$=10	4.7 (5.5)	1.6 (2.1)	4.8 (5.7)	9.9 (8.0)	3.0 (3.3)	10 (8.6)	4.1	1.9	4.0	1.3 (1.2)	2.0 (2.4)
L:CC	2.9 (13)	2.6 (12)	2.7 (12)	29 (200)	30 (200)	26 (170)	2.5	1.6	2.6	1.6 (0.9)	1.8 (2.5)
L:PE	6.9 (60)	7.0 (61)	6.3 (51)	80 (1200)	81 (1200)	63 (940)	4.3	4.0	5.0	2.0 (1.1)	2.2 (2.2)
L:70∘	7.2 (8.5)	2.1 (2.4)	4.7 (5.7)	15 (11)	3.4 (3.1)	8.5 (7.4)	4.8	2.2	3.0	1.2	2.2
L:20∘	2.8 (3.1)	1.5 (1.9)	7.1 (7.7)	4.4 (4.2)	2.7 (2.9)	10.7 (9.5)	3.0	1.7	6.8	1.3	1.9 (2.5)
L:$\mu$=1e-3	1.1	0.36	1.2	3.0 (2.3)	1.2	3.2 (2.5)	0.88	0.28	0.95	0.56	2.4 (3.2)
L:$\mu$=0.1	14 (17)	3.4 (5.6)	14 (17)	15 (24)	7.4 (24)	15 (22)	12	3.5	12	1.7	1.7
HII-far:
S	1.8e-3	1.8e-3	2.0e-3	4.7e-5	1.4e-5	2.1e-5	1.4e-3	1.4e-3	1.4e-3	4.3e-4	0.63
M	2.9e-3	2.7e-3	4.5e-3	1.5e-3	1.6e-3	3.5e-3	2.4e-3	1.9e-3	2.3e-3	7.4e-4	1.3 (2.5)
L	0.55	0.20	0.58	1.1	0.90	1.0	0.43	0.19	0.44	0.28	1.5 (2.2)
WIM:
S	9.8e-4	9.8e-4	1.1e-3	6.4e-4	6.4e-4	9.4e-6	9.7e-4	9.7e-4	1.1e-3	7.2e-4	0.49
S:LoV	9.8e-4	9.8e-4	1.1e-3	7.0e-4	7.0e-4	4.3e-6	9.7e-4	9.7e-4	8.6e-4	7.2e-4	0.42
M	9.9e-4	9.8e-4	1.1e-3	1.4e-3	1.4e-3	1.2e-4	9.7e-4	9.7e-4	9.4e-4	7.2e-4	0.36
M:LoV	1.0e-3	1.0e-3	1.1e-3	6.0e-4	6.0e-4	2.0e-4	9.7e-4	9.7e-4	8.5e-4	7.3e-4	0.25
L	0.032	0.021	0.088	5.3	5.3	0.085	3.6e-3	3.6e-3	3.9e-3	3.5e-3	0.37
L:LoV	1.5e-3	1.3e-3	2.3e-3	1.1e-3	7.3e-3	1.2e-3	1.1e-3	1.1e-3	8.0e-4	8.0e-4	0.25
Corona:
S	0.011	0.011	0.013	0.012	0.012	3.5e-5	0.011	0.011	0.013	0.060	0.63
M	0.012	0.012	0.014	0.62	0.63	0.010	0.011	0.011	0.022	0.060	0.86
L	0.27	0.24	0.086	24	24	3.1	0.076	0.069	0.15	0.066	0.57
L:$\tau$=100	0.14	0.05	0.21	0.17	0.13	0.21	0.023	0.023	0.013	0.085	0.6 (0.8)
CGM:
Default	0.83	0.82	0.63	23	23	22	0.040	0.040	0.042	0.34	0.28
$\gamma$=5/3+PE	2.5	2.5	2.0	63	63	63	0.14	0.14	0.13	0.56	0.25
${\tau }_{\mathrm{low}}$	0.031	0.031	0.036	2.0	2.0	1.9	0.024	0.021	0.019	0.021	0.27
${\tau }_{\mathrm{low}}$+PE	0.047	0.046	0.057	2.9	2.9	2.9	0.035	0.032	0.023	0.025	0.27

Equations

$$\frac{\mathrm{d}{\mathbf{v}}_d}{\mathrm{dt}}$$

$${\rho }_g\left(\frac{\partial }{\partial t}+{\mathbf{u}}_g\cdot \nabla \right){\mathbf{u}}_g=$$

$$-\int {\mathrm{d}}^3{\mathbf{v}}_d{f}_d(\mathbf{x},{\mathbf{v}}_d){\mathbf{a}}_{\mathrm{gas}\text{mathchar}45dust}({\mathbf{v}}_d,...),$$

$$P$$

$$c_s^2$$

$$t_s$$

$$t_L$$

$$X_0\equiv ⟨X({\rho }_g^0,{\mathbf{w}}_s^0,...,t=0)⟩$$

$$t_s^{-1}$$

$$Z_{\mathrm{grain}}$$

$$P(q)\propto \exp (-|q/q_0{|}^{\gamma })$$

$${\dot{u}}_0(v/v_{\ast }{)}^{\xi }{\partial }^{2}P(v)/\partial v^2\sim -t_{\mathrm{damp}}^{-1}(v/v_{\ast }{)}^{\psi }P(v)$$

$$\gamma$$

$$P_M\propto \rho P_V$$

$$S$$

$$S_M$$

$${\rho }_0$$

$$C$$

$$⟨{\rho }_g⟩\delta {\mathbf{u}}_g^2/2\sim \delta {\mathbf{B}}^2/8\pi$$

$$[\delta \ln ({\rho }_g/⟨{\rho }_g⟩){]}^2=\ln [1+(b|\delta {\mathbf{u}}_g|/c_s{)}^2],$$

$$|\delta \mathbf{B}|/|\mathbf{B}|\sim |\delta {\mathbf{u}}_g|/v_A$$

$$|\delta {\mathbf{v}}_d|\sim |\delta {\mathbf{u}}_g|,$$

$$\frac{(\mathbf{B}\cdot \nabla )\mathbf{B}}{4\pi }\sim \frac{|\mathbf{B}||\delta \mathbf{B}|}{4\pi L_{\mathrm{box}}}\sim \frac{{\rho }_d{\mathbf{w}}_s}{t_s}\approx {\rho }_d\mathbf{a}$$

$$\hat{\mathbf{k}}\approx \hat{\mathbf{B}},$$

$$|\delta \mathbf{B}|/|\mathbf{B}|\sim |\delta {\mathbf{u}}_g|/v_A.$$

$$|{\mathbf{v}}_d|\sim |\mathbf{a}|t_s\sim |{\mathbf{w}}_s|,$$

Full text

1 Introduction

Almost all astrophysical, planetary, and atmospheric fluids are laden with grains of dust, which play a central role in many astrophysical processes including in planet and star formation; in the attenuation and extinction of observed light; in cool-star, brown dwarf, and planetary evolution; in atmospheric dynamics; in astro-chemistry; in feedback and outflow-launching from star-forming regions, cool stars, and active galactic nuclei (AGN); and in inter-stellar cooling and heating (see Draine, 2003; Dorschner, 2003; Apai & Lauretta, 2010, for reviews). The dynamical interactions between dust and gas therefore are of fundamental importance in astrophysics.

Recently, Squire & Hopkins (2018b) showed that dust-gas mixtures are unstable to a broad class of instabilities, which they referred to as “Resonant Drag Instabilities” (RDIs). The Squire & Hopkins (2018b) instabilities manifest whenever a fluid, gas, or plasma system contains dust streaming with non-zero drift velocity ${\bf w}_{s}\equiv{\bf v}_{d}-{\bf u}_{g}$ relative to the gas (where ${\bf v}_{d}$ and ${\bf u}_{g}$ are the dust and gas velocities, respectively). Although a broad range of wavelengths are unstable, the resonances which produce the most rapidly-growing instabilities occur when the natural frequency of a linear gas mode (e.g. a sound wave, MHD wave, or epicyclic oscillation) matches the natural frequency of a dust mode (e.g. advection, with frequency ${\bf w}_{s}\cdot{\bf k}$, or gyro oscillations). Every such pair of modes produces an instability, with a unique growth rate, resonance, and linear mode structure. Since dust-gas drift may be caused by many external forces, such as radiative absorption or scattering by dust or gas, gravity in quasi-hydrostatic systems, centrifugal or coriolis forces in rotating systems, or large-scale hydrodynamic or pressure forces, these instabilities will develop in a range of astrophysical environments.

In a series of papers, Hopkins & Squire (2018b), Squire & Hopkins (2018a), and Hopkins & Squire (2018a) analytically explored various examples of these instabilities in some astrophysical systems. Hopkins & Squire (2018a) (hereafter Paper I) focused on the case of instabilities involving charged dust in magnetized gas, relevant in the warm interstellar medium (ISM), circum and inter-galactic medium (CGM/IGM), HII regions, supernovae (SNe) ejecta and remnants, the Solar and stellar coronae, cool-star winds, AGN outflows and obscuring torii, and giant molecular clouds (GMCs). They showed that a variety of instabilities appear with different properties and growth rate scalings, even in the case of a homogeneous gas obeying ideal MHD (a good approximation in most of these regimes), with a single group of grains interacting via drag and Lorentz forces. However, their analysis was restricted to analytic, linear perturbation theory. Moseley et al. (2018) presented simulations of the un-magnetized and un-charged instabilities in the non-linear regime. Seligman et al. (2018) presented a case study of one example in the magnetized regime, and found that the introduction of a magnetic field produced novel dust behaviors and outcomes in both the linear and non-linear regimes of the instability. That first study necessarily neglected much of the large parameter space.

In this paper, we present a large survey of $\sim 40$ simulations111Animations and additional visualizations of the simulations here are available at \hrefhttp://www.tapir.caltech.edu/ phopkins/Site/animations/dust-and-gas-in-astrophysic\urlhttp://www.tapir.caltech.edu/ phopkins/Site/animations/dust-and-gas-in-astrophysic that explore the non-linear regime of these instabilities in a representative range of the astrophysically relevant parameter space for charged dust in magnetized gas. These idealized experiments inform our understanding of the mechanisms responsible for the growth and saturation of the instabilities, the non-linear structure of the dust and gas, and the potential theoretical and observational ramifications. They are complex because the instabilities depend on six dimensionless parameters, and as shown in Paper I, at any given wavenumber ${\bf k}$, the linear dispersion relation typically features $\sim 3-7$ unstable modes (each of which has growth rates that depend on the mode angle $\hat{\bf k}$, at a given $|{\bf k}|$). This inherent complexity further underscores the necessity for numerical simulations that explore different non-linear regimes. We show that a diverse variety of behaviors arise, depending both on the physical parameters of the system and the spatial scales studied, all of which may have important astrophysical consequences.

This paper is organized as follows. § 2 presents our methodology, and § 3 discusses the parameter space surveyed (see also Table 1 and Fig. 1). § 4 presents several results from the simulations (e.g. morphologies, saturated fluctuation amplitudes and PDFs). § 5 discusses these results in more detail and compares them to theoretical expectations, attempting to identify classes of saturation mechanisms. We summarize and conclude in § 7.

2 Methods & Simulation Setup

2.1 Numerical Methods & Equations Solved

The numerical methods adopted here have been described in detail in Moseley et al. (2018); Seligman et al. (2018), and we briefly summarize them here. Our simulations were run with the code GIZMO (Hopkins, 2015),222A public version of the code, including all methods used in this paper, is available at \hrefhttp://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html\urlhttp://www.tapir.caltech.edu/ phopkins/Site/GIZMO.html using the Lagrangian “meshless finite volume” (MFV) method for MHD, which has been extensively tested on problems involving multi-fluid MHD instabilities, MRI, shock-capturing, and more (Hopkins & Raives, 2016; Hopkins, 2016b, 2017; Su et al., 2017). Grains are integrated using the “super-particle” method (see, e.g. Carballido et al., 2008; Johansen et al., 2009; Bai & Stone, 2010b; Pan et al., 2011), whereby the motion of each dust “particle” in the simulation follows Eq. (1) below, but each represents an ensemble of dust grains with similar size, mass and charge (denoted $\epsilon_{\rm grain}$, $m_{\rm grain}$, $q_{\rm grain}$, respectively). The numerical methods used for the integration are described and tested in Hopkins & Lee (2016); Lee et al. (2017). The back-reaction is accounted for as in Moseley et al. (2018) (see App. B) and the the Lorentz forces is evolved using a Boris integrator.

Each individual grain (dust super-particle) in the code obeys

[TABLE]

where ${\bf a}_{\rm ext,\,dust}$ is a constant external acceleration, $t_{s}$ is the drag coefficient or “stopping time” and $t_{L}$ the gyro or Larmor time. The gas obeys the ideal MHD equations – the standard advection equation, $\partial\rho_{g}/\partial t=-\nabla\cdot(\rho_{g}\,{\bf u}_{g})$, for the gas density $\rho_{g}$, and the standard induction equation, $\partial{\bf B}/\partial t=\nabla\times({\bf u}_{g}\times{\bf B})$, for the magnetic field ${\bf B}$ – with the addition of a back-reaction force from the grains in the momentum equation. In particular, whenever drag or Lorentz forces exert a force $m_{\rm grain}\,\mathrm{d}{\bf v}_{d}/\mathrm{dt}$ on a grain within a given gas cell, an equal-but-opposite force is applied to the gas. This is treated as a usual momentum flux within GIZMO, which numerically guarantees exact force balance and total momentum conservation. The gas momentum equation reads

[TABLE]

where $f_{d}({\bf x},\,{\bf v}_{d})$ is the phase-space density distribution of dust (i.e. differential mass of grains per element $\mathrm{d}^{3}{\bf x}\,\mathrm{d}^{3}{\bf v}_{d}$) and ${\bf a}_{\rm ext,\,gas}$ is an external gas acceleration that we set to zero in the simulations here. The gas obeys an exactly polytropic equation of state with thermal pressure $P$ and sound speed $c_{s}$:

[TABLE]

In our default simulations we assume Epstein drag, which can be approximated to very high accuracy with the expression (valid for both sub and super-sonic drift)

[TABLE]

where $\bar{\rho}_{\rm grain}^{\,i}$ and $\epsilon_{\rm grain}$ are the internal grain density and radius, respectively. The Larmor time is:

[TABLE]

where $m_{\rm grain}$ and $q_{\rm grain}=Z_{\rm grain}\,e$ are the grain mass and charge.

2.2 Initial Conditions

We initialize a periodic, cubic box of side-length $L_{\rm box}$ with uniform gas density $\rho_{g}^{0}\equiv M_{\rm gas,\,box}/L_{\rm box}^{3}$ and dust density $\rho_{d}^{0}\equiv\mu\,\rho_{g}^{0}$ (dust-to-gas ratio $\mu$), gas velocity ${\bf u}_{g}^{0}=0$, and dust drift ${\bf w}_{s}^{0}=|{\bf a}|\,t_{s}^{0}\,(1+\mu)^{-1}\,(1+\tau^{2})^{-1}\,[\hat{{\bf a}}-\tau\,(\hat{{\bf a}}\times\hat{\bf B}_{0})+\tau^{2}\,(\hat{{\bf a}}\cdot\hat{\bf B}_{0})\,\hat{\bf B}_{0}]$ (see Paper I, section 3.1). Here

[TABLE]

is the initial homogeneous value of some variable $X$, ${\bf a}\equiv{\bf a}_{\rm ext,\,dust}-{\bf a}_{\rm ext,\,gas}$ is the difference between dust and gas accelerations, and $\tau\equiv t_{s}^{0}/t_{L}^{0}$ parameterizes the grains’ magnetization. The homogenous, steady-state equilibrium solution preserves this quasi-equilibrium while the whole box uniformly accelerates333As shown explicitly in Hopkins & Squire (2018b) (see § 2.1 and Appendix B therein) and Paper I, the dynamics of this local (un-stratified) problem are manifestly identical in the stationary and free-falling or uniformly-accelerating frames moving with the homogeneous solution, or (by extension) if we add an equal-and-opposite mean acceleration on the gas such that the mean acceleration of the entire system vanishes. The problem is trivially invariant to any uniform velocity boost. We therefore will perform all our analysis in the co-moving (free-falling) frame. We have also verified (for numerical testing) that our results are identical up to machine error (as they should be given our Lagrangian code) if we instead add an explicit uniform acceleration and/or boost to gas and dust to ensure the homogenous ${\bf u}_{g}(t)=\bf{0}$ in the lab frame. with ${\bf u}_{g}({\bf x},\,t)=\langle{\bf u}_{g}(t)\rangle={\bf a}_{\rm ext,\,gas}\,t+{\bf a}\,\mu\,t/(1+\mu)$.

We can make the coupled dust-gas equations dimensionless by working in units of the equilibrium sound speed $c_{s}^{0}$, gas density $\rho_{g}^{0}$, and “weighted grain size” $\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}/\rho_{g}^{0}$. Then, for a given equation-of-state, the dynamics of the problem (at infinite numerical resolution) are entirely determined by six dimensionless parameters: (1) the acceleration $\tilde{a}\equiv|{\bf a}|\,\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}/((c_{s}^{0})^{2}\,\rho_{g}^{0})$, (2) the box size or grain “size parameter” $\tilde{\alpha}\equiv\bar{\rho}_{\rm grain}^{\,i}\,\epsilon_{\rm grain}/\rho_{g}^{0}\,L_{\rm box}$, (3) the grain “charge parameter” $\tilde{\phi}\equiv 3\,Z_{\rm grain}^{0}\,e/(4\pi\,c\,\epsilon_{\rm grain}^{2}\,(\rho_{g}^{0})^{1/2})$, (4) the dust-to-gas ratio $\mu\equiv\rho_{d}^{0}/\rho_{g}^{0}$, (5) the plasma $\beta\equiv P_{0}/(|{\bf B}_{0}|^{2}/8\pi)$, and (6) the angle $|\cos{\theta_{{\bf B}{\bf a}}}|\equiv|\hat{\bf B}_{0}\cdot\hat{{\bf a}}|$ between the initial field direction $\hat{\bf B}_{0}$ and $\hat{{\bf a}}$. Note that in our linear theory perturbation analysis we chose to work with a different, but mathematically equivalent, set of dimensionless variables: (1) $|{\bf w}_{s}^{0}|/c_{s}^{0}$, (2) $L_{\rm box}/c_{s}^{0}\,t_{s}^{0}$, (3) $\tau$, (4) $\mu$, (5) $\beta$, (6) $\theta_{{\bf B}{\bf a}}$.

Our default simulations adopt $N_{\rm gas}=128^{3}$ gas resolution elements and an equal number of dust elements, $\mu=0.01$ (the ISM mean), and an isothermal ($\gamma=1$) equation of state (appropriate for most ISM/CGM/HII region conditions of interest). But we vary all of this below. For simplicity, we assume throughout that grains are all of the same size and charge, and that the grain charge $q_{\rm grain}$ is fixed during the simulation (as appropriate for large grains in isothermal gas) with grain Larmor time $t_{L}\equiv m_{\rm grain}\,c/|q_{\rm grain}\,{\bf B}|$. The “non-default” simulations discussed in § 4.2 relax some of these restrictions, exploring different equations of states, dust-to-gas ratios, resolutions, and allowing $q_{\rm grain}$ to vary with local gas parameters.

For convenience, throughout we adopt the Cartesian ($xyz$) axis convention with $\hat{\bf z}\propto{\bf B}_{0}$, $\hat{\bf x}\propto{\bf a}_{\bot}$ (i.e. the $x-z$ plane is defined to contain ${\bf a}$, so ${\bf a}={\bf a}_{\bot}+{\bf a}_{\|}={\rm a}_{\bot}\hat{\bf x}+{\rm a}_{\|}\hat{\bf z}$), and $\hat{\bf y}\propto{\bf B}_{0}\times{\bf a}$ (the mutually-perpendicular direction). In our 3D visualizations, the width/depth/height dimensions correspond to $x/y/z$.

3 Parameter Space Explored

Because the possible parameter space is enormous (the six dimensions above, plus the choice of equation-of-state, drag law, and charge law), we do not attempt to survey it systematically, but instead choose several unique parameter combinations motivated by values one might expect in different astrophysical systems. For more information, see § 9 of Paper I, which discusses each of these physical regimes extensively (Fig. 1 is adapted from this). The baseline parameters for our “default” simulation set are given in Table 1 and illustrated in Fig. 1. We give a brief description of each system in the following paragraphs.

The parameters HII-near and HII-far (from Fig. 7 in Paper I) correspond to plausible parameters in a massive HII region at two different radii from the star(s). Specifically, HII-near corresponds to parameters expected for HII regions with $\sim 0.1\,\mu{\rm m}$ grains at a distance $r\sim 0.1\,$pc from an OV-star or group or stars with luminosity $\sim 10^{6}\,L_{\sun}$ (providing the radiation pressure on the grains), local gas density $\sim 10^{4}\,{\rm cm^{-3}}$, temperature $\sim 10^{4}\,$K, plasma $\beta\sim 10$, $t_{s}$ calculated including both Epstein and Coulomb drag terms, and the grain charge calculated including photo-electric and collisional charging in that radiative environment (and accounting for saturation of the grain charge). HII-far takes the same system, and re-calculates all properties assuming a distance $r\sim 1\,$pc from the star (assuming the gas density falls $\propto r^{-2}$ and $\beta$ is constant). The important difference, for our purposes, is that the equilibrium drift velocity ${\bf w}_{s}$ is super-sonic in HII-near (where the radiation field is stronger), and sub-sonic in HII-far.

Likewise AGB is chosen to represent grains near the base (at $r\sim 100\,R_{\odot}$) of a dust-driven wind from a cool giant star, with a steady-state wind mass-loss rate $\dot{M}_{w}\sim 10^{-6}\,M_{\odot}\,{\rm yr^{-1}}$, wind velocity $v_{w}\sim 10\,{\rm km\,s^{-1}}$, temperature $T\sim 2000\,$K, $\beta\sim 1$, stellar luminosity $\sim 10^{3}\,L_{\odot}$ providing the grain acceleration, and $\sim 0.1\,\mu\,{\rm m}$ grains with a similar calculation of the charge and drag parameters. The distinguishing feature of this case is that the high gas densities ($\sim 10^{12}\,{\rm cm^{-3}}$) mean drag (collisional) grain coupling strongly dominates over Lorentz forces, so $\tau\ll 1$.444 As shown in Paper I, at low $\tau\ll 1$, the gyro RDIs are formally present, but are generally stabilized, except at either very large $k$ or very specific angles where they become degenerate with the MHD-wave RDIs. For this reason we do not plot them in Fig. 1 for the AGB case.

WIM represents a random patch of the diffuse warm interstellar medium, with $\beta\sim 1$, $T\sim 10^{4}\,$K, gas density $\sim 1\,{\rm cm^{-3}}$. We assume that the radiation energy density accelerating grains is comparable to the thermal energy density. Here $\tau\gg 1$ because the gas is diffuse, and the drift is highly sub-sonic.

Corona has parameters similar to those expected near the base of the solar Corona: $r\sim R_{\odot}$, $L\sim L_{\odot}$ (with both gravity and radiation contributing to the drift), $n\sim 10^{8}\,{\rm cm^{-3}}$, $\beta\sim 0.001$, $T\sim 10^{6}\,$K, for $\sim 0.1\,\mu{\rm m}$ grains. Here $\beta\ll 1$ makes this regime distinct; this also means the drift is super-sonic but sub-Alfvénic, and $\tau\gg 1$.

CGM represents parameters that could be present in the circum-galactic medium at $r\sim 100\,$kpc from a bright quasar with $L\sim 10^{13}\,L_{\odot}$, and $n\sim 10^{-2}\,{\rm cm^{-3}}$, $T\sim 10^{6}\,$K, $\beta\sim 1000$. The low density means $\tau\sim 10^{3}-10^{7}$ is large, while the high luminosity provides a super-Alfvénic equilibrium drift, producing a distinct mode structure.

Example is not chosen to match a particular system, but is the case examined in Seligman et al. (2018), including a resolution study. It lies between HII-near and WIM, so we include it here for comparison.

For each of these parameter sets, Fig. 2 shows the linear growth rates of the instabilities as a function of wavelength. For a given $k=|{\bf k}|$, these are calculated by choosing a mode angle $\hat{\bf k}$, which is either parallel to the drift or magnetic field, or satisfies one of the various resonant conditions at which a natural dust oscillation frequency matches one of the gas (e.g.  dust advection and Alfvén waves). For each class of resonance, we plot the maximum growth rate for each $\hat{\bf k}$. More information, including how to calculate the resonant angles, is given in Paper I.

Although all scales are unstable, different modes have growth rates that depend on the scale differently. In the physical systems explored above, the dynamic range between the largest global scales (where our local box treatment would be inappropriate) and the smallest scales where the instability operates is so large that it is impossible to resolve in a single box. Therefore, we construct a series of boxes for most cases here,555We only consider a single box-size in Example, as the results resemble HII-far, and CGM, as the mode structure does not change over a large range of wavelengths. each of which resolves a different range of wavelengths.

The parameters presented above are plausible, but will vary within and between different astrophysical regimes. Parameters that depend on highly uncertain grain chemistry or physical structure, such as $\tau$ (which depends on the grain charge), are particularly uncertain. Other quantities held fixed in our study – e.g. the gas equation-of-state, or the dust-to-gas ratio – might vary between regimes. For this reason we consider a number of variations in gas and dust physics compared to the “default” simulations. These are noted in Table 2 and discussed in detail in § 4.2.

The simulations presented here are scale free, and can represent any system that has the same dimensionless parameters. They are not strictly tied to the one specific physical system described above, but rather provide a well-motivated starting point for our study. For example, as discussed in Paper I, at different stages in the evolution of a supernova remnant (SNR), the SNR will pass through extended periods with parameters broadly similar to the AGB, Corona, and HII-near regimes above (although details will differ). Similarly, many regions of GMCs and the obscuring dusty torus around AGN will feature parameters resembling the AGB case. Some further intuition can be gained from Fig. 1 (although note that there are other important parameters that cannot be easily shown on a two-dimensional plot)

We broadly classify the behavior of the saturated simulations into two regimes. The first is the “clumped” regime, where the dust becomes strongly concentrated; this occurs, for example, in HII-near (Fig. 4), HII-far (Fig. 5), Example (Fig. 16), and AGB (Fig. 6). The second is the “disperse” regime, where the dust is expelled from certain regions at high velocity but remains relatively homogeneous; this occurs in (some scales in) WIM (Fig. 7), CGM (Fig. 13) and Corona (Fig. 8).

4 Results

4.1 Default Simulation Set

We evolve each box well into its non-linear growth phase. As a first example, Fig. 3 shows visualizations666The gas+dust visualization is constructed by interpolating the continuous gas properties onto the $x/y/z$ axis surfaces, and taking all dust “super-particles” that lie within a thin slice of width equal to roughly the median inter-particle separation and projecting them onto the surface. The dust-density visualization is constructed by first calculating the dust density in the vicinity of each dust particle element using a kernel density estimator as in Moseley et al. (2018), then plotting each dust particle in the projected 3D space, color-coded by the density (with a constant transparency). The range on the color-bars is scaled to include some fraction (typically $>90\%$) of the plotted elements, in order to show contrast, but there are always some elements with higher/lower values. of the dust and gas in the large HII-near simulation box, at a time well into the non-linear phase of the evolution. Large fluctuations and coherent structure are visible in all gas quantities and in the dust.

Figures 4, 5, 6, 7, and 8 show the morphologies of gas and dust in the default HII-near, HII-far, AGB, WIM, and Corona sets, respectively. We show gas and dust properties during two different simulation times, corresponding to the early non-linear and saturated phases of evolution. For each parameter set, we show the boxes of different sizes in parallel columns. Sequential simulations from left-to-right are $\sim 300\times$ larger in size scale; an entire box at the left is approximately $\sim 1/2$ of a pixel/element in the box to its immediate right. We only show the gas quantities that exhibit the most obvious morphological structure. Fig. 11, 13, 12, and 19 show the same for the Example, CGM and HII-near sets. Parallel columns now compare compare physical variations (e.g. changing the equation of state of the gas, or dust charge scaling). These are discussed further in § 11.

Figure 9 show the magnitude of fluctuations in dust velocity versus time for each simulation. All dust fluctuations are plotted in units of the expected linear-theory growth timescales, which allows us to compare all of the simulations on the same time axis. Since there is no single growth time in a simulation (see Fig. 2), we consider the maximum growth rate of modes at the box scale.777Specifically we find the minimum “box scale growth time” $t_{\rm grow}(L_{\rm box})\equiv 1/\Im{(\omega_{\rm max})}$ where $\omega_{\rm max}$ corresponds to the mode with $|{\bf k}|=2\pi/L_{\rm box}$ that has the largest positive value of $\Im(\omega)$ (marginalizing over direction $\hat{\bf k}$). In most cases, modes near the resolution scale ($L_{\rm box}/N$) grow $\sim 10\times$ faster than those at the box scale,888This estimate is not always accurate and depends on details of the dispersion relation. E.g. in some cases, a different mode starts to dominate at some mid-range scale. and so the initial rapid growth is dominated by these small-scale modes. However, these modes generally saturate at a lower amplitude, so the later growth is then dominated by the box-scale modes, at rates in agreement with linear theory. This is surprising given the ambiguity of defining the total growth rate, and the fact that small-scale modes have already become non-linear. All of the boxes eventually reach saturation, with the fluctuations in all properties in a statistical steady-state. The cases in Fig. 9 that continue off the plot have been evolved to longer times to ensure that they are approximately saturated.999We have verified that all simulation quantities within a given simulation saturate on a similar timescale, as shown in Fig. 9. Explicit demonstrations of this, as well as plots of the full PDFs of the salient quantities vs. time, for representative simulations, are shown in both Moseley et al. (2018) and Seligman et al. (2018).

In Table 2, we provide the magnitude of the time-averaged saturated fluctuations101010We define the fluctuation $\delta X$ in quantity $X$ as the rms ($1\sigma$) deviation, i.e. $\delta X\equiv\langle(X-\langle X\rangle)^{2}\rangle^{1/2}$, where $\langle U\rangle\equiv(\int U\,\varpi\,d^{3}{\bf x})/(\int\varpi\,d^{3}{\bf x})$ and $\varpi({\bf x})$ is a weight. For ${\bf u}_{g}$, ${\bf v}_{d}$, and ${\bf B}$, it is most physical to relate the fluctuations to the energy in each component. So, for ${\bf u}_{g}$ we use the gas-mass-weighted average ($\varpi=\rho_{g}$, such that the total kinetic energy of fluctuations is $(1/2)\,M_{\rm gas}^{\rm box}\,|\delta{\bf u}_{g}|^{2}$), for ${\bf v}_{d}$ we likewise use the dust-mass-weighted average ($\varpi=\rho_{d}$, so dust kinetic energy is $(\mu/2)\,M_{\rm gas}^{\rm box}\,|\delta{\bf v}_{d}|^{2}$), and for ${\bf B}$ we use the volume average ($\varpi=1$, so magnetic energy is $V_{\rm box}\,|\delta{\bf B}|^{2}/8\pi$). For $\rho_{g}$, $\rho_{d}$, and $\rho_{d}/\rho_{g}$ we quote the volume-averaged fluctuations ($\varpi=1$). The differences between mass and volume weighting are discussed further below. These are all defined in the center-of-mass (i.e. co-moving or free-accelerating) frame. in each component of ${\bf u}_{g}$, ${\bf v}_{d}$, ${\bf B}$, $\rho_{g}$, $\rho_{d}$. Figure 10 plots these statistically averaged quantities against one another in various forms, and against some expectations from quasi-linear theory for the saturated state. This is discussed further in § 5 below.

Figures 11, 13, 14, 12, 15, 19, 20, consider further comparisons of different physics and parameter variations, as discussed below.

Figures 21, 22, 23, 24, 25, 26, 27 examine various statistics of each run in more detail, plotting the probability distribution functions (PDFs) of fluctuations of different quantities in the saturated state.111111Given our Lagrangian numerical method, at our default simulation resolution, dust under-densities $\rho_{d}/\rho_{d}^{0}\ll 10^{-6}$ cannot be resolved, but these are much smaller than any plotted in Figures 21-27 and much smaller than the typical fluctuations in Table 2. There is no formal upper limit to the maximum resolved concentration, but the dust becomes increasingly over-resolved relative to the gas at very large $\rho_{d}/\rho_{g}\gg 10^{6}$ – however this is precisely the regime where we expect other physics to dominate, discussed below.

4.2 Different Equations-of-State, Drag Laws, and Grain Charge Scalings

In this subsection, we discuss the effects of additional physical variations in the simulations. We discuss the physical applicability of different drag laws in different regimes (expanding the discussion from Paper I), and explore different choices. We outline three idealized, physically motivated regimes of the grain charge scaling, and compare the effects of each of these in the simulations. We justify our usage of an isothermal equation of state for the majority of our simulations, and discuss astrophysical situations where different equations of state are appropriate. We also investigate how changes to the magnetic field, drift geometries, and additional grain parameters alters the non-linear evolution of the instabilities.

4.2.1 Drag law

The scalings and physical applicability of different drag laws are discussed in detail in Paper I. There it is shown that for any regime where the MHD RDIs are important, the drag is dominated by either Epstein or Coulomb drag (e.g. Stokes drag, important only when the grain is larger than the gas mean free path and Reynolds numbers are less than unity, is never relevant for the physical grain sizes in magnetized media considered here). Our default simulations adopt the Epstein drag scalings, which, for fixed grain properties, give

[TABLE]

($a_{E}=9\pi\,\gamma/128$). Likewise, Coulomb drag gives

[TABLE]

where $a_{C}=(2\gamma^{3}/9\pi)^{3}$ and $f(T,U)=f_{\rm ion}\,(z_{i}\,U/T)^{2}\ln{\Lambda}$ is a function of the temperature $T$, ionization fraction $f_{\rm ion}$, grain electrostatic potential $U$, gas ion charge $z_{i}$, and a Coulomb logarithm $\ln{\Lambda}$.

It is difficult to evaluate Coulomb drag without specifying the physical system (to calculate $f_{\rm ion}$, $T$, $\Lambda$, etc.), whereas Epstein drag is determined by scale-free hydrodynamic quantities. Epstein drag dominates at high $w_{s}/c_{s}$, scaling as $\sim w_{s}/c_{s}$, while Coulomb drag becomes weaker, scaling as $(w_{s}/c_{s})^{-3}$. Therefore, Epstein drag always dominates in super-sonic cases (where $w_{s}/c_{s}\gtrsim 1$ or $|\delta{\bf v}_{d}|/c_{s}\gtrsim 1$). Moreover, Epstein drag often dominates in sub-sonic cases if the ionization fraction is small or if the grains reach charge saturation (see Paper I). We therefore expect Epstein drag to be a good approximation in many cases. Even in sub-sonic regimes where Coulomb drag may dominate, if the temperature and grain potential $U$ are fixed (as they are in our default runs), then Epstein and Coulomb drag scale identically (as $t_{s}^{-1}\propto\rho$). Because the normalization of $t_{s}$ is arbitrary (it simply goes into the dimensionless parameter $c_{s}\,t_{s}/L_{\rm box}$), our predictions therefore apply equally well to either drag law in this sub-sonic regime, provided one re-scales the simulations dimensionless parameters appropriately.

4.2.2 Grain charge

Scalings of grain charge remain uncertain and depend on a variety of environmental and local factors (see Paper I and, e.g. Draine & Sutin, 1987; Weingartner & Draine, 2001b, a; Tielens, 2005). However, if we assume that external environmental properties are fixed (e.g. a UV radiation field) and grain material properties are fixed (e.g. compositions, sizes), then three regimes emerge, which we define as:

[TABLE]

(a) “Fixed Charge”: if the grain charge is saturated (maximal/minimal), or the charge-aggregation/equilibration times are long compared to the timescales we evolve, then the charge is approximately fixed. (b) “Collisional Charging”: if the charging is dominated by collisional processes, is sufficiently fast, and is un-saturated, it depends solely on the temperature as $Z_{\rm grain}\propto T$. (c) “Photo-Electric Charging”: if the charging is dominated by photo-electric processes (and again is also un-saturated and sufficiently rapid), then $Z_{\rm grain}\propto T^{1/2}\,\rho_{g}^{-1}$. Our default simulations assume regime (a), i.e. fixed charge, which is common for large grains in well-ionized environments.121212We assume throughout that we have large grains with $|\langle Z_{\rm grain}\rangle|\gg 1$, so charge can be treated as a continuous, smooth quantity. In an isothermal gas, our default case, regime (b), collisional charging, is identical (because $Z_{\rm grain}$ depends only on $T$ for specified grain properties). Regime (c), photo-electric charging, is only different in the isothermal case if the gas is highly-compressible so that fluctuations in $\rho$ are large.

In order to explore the non-linear effects of Regime (c), we re-run HII-near L, but adopt the photo-electric scaling $t_{L}^{-1}\propto Z_{\rm grain}\propto T^{1/2}\,\rho_{g}^{-1}$ with the same initial homogeneous value of $Z_{\rm grain}$ such that the constants in Table 1 are identical (i.e. only the scaling is changed). We also compare a run with $t_{L}^{-1}\propto T$ corresponding to collisional charging, adopting $\gamma=5/3$ so that $T$ will actually vary. These HII-near L comparisons are shown in Fig. 12. It is evident that the photo-electric charging law produces dust structures with distinct topologies. The dust is more coherent, and the gas is dragged along into regions with higher dust densities and higher velocities.

Given the (significant) physical uncertainties in predicting grain charge, we also explore the effects of re-scaling the assumed charge from our default assumptions. The results of this is shown in the left column of Fig. 12 for HII-near, the right column of Fig. 13 for CGM, and the right column of Fig. 14 for Corona. We find that increasing the charge (holding all else equal) produces faster-growing, more violent instabilities.

4.2.3 Gas equation of state

Our default simulations adopt an isothermal equation of state ($\gamma=1$). This is a good approximation to many astrophysical cases of interest (e.g. the WIM, HII regions, the Solar corona), where cooling is efficient. Moreover, the majority of our runs are only weakly-compressible, so changing $\gamma$ (with all other parameters fixed) only has a weak effect. In Figure 11, we present a comparison of Example with $\gamma=5/3$ (with the initial pressure and all other parameters identical to the “default” run). It is apparent that changing the equation of state in the weakly compressible limit has very little effect on the non-linear evolution of the instability. This near-independence of the instabilities from changes to the equation of state also suggests that different charge scalings or drag scalings (i.e. a different dependence on $T$, $\rho$, etc., discussed above) will produce similar results to our default assumption (constant charge + Epstein drag), because gas properties do not vary dramatically regardless of $\gamma$. On smaller scales than examined here, different equation of states are more appropriate because the growth time of the instabilities may become short compared to the cooling time.

The few of our runs that are more strongly compressible in the gas all have $|\delta{\bf v}_{d}|/c_{s}\gg 1$ (as they must, since compressibility requires $|\delta{\bf u}_{g}|/c_{s}\gg 1$ and typically $|\delta{\bf u}_{g}|\lesssim|\delta{\bf v}_{d}|$). This implies that Epstein drag should better approximate dust drag forces (see § 4.2.1). It is difficult to identify any parameter regime where Coulomb drag both dominates over Epstein drag, and has the potential to significantly alter our findings (due to its temperature or charge dependence in a non-isothermal EOS). Of our compressible simulations, we note that most represent parameters motivated by the very regions we expect to have $\gamma\approx 1$ on the scales of interest. The only case where one might expect $\gamma\neq 1$ among this set is run CGM, representative of the CGM around a bright source (where cooling is inefficient so $\gamma\approx 5/3$ is more appropriate). We therefore re-run this simulation with $\gamma=5/3$ and Epstein drag (see Fig. 13). In this case the physically realistic charge scaling is not obvious, although it is probably saturated (see Paper I). For the sake of exploring different regimes we assume un-saturated photo-electric charging, so $t_{L}^{-1}\propto Z\propto T^{1/2}/\rho\propto\rho^{(\gamma-3)/2}$. The resulting differences with the default parameters are quite small (cf. left and middle panels in Fig. 13)

4.2.4 Other parameter variations: Magnetic field orientation and strength, Drift geometry and speed, and Dust-to-gas ratios

Aside from these variations in drag-law, gas properties, and grain charge, we consider variations of other parameters compared to the default simulations in Figs. 11-19.

We examine how changes to the magnetic and drift geometries and magnitudes effect the evolution. We decrease $\beta$ by a factor of $2$ (from $\beta=2$ to $\beta=1$) in Example, keeping all other parameters fixed. This has only minor effects on the qualitative saturation behavior, but it does shift the resonant angles, effectively rotating the box (see the right column of Fig. 11). We compare WIM S/M/L with very low drift velocity $|{\bf w}_{s}|/c_{s}\sim 0.01$ ($5\times$ lower than our default case). This makes the growth rates of some modes prohibitively long, requiring too long an evolution to capture; nonetheless, we see many similar qualitative features emerge as in the default case, but with weaker saturated amplitudes (see Fig. 14). In Fig. 15, we re-run HII-near L varying the angle between ${\bf B}$ and ${\bf a}$. This does not qualitatively change the behavior, but does shift the angle of the streams of dust in the saturated state (which tend to align with ${\bf a}$), and produces slightly faster growth rates as ${\bf a}$ becomes more orthogonal to ${\bf B}$ (this is predicted by linear theory, but is not obvious, as ${\bf w}_{s}$ actually becomes smaller).

We also investigate how changes to the dust properties effect the evolution. As shown in Seligman et al. (2018), removing dust charge in our Example box gives a linear growth rate approximately 2-3 orders of magnitude lower. We show the result of simulating this at similar times to our charged cases in Fig. 19, to show that no instability has grown. This also makes it clear that our numerical scheme does not introduce any artificial dust clumping. Finally, we simulate Example and HII-near with lower and higher dust-to-gas ratios ($\mu=0.001$ and $\mu=0.1$), which are shown in Figs. 16 and 18). The lower-$\mu$ cases are the most interesting: as expected, the growth rate decreases by a factor $\sim 3$ (see Fig. 17) and the gas density and turbulent velocity fluctuations decrease since the forcing from the dust on the gas is weaker. Surprisingly however, lowering $\mu$ appears to enhance the strength of the saturated dust density fluctuations. This may occur because the weaker gas turbulence is less efficient at disrupting dense dust structures.

4.3 Additional Numerical Tests

In addition to the generic code validation tests described in § 2, previous papers have tested the numerical methods here (Carballido et al., 2008; Johansen et al., 2009; Bai & Stone, 2010b; Pan et al., 2011; Hopkins & Raives, 2016; Hopkins, 2016b, 2017; Su et al., 2017; Hopkins & Lee, 2016; Lee et al., 2017; Moseley et al., 2018; Seligman et al., 2018). For example Hopkins & Lee (2016) show that the “finite-sampling” effects in super-particle methods, which introduce some shot noise in the particle densities and divergence between particle trajectories and gas (at the integration error level) in the perfectly-coupled limit (see Price & Federrath, 2010; Ayliffe et al., 2012; Genel et al., 2013; Tricco et al., 2017), introduce $\sim 0.01-0.05\,$dex scatter in the dust density in a super-sonically turbulent medium ($\delta{\bf u}_{g}/c_{s}\gg 1$). In both Hopkins & Lee (2016) and Moseley et al. (2018) we show in a variety of tests that even the “worst-case” scatter of this nature is completely negligible for our predictions,131313We note that Tricco et al. (2017) contains an error in their comparison of “one-fluid” dust treatments (which cannot represent multi-valued dust velocity distribution functions, essential for the present study; Cadiou et al. 2019; Lebreuilly et al. 2019) to the “tracer particle” ($\mu=0$) method used in Hopkins & Lee (2016). The discussion in Tricco et al. (2017) incorrectly compares to simulations in Hopkins & Lee (2016) which adopted different dust parameters ($\tilde{\alpha}$), a different dust drag law, and did not use the same weighting for the PDFs/dispersions measured. If one correctly compares the simulations from Hopkins & Lee (2016) which have similar $\tilde{\alpha}=0.001-0.01$ to the cases studied in Tricco et al. (2017), and weights the PDFs identically, then the agreement is reasonably good, and in fact the simulations in Hopkins & Lee (2016) predict slightly smaller fluctuations in the dust-to-gas ratio $\rho_{d}/\rho_{g}$ than those in Tricco et al. (2017), the opposite of their conclusion. and that in sub-sonically turbulent (or laminar) media this drops to $\ll 0.01\,$dex scatter. In either case this is much smaller than the saturated dispersions in $\rho_{d}$ here.

As in Moseley et al. (2018) and Seligman et al. (2018), we have re-run variants of several of our standard simulations with varied numerical choices, including: (1) a different hydrodynamic solver (for Example, WIM L, HII-far S/M/L), we use the “meshless finite mass” (MFM) method from Hopkins (2015); (2) a different scheme for calculating gradients and reconstruction of the magnetic field at cell faces (for Example, Corona S), specifically the “constrained gradient” MHD scheme in Hopkins (2016b); (3) using a naive (non-manifestly energy-conserving) explicit leapfrog integrator instead of the Boris integrator for the Lorentz forces on grains (in Example, HII-near S/M/L Corona S/M/L, WIM M/L, CGM); and (4) using different initial conditions, namely glass-like instead of lattice initial particle configurations (for Example). Choices (2), and (4) have no significant effect on any results we explore. We find that choice (1), the MFM method, moderately suppresses the initial maximum growth rates of the modes in the box, owing to its slightly less-accurate integration of highly-subsonic flows near the grid scale (effectively, the fastest-growing resolved modes in the simulation correspond to those at $\sim 4-6$ times the grid scale, instead of $\sim 2-3$ times the grid scale). Choice (3) has negligible effects on boxes Example and HII-near S/M/L, where the dust is strongly clumped and $\tau$ is not extremely large. However for boxes Corona, WIM, and CGM, where $\tau\geq 100$ is larger and the systems often end up in the “disperse” mode, we find that it becomes important to use an integrator (like our default Boris scheme) that manifestly conserves Lorentz orbits. Otherwise, the dust tends to drift and spiral outwards and gain energy artificially, an effect that is well-known in standard particle-in-cell codes.

We have also run a number of resolution tests. Our default runs here use $2\times 128^{3}$ elements ($128^{3}$ each in gas and dust). All of our runs here have also been run at lower resolution ($2\times 32^{3}$ and $2\times 64^{3}$), and one at higher ($2\times 256^{3}$). Figure 20 compares the saturated state in this resolution series. The qualitative behavior is similar across this range of resolutions. However, without some sort of physical isotropic viscosity, or non-fluid effects (e.g. the ion gyro radii), these instabilities are linearly unstable at all wavelengths, with growth rates that increase with $k$. Thus, increasing resolution and keeping all else fixed will always resolve new, faster-growing modes at smaller scales, so it is not possible to undertake a formal convergence analysis. A more formal resolution study, as well as a number of additional numerical validation tests and more extensive discussion, are presented in Moseley et al. (2018).141414Note that in Moseley et al. (2018), where we studied the pure acoustic RDI, we noted that the high-$k$ modes (with $|{\bf w}_{s}|\gg c_{s}$) were difficult to resolve in some cases because the growth rates become sharply-peaked around a narrow resonant angle of width $\delta\theta\propto\mu^{2/3}\,k^{-1/3}\,|{\bf w}_{s}|^{-1}$. As shown in Paper I, the MHD modes have a much more complex but also broader resonant structure, especially at the lower ${\bf w}_{s}$ and higher $\tau$ values studied here. This contributes to making our results significantly less sensitive to resolution. For example, Fig. B3 in Moseley et al. (2018) compares the PDFs of the same properties studied here (e.g. dust-to-gas ratio) as a function of the relative number of dust and gas elements (varying $N_{\rm dust}/N_{\rm gas}$ by a factor of $64$) – as expected, the resulting PDFs are essentially identical up e.g. Poisson sampling in the tails of the distribution. We have verified this with our runs here as well. As shown in Hopkins & Lee (2016); Lee et al. (2017); Moseley et al. (2018), these PDFs are also robust to the details of kernel density estimators for dust or gas.

The saturated amplitudes of the velocity and density fluctuations are relatively insensitive to resolution in run Example. This suggests that most of the power is dominated by the well-resolved largest modes in the box, rather than the fastest-growing but smallest-scale modes. This is consistent with our power-spectrum analysis in Seligman et al. (2018), where we showed explicitly in run Example that most of the power in the saturated state was in modes a factor $\sim 2$ around the box scale. This may not be true in some other runs, e.g. HII-far L, where the dust clumps on small scales.

4.4 Distribution Functions

Figures 21-27 examine the probability distribution functions (PDFs) of fluctuations for the seven simulations considered in this paper. We see that the large diversity of behavior in morphology, saturated amplitude, and anisotropy is mirrored in the diverse array of PDF shapes. The PDFs are often highly non-Gaussian; this comes (broadly speaking) in two “forms.”

First, some PDFs exhibit substantial sub-structure, e.g. AGB-M or HII-far in $\delta{\bf u}_{g}$ or $\delta{\bf B}$, which have multiple bumps (sharp inflection points/changes in curvature) or even second peaks. This is directly related to coherent, large-scale morphological structures in Figs. 4-7. Each such “peak” corresponds to “patches” (sub-volumes of the parent box) which have either much larger or smaller dust density within them, and are have essentially de-coupled from one another, evolving non-linearly almost independently.

Second, most of the “smooth” PDFs have non-Gaussian tails, which are most commonly exponential or “stretched exponential” in form: i.e.

[TABLE]

for some $q$ and $\gamma$. For example, the $\delta{\bf B}$ PDFs in almost all cases have this form, as do the $\delta{\bf u}_{g}$ PDFs in Corona and WIM S/M and HII-far S/M. PDFs with exponential or stretched-exponential tails are common in certain types of gas turbulence, velocity distributions of granular gases, and passive scalar concentrations in sub-sonic incompressible turbulence (Ruiz-Chavarria et al., 1996; Yakhot, 1997; Ben-Naim & Krapivsky, 2000; Antal et al., 2002; Ernst & Brito, 2002; Kohlstedt et al., 2005; Aranson & Tsimring, 2006; Monchaux et al., 2010; Hopkins, 2013b; Colbrook et al., 2017). This generically arises from a competition between driving and dissipation. Consider the distribution of velocities in a statistically homogeneous system (so $P({\bf u}_{g},\,{\bf x})=P({\bf u}_{g})$), where the velocities are driven by an uncorrelated stochastic process with the specific energy injection rate $\sim\dot{u}_{0}\,(v/v_{\ast})^{\xi}$ (the effective “diffusion coefficient” in velocity space), and damped151515Note this damping corresponds to any process damping the fluctuations, not necessarily the dust drag. with rate $t_{\rm damp}^{-1}\,(v/v_{\ast})^{\psi}$. In steady-state, if these driving and damping processes dominate, the PDF obeys

[TABLE]

for each component, the solutions of which obey $P(v)\propto\exp{(-|v/v_{\ast}|^{\gamma})}$ at large $v$ with161616For $\psi-\xi<-2$, the solutions become peaked functions $P(v)$ which asymptote to constant $P>0$ as $v\rightarrow\infty$.

[TABLE]

So, for the simple case of white-noise (Brownian) driving and constant damping, or any case where driving and damping depend on $v$ in the same manner ($\psi\approx\xi$), the tails are exponential. The characteristic width of the PDF is simply given by $v_{\ast}^{2}\sim\dot{u}_{0}\,t_{0}$; i.e. the energy injected in a damping/dissipation time.

While it is easy to qualitatively understand the range of PDF tails in this manner, such arguments fall far short of a predictive model. In other words, it is not possible (given the arguments here alone) to predict the PDFs and structure functions purely from the various simulation parameters. For example, it is not a priori obvious what the relevant driving and damping rates should be in the saturated regime. If turbulence dominates both the non-linear forcing and the damping (e.g. eddies shearing apart growing modes, as we argue sets the saturation amplitude of some modes below) then both injection and damping times might scale with eddy turnover times, as in e.g. the standard theory of supersonic turbulent density fluctuations (Vazquez-Semadeni, 1994; Hopkins, 2013a, b; Squire & Hopkins, 2017). But driving could also arise directly from mode growth, or from non-linear parasitic modes, while damping could also stem from drag (acting on dust velocity fluctuations) or sound waves (for gas pressure fluctuations) – the dominant terms do not have to be the same for each type of fluctuation.

Some of the PDFs exhibit strong asymmetries, with a much stronger “tail” in one direction (e.g. ${\bf u}_{g}$ and ${\bf v}_{d}$ in AGB L/XL, $\rho_{d}$ in HII-far M/L or CGM). In cases like AGB or HII-far where the tails extend to both larger ${\bf u}_{g}$ in the direction of acceleration ${\bf a}$ and larger $\rho_{d}$, this relates to the fact that regions with non-linearly larger local $\rho_{d}/\rho_{g}$ experience faster local growth of the instabilities, and more efficient acceleration of the coupled dust-gas mixture (because the gas acceleration scales as $\sim\mu\,{\bf a}$). Cases with a large tail in $\rho_{d}$ towards smaller $\rho_{d}$ (CGM, WIM) likely arise because dust is locally expelled from small pockets that are local vorticity maxima. This is essentially the generic mechanism studied in Hopkins (2016a) and is well-known in terrestrial particulate “preferential concentration” studies (Squires & Eaton, 1991; Fessler et al., 1994; Rouson & Eaton, 2001; Gualtieri et al., 2009; Monchaux et al., 2010). Essentially, grains are centrifugally “flung out” of high-vorticity regions, and concentration is a side effect as grains collect in-between. Finally, there are examples that are highly non-Gaussian but do not neatly fit into any of the classes above (e.g. ${\bf v}_{d}$ in CGM or WIM L).

We also compare the volume-weighted ($P_{V}$, the probability that a given random volume has some value) and mass-weighted (dust weighted by dust mass, gas by gas mass; $P_{M}$) statistics. If the density variations are small, these must be similar. Indeed for the gas velocity and magnetic field (${\bf u}_{g}$ and ${\bf B}$) the two statistics rarely differ dramatically, as the gas density fluctuations are small. The few cases (e.g. HII-near L) where some differences appear are those with the most dramatic gas density fluctuations, and they still do not differ qualitatively. Even for the dust, where $\rho_{d}$ varies strongly, the differences are rarely dramatic in ${\bf v}_{d}$ (an exception is HII-near L, where dust in the high-$\rho_{d}$ filaments is coherently moving relative to the low-$\rho_{d}$ “background,” giving rise to larger dispersions when volume-weighted).

In the density PDFs, volume and mass weighting makes more difference. Some of this is by definition: because $\rho=d{\rm Mass}/d{\rm Volume}$,

[TABLE]

if both are measured in differentially small regions. So $P_{M}(\rho)$ will always be biased to higher $\rho$ than $P_{V}(\rho)$. Let us define the mass and volume-weighted averages $\langle(\cdot)\rangle_{M}$ and $\langle(\cdot)\rangle_{V}$ respectively, the volume-weighted variance $S$ in the log-density field $\chi\equiv\ln{\rho}$, and the usual “clumping factor” $C$:

[TABLE]

where $\rho_{0}=\langle\rho\rangle_{V}$ follows from our definitions. If $P_{V}(\rho)$ is exactly log-normal, then mass-conservation implies: (1) the volume-weighted mean/median/mode of $\chi$ is $\langle\ln(\rho/\rho_{0})\rangle_{V}=-S_{V}/2$, (2) $P_{M}$ is also lognormal with $S_{M}=S_{V}$ and mass-weighted $\langle\ln(\rho/\rho_{0})\rangle_{M}=+S_{V}/2$, (3) $C=\exp{(S_{V})}$.

One way that we can quantify the deviations from log-normal in these PDFs is via the extent to which $S_{M}\neq S_{V}$ or $C\neq\exp{(S_{V})}$ (this is commonly used as a diagnostic in supersonic turbulence studies). This can be read off from the values of $S_{V}=(\delta\ln{[\rho/\rho_{0}]})^{2}$ and $C$ in Table 2. In several cases (e.g. AGB) the results are consistent with log-normality ($C\approx\exp{(S_{V})}$). But in some very strongly-clumped cases we have, for the dust, $S_{M}>S_{V}$ and $C\gg\exp{(S_{V})}$ (e.g. HII-near and HII-far M,L have $S_{V}\sim 2$ and $S_{M}\sim 4-6$). Specifically, given the measured $\delta\ln{(\rho_{d}/\langle\rho_{d}\rangle)}=S_{V}^{1/2}$ (for the dust), and assuming log-normality ($C=\exp{(S_{V})}$) would predict $C\sim 5-10$ in these cases. In contrast the measured $C$ is $C\sim 100-300$ in HII-near L and HII-far L, and $C\sim 1500$ in HII-far M. Conversely, the PDFs with tails towards low $\rho_{d}$ have smaller $C$ than would be predicted from $\delta\ln{(\rho_{d}/\langle\rho_{d}\rangle)}$; e.g. the CGM and WIM runs have $C\sim 1.0-1.1$. This owes directly to the large, asymmetric tails in $\rho_{d}$ visible in Figs. 21-27.

Another, perhaps simpler, way of emphasizing this is to note that the most extreme cases (e.g. HII-far M, or HII-near L and Example with lower $\mu=0.001$) have extreme tails, where $\gtrsim 0.1\%$ of the dust mass lies at densities $\gtrsim 10^{7}$ times larger than the mean (reaching as high as $\sim 10^{9}-10^{12}$). This corresponds to values $\gtrsim 15$ standard deviations from the volume-weighted median of the density PDF.

5 Discussion

From the analysis in § 4, we can identify a number of important conclusions about both the linear and early non-linear phases, as well as the saturated and late non-linear phases.

5.1 Linear and Early Growth Phases

In all cases, the instabilities grow and produce strongly non-linear properties in the dust (see Figs. 4-13). This is not surprising, since all cases here are linearly unstable (Fig. 2). However, it does show that non-linear growth can occur even when the growth timescale is shorter than the dust “stopping time” or Larmor time (i.e. $\Im{(\omega\,t_{s})}>1$ or $\Im{(\omega\,t_{s})}>\tau$, in Fig. 2). This implies that the instabilities grow faster than the timescale for the system to reach a “new” equilibria, which is not a trivial conclusion. 2. 2.

Generally, the instabilities exhibit the most rapid initial growth at the smallest scales. All wavelengths here are linearly unstable, with growth rates generally increasing at smaller wavelengths (Fig. 2), because we do not include explicit dissipation (e.g. viscosity).171717In fact (see Paper I) some of instabilities (e.g. the Alfvén wave and gyro resonances) have growth rates that can continue to rise with decreasing wavelength even below the field-parallel viscous scale, potentially down to the ion gyro radii or even further. Thus, within each parameter set, the smaller boxes evolve more rapidly in fixed physical units (e.g. $t_{s}^{0}$). Moreover within each box, in most cases we see small-scale modes saturate first, with the scale of structures increasing until box-scale modes saturate (compare e.g. early and later times in Fig. 4, 5, 11). 3. 3.

The simulations broadly reproduce the expected growth rates from linear theory in the early linear growth phase, as depicted in Fig. 9. We demonstrate and discuss this in more detail for a few representative acoustic RDIs in Moseley et al. (2018) (Fig. 2 therein) and our Example case in Seligman et al. (2018) (Figs. 4-5). For modes that are unstable on all scales (e.g. the Alfvén MHD-wave RDI), the early growth rate appears to match the linear theory prediction for wavenumbers $k\sim(1-3)/\Delta x_{0}$ (where $\Delta x_{0}=L_{\rm box}/N_{\rm 1D}$ is the initial inter-element grid spacing; again see also Moseley et al. 2018, Fig. 2). This suggests that the growth rate can be recovered even if the full wavelengths are resolved by just $\sim 2-6$ inter-particle spacings.181818For rigorous demonstration of this in idealized test problems, as well as formal convergence studies, see Moseley et al. (2018), Appendices A-B. This is also consistent with the results using different codes and numerical methods, e.g. Johansen & Youdin (2007). At later times, the box-scale modes take over as the strongest growing modes, until saturation. 4. 4.

The presence of magnetic fields and grain charge, along with associated Lorentz forces, play a critical role in the linear development of the instability. For several cases studied here, if we artificially remove the dust charge (Lorentz forces), or magnetic fields, the instabilities become stable or orders-of-magnitude more slowly-growing. Fig. 19 shows that no structure develops at late times in the Example run with zero grain charge. Moreover, increasing the grain charge ($q_{\rm grain}$ and $\tau$) produces faster-growing instabilities with more violent saturated dust behavior, in both the “clumped” and “dispersed” regimes (see Figs. 12, and 13). This is despite the fact that Lorentz forces decrease the “equilibrium” dust drift velocities, and magnetic fields increase the pressure support of the gas, which would naively appear to be “stabilizing” effects. But these terms also introduce a variety of new dust and gas modes (e.g. Alfvén and slow waves, dust gyro motion) which in turn dramatically increase the number of accessible “resonances” for the instabilities (as well as introducing new energy sources for the instabilities). 5. 5.

The “resonances” where linear growth rates are maximized are sufficiently vigorous that they can often be identified well into the non-linear evolution. These maximal resonances occur where the “natural frequency” of advection ${\bf w}_{s}\cdot{\bf k}$ or gyro motion $t_{L}^{-1}$ in the dust matches the “natural frequency” of an Alfvén or magnetosonic wave in the gas. Even in the non-linear phases of evolution, these resonances manifest as particular angles or wavelengths of the structures that form, as can be seen in Figs. 16-17.

These results are all, to some extent, predicted from the linear theory in Paper I, but we both verify the linear theory results and confirm that these conclusions persist even well into non-linear evolution.

5.2 Saturation & Late Non-Linear Phases: Generic Conclusions

As noted above, the saturated states exhibit some qualitatively different behaviors, but there are some generic conclusions that apply to all of our runs. We discuss these first, before exploring the physics that is distinct in different saturated states.

All of the systems saturate in a turbulent quasi-steady state. This is evident in Fig. 9, where the velocity dispersion has clearly reached saturation, although in some cases certain fluctuations continue to grow very slowly. The finite velocity dispersion of the dust grains in the saturated state does not shut down fluctuations, even in the cases that reach a nearly isotropic grain velocity distribution function (e.g. Fig. 27). In other words, the turbulence reaches a saturated steady state, rather than the instability quenching itself. 2. 2.

The predicted structure here is qualitatively completely different from that formed in “passive” dust experiments, in which the dust moves as a tracer particle in externally driven turbulence (neglecting the forces from dust on gas, i.e. the momentum-conserving “back-reaction” terms). In other words, “passive” dust cannot generate these instabilities or structures – indeed, the RDIs do not exist with “passive” dust. Some of the most noticeable differences are found in the large-scale dust morphology; there is a much stronger prevalence of filaments and sheets here, compared to no strong anisotropy present in “passive dust” simulations. Also in “passive dust” simulations, the PDFs of dust density do not have the same shape or qualitative scalings as those presented here (compare Hopkins & Lee 2016; Lee et al. 2017, or the discussion in Moseley et al. 2018). In most cases we study, the dust density fluctuations are vastly larger in our simulations (with “active” dust which can drive the RDIs) as compared to an analogous “passive-dust” case. For example, with $\tilde{\alpha}\sim 0.001$, similar to some of our most strongly-unstable and dust-clumped cases in HII-near, HII-far, and AGB L, Lee et al. 2017 typically found $\delta\ln{(\rho_{d}/\rho_{g})}\sim 0.01$ for “passive” grains (orders-of-magnitude smaller than our result here). Perhaps most importantly, the scaling of the strength of dust clustering with $\tilde{\alpha}$ ($L_{\rm box}$) or $\tilde{\phi}$ ($\tau$) in “passive dust” studies is, in many cases, almost opposite those here (e.g. Lee et al. 2017 found fluctuations in dust-to-gas ratio with “passive” dust were only strong at $\tilde{\alpha}\gg 1$). 3. 3.

Details of the gas equation-of-state, the functional form of the drag law (Coulomb+Epstein or just Epstein), or the grain charge scaling (dependence on local temperature and density) do not qualitatively alter our conclusions, although they certainly have quantitative effects (Table 2, and Figs. 11, 13, and 12). Larger dust charge generally produces more violent saturation (Figs. 12 and 13). The non-linear behavior of the instabilities does not depend sensitively on particular alignments or anti-alignments between acceleration and magnetic field directions (except insofar as the resonant angles change; e.g. Fig. 15), and in fact, cases where the two are more strongly anti-aligned can even grow faster, despite weaker grain drift. Likewise, modest variation in parameters like the equilibrium grain drift velocity (relative dust-gas acceleration) or magnetic $\beta$ do not qualitatively alter the behavior or character of saturation (Fig. 11 & 14). Lowering the dust-to-gas ratio produces slower initial growth and weaker gas turbulence, as expected. However, surprisingly, it can produce non-linear clustering in the dust that is as strong, or even stronger, than higher-$\mu$ cases (Figs. 16, 18, 28, 29). 4. 4.

Most systems are driven towards approximate equipartition between gas velocity and magnetic field fluctuations. This agreement – i.e.

[TABLE]

is at the order-of-magnitude level, as shown in Fig. 10.191919There are a few notable exceptions with $\langle\rho_{g}\rangle\,\delta{\bf u}_{g}^{2}/2\gg\delta{\bf B}^{2}/8\pi$. Two are WIM-L and Corona-L:$\tau$=100, although the variant WIM-L:LoV (“low drift velocity”) and “default” (higher-$\tau$) Corona-L have $\langle\rho_{g}\rangle\,\delta{\bf u}_{g}^{2}/2\approx\delta{\bf B}^{2}/8\pi$ to within $\sim 10\%$. In both of these exceptions, the gas moves nearly incompressibly and two-dimensionally, so the ${\bf B}$ fields are moved in the $xy$ plane but not compressed, generating negligible ${\bf B}$ fluctuations (see Fig. 7). Our default (high-$\tau$) CGM boxes also exhibit low $\delta{\bf B}$, though here it may be because the fluctuations are dominated by small-$k$, random gas motions which do not cause an effective coherent dynamo. This result is independent of the initial $\beta$ (from $\beta_{0}\sim 0.001-1000$). In some cases this involves strong amplification of ${\bf B}$ fields (e.g. from $\beta\sim 1$ to $\sim 0.01$ in AGB-XL). Because of this, while the instabilities can drive highly super-sonic turbulence in some cases, it is usually trans-Alfvénic. 5. 5.

All of the instabilities examined saturate with sustained gas turbulence. To rough order-of-magnitude, saturation often occurs when the eddy turnover timescale on the box scale becomes shorter than the box-scale linear growth timescale i.e. $|\delta{\bf u}_{g}|/L_{\rm box}\sim\Im{(\omega[L_{\rm box}])}$. However, as discussed below, for both very small boxes and some magnetically-dominated boxes, other criteria (e.g. equipartition between magnetic tension and driving by dust) may instead set the saturation amplitude (see Fig. 10). In any of these cases, the strength of the saturated gas turbulence increases with dust-to-gas ratio and box size/wavelength (see Fig. 10, Table 2). This can be understood physically, since the forcing of the dust onto the gas becomes stronger relative to pressure and magnetic forces. Provided some gas velocity fluctuations, the gas density fluctuations roughly follow the usual relation for pure isothermal MHD turbulence,

[TABLE]

but with substantial variation in the “compressibility” $b\sim 0.2-1$ (see Fig. 10). ${\bf B}$ can be related to ${\bf u}_{g}$ as described above. Anisotropy in the gas properties can usually be understood as a direct reflection of the anisotropy in the fastest-growing linear modes at the box scale (see below). 6. 6.

The dust saturation is ubiquitously more complex than the gas saturation. In some cases, the dust exhibits extremely strong “clumping” or clustering, with a wide range of distinct morphologies and topologies (e.g. differently-connected sheets, filaments, or clumps). In the most extreme cases simulated here, the dust over-densities reach magnitudes of $\sim 10^{9-12}$ times the mean, as seen in Figs. 21-29! In other cases, the dust is “dispersed” throughout the box, with nearly-isotropic, large, velocity dispersions. Qualitatively, the anisotropy of ${\bf v}_{d}$ and relation between ${\bf v}_{d}$ and ${\bf u}_{g}$ reflect those of the fastest-growing linear modes at the box scale (see below). Generically, on “intermediate” and “large” scales $k\lesssim 1/(\mu\,c_{s}\,t_{s})$, we expect and see $|\delta{\bf v}_{d}|\gtrsim|\delta{\bf u}_{g}|$, while on “small” scales $k\gtrsim 1/(\mu\,c_{s}\,t_{s})$, $|\delta{\bf v}_{d}|\lesssim|\delta{\bf u}_{g}|$ (see Fig. 10).202020This follows from the linear RDI behavior and can be understood from a local-balance-type argument from the equations for “forcing” the dust via gas. Dimensionally, a linear perturbation $\delta{\bf v}_{d}$ should have $\omega\,\delta{\bf v}_{d}\sim\delta{\bf u}_{g}/t_{0}$ where $t_{0}\approx t_{s}$ if drag dominates, or $t_{0}\approx t_{L}$ if Lorentz forces dominate, so $|\delta{\bf v}_{d}|\sim|\delta{\bf u}_{g}|/(\omega\,t_{0})$. But generically for the MHD-wave RDI-type modes, $\omega\sim t_{0}^{-1}\,(\mu\,k\,c_{s}\,t_{0})^{\nu}$ with $\nu\sim 1/3-2/3$ depending on wavelength, so the scaling switches from $|\delta{\bf v}_{d}|\gtrsim|\delta{\bf u}_{g}|$ to $|\delta{\bf v}_{d}|\lesssim|\delta{\bf u}_{g}|$ around $k\sim 1/(\mu\,c_{s}\,t_{0})$. On the very largest scales, the dust density fluctuations $\delta\rho_{d}$ are comparable to gas density fluctuations $\delta\rho_{g}$. But, while $\delta\rho_{g}$ decreases with scale, $\delta\rho_{d}$ does not, because there is no internal pressure resisting compressions. In fact, some of the most extreme dust-density fluctuations appear when the gas is nearly incompressible, and they actually become stronger at lower dust-to-gas ratios (Figs. 28-29). 7. 7.

The statistics of both dust and gas fluctuations are often highly non-Gaussian (Figs. 21-29), with exponential or “stretched exponential” tails and, in some cases, coherent sub-structure. This is generally associated with strong intermittency and stochastic driving in dissipative systems (§ 4.4). These strong deviations from Gaussianity mean, for example, that the mass-weighted dust density fluctuations can deviate substantially from volume-weighted fluctuations, and in some cases a significant fraction of the dust mass ($>0.1\%$) can reside at values $\gtrsim 15-20\,$ standard deviations from the median.

5.3 Saturation: “Clumped” States

Although it is clear here that the saturated states are diverse and occupy a continuum of properties, we attempt to classify them into two very broad “regimes”, based on their morphology and resemblance to intuition from linear theory. First, we note that despite their obvious differences, boxes Example, AGB, HII-near and HII-far have several qualitative properties in common. These runs all have $\beta\gtrsim 1$, and $\tau\lesssim 100$, a value that is not too large. They all share a defining feature, that the dust is strongly “clumped” and remains highly anisotropic even in saturation. Prominent clumps, plumes, filaments, and sheets appear, even when the gas is only weakly compressible.

In these, the “medium” and “large” boxes (Example, AGB M/L/XL, HII-near M/L, HII-far M/L) saturate with significant anisotropy or bias in $\delta{\bf u}_{g}$ along the direction of the acceleration ${\bf a}$ (as opposed to e.g. ${\bf w}_{s}$ or ${\bf B}$). The components in the perpendicular direction are not negligible and the strength of the anisotropy varies, owing to mixing from the Lorentz forces. These runs also generally have

[TABLE]

(kinetic/magnetic energies similar) with anisotropies oriented in the same plane(s). Moreover

[TABLE]

with $\delta{\bf v}_{d}$ typically slightly larger, but not by more than a factor of a few. In the largest boxes, $\delta\rho_{d}/\rho_{d}\sim\delta\rho_{g}/\rho_{g}$, (Fig. 10), and the PDFs become increasingly Gaussian/lognormal (Figs. 22 and 23) especially at high-$\mu$ (Fig. 18). In the intermediate-size-scale (M) boxes $\delta\rho_{d}/\rho_{d}>\delta\rho_{g}/\rho_{g}$.

These behaviors can all, remarkably, be predicted (at least qualitatively) by the linear properties of the fastest-growing modes at the box scale. These predictions are discussed in detail in Hopkins & Squire (2018b) (Fig. 2) and Paper I (§ 4-5), and we briefly summarize them here. If the relevant modes at the box scale are the MHD-wave RDI modes or the aligned modes ($\hat{\bf k}=\hat{{\bf w}}_{s}$, also called “pressure-free” or quasi-sound/drift modes), then for modest magnetization ($\beta^{-1}$ and $\tau$ not too large) most of the insight can be gained from considering the much simpler pure hydro case (see Paper I for further discussion). At $k\lesssim\mu/c_{s}\,t_{s}$ (AGB L/XL and HII-near L, where $L_{\rm box}/c_{s}\,t_{s}\gtrsim 1000$), the aligned “pressure-free” mode dominates, where internal pressure effects of the gas are weak compared to the bulk force from dust on gas. In this mode $\delta\rho_{d}/\langle\rho_{d}\rangle\approx\delta\rho_{g}/\langle\rho_{g}\rangle$ and $\delta{\bf v}_{d}\approx\delta{\bf u}_{g}$ (i.e. dust and gas fluctuate together; see Fig. 10), with $\delta{\bf u}_{g}\propto\hat{\bf k}$ (fluctuations are longitudinal) and maximum growth rates at $\hat{\bf k}=\hat{{\bf w}}_{s}$. Because of the weak pressure effects, ${\bf B}$ is driven passively by the velocity fluctuations so $|\delta{{\bf B}}|/|{\bf B}|\sim|\delta{\bf u}_{g}|/v_{A}$ (with $\delta{{\bf B}}$ orthogonal to $\delta{\bf u}_{g}$ in the ${\bf B}-\delta{\bf u}_{g}$ plane).212121For example, if we assume $\delta{\bf u}_{g}\propto{\bf a}\propto(\sin{\theta_{{\bf B}{\bf a}}},\,0,\,\cos{\theta_{{\bf B}{\bf a}}})$ in $(\hat{x},\,\hat{y},\,\hat{z})$, then since $\hat{\bf B}=\hat{z}$, for linear perturbations $\delta{\bf B}\propto(-\cos{\theta_{{\bf B}{\bf a}}},\,0,\,\sin{\theta_{{\bf B}{\bf a}}})$. So for our runs with $|{\bf B}_{0}\cdot{\bf a}|=\cos{\theta_{{\bf B}{\bf a}}}=1/\sqrt{2}$, this gives $\delta{\bf u}_{g}\propto(1/\sqrt{2},\,0,\,1/\sqrt{2})$ and $\delta{\bf B}\propto(-1/\sqrt{2},\,0,\,1/\sqrt{2})$ (i.e. the absolute magnitude of the anisotropy is the same in each direction for $\delta{\bf u}_{g}$ and $\delta{\bf B}$). For our runs with $|{\bf B}_{0}\cdot{\bf a}|=0.05$, this gives $\delta{\bf u}_{g}\propto(0.99,\,0,\,0.05)$ and $\delta{\bf B}\propto(-0.05,\,0,\,0.99)$ (so the dominant direction for $\delta{\bf u}_{g}$ is $\hat{x}$, while that for $\delta{\bf B}$ is $\hat{z}$). These compare well to the results in Table 2. Note that initially, $\hat{{\bf w}}_{s}$ is not aligned with ${\bf a}$ (for non-zero $\tau$), so the modes produce the “sheets” of overdensity in dust perpendicular to $\hat{{\bf w}}_{s}$. However, because the pressure effects are weak in these modes, the non-linear forcing from ${\bf a}$ tends to overwhelm competing forces like magnetic tension, and push the system to drift in the $\hat{{\bf a}}$ direction (giving $\delta{\bf u}_{g}\propto{\bf a}$).

In the “intermediate” boxes (Example, AGB, HII-near, and HII-far M) the “mid-$k$” MHD-wave modes dominate at $\mu\lesssim k\,c_{s}\,t_{s}\lesssim\mu^{-1}$. Again the linear modes have $\delta{\bf v}_{d}\sim\delta{\bf u}_{g}$, $|\delta{{\bf B}}|/|{\bf B}|\sim|\delta{\bf u}_{g}|/v_{A}$. At mid-$k$, the initially fastest-growing modes approximately satisfy $\hat{\bf k}\bot{\bf B}$ if $|{\bf w}_{s}|\ll c_{s}$ (the Alfvén or slow RDIs) or $\hat{\bf k}\bot\hat{{\bf w}}_{s}$ if $|{\bf w}_{s}|\gg c_{s}$ (the fast RDI). This produces the perpendicular sheets and filaments extended along $\hat{{\bf w}}_{s}$, which are seen at early times. It also explains the observed anisotropies, although these are weaker because the linear modes have a mix of components in each direction. Perhaps most notably, these boxes have $\delta\rho_{d}/\langle\rho_{d}\rangle\gg\delta\rho_{g}/\langle\rho_{g}\rangle$, which, as shown in Paper I, are likely related to the linear modes, which satisfy $\delta\rho_{d}\sim\delta\rho_{g}\,\Im{(\omega\,t_{s})}$, i.e. $\delta\rho_{d}/\langle\rho_{d}\rangle\approx(\Im{(\omega\,t_{s})}/\mu)\,\delta\rho_{g}/\langle\rho_{g}\rangle$. Because $\Im{(\omega\,t_{s})}\sim(k\,c_{s}\,t_{s})^{1/2}$ in this mode, the relative strength of $\delta\rho_{d}/\delta\rho_{g}$ increases at smaller scales and smaller $\mu$, consistent with our experiments (Figs. 10 and 28-29). While this provides a reasonable qualitative explanation for the observed trends, we do caution that the magnitude of the saturated $\delta\rho_{d}/\delta\rho_{g}$ is often significantly larger than that predicted by linear theory.

In the smaller boxes (AGB S, HII-near S, HII-far S), the “high-$k$” MHD-wave modes dominate, with $k\gtrsim 1/(\mu\,c_{s}\,t_{s})$. The fastest-growing mode directions are the same as the “mid-$k$” modes (Paper I), with a similar anisotropy (here $\delta{\bf v}_{d}\propto\delta{\bf u}_{g}\propto\hat{\bf k}$ to leading order, giving anisotropy in the $xy$ plane for $\delta{\bf u}_{g}$, and in $z$ for $\delta{\bf B}$). Again $|\delta{{\bf B}}|/|{\bf B}|\sim|\delta{\bf u}_{g}|/v_{A}$ in the linear mode and saturated turbulence (Table 2). As in the mid-$k$ modes, we have $\delta\rho_{d}\sim\delta\rho_{g}\,\Im{(\omega\,t_{s})}$, but now with $\Im{(\omega\,t_{s})}\sim(k\,c_{s}\,t_{s})^{1/3}\gtrsim 1$, so the ratio $\delta\rho_{d}/\delta\rho_{g}$ continues to rise with smaller $\mu$ or at smaller $k$ (making $\delta\rho_{d}$ weakly-dependent on box size, with a small decrease to smaller $L_{\rm box}$). This occurs even though $\delta\rho_{g}\rightarrow 0$ (see Fig. 10). One notable difference from the mid-$k$ modes, however, is that the linear perturbations feature $|\delta{\bf v}_{d}|\sim|\delta{\bf u}_{g}|/\Im{(\omega\,t_{s})}\ll|\delta{\bf u}_{g}|$. This feature is also seen in the saturated turbulence (Fig. 10).

We stress that despite their common elements, there are important differences across these clumped boxes, beyond just the magnitude of the effects. The morphology, topology, and even dimensionality (e.g. clumps, filaments, sheets) of the dust structures varies and depends on a complicated mix of both the global parameters (e.g. $|{\bf w}_{s}|/c_{s}$, $\tau$, $\beta$, etc.), as well as scale, owing to the complex superposition of different modes. Box AGB, with initial $\tau\ll 1$, is closest to the pure-hydrodynamic cases studied in Moseley et al. (2018). As a result it saturates in primarily compressible, supersonic magnetosonic turbulence, with the saturation amplitudes for ${\bf u}_{g}$ in boxes M/L/XL well-predicted by the eddy turnover time argument (tested in detail therein), and $\delta{\bf B}$ following from $\delta{\bf u}_{g}$. Box Example, with higher $\tau\sim 30$, saturates in primarily incompressible MHD turbulence (see Seligman et al., 2018). In this case, the saturation amplitude of $\delta{\bf B}$ (especially its variations with $\mu$) is more accurately predicted by assuming force balance between forcing from dust and magnetic tension of the dominant (box-scale) modes, with $\delta{\bf u}_{g}$ following from $\delta{\bf B}$, i.e.:

[TABLE]

see Fig. 10 and Seligman et al. 2018).

Although it is beyond the scope of this work to study cases where the external drift driving these instabilities initially is time-variable, it is worth noting that even if that drift were somehow “turned off” (which should allow the induced turbulence to decay), there is no obvious mechanism to disperse the dust-to-gas fluctuations formed. Also, in environments with some externally driven turbulence, it would be interesting to explore whether the net effect of this turbulence is to enhance the dust-to-gas ratio fluctuations (as occurs in the absence of RDIs; see Hopkins & Lee 2016; Lee et al. 2017) or to limit the saturation of the RDI-induced clumping.

5.4 Saturation: “Dispersed” or “Granular” States

We refer to the second regime as “Dispersed” or “Granular”, because the dust is generally more dispersive in these runs (boxes CGM, WIM, Corona). It appears that the transition between the two regimes occurs as $\tau$ becomes very large, specifically $\tau\gtrsim 100$. In this regime, the saturated states of the instabilities begins to differ from the description above, and the dust has more isotropic velocity dispersion and notably smaller density fluctuations (especially at the high-$\rho_{d}$ end, which is suppressed relative to low-$\rho_{d}$; see Figs. 25-27).

Many of the saturated properties are consistent with the dominant linear modes, as observed in the previous regime. Unlike the “clustered” boxes, which are dominated by a combination of the low-$k$ “pressure free” (and quasi-sound/drift) modes and mid/high-$k$ MHD-wave modes, at sufficiently high-$\tau$ the instabilities become increasingly confined along ${\bf B}$. Boxes CGM, WIM, Corona are dominated by a combination of the strong ${\bf B}$-aligned “cosmic-ray-like” instabilities (see Paper I for details), together with the related gyro RDIs in WIM S and Corona S/M (Fig. 2). In linear theory, the fastest-growing modes (in both cases) have

[TABLE]

with ${\bf B}$ field fluctuations transverse ($\delta{\bf B}$ preferentially in the $xy$ plane, similar to an Alfvén wave), with222222As discussed in Fig. 7, WIM L (the default run, with larger drift velocity) is the one notable exception with $|\delta{\bf B}|/|{\bf B}|\ll|\delta{\bf u}_{g}|/v_{A}$.

[TABLE]

Like the “clumped” case the gas fluctuations $\delta\rho_{g}$ can be related to $\delta{\bf u}_{g}$ with the usual MHD turbulence scalings. However, the anisotropy is often weak, because (1) the overall turbulence is isotropized and (2) the linear modes have components in all directions. Also like the “clumped” case, the intermediate/large-scale boxes (CGM M/L, Corona M/L, WIM L) have a saturation amplitude of the gas that is reasonably well explained by equating eddy turnover and growth timescales ($\delta{\bf u}_{g}\sim\Im{(\omega[L_{\rm box}]\,L_{\rm box})}$),232323The notable apparent exceptions in Fig. 10 are the CGM low-$\tau_{\rm low}$ runs. However, Fig. 13 show this initially grows vigorously at high-$k$ ($\sim 100/L_{\rm box}$) and in fact it reaches a large $|\delta{\bf u}_{g}|\sim c_{s}$, before the “sheets” break up and disperse the dust suppressing growth of larger-scale modes, and $\delta{\bf u}_{g}$ actually decays somewhat before reaching its equilibrium value. If we use the higher $k$ at which the largest rapidly growing modes are present, and the larger ${\bf u}_{g}$, before the isotropized dust orbits lead to less coherent turbulent motions in the gas, then these runs are plausibly consistent with the “eddy turnover time” saturation scaling. while the small-scale242424As discussed in Paper I, when $\tau\gg 1$, the dividing line between “small” and “intermediate” scales is not simply $k\gtrsim 1/(\mu\,c_{s}\,t_{s})$, but can become a rather complicated expression of $\tau$, $\beta$, etc. Here they can be effectively defined by the presence in Fig. 2 of modes with the “high-$k$” asymptotic scalings. boxes (Corona S, WIM S/M) have a saturation amplitude $\delta{\bf u}_{g}$ that is better explained by the same “high-$k$” scaling as the “clumped” cases above (see final paragraph of § 5.3).

In the linear modes of these high-$\tau$ cases, the ${\bf v}_{d}$ perturbation is approximately a gyro orbit, i.e. preferentially equal power in the $xy$ direction. The scaling of $\delta{\bf v}_{d}/\delta{\bf u}_{g}$ is similar to the “mid-$k$” and “high-$k$” MHD-wave cases discussed above (i.e. substantially smaller in the high-$k$ limit) but enhanced by a factor between $\sim\mu^{-1/2}$ and $\mu^{-1}$ (depending on wavelength in the out-of-resonance gyro or cosmic-ray like mode; see Paper I, § 4). This is directly evident in Fig. 10, which shows that $|\delta{\bf v}_{d}|/|\delta{\bf u}_{g}|\gtrsim 1$ in many of these strongly magnetized cases. Finally, $\delta\rho_{d}$ also scales qualitatively like the mid/high-$k$ MHD-wave modes, in that it is weakly dependent on $L_{\rm box}$ or $k$ (while $\delta\rho_{g}$ decreases at lower $L_{\rm box}$). However, in both the aligned cosmic-ray-like and gyro modes, the Lorentz motion is (to leading order) incompressible, with the dust density fluctuations suppressed by a factor $\sim\mu^{1/2}$ and gas density fluctuations by a factor $\sim\mu$ (see Paper I, § 6.4). This suggests that $\delta\rho_{d}\sim\exp{(\delta\ln{\rho_{d}})}-1$ should be an order-of-magnitude lower compared to similar “clumped” runs. This intuition provides a surprisingly good fit to the difference between “clumped” and “dispersed” runs in Fig. 10.

This regime is analogously very broad, and there is no single, transcendent behavior that defines it. In Corona M/L and CGM, the aligned modes initially produce rather thick “sheets” in the $xy$ plane perpendicular to ${\bf B}$. These form as dust particles move slowly relative to each other along the field lines, and collapse into thin sheets, with an increase in $\delta\rho_{d}$. However, once the sheet becomes thin, the acceleration on dust in direction $\hat{{\bf a}}\neq\hat{{\bf B}}$ pushes the dust with a component transverse to ${\bf B}$ only at one point along ${\bf B}$. This excites gyro motion of dust about ${\bf B}$, but also drags the field and “bends” ${\bf B}$ locally (as opposed to simply pushing the entire field line uniformly, as in the initial state), generating a magnetic tension and exciting Alfvén waves. That, in turn, can re-orient the gyro motion (bending or “dispersing” the sheet). How “isotropized” the dust is – and, correspondingly, how uniformly the dust is spread – depends on how easily the field can be bent. Thus in CGM, with high-$\beta$, the fields and corresponding ${\bf v}_{d}$ can be fully isotropized; in contrast, in Corona-M, the low-$\beta$ and small scales mean the energy in the dust cannot fully re-orient the fields, and the dust motion remains primarily in the $xy$ plane. The maximum dust velocity dispersion is set by equating the “pumping” of the gyro motion (acceleration ${\bf a}$) with damping by drag, which just gives

[TABLE]

with isotropic ${\bf w}_{s}$.

In contrast, in Corona S, and WIM S/M, the dust is collected in dense “lines,” or, more precisely, closed vertical sheets or “tubes”, oriented along the ${\bf B}$ direction. The dominant modes are gyro modes (not the aligned modes), so the “horizontal sheets” discussed above do not form.252525We have checked that there is no apparent correlation between the presence of these tubes and whether or not the box is large enough that the dust “wrapping” time ($L_{\rm box}/{\bf w}_{s}$) is shorter or longer than the mode growth time, which might artificially contribute to such structures. We have also checked for any dependence on whether the dust gyro radii are resolved or unresolved in the gas cells (because of our super-particle approach, dust orbits are always resolved). This produces the “granular” appearance of $xy$ slices. Without the “sheets” to bend the fields strongly at individual points, the field lines move coherently in aligned “columns” or “tubes.” Therefore, the turbulence is essentially two-dimensional, and because the scales are small the gas is weakly-compressible. These features are also evident in box WIM-L, although this also has some aligned modes, giving it a mix of properties.

It is also worth noting that in some boxes (e.g. Corona-S), the modes do not appear to reach the grid scale, even after the simulation has been evolved for longer than the box-scale mode growth time. This tends to occur when there are strongly-growing, dominant gyro modes on somewhat smaller scales which can retain their dominance, even at late times. There is some hint that these may continue growing, along with box-scale modes, more slowly (e.g. linearly in time). Unfortunately it is computationally prohibitive to run these boxes to arbitrarily long times.

We emphasize that there does not appear to be a sharp “threshold” where behavior changes between “clumped” and “dispersed” modes. Obviously, some “clumped” runs feature gyro modes and closer-to-isotropic grain velocity dispersions, while some “disperse” runs here (like Corona-M) retain coherent dense grain structures well into their non-linear evolution. Rather, there is a spectrum of different behavior in different regimes, and different parameters involve a different mixture of these.

6 Other Physics (Not Included Here)

The simulations here are intentionally idealized, designed to study the physics and non-linear behaviors of the instabilities identified in Paper I. Of course, in different specific physical applications, there are an infinite variety of additional physics which could also be important, some of which may modify the instabilities themselves. Some of these cases will be studied in future work modeling observables in specific physical systems (e.g. AGB outflows and dust extinction/density estimation; Steinwandel et al., in prep.), or analytic studies focused on the linear instabilities in more complex systems (e.g. three-fluid radiation-MHD RDIs, Squire et al., in prep.). Other examples (viscosity, non-ideal MHD, stratification, differential rotation) are discussed in detail in Hopkins & Squire (2018b), Squire & Hopkins (2018a), and Hopkins & Squire (2018a). But we briefly review some here for context.

Small-scale transport processes (e.g. viscosity, conductivity) are discussed in Paper I. For the physical examples which motivate our parameter study in Table 1, we specifically chose the parameters so that our smallest resolvable scales in the smallest (“S”) boxes correspond to the scale where viscous effects could begin to become significant ($k^{2}\gtrsim|\omega|/\nu$, where $\nu$ is the kinematic viscosity), so that all resolved scales in all our boxes correspond to the “inertial range” where viscous effects are weak (see § 9 in Paper I, where these scales are estimated). But we also note, as shown explicitly in Paper I (e.g. Fig. 7 there and Squire & Hopkins 2018a), that most of the instabilities here are undamped below the viscous scale (down to at least the ion gyro-radii), like e.g. the cosmic ray resonant and non-resonant instabilities (Wentzel, 1968; Bell, 2004). We also choose motivating parameters where non-ideal MHD effects should be negligible (these terms generally require an ionization fraction $f_{\rm ion}\ll 10^{-8}$ to be significant for the modes studied here; see § 8.3 in Paper I), but note that even in overwhelmingly neutral gas (e.g. proto-planetary disks), ambipolar diffusion and the Hall effect do not actually damp many of the instabilities studied here (but instead modify them into additional branches; see e.g. Squire & Hopkins 2018a, Figs. 7-8). Other small-scale effects such as the current carried by dust grains (§ 8.3 in Paper I), finite inter-grain separation effects, or degeneracy pressure are completely negligible in any regime motivating our study here, or any regime where ideal MHD could actually apply (e.g. the dust current is sub-dominant by a factor $\sim 10^{10}\,f_{\rm ion}\,(T/10^{4}\,{\rm K})$ in the induction equation).

On the largest spatial scales, global effects such as stratification, gradients in the external forces driving drift, rotation and shear make our local periodic-box approximation invalid. We have therefore specifically chosen parameters such that the largest scales in our largest boxes (“L/XL”) in Table 1 correspond to these scales. All of these global effects introduce new RDIs such as the Brunt-Väisälä, vertical settling, and epicyclic RDIs (see e.g. Appendix C in Hopkins & Squire 2018b or § 5 & 6 in Squire & Hopkins 2018a).

There are also various effects which directly involve dust-gas interactions and operate on relatively long timescales such as gas cooling; dust accretion, mantle formation, and chemical enrichment (generally relevant in cold or very dense gas like AGB outflows); or sputtering of dust in hot gas (relevant in gas with $T\gtrsim 10^{6}\,{\rm K}\,(\mu/0.01)$). These can modify the equation-of-state of the gas, dust-to-gas ratios or dust size/charge, but generally operate on timescales much longer than (or at most comparable to) the RDI growth times. It is therefore unlikely that these strongly modify the RDIs directly: instead a medium with e.g. efficient cooling or sputtering would simply shift the parameters of interest here, having effectively lower temperature or dust-to-gas ratios, respectively. But it is likely that these processes are themselves strongly modified by the RDIs, as they depend sensitively on clumping factors of the dust and gas and their (grain-size-dependent) cross-correlation, a subject worthy of future study.

One effect which could modify the instabilities here is externally-driven turbulence. Although we clearly show the RDIs persist in the presence of self-excited turbulence, it is possible that if some other process drives gas turbulence much more strongly (on a given scale) it could significantly modify the RDIs. This could potentially set the saturation amplitude by e.g. shearing apart different modes if the eddy turnover times become much shorter than RDI growth times (see Squire & Hopkins 2018a, § 8), although it is also well-known from both theory and laboratory experiments that even test-particle grains ($\mu=0$, so no RDIs are present) are actually strongly clustered on small-scales in externally-driven turbulence (see e.g. Squires & Eaton, 1991; Fessler et al., 1994; Bracco et al., 1999; Cuzzi et al., 2001; Rouson & Eaton, 2001; Monchaux et al., 2010; Pan et al., 2011; Hopkins & Lee, 2016). This is most likely to be important in physical systems motivating boxes like WIM or CGM:$\tau_{\rm low}$, where strength of the saturated gas turbulence from the RDIs alone (Table 2) is significantly weaker than the level of turbulence observed.

A variety of physical effects will also become important specifically when the local dust-to-gas ratio $\mu$ becomes very large. The most obvious is grain-grain collisions, as the ratio of stopping to collision time scales as $t_{\rm stop}/t_{\rm collision}\sim\mu\,(\delta v_{\rm dd}/c_{s})$ (where $\delta v_{\rm dd}$ is the rms relative grain-grain velocity difference on infinitesimally small spatial scales). Obviously this will become important at the very high $\mu$ reached by some of our simulations on small scales, but how it modifies the distribution of $\mu$ itself and subsequent RDIs depends on the outcome of said collisions (which can produce grain growth/coagulation, elastic or inelastic scattering/bouncing, or grain shattering, depending sensitively on the local conditions and grain chemistry). In some cases, the self-gravity of a cluster of grains can become important before collisions (though this is more likely important in proto-planetary disks, it could occur under exceptional ISM conditions, see e.g. Hopkins 2014). And if radiation pressure drives the dust drift, then the dust can become self-shielding if the dust-to-gas ratio within a patch of size $\sim\lambda$ exceeds $\mu(\lambda)\gg\tilde{\alpha}\,(\lambda/L_{\rm box})^{-1}\,(1+\lambda_{\rm rad}/\epsilon_{\rm grain})$ (where $Q\sim(1+\lambda_{\rm rad}/\epsilon_{\rm grain})^{-1}$ crudely approximates the grain absorption efficiency for incident radiation with wavelength $\lambda_{\rm rad}$). In that regime the the radiation can be multiply-scattered, leading to more complicated three-fluid (radiation, dust, gas) RDIs. The very high local $\mu$ reached here in the absence of these effects should be considered motivation for studying these different regimes.

7 Summary & Conclusions

We have presented the first simulation parameter study of the non-linear regime of the Squire & Hopkins (2018b) “resonant drag instabilities” of charged dust grains in magnetized gas. Because the parameter space of the instabilities is large and this is a first study, we focused on several sets of initial conditions that are broadly representative of different astrophysical regimes. In all cases studied, the instabilities produce highly non-linear structures and fluctuations, often including strong turbulence and magnetic field amplification. Strong anisotropy and non-linear features appear in the dust, including orbit-crossing, density fluctuations, and complicated velocity distribution functions. This necessitates numerical simulations using dust particle methods which can follow the full non-linear dust velocity distribution function, as opposed to just the fluid limit.

Even within our small survey, the simulations exhibit a diverse array of behaviors. In Paper I we demonstrated using linear perturbation theory that a homogeneous MHD gas coupled to dust via Lorentz forces and drag exhibits around $10$ different “instability families”. These families include the Alfvén, slow, and fast magnetosonic MHD-wave RDIs; the three corresponding gyro RDIs; the “pressure-free” and related “quasi-drift” and “quasi-sound” modes; and the “cosmic ray streaming”-type modes. These families all have different linear growth rates and mode structures, but often overlap and occur within the same system. Without introducing additional physics or constraints (e.g. Braginskii viscosity, which suppresses the growth rate of the magnetosonic modes, but not the Alfvén modes), it is generally impossible to construct a simulation that isolates a single instability family. Our range of initial conditions were chosen both to be representative of different astrophysical regimes and also to exhibit different dominant, fastest-growing instability families. The resultant non-linear evolution yields a remarkable diversity of outcomes.

We show that over the course of the simulations, the instabilities become violently non-linear on all scales. Their non-linear outcomes can result in dust being highly concentrated (the “clumped” regime), or dispersed with large isotropic velocity dispersions (the “disperse” regime). In the “clumped” case, there are a wide range of morphologies, topologies, and dimensionalities of the clumped structures, depending on the parameters of the system, with dust in e.g. multiply-connected sheets, filaments, or point-like clumps. In the most extreme cases, the dust can reach enormous overdensities in these idealized tests ($>0.1\%$ of the dust mass at $\gg 10^{9}$ times the mean dust density, with volume-averaged “clumping factors” $>10^{4}$). The dust clumping does not depend systematically on the spatial scale or the compressibility of the gas: clumping can be stronger on small scales in nearly incompressible gas than on large scales in highly compressible cases. Surprisingly, the clumping is stronger at lower dust-to-gas ratios, although the growth times of the instabilities are longer.

In contrast, in the “disperse” cases, the dust can be accelerated to highly super-sonic isotropic velocity dispersions (even undergoing first-order Fermi acceleration with stretched-exponential velocity “tails” that would reach relativistic velocities in some physical systems) and dispersed nearly-uniformly over the box (clumping factors as small as $\sim 1.02$). These cases are akin to well-studied cosmic ray instabilities that self-excite diffusive behavior (i.e. self-generating dust diffusion). The growth times of these instabilities can be shorter than either the dust drag/stopping or gyro times, and are often extremely short relative to other timescales in the gas (e.g. the sound-crossing or dynamical times, for physical systems of interest). This will have a huge range of important physical ramifications for essentially all regimes where dust is present.

The dust drives anisotropic turbulence in the gas, whose properties may be qualitatively understood via heuristic quasi-linear theory. The gas turbulence is stronger and more compressible if the box scale is larger, or if the dust-to-gas ratio is larger. Moreover, the gas velocity and density fluctuations are correlated in approximately the same manner observed for pure-gas MHD turbulence. In the saturated state, the instabilities tend to produce equipartition between gas velocity and magnetic field fluctuations. The latter means that systems can have the magnetic fields strongly amplified by the instabilities (e.g. decreasing $\beta$ by factors of $\sim 10^{4}$, in the most extreme cases here).

The gas turbulence can be qualitatively magnetosonic or Alfvénic, but in addition to the dust properties, it has many characteristics that are unique. PDFs of various fluctuations (${\bf u}_{g}$, ${\bf B}$, $\rho_{g}$, ${\bf v}_{d}$, $\rho_{d}$) are typically highly non-Gaussian, with exponential or stretched-exponential tails indicative of strong intermittency and stochastic driving in a highly-dissipative system. These large tails mean that sometimes the mass and volume-weighted statistics can yield quite different results, and the tails can be many orders-of-magnitude more populated than Gaussian (e.g. many of our systems have $>0.1\%$ of their dust or gas mass at $\gtrsim 10-20$ “standard deviations” in some property).

The parameter space here is highly multi-dimensional, so it is difficult to make conclusions that can be robustly applied to all regimes. However it appears that choices such as the detailed equation-of-state, form of the dust charge scaling with ambient gas properties, exact magnetic $\beta$, initial direction of field orientiation (relative to the relative dust-gas acceleration), or exact dust-to-gas ratio do not dramatically (qualitatively) change the character of the solutions. As anticipated from linear theory in Paper I, the most important parameters that determine the qualitative behavior appear to be the physical scale and ratio of magnetic (Lorentz) to drag (aerodynamic+Coulomb) forces on dust.

Future work will be necessary to investigate how these instabilities will manifest in physical systems. We stress that the models here are intentionally idealized: we follow a single grain species subject to a constant external differential dust/gas acceleration in periodic, initially homogeneous gas boxes, with an ideal equation-of-state. Therefore our model names (e.g. AGB) should not be taken literally – these are not intended to be realistic physical simulations of those systems. Rather, the names are chosen reflect an example of a system where the key dimensionless parameters for these instabilities ($\tau$, $\beta$, $|{\bf w}_{s}|/c_{s}$) are similar to the box simulated. Nor is it obvious, yet, how the simulations here directly translate to observational consequences: even if the behavior in fully “realistic” systems is similar to that predicted here, this will clearly have consequences for dust extinction curves and emission, cooling and dust chemistry, and many other areas, but the magnitude or even the sign of these effects depends sensitively on the exact observations considered, as well as physical conditions (e.g. chemical conditions, optical depth, additional radiative transfer effects) beyond those modeled here. The physical systems simply provide helpful motivation for our survey and in future work we will explore more realistic scenarios.

This investigation has begun to elucidate the broad non-linear behaviors of the resonant drag instabilities in magnetized gas. However, it raises more questions than it answers, some of which include:

What are the effects of a broad spectrum of grain sizes and charges? Here we intentionally simulate a single species of grains, so that the growth rates and dominant modes can be clearly defined and studied. However in almost all astrophysical situations there will be a wide range of grain sizes. In some circumstances the largest grains (which tend to dominate the dust mass) will dominate the dynamics, while in others there is an intricate mix of which grains dominate which terms in the relevant equations (for a more detailed discussion, see Paper I). Further simulations will be required to understand when different grain sizes are effectively independent, and when they will have strong non-linear interactions via the gas (see e.g. Bai & Stone, 2010a, for examples in the streaming RDI). 2. 2.

Is there a meaningful way to define convergence, and incorporate the effects of all relevant size scales in a single simulation? The dynamic range over which the instabilities are present with interesting growth rates is enormous, and far larger than can be simulated at present. Because we find that the gas turbulence has most of its power at the largest (driving) scales, there is hope that, like MHD turbulence, certain bulk properties (e.g. the bulk power in turbulence, the dissipation rate, and most of its effects) can be “converged” even if the Kolmogorov scale is unresolved. It is unclear if a meaningful convergence criteria can even be defined for the dust, since some properties, like the dust density fluctuations, are not uniquely dominated by the driving-scale modes. 3. 3.

What is the effect of stratification or time-variability on the instabilities? As the size scale of the simulations increases, the physical system would become stratified and non-uniform in space or time. We have shown that stratification and other large-scale terms such as shear, rotation, and differential acceleration (e.g. Coriolis forces) can all introduce additional instabilities, some of which have faster growth rates on large scales than those here (see Squire & Hopkins, 2018a). So, this could introduce unique and important phenomena. On large time-scales, the “forcing” terms driving dust drift (hence the instabilities) could fluctuate: if this occurs rapidly compared to growth timescale it may produce unique phenomena as well. It would be particularly interesting to study the analog of “decaying turbulence” when external acceleration/forcing is suddenly “turned off” once the simulations have reached saturation. 4. 4.

What is the nature of the non-linear gas turbulence and its intermittency? Any of these simulations provide ample opportunity to study the nature of the turbulence. We have only explored simple, “zeroth order” diagnostics, but more detailed studies of the Eulerian and Lagrangian power spectra and structure functions will provide substantial insight into the structure of the turbulence. The highly non-Gaussian behavior we see suggests that the turbulence may be quite different from a simple Kolmogorov-type picture (although of course pure-MHD turbulence is already substantially non-Gaussian in some measures). 5. 5.

Can we develop a predictive theory for the dust clumping and turbulence? Here we are able to develop some qualitative insights into when dust should exhibit clustering. Further, from quasi-linear theory, we can heuristically understand why the dust density fluctuations become stronger at lower dust-to-gas ratio, and are weakly-dependent on the spatial scale or compressibility of the gas. However, this is far from a predictive or quantitative theory. Most past work on the non-linear saturation of particle clustering in turbulence assumed “passive” grains (i.e. neglected the forces of dust on gas). However, in the “passive” case these instabilities do not exist, and the predicted saturated dust clustering (neglecting back-reaction terms) can be very different from what we find here. It is clearly important to develop new theoretical models which can explain the strength and nature of the observed dust structures. 6. 6.

What are the important regimes of parameter space yet to explore? Our parameter survey is far from complete. We have chosen a few examples which are interesting and plausibly motivated, but there are many variations possible within similar physical systems. Many physical systems with dust have fast-growing RDIs, but with parameters quite different to those studied here (e.g. AGN and their outflows, dense GMCs, supernova remnants, the ISM of primordial neutral galaxies, planetary atmospheres, proto-stellar disks; see Paper I for further discussion). 7. 7.

How does the introduction of additional grain physics modify the predictions here? As noted above, in the non-linear regime, some of our simulations reach dust concentrations $\gg 10^{9}$ times the mean density. Obviously, other physics will become important long before these densities are reached, including e.g. dust self-shielding (if radiation drives the differential acceleration), self-gravity of the grains, grain collisions, and dust current. Some of these could suppress clustering, but others could make it stronger. The high densities and low local relative grain-grain velocities in these regimes may make them ideal for coagulation. But much of the volume being dust-poor might suppress dust growth via accretion of ions from the gas. In the disperse regime, some of our simulations produce isotropic dust velocities many times the sound speed, which may introduce sputtering and shattering in grain collisions.

Finally, we stress that we have explored only the ideal MHD case of the Squire & Hopkins (2018b) instabilities. A host of other instabilities exist which appear when other physics are present (e.g. external or self gravity, stratification, centrifugal or Coriolis forces, non-ideal MHD, kinetic MHD, strong coupling of multiply-scattered radiation, etc.). In future work, we hope to explore these cases in more detail, together with some of the questions above.

Acknowledgments

We thank Alexander Kaurov and Ulrich Steinwandel for a number of enlightening discussions and useful comments. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214. Support for J.S. was provided by Marsden Fund grant UOO1727 and a Rutherford Discovery Fellowship, managed through the Royal Society Te Apārangi. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.