The fundamentally different dynamics of dust and gas in molecular clouds

Abstract

We study the behaviour of large dust grains in turbulent molecular clouds (MCs). In primarily neutral regions, dust grains move as aerodynamic particles, not necessarily with the gas. We therefore directly simulate, for the first time, the behaviour of aerodynamic grains in highly supersonic, magnetohydrodynamic turbulence typical of MCs. We show that, under these conditions, grains with sizes a ≳ 0.01 micron exhibit dramatic (exceeding factor ∼1000) fluctuations in the local dust-to-gas ratio (implying large small-scale variations in abundances, dust cooling rates, and dynamics). The dust can form highly filamentary structures (which would be observed in both dust emission and extinction), which can be much thinner than the characteristic width of gas filaments. Sometimes, the dust and gas filaments are not even in the same location. The ‘clumping factor’ |$\langle n_{\rm dust}^{2} \rangle / \langle n_{\rm dust} \rangle ^{2}$| of the dust (critical for dust growth/coagulation/shattering) can reach ∼100, for grains in the ideal size range. The dust clustering is maximized around scales ∼ 0.2 pc (a/μm) (n_gas/100 cm^{− 3})^{− 1}, and is ‘averaged out’ on larger scales. However, because the density varies widely in supersonic turbulence, the dynamic range of scales (and interesting grain sizes) for these fluctuations is much broader than in the subsonic case. Our results are applicable to MCs of essentially all sizes and densities, but we note how Lorentz forces and other physics (neglected here) may change them in some regimes. We discuss the potentially dramatic consequences for star formation, dust growth and destruction, and dust-based observations of MCs.

accretion, accretion discs, instabilities, turbulence, planets and satellites: formation, galaxies: formation, cosmology: theory

1 INTRODUCTION

Dust is ubiquitous in astrophysics, and critical to understanding phenomena as diverse as star and planet formation, feedback from stars in galaxy formation, and the origin and fate of certain heavy elements. Even if it is only to correct foreground contamination or extinction, understanding the dust size distribution and dust-to-gas ratio, and any possible variations (hence variations in the extinction curve, for example) is necessary to almost every area of astronomy.

Despite this, there has been little theoretical work to understand the dynamics of dust as aerodynamic particles in the cold interstellar medium (ISM). For example, a critical process, which could produce fundamentally new phenomena, is the inevitable fluctuation of large dust grain densities in a turbulent medium (so-called turbulent concentration). It is well known that in a primarily neutral, dense gas, massive dust grains (which contain a large fraction of all the ISM metals) behave as aerodynamic particles (the dominant force is drag from collisions with atoms/molecules). As such, they can, under the right conditions, decouple from the gas, and can clump or disperse independent from gas density fluctuations.

Much attention has, in fact, been paid to the question of grain density fluctuations (arising from this mechanism and others) in protoplanetary discs. When stirred by turbulence or trapped in various instabilities, the number density of grains can fluctuate by orders of magnitude relative to the gas. This has been seen now in a wide variety of situations, including or excluding grain collisions, in magnetized and non-magnetized discs, and in turbulence driven by self-exciting (‘streaming’) instabilities, gravitational instabilities, the magneto-rotational instability, convection, and Kelvin–Helmholtz instabilities (see e.g. Bracco et al. 1999; Cuzzi et al. 2001; Youdin & Goodman 2005; Johansen & Youdin 2007; Carballido, Stone & Turner 2008; Bai & Stone 2010a,b; Pan et al. 2011; Dittrich, Klahr & Johansen 2013; Jalali 2013; Hopkins 2014a). In the terrestrial turbulence literature as well, ‘preferential concentration’ of aerodynamic particles is well studied with both laboratory experiments (Squires & Eaton 1991; Fessler, Kulick & Eaton 1994; Rouson & Eaton 2001; Gualtieri, Picano & Casciola 2009; Monchaux, Bourgoin & Cartellier 2010) and numerical experiments (Cuzzi et al. 2001; Yoshimoto & Goto 2007; Hogan & Cuzzi 2007; Bec et al. 2009; Pan et al. 2011; Monchaux, Bourgoin & Cartellier 2012) demonstrating that gas is unstable to the growth of large-amplitude inhomogeneities in the grain density.

Many authors have pointed out that the relevant phenomena appear to be scale free, if the dust grains are sufficiently large so that their ‘stopping’ (friction or drag) time-scale t_s corresponds to an eddy turnover time t_e for eddies which lie within the inertial range of turbulence (Cuzzi et al. 2001; Hogan & Cuzzi 2007; Bec et al. 2009; Olla 2010; Hopkins 2016). Qualitatively, if t_s ≪ t_e, grains are well coupled to gas, so should move with the flow (although they may still ‘settle’ and exhibit non-trivial dynamics); if t_s ≫ t_e, grains are effectively decoupled from the local gas flow; when t_s ∼ t_e, grains can be ‘flung out’ of regions of high vorticity by centrifugal forces, and collect in regions of high strain (Yoshimoto & Goto 2007; Bec et al. 2008; Wilkinson, Mehlig & Gustavsson 2010; Gustavsson et al. 2012).

In principle, the same mechanisms which would generate large grain density fluctuations in a protoplanetary disc with ∼1–100 cm boulders could operate on micron-sized dust in a giant molecular cloud (GMC). In fact, Hopkins (2014b) applied the analytic scalings derived as a function of the dimensionless ratio of grain stopping time to dynamical time in protoplanetary discs to GMCs, and argued that this implied micron-sized dust could cluster strongly on scales large enough to alter stellar abundances (potentially dramatically, in the most extreme cases).

However, essentially all of the numerical and experimental studies to date have focused on sub-sonic, either incompressible or weakly compressible gas with inefficient cooling (adiabatic), usually without magnetic fields. In the cold ISM, on the other hand, the gas is rapidly cooling (effectively close to isothermal), magnetized, highly compressible (with density fluctuations of factors of thousands), and supersonically turbulent with Mach numbers ≳ 10 on the scales of large clouds. Moreover, in incompressible gas, the velocity field is divergence free, so the vorticity and strain dominate grain aggregation; in highly supersonic turbulence, the density field is a network of filamentary shocks and rarefactions, each of which can trap or disperse dust (Booth, Sijacki & Clarke 2015). And the ‘stopping time’ t_s in the supersonic case is not a constant (as it is for grains of a given size in a sub-sonic medium), but depends on the local density and dust-gas relative velocity. It is not obvious that the dynamics should be even qualitatively similar in these cases.

In this paper, we therefore for the first time explore the dynamics of dust grains in a neutral, highly compressible, supersonically turbulent medium, under conditions which resemble observed atomic/molecular clouds, clumps, and cores. We show that a wide range of dust grain sizes show dramatic clustering effects, in many cases more than would be expected from simply scaling up the sub-sonic case, and discuss the implications for different observations and theoretical models.

2 METHODS

2.1 Magnetohydrodynamics and turbulent driving

Our simulations use Gizmo (Hopkins 2015),¹ a mesh-free, Lagrangian finite-volume Godunov code designed to capture advantages of both grid-based and smoothed-particle hydrodynamics (SPH) methods, built on the gravity solver and domain decomposition algorithms of Gadget-3 (Springel 2005). In Hopkins (2015) and Hopkins & Raives (2015) we consider extensive surveys of test problems in both hydrodynamics and MHD, and demonstrate accuracy and convergence in good agreement with well-studied regular-mesh finite-volume Godunov methods and moving-mesh codes (e.g. Athena & Arepo; Stone et al. 2008; Springel 2010). We run Gizmo in its Meshless-Finite Mass (MFM) mode but have verified that Meshless Finite-Volume (MFV) mode produces nearly identical results (as expected from the previous studies). Note that in Hopkins (2015) and, Hopkins & Raives (2015), we demonstrate excellent agreement between Gizmo and high-resolution, state-of-the-art moving mesh and grid-based codes for simulations of both supersonic and sub-sonic MHD turbulence.

The turbulent driving routines follow Bauer & Springel (2012). Briefly, a periodic box is stirred via the usual method in e.g. Schmidt, Federrath & Klessen (2008), Federrath, Klessen & Schmidt (2008), Price & Federrath (2010), where a small range of modes corresponding to wavelengths between 1/2and1 times the box size are driven in Fourier space as an Ornstein–Uhlenbeck process, with the compressive part of the acceleration projected out via a Helmholtz decomposition in Fourier space so that the driving is an adjustable mix of compressible modes and incompressible/solenoidal modes. Our specific implementation has been verified for both hydro and MHD cases in Hopkins (2013a) and Hopkins (2015), Hopkins & Raives (2015). Our parameter choices for the driving generally follow Bauer & Springel (2012), table 4, but we specify the relevant Mach numbers below. We initialize a uniform seed field |${\boldsymbol {B}} = B_{0}\hat{z}$|⁠, which determines the saturated mean-field strength.

2.2 Dust dynamics

We follow Carballido et al. (2008), Hogan, Cuzzi & Dobrovolskis (1999), Johansen & Youdin (2007), Johansen, Youdin & Mac Low (2009), Bai & Stone (2010a), Pan et al. (2011) and model the dust via a collection of ‘superparticles,’ each one of which represents an ensemble of grains of a fixed size, whose trajectories are integrated on-the-fly through the fluid. This is essentially a Monte Carlo ‘tracer particle’ approach.

Draine & Salpeter (1979a) show that grains obey the following equation of motion:

\begin{eqnarray} \frac{{\rm d}{\boldsymbol {u}}_{d}}{{\rm d}t} &=& -\frac{{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}}{t_{\rm {s}}} \end{eqnarray}

(1)

\begin{eqnarray} t_{\rm {s}} &\equiv & \frac{\pi ^{1/2}}{2\sqrt{2}}\,\left( \frac{\bar{\rho }_{d}\,a_{d}}{c_{\rm {s}}\,\rho _{\rm gas}} \right)\,\left( 1 + \left|\frac{3\pi ^{1/2}}{8}\,\frac{{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}}{c_{\rm {s}}} \right|^{2}\right)^{-1/2} \end{eqnarray}

(2)

where |${\boldsymbol {u}}_{d}$| is the grain velocity, d/dt is a Lagrangian derivative, c_s and ρ_gas the isothermal sound speed and density of the gas, |$\bar{\rho }_{d}\approx 2.4\,{\rm g\,cm^{-3}}$| is the internal (material) grain density (Draine 2003), and a_d ∼ 0.001-1 μm is the grain radius.²

Note that in the sub-sonic limit, this becomes the well-studied Stokes expression with constant ‘stopping time’ t_s. In the supersonic case, c_s and ρ_gas depend on position (since the gas is compressible), and the term in |${\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}$| can be important. This replaces the sound speed with the ‘total gas-dust’ velocity |${\sim } \sqrt{c_{\rm {s}}^{2} + |{\boldsymbol {u}}_{d} - {\boldsymbol {u}}_{\rm gas}|^{2}}$|⁠.

We solve this equation by kernel-interpolating the quantities c_s, ρ_gas, and |${\boldsymbol {u}}_{\rm gas}$| and their derivatives from gas particle positions (where they are determined by the MHD solver) to the grain particle position, i.e. |$\rho _{\rm gas}({\boldsymbol {x}}_{d,\,i}) = \sum \,W({\boldsymbol {x}}_{{\rm gas},\,j} - {\boldsymbol {x}}_{d,\,i},h_{i})\,\rho _{\rm gas}({\boldsymbol {x}}_{{\rm gas},\,j})$| where W is normalized so that |$1=\sum \,W({\boldsymbol {x}}_{{\rm gas},\,j} - {\boldsymbol {x}}_{d,\,i},\,h_{i})$|⁠. The functional form of W and kernel size h are identical to those used for hydrodynamic operations (see Hopkins 2015), so the interpolation is numerically stable and consistent. Equation (1) is then solved exactly over half-timesteps Δt/2, assuming the interpolated quantities vary linearly in time and space, and we use this to determine the mean numerical acceleration over the corresponding interval |$\bar{{\boldsymbol {a}}}_{i}(t,\,t+\Delta t/2) = [{\boldsymbol {u}}_{d}^{\rm exact}(t+\Delta t/2)-{\boldsymbol {u}}_{d}(t)]/[\Delta t/2]$|⁠. This maintains good behaviour even in the limit t_s ≪ Δt. The particle trajectories are then integrated with a semi-implicit leapfrog scheme (this is already well tested in our code for collisionless particles in e.g. cosmological simulations). To ensure numerical stability, grain particles obey the usual timestep limits for all collisionless particles (e.g. |${\rm d}t < {\rm MIN}[0.1\,(h/|{\rm d}{\boldsymbol {u}_{d}}/{\rm d}t|)^{1/2},\,0.2\,|\nabla \cdot {\boldsymbol {u}}_{d}|^{-1}]$|⁠), plus a Courant criterion given by the minimum of the timestep of neighbour gas particles or |$0.05\,h/\sqrt{c_{\rm {s}}^{2}+|{\boldsymbol {u}}_{\rm gas}-{\boldsymbol {u}}_{d}|^{2}}$|⁠.

2.3 Units

We adopt an isothermal equation of state (γ = 1) for the gas; this is a reasonable assumption for molecular clouds over the density and temperature range of interest here, and enables more direct comparison with previous studies of supersonic turbulence and star formation (see e.g. Li, Mac Low & Klessen 2005; Krumholz, Klein & McKee 2007; Federrath et al. 2008; Schmidt et al. 2009; Kritsuk et al. 2011; Molina et al. 2012; Konstandin et al. 2012b; Hopkins 2013b).

With this assumption, the inviscid ideal MHD equations are inherently scale-free. The box length L_box, mean gas density |$\langle \rho _{\rm gas} \rangle \equiv M_{\rm gas} / L_{\rm box}^{3}$|⁠, and sound speed c_s can therefore be freely rescaled to any physical values. Physically, this means our results are entirely determined by three dimensionless numbers: the box-averaged Mach number |$\mathcal {M} \equiv \langle |{\boldsymbol {u}}_{\rm gas}|^{2} \rangle ^{1/2} / c_{\rm {s}}$|⁠, the (coherent) magnetic mean-field strength |$|\langle {\boldsymbol {B}} \rangle | / (\rho ^{1/2}\,c_{\rm {s}})$|⁠, and the ‘grain size parameter’ α:

\begin{eqnarray} \alpha \equiv \frac{\bar{\rho }_{d}\,a_{d}}{\langle \rho _{\rm gas} \rangle \,L_{\rm box}} \end{eqnarray}

(3)

(i.e. the value of |$\bar{\rho }_{d}\,a_{d}$| in code units; note that only this combination of grain size and density appears in the equations, so a_d and |$\bar{\rho }_{d}$| are formally degenerate).

A grain parameter α translates to physical grain size a_d as:

\begin{eqnarray} a_{d}&=&0.4\,\alpha \,\mu {\rm m}\,\left( \frac{L_{\rm box}}{10\,{\rm pc}} \right)\,\left( \frac{\langle n_{\rm gas} \rangle }{10\,{\rm cm^{-3}}} \right)\,\left( \frac{\bar{\rho }_{d}}{2.4\,{\rm g\,cm^{-3}}} \right)^{-1} \end{eqnarray}

(4)

To aid in rescaling to physical units, if we assume the simulations sample Milky Way-like GMCs that lie on the observed linewidth–size relation (⁠|$\mathcal {M} \sim (R/R_{\rm sonic})^{1/2}$| with R_sonic ∼ 0.1 pc being the sonic length), and size–mass relation (M_GMC ∝ R_GMC, or 〈Σ_GMC〉 ∼ 300 M_⊙ pc⁻² ∼ constant), and that our boxes sample ‘typical’ sub-regions of the clouds, we obtain:

\begin{eqnarray} L_{\rm box}&&\sim 10\,{\rm pc}\,\left(\frac{\mathcal {M}_{\rm box}}{10}\right)^{2}\,\left( \frac{R_{\rm sonic}}{0.1\,{\rm pc}} \right) \end{eqnarray}

(5)

\begin{eqnarray} \langle n_{\rm gas} \rangle &\equiv & \frac{\langle \rho _{\rm gas} \rangle }{\mu \,m_{p}} \sim 10\,{\rm cm^{-3}}\,\left( \frac{\langle \Sigma _{\rm GMC}\rangle }{100\,{\rm M}_{{\odot }}\,{\rm pc^{-2}}} \right)\,\left( \frac{R_{\rm GMC}}{100\,{\rm pc}} \right)^{-1} \end{eqnarray}

(6)

where we assume a mean molecular weight μ ≈ 2.3 to convert between ρ_gas and n_gas.

Moreover, to the extent that turbulence is (approximately) self-similar, we can think of our lower Mach number simulations as sampling ‘sub-volumes’ of our higher Mach number simulations (similar to how the effective box size scales with Mach number above, if we assume a linewidth–size relation). At infinite resolution, a sufficiently large box of high |$\mathcal {M}$| should contain all possible realizations of smaller-|$\mathcal {M}$| sub-regions. Moreover, in the limit where the initial mean-field is weak (⁠|$|\langle {\boldsymbol {B}} \rangle | / (\rho ^{1/2}\,c_{\rm {s}}) \lesssim 1$|⁠), the magnetic field dynamics are dominated by those produced by the turbulent dynamo itself (⁠|$\langle |{\boldsymbol {B}}| \rangle \gg |\langle {\boldsymbol {B}} \rangle |$|⁠) and are just a function of |$\mathcal {M}$|⁠, and the mean-field value is irrelevant. This is (by construction) the case in many of our simulations, although we also consider strong mean-field cases and show they have weak effects on grain clustering.

So in a sense, there is really one dominant dimensionless parameter (α) which specifies the physics. But because we are limited by computational cost, it is more convenient to run separate boxes of different |$\mathcal {M}$|⁠; given the linewidth–size relation above, resolving the same smallest physical scale as in one of our 256³, |$\mathcal {M}=2$| boxes in a |$\mathcal {M}=10$| box would require a ∼6400³ (trillion-particle) simulation!

Because we do not explicitly include the ‘back-reaction’ of dust grains on gas, the absolute value of the dust-phase metallicity (or equivalently, the mean dust abundance) |$Z_{d} \equiv \rho _{\rm dust} / \rho _{\rm gas} = (4\pi /3\,\bar{\rho }_{d}\,a_{d}^{3}\,n_{\rm dust}) / (\mu \,m_{p}\,n_{\rm gas})$| does not enter our equations. The simulations predict relative fluctuations in Z_d but can be freely rescaled to any mean 〈Z_d〉, modulo caveats below.

2.4 Neglected dust physics

We neglect several processes in this study:

• Dust–dust collisions: the mean free path to these is |${\sim }(n_{\rm dust}\,\sigma _{d})^{-1}\sim \alpha \,Z_{d}^{-1}\,L_{\rm box}$| (see Section 2.3 for definitions); stopping/deceleration lengths in the gas are |${\sim } \alpha \,(1+\mathcal {M}^{2})^{-1/2}\,L_{\rm box}$|⁠, so dust-dust collisions are sub-dominant by a factor |${\sim } Z_{d}/(1+\mathcal {M}^{2})^{1/2} \ll 1$|⁠. This does not mean dust collisions are uninteresting, as they can play a key role in modifying the dust size distribution over time. But while the dust dynamics we study may critically alter dust collisions, dust collisions do not (usually) significantly alter the dust dynamics.

• Destruction/creation: we study dust dynamics in cold cloud regions, so we do not consider sources of new dust and/or destruction by shocks/sputtering from SNe and stellar winds, although these can occur inside GMCs later in the cloud lifetime. Turbulent shocks inside the cold cloud regions do not reach sufficient temperatures to destroy grains (Draine & Salpeter 1979b).

• Coulomb forces: following Draine & Salpeter (1979a), for low-temperature (T ≲ 10⁵ K) gas we expect the ratio of Coulomb forces to collisional drag (for grains in the size range of interest here) to be ∼10 f_ion, where f_ion ≲ 10⁻⁷ is the ionized fraction of gas in GMCs, so this is negligible.

• Radiation pressure: near massive stars, this can dominate grain dynamics, but not in random portions of the cloud. Assuming geometric absorption, the ratio of radiation pressure to drag forces at a distance r_* from an O-star (luminosity L_*) is |${\sim } 0.5\,(L_{\ast }/10^{4}\,{\rm L}_{\odot })\,(n_{\rm gas}/10\,{\rm cm^{-3}})^{-1}\,(r_{\ast }/{\rm pc})^{-2}\,(|{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}|/$| 10 km s^{− 1})^{− 2}. The radius where radiation pressure dominates is within the Stromgren sphere (i.e. fundamentally different conditions from what we simulate) for essentially all stars.

• Lorentz forces: the acceleration from the Lorentz force is |${\rm d}{\boldsymbol {u}}_{d}/{\rm d}t = (z_{d}\,e/m_{d}\,c)\,({\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{g})\times {\boldsymbol {B}}$| (where z_d e is the grain charge, m_d its mass, |${\boldsymbol {B}}$| the local magnetic field, and c the speed of light); if we combine this with the expectation value of 〈z_d〉 calculated by Draine & Sutin (1987) for grains in cold molecular clouds (as a function of grain size and temperature), the ratio of Lorentz to collisional forces becomes |${\sim } 0.1\,B_{\mu {\rm G}}\,\langle z_{d} \rangle \,a_{\mu {\rm m}}^{-2}\,n_{10\,{\rm cm^{-3}}}^{-1}\,T_{100\,{\rm K}}^{-1/2}\,x^{-1/2} \sim B_{\mu {\rm G}}\,T_{100\,{\rm K}}^{1/2}\,a_{\mu {\rm m}}^{-1}\,n_{10\,{\rm cm^{-3}}}^{-1}\,x^{-1/2}$|⁠, where |$x\equiv 1+9\pi \,|{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{g}|^{2}/64\,c_{\rm {s}}^{2}$|⁠. If magnetic fields follow the expected equipartition of the supersonic turbulent dynamo (E_B ∼ 0.05 E_K; see Kritsuk et al. 2011; Federrath et al. 2014), then we can further simplify and obtain |${\sim } 0.1\,T_{30\,{\rm K}}\,a_{\mu {\rm m}}^{-1}\,n_{10\,{\rm cm^{-3}}}^{-1/2}$|⁠. Therefore, over much of the physical parameter space of interest, Lorentz forces are sub-dominant to collisional drag. However, they are by no means negligible and clearly could be dominant in some regimes (see e.g. Yan, Lazarian & Draine 2004). In future work, we will extend our simulations with explicitly coupled Lorentz force equations, but these are complicated to include with drag in a numerically stable manner and depend on some model for grain charging, so we will neglect them for now (with the appropriate caveats).

• Back-reaction: we neglect the loss of momentum from the gas to the grains, which scales with the dust-to-gas mass ratio Z_d. Unlike the protoplanetary disc case, where grain concentrations might reach Z_d ≫ 100, we have 〈Z_d〉 ∼ 0.01 ≪ 1; so this is usually negligible. It is always negligible for sufficiently low-metallicity clouds. However, the maximum dust concentrations we simulate do correspond to Z_d ≫ 1 if the cloud has solar metallicity, so in future work we will also consider this in more detail.

3 RESULTS

3.1 Qualitative behaviours: critical thresholds for different phenomena

Figs 1 and 2 show images of representative times during some of our simulations. It is clear that the dust and gas dynamics differ, sometimes dramatically.

Figure 1.

Image of a very high-resolution (1024²) 2D simulation of aerodynamic grains in supersonic MHD turbulence, here with grain size parameter |$\alpha \equiv (\bar{\rho }_{d}\,a_{d}) / (\langle \rho _{\rm gas} \rangle \,L_{\rm box}) =0.01$| (see equation 4 for how this relates to physical grain sizes a ∼ 0.001-1 μm) and rms Mach number |$\mathcal {M}\sim 5$|⁠. Time shown is after the rms Mach numbers and magnetic energy reach steady state. We show the full simulation box (left) and zoom-in of a dense region (right). Colour shows gas density relative to the mean (n_gas/〈n_gas〉), on a logarithmic scale (see colour bar); black points show the dust (super)particles. As expected, gas forms a filamentary network of shocks and rarefactions. Dust loosely traces the same on large scales, but with much more detailed small-scale structure. Much of the dust lies along razor-thin filaments, some of which are not associated with a gas filament.

Figure 2.

Images (as Fig. 1) of our standard 3D, 256³ MHD simulations (each shows a thin slice through z = 0), here with Mach number |$\mathcal {M}\sim 10$|⁠. Time is after run each reaches steady state; colour shows gas density; black points show dust. Each image shows a different dust size α = 0.01, 0.1, 1 (top-to-bottom). For small α ≲ 0.1, dust traces gas roughly on large scales, but shows very thin, small-scale filaments which can be overdense relative to the gas (seen in Fig. 1). Intermediate α ∼ 0.1 exhibit more dramatic large-scale dust clumping. Large α ∼ 1 dust is only weakly coupled to the gas, and remains at approximately the mean density everywhere.

Since here the dust does not act on gas, the gas dynamics are identical to those expected for supersonic MHD turbulence. We confirm (for detailed analysis see Hopkins 2015) that the highly supersonic cases here (⁠|$\mathcal {M}\sim 10$|⁠) develop a velocity scaling similar to the linewidth–size relation observed, with rms velocities on scale λ of |$\langle \mathcal {M}^{2}(\lambda ) \rangle \sim \mathcal {M}(L_{\rm box})\,(\lambda /L_{\rm box})^{1/2}$|⁠, down to a sonic length |$R_{\rm sonic}=\lambda (\mathcal {M}=1) \sim L_{\rm box}\,\mathcal {M}(L)^{-2}$|⁠, also as expected from previous work (Scalo et al. 1998; Schmidt et al. 2009; Konstandin et al. 2012a). Below this scale, the turbulence becomes sub-sonic and we expect a Kolmogorov (1941) cascade, |$\mathcal {M}(\lambda <R_{\rm sonic}) \sim (\lambda /R_{\rm sonic})^{1/3}$|⁠. The gas forms a filamentary network of shocks and rarefactions; the characteristic width of filaments and dense structures is of order R_sonic (see Vazquez-Semadeni 1994; Padoan, Nordlund & Jones 1997; Klessen 2000; Kritsuk et al. 2007). In our simulations with initially weak mean-field strengths, the initial magnetic fields grow exponentially until saturating with magnetic energy ∼5 per cent of the kinetic energy (also as expected; Schekochihin et al. 2004; Brandenburg & Subramanian 2005; Federrath et al. 2014). As expected, for our driving routines, after a few dynamical times, the simulations reach a steady state, and the statistics of turbulent velocity fluctuations and dust dynamics do not evolve.

Note that the free-streaming length of a dust grain (relative to gas) can be estimated as |$L_{\rm stream} \sim \langle |{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}| \rangle \,t_{\rm {s}}$|⁠, where we expect (and confirm in our simulations) that the typical relative velocity |$\langle |{\boldsymbol {u}}_{d}-{\boldsymbol {u}}_{\rm gas}| \rangle$| corresponds to the ‘eddy velocity’ of turbulent modes on a scale ∼L_stream (much smaller modes do not strongly perturb the dust, and much larger coherent modes simply entrain both dust and gas together; see Voelk et al. 1980; Ormel & Cuzzi 2007; Pan & Padoan 2010). If we combine this with the expressions for t_s and the velocity scalings above, we can approximate the solution as:

\begin{eqnarray} {L_{\rm stream}} &\sim& \left\lbrace \begin{array}{@{}l@{\quad }l@{}} \displaystyle \frac{\bar{\rho }_{d}\,a}{\rho _{\rm gas}} = \alpha \,\left(\frac{n_{\rm gas}}{\langle n_{\rm gas} \rangle }\right)^{-1}\,L_{\rm box} \ \ \ \ \ { (L_{\rm stream} > R_{\rm sonic}) } \\ {\displaystyle R_{\rm sonic}^{-1/2}\,\left(\frac{\bar{\rho }_{d}\,a}{\rho _{\rm gas}}\right)^{3/2} = \alpha ^{3/2}\,\mathcal {M}\,\left(\frac{n_{\rm gas}}{\langle n_{\rm gas} \rangle }\right)^{-3/2}\,L_{\rm box}}\\ {\qquad\qquad\qquad\qquad\qquad\qquad\qquad (L_{\rm stream} < R_{\rm sonic}) } \end{array}\right. \end{eqnarray}

(7)

where the behaviour differs depending on whether L_stream is above or below the sonic length (the motion is super or sub-sonic). Equating the two, we arrive at the critical α above which L_stream > R_sonic

\begin{eqnarray} \alpha (L_{\rm stream} = R_{\rm sonic}) &=& \mathcal {M}^{-2}\,\left(\frac{n_{\rm gas}}{\langle n_{\rm gas} \rangle }\right)^{-1} \end{eqnarray}

(8)

Consider the following limits:

• α ≳ 1: In isothermal strong shocks the maximum density enhancement is |${\sim } \mathcal {M}^{2}$|⁠; therefore, if α ≳ 1, we expect L_stream ≳ L_sonic even at these large post-shock densities, so dust can stream through even the densest structures assuming they have sizes of order the sonic length. The dust is therefore always near its mean density (little variance in n_dust) while the gas dynamics proceed as usual (large variance in n_gas). This is the weakly coupled limit.

• |$1/\mathcal {M}^{2} \lesssim \alpha \lesssim 1$|⁠: The dust cannot ‘break out’ of the most dense structures. However, it can cluster on scales larger than the sonic length at the mean density. Therefore dust clustering is imprinted on the medium at intermediate densities and then turbulent compressions ‘trap’ the fluctuations in dense regions.

  In very low-density regions where n_gas ≲ α 〈n_gas〉, L_stream → L_box, and the dynamics resemble the weakly coupled case (although the low-density regions are precisely those where our neglect of Lorentz forces on dust is a poor approximation, and these may recouple the dust and gas).

  In high density regions where |$n_{\rm gas} \gtrsim \alpha \,\mathcal {M}^{2}\,\langle n_{\rm gas} \rangle$|⁠, the clustering scale drops below the sonic length, and the dynamics begin to resemble the sub-sonic case. In this limit, grains can be efficiently expelled from regions of high vorticity and trapped along lines of high strain, leading to their alignment in narrow filamentary structures (Cuzzi et al. 2001; Rouson & Eaton 2001; Bec et al. 2009; Pan et al. 2011; Monchaux et al. 2012). Assuming the Kolmogorov scale is arbitrarily small, in this limit the maximum clustering amplitude becomes self-similar, because all grains ‘see’ eddies which lie in an effectively infinite inertial range and are resonant with their own streaming time-scale (so the clustering dynamics for grains of different sizes are simply rescaled in size; Hogan & Cuzzi 2007; Yoshimoto & Goto 2007; Bec et al. 2008). Detailed discussion of this limit can be found in Hopkins (2016); both experiments and numerical simulations of sub-sonic turbulence suggest that, provided infinite resolution, the small-scale dispersion in the dust-to-gas ratio converges to ∼0.4–0.5 dex.

• |$\alpha \ll 1/\mathcal {M}^{2}$|⁠: The dust is strongly coupled down to below the sonic scale, for all densities equal to or larger than the mean. This case is in the sub-sonic limit above over most of the domain. Since the dust is well coupled at the mean density, gas being compressed into structures of order the sonic scale has approximately the mean dust-to-gas ratio, with little scatter. So although the sub-sonic processes described above can operate, the variance we expect at intermediate to high densities is smaller because there are no ‘seed’ fluctuations in the dust-to-gas ratio in the diffuse medium.

Each of these behaviours are illustrated in Figs 1 and 2. Table 1 gives a complete list of the simulations, and summarizes some of their salient properties, which we will discuss below.

Table 1.

3D simulations.

Grain	Mach	Alfvén	Mean	Clumping	Clumping
size	number	Mach	field	factor	factor
α	|$\mathcal {M}$|	|$\mathcal {M}_{A}$|	|$| \langle {\boldsymbol {B}} \rangle |$|	C	(dense gas)
0.001	10	6	0.2	32	42
0.01	10	6	0.2	57	78
0.005–0.015	10	6	0.2	49	67
0.03	10	6	0.2	73	110
0.10	10	6	0.2	39	68
1.0	10	6	0.2	4.9	6.2
0.01	2	5	0.05	21	28
0.1	2	5	0.05	23	33
1.0	2	5	0.05	4.7	5.6
0.03	10	∞	0	52	74
0.03	10	6	0.5	73	110
0.03	10	2.5	4.3	75	120
0.03	10	0.4	27	55	93

Notes. Parameters describing the simulations analysed in the text: Each has unit box length, sound speed, and mean density in code units (see Section 2.3 for physical units), and is a 3D MHD simulation with driven isothermal turbulence assuming a natural mix of compressive and solenoidal modes. All runs have 256³ gas and 2 × 256³ dust particles.

(1) α: Grain size parameter (see equation (4) to relate this to physical grain size). Most runs adopt a single α; however we include one run adopting a uniform logarithmic (random) distribution of α between 0.005–0.015, to test the effects of non-uniform sizes.

(2) |$\mathcal {M}$|⁠: Steady-state turbulent rms Mach number on the box scale (set by the driving routines).

(3) |$\mathcal {M}_{A}$|⁠: Steady-state volume-averaged Alfvén Mach number |$\mathcal {M}_{A} = \langle |{\boldsymbol {v}}| / v_{A} \rangle$|⁠, where the Alfvén speed |$v_{A} = |{\boldsymbol {B}}|/{\rho ^{1/2}}$|⁠.

(4) |$| \langle {\boldsymbol {B}} \rangle |$|⁠: Time-averaged mean-field strength (in code units): |$| \langle {\boldsymbol {B}} \rangle |_{\rm phys} \sim 1.3\,\mu G\,| \langle {\boldsymbol {B}} \rangle |_{\rm code}\,\sqrt{(\langle n_{\rm gas} \rangle /10\,{\rm cm^{-3}})\,(T/100\,{\rm K})}$|⁠. Note that rms field strengths can be much larger.

(5) C: Time-averaged dust clumping factor |$C \equiv \langle n_{\rm dust}^{2} \rangle / \langle n_{\rm dust} \rangle ^{2}$|⁠.

(6) C (dense gas): Dust clumping factor measured only in the dense (n_gas > 〈n_gas〉) gas.

3.2 The distribution of dust and gas densities

We now consider these effects more quantitatively. In Figs 3 and 4, we measure the density of both dust and gas around every dust particle in the simulation, averaged on a smoothing scale h_min,³ and plot the normalized 2D histogram of all points in the n_gas − n_dust plane. Since, as noted above, the turbulence becomes steady state after ∼1 crossing time, the results at any individual time output are statistically identical; we therefore simply combine all outputs after the first few crossing times to reduce the sampling noise.

Figure 3.

Distribution of dust and gas densities in the 3D, Mach number |$\mathcal {M}\sim 10$| (Alfvén |$\mathcal {M}_{A}\sim 6$|⁠) simulations from Fig. 2. We measure dust and gas density around each dust particle, at all times after the system reaches steady state. We plot isodensity contours at fixed probability density levels dP/d log n_gas d log n_dust = 10⁻¹, 10⁻², 10⁻⁴, 10⁻⁷ (black, green, blue, red, respectively). Each column shows a different grain size parameter α = 0.001, 0.01, 0.03, 0.1, 1 ( left-to-right). Each row smooths the density around each particle on a progressively larger scales corresponding to a spherical radius h_min = 0 (no smoothing, top), h_min = 0.001 L_box ( middle), and h_min = 0.01 L_box (approximately the sonic length R_sonic; bottom). We show n_dust = 〈n_dust〉 (dust at constant density) and n_dust = (〈n_dust〉/〈n_gas〉) n_gas (perfect dust-gas coupling, i.e. δ = 1) as dotted lines. At high α ≳ 1, the dust is decoupled from the gas and remains close to 〈n_dust〉 independent of n_gas. At low α ≪ 0.001, the dust is tightly coupled to the gas, with significant scatter on very small scales that is quickly averaged-out on larger scales. At intermediate α ∼ 0.01–0.1, the dust decouples from gas at low gas densities, and traces it on average at high densities, but with large fluctuations in n_dust/n_gas even at the highest gas densities. These fluctuations are on larger scales for the larger α, and for α ≳ 0.01 are only weakly averaged-down as we smooth over scales as large as ∼R_sonic. For α = 0.01, we also compare (thin lines) one case where the grains have not a single size but a random distribution of sizes between α = 0.005–0.015; because grains in the same locations experience different drag, this slightly reduces the maximum dust clustering, but the effect is weak.

Figure 4.

Distribution of dust and gas densities in the 3D, Mach number |$\mathcal {M}\sim 2$| simulations; style is identical to Fig. 3. As expected, the lower Mach numbers produce much smaller gas-density fluctuations, so the dynamic range sampled is significantly smaller than Fig. 3. At fixed n_gas within the range probed, the dynamics are similar, and the dust-to-gas ratio fluctuates by a similar (albeit slightly smaller) amount, compared to |$\mathcal {M}=10$|⁠.

As expected from our arguments above, at sufficiently low n_gas and large α, the dust decouples from the gas, residing at the mean dust density independent of the gas density, with small fluctuations. At higher densities, the dust tracks the gas on average, 〈n_dust(n_gas)〉 ∝ n_gas, but with obvious scatter. Note the scatter is much larger than the Poisson noise and is numerically converged (see Appendix B).

Because the dust does not alter the gas dynamics in our simulations, the bivariate distribution P(n_gas, n_dust) is separable into P(n_gas) P(n_dust | n_gas). The distribution of gas density P(n_gas) in isothermal MHD turbulence is well studied, and we confirm the usual lognormal form (Passot, Pouquet & Woodward 1988; Vazquez-Semadeni 1994; Konstandin et al. 2012b; Molina et al. 2012; Federrath & Banerjee 2015) with subtle non-lognormal deviations consistent with those predicted in Hopkins (2013b) and confirmed in Federrath (2013). The interesting behaviour we wish to study here is encapsulated in the non-universal dust-to-gas ratio P(n_dust | n_gas).

Fig. 5 therefore collapses our 2D distribution functions into the distribution of dust-to-gas ratio, which we define for convenience relative to the mean in the box:

\begin{eqnarray} \delta &\equiv & \frac{ n_{\rm dust}/ n_{\rm gas}}{\langle n_{\rm dust} \rangle / \langle n_{\rm gas} \rangle } \end{eqnarray}

(9)

We find in every case a broad distribution, with a lognormal ‘core’ and low-δ behaviour, and a power-law-like tail at the highest δ.

Figure 5.

Distribution of dust-to-gas ratios. We calculate δ (dust-to-gas ratio relative to mean) around each dust particle (at any gas density) in the distribution functions from Figs 3 and 4, and plot the resulting PDF dP/d log (δ) (time-averaged over the simulation). Top: Mach number |$\mathcal {M}\sim 10$| cases from Fig. 3 (different α as labelled). As in the previous figure we also show the case with α = 0.005–0.015 as the thin line in the same style as the α = 0.01 case; the difference owing to a distribution of grain sizes of modest width is small. Solid lines show the case with h_min = 0; dotted lines show h_min = 0.01 L_box. For h_min > 0, the mean of the distributions is shifted by averaging closer to δ ∼ 1, but the scatter is only reduced by a modest amount. The core of each distribution is approximately lognormal, but with a large ‘tail’ towards high-δ; this primarily arises in the low-density regions where the dust and gas decouple (see Fig. 6). Bottom: same, for the |$\mathcal {M}=2$| cases. The scatter is reduced, primarily because the lower |$\mathcal {M}$| means the dynamic range of gas densities is much smaller; this most noticeably suppresses the ‘tail’ coming from very low n_gas.

The origin of this behaviour is demonstrated in Fig. 6, which plots δ at fixed n_gas; this is equivalent to P(n_dust | n_gas). At a given gas density, δ (or n_dust) is distributed approximately lognormally. The dispersion σ of the lognormal depends relatively weakly on n_gas, while the mean 〈ln (δ)〉_dust shifts systematically. At low gas densities (n_gas ≲ α 〈n_gas〉) there is some excess at high-δ with respect to a lognormal fit, but most of the deviation from a lognormal in Fig. 5 arises because it represents an integral over the different lognormal P(n_dust | n_gas) with different mean values.

Figure 6.

Distribution of dust-to-gas ratios (δ) at different gas densities n_gas (colours as labelled), for a single representative simulation (α = 0.03, |$\mathcal {M}=10$|⁠, |$\mathcal {M}_{A}=6$|⁠). Top: distribution measured around each dust particle. At each n_gas, the δ distribution is roughly lognormal (dashed line shows a lognormal with the same mean and variance), albeit with wider tails; this motivates our fitting function (equation 10). The variance and mean depend on n_gas; the high-δ tail in Fig. 5 comes from the weakly coupled limit (low n_gas), where n_dust ∼ 〈n_dust〉 so 〈δ〉 ≫ 1. At high n_gas, 〈δ〉 ∼ 1 but with large scatter. Middle: as top, except we sample δ around random points in space (the volume-weighted PDF, instead of the dust-mass weighted PDF; see Section 3.2.1). The relative normalizations and median-δ values shift because low-density regions (with large volume-filling factors) are preferentially sampled, but the scatter is similar. Bottom: same, except we sample δ around random gas elements, (the gas-mass weighted PDF). Median δ values shift closer to ∼1, but the scatter is similar.

A physical argument for why the dust-to-gas ratio should be distributed lognormally in the high-density limit is presented in Hopkins (2016). Essentially, each encounter between a Lagrangian ‘parcel’ of grains and a vorticity/strain structure or eddy in the turbulence imparts an essentially random multiplicative factor on the local grain density (the factor depends on the magnitude of the vorticity/strain, and orientation of the eddy relative to the grain velocity). Integrating over time and structures of a wide range of sizes, this random multiplicative process produces a quasi-lognormal distribution. Unlike the gas density fluctuations (where lognormal behaviour is specific to isothermal gas) this expectation is independent of the gas equation of state, and has been seen in sub-sonic experiments with strictly adiabatic, incompressible gas (Hogan et al. 1999; Bec, Cencini & Hillerbrand 2007; Pan et al. 2011).

Predicting the magnitude of the fluctuations σ is more challenging. Again, some analytic arguments are presented in Hopkins (2016): they show that a single encounter with a ‘resonant’ structure with eddy crossing time ∼t_stop and coherence length ∼L_stream leads to an approximate factor ∼2 (0.3 dex) multiplicative effect on the dust-to-gas ratio. Larger/smaller structures produce weaker effects. Broadly similar multiplicative effects occur when dust grains pass through a shock or rarefaction with gas velocity gradient |$|\nabla \cdot {\boldsymbol {v}}| \sim 1/t_{\rm stop}$| (see e.g. Booth et al. 2015; Lorén-Aguilar & Bate 2015). Assuming the dispersion is dominated by a couple such encounters per global grain-crossing time, values σ ∼ 0.3–0.6 dex (seen in Fig. 6) are plausible. However, the magnitude of this effect depends on the geometry and filling factor of structures in the turbulence, so more detailed analytic models are needed.

3.2.1 Approximate fitting functions

Based on the above, the bivariate distribution of dust and gas densities (around a random dust particle) in Figs 3 and 4 can be approximated by

\begin{eqnarray} &&{\frac{{\rm d}P_{\rm dust}}{{\rm d}\ln {n_{\rm gas}}\,{\rm d}\ln {n_{\rm dust}}} = \frac{1}{2\pi \sqrt{S_{\rm gas}\,S_{\rm dust}}}\, \exp {\left[ -\left(\frac{\Delta _{\rm gas}^{2}}{2\,S_{\rm gas}} + \frac{\Delta _{\rm dust}^{2}}{2\,S_{\rm dust}}\right) \right]}}\nonumber \\ &&\!\!\!\!S_{\rm dust} = S_{\rm dust}(n_{\rm gas},\,h_{\rm min}) = \left( \sigma _{\log _{10}[\delta (n_{\rm gas})]}\,\ln {10} \right)^{2}\nonumber \\ &&\!\!\!\!\Delta _{\rm dust} \equiv \ln {(\delta )} - \langle \ln {\delta (n_{\rm gas})} \rangle _{\rm dust}\nonumber \\ &&\!\!\!\!\Delta _{\rm gas} \equiv \ln {\left(\frac{n_{\rm gas}}{\langle n_{\rm gas}\rangle }\right)} + \frac{S_{\rm gas}}{2} \end{eqnarray}

(10)

where |$S_{\rm gas} = S_{\rm gas}(\mathcal {M})$| is a constant for a given simulation (given by the usual relations in supersonic isothermal turbulence; see Appendix A), while S_dust and 〈ln δ(n_gas)〉_dust are the variance and mean in ln δ at a given n_gas.

Fig. 7 plots the variables in equation (10), as a function of n_gas. Specifically we measure the logarithmic mean 〈log₁₀(δ)〉 = 〈ln (δ)〉/ln 10 and dispersion |$\sigma _{\log _{10}[\delta (n_{\rm gas})]}$| in δ(n_gas).

Figure 7.

Top: 1σ logarithmic dispersion in the dust-to-gas ratio log₁₀(δ) (or equivalently dust density log₁₀(n_dust)) at fixed gas density n_gas (see equation 10). We show the simulations from Fig. 3 with |$\mathcal {M}\sim 10$| (left) and those from Fig. 4 with |$\mathcal {M}\sim 2$| ( right). For both, we show the results smoothed within a radius h_min = 0 (solid) and h_min = 0.01 L_box ( dotted). The variance is generally maximized for α ∼ 0.01–0.1, but is generally only weakly dependent on α and n_gas. Smoothing on larger scales (h_min ∼ 0.01 L_box ∼ R_sonic) decreases the dispersion, but only by a modest ∼0.1 dex. Bottom: Logarithmic mean 〈log₁₀(δ)〉 in the dust-to-gas ratio relative to mean (δ) as a function of gas density. For α ≲ 1, at low densities (n_gas ≲ α〈n_gas〉), the dust decouples from the gas so the mean |$\langle \delta \rangle \propto n_{\rm gas}^{-1}$|⁠; at high densities the dust and gas couple more tightly so 〈δ〉 → 1. For α ≳ 1 the dust streams through sonic-length structures and remains relatively poorly coupled until much higher densities (see Section 3.1).

Figure 8.

Image of dust and gas (as Fig. 2) in simulations with |$\mathcal {M}\sim 10$| and α = 0.03 but varied magnetic mean-field strength: no magnetic field (⁠|$\mathcal {M}_{A}=\infty$|⁠; left) or very strong mean-fields (⁠|$\mathcal {M}_{A}=0.4$|⁠; right). Top: x − z projection, where |$\hat{z}$| is the mean-field direction. Bottom: x − y projection. The strong-field case clearly produces global anisotropy and more filamentary structure; it also slightly enhances the segregation of dust and gas.

The dispersion/variance in log δ is maximized (at most n_gas) for α ∼ 0.01–0.1. It drops for α ≲ 0.001 (when the grains are strongly coupled so track gas closely), and α ≳ 1 (when the grains are weakly coupled, so stay near 〈n_dust〉). In general the dispersion decreases weakly with higher n_gas as grains become more tightly coupled (except α ≳ 1, where the grains are very weakly coupled so only begin to cluster at the highest densities). But in all cases the dependence of |$\sigma _{\log _{10}[\delta (n_{\rm gas})]}$| on n_gas is weak, with typical values ∼0.3–0.6 dex (S_dust ∼ 0.5–2). The mean 〈log₁₀(δ)〉 shows a clear transition: at low n_gas, in the weakly coupled regime, dust resides near 〈n_dust〉 so |$\langle \delta \rangle \sim n_{\rm gas}^{-1}$|⁠, while at high n_gas, in the tightly coupled regime, 〈δ〉 → 1.

Because our default calculation measures the properties around each dust particle, this is the dust-mass-weighted probability distribution function (PDF) dP_dust. We could, instead, uniformly sample the volume (giving the volume-weighted PDF dP_vol) or the gas mass dP_gas. The differences between each are discussed in detail in Appendix A; for the pure point-wise PDF (h_min → 0), they are all trivially related: dP_gas ∝ n_gas dP_vol ∝ (n_gas/n_dust) dP_dust. As expected, then, volume-weighted PDF shifts the weight towards low-density regions, which occupy a larger volume, while the gas-mass weighted PDF shifts the weight towards higher gas densities and gas-to-dust ratios, bringing the average δ closer to unity. However, it is easy to show for equation (10) that these transformations preserve the lognormal shape of the PDF, and do not change the variance S_dust or S_gas. They simply shift the mean values of Δ_dust and Δ_gas.⁴

3.2.2 Dependence on smoothing scale and Mach number

Note in Figs 3 and 4, we consider three values of h_min. First, h_min = 0, i.e. the fluctuations measured on the smallest possible (resolution) scale. Then, h_min = 0.001 L_box and h_min = 0.01 L_box. These give the average density calculated around each point, averaged within a finite volume-averaging spherical radius equal to h_min. Recall, for our |$\mathcal {M}=10$| simulations, h_min = 0.01 L_box corresponds to the sonic length R_sonic. As must occur, the variation in n_dust/n_gas decreases as h_min increases (obviously, in the limit h_min → L_box, all points collapse to exactly the box-averaged 〈n_dust〉 and 〈n_gas〉). For smaller values of α, the variance in n_dust at fixed n_gas decreases more rapidly as we increase h_min – this is because, as noted above, the physical scale of the clustering is smaller for smaller α (all else being equal), so by smoothing on a fixed scale h_min we are averaging-out more of the small-scale variations. Still, for most of our simulations, significant variation persists even with h_min = 0.01.

At low Mach numbers (our |$\mathcal {M}=2$| suite), we robustly find that the variance in log δ at high n_gas/〈n_gas〉 is lower for a given α. This is not surprising. For one, the typical magnitude of gas density fluctuations is smaller in these cases, which in turn generates smaller ‘seed’ fluctuations for dust. More importantly, recall, in the sub-sonic limit, |$L_{\rm stream} \propto \mathcal {M}$| at fixed α and n_gas/〈n_gas〉, while |$\alpha _{\rm {s}} \propto \mathcal {M}^{-2}$| – in other words, the free-streaming length of the dust relative to the box size, and relative to the sonic length (⁠|$L_{\rm stream}/L \propto \mathcal {M}^{3}$|⁠), are lower by factors of ∼5 and ∼125, respectively, in our |$\mathcal {M}\sim 2$| cases. Another way to think of this is to simply note that, in physical units (at fixed α and assuming clouds lie on the local linewidth–size relation), our |$\mathcal {M}\sim 2$| runs correspond to L_box ∼ 0.4 pc instead of ∼10 pc at |$\mathcal {M}=10$|⁠, and correspondingly the physical grain size a_d is a factor ∼25 smaller. So it is not, in fact, physical to think of our |$\mathcal {M}=2$| runs as ‘the same’ physical setup as |$\mathcal {M}=10$| with simply a different Mach number – they are measuring grain fluctuations corresponding to different physical grain sizes on much smaller physical scales. Given this, the remarkable fact is how similar the trends are, implying that the dependence of the dust-to-gas fluctuations on scale, turbulent properties, and grain size are surprisingly weak.

3.2.3 Dependence on mean-field strength

Thus far, we have focused on simulations with small mean (coherent box-scale) magnetic fields. This simplifies our study as the saturated field depends only on the sonic Mach number |$\mathcal {M}$|⁠. However, this does not mean the fields are negligible: in our |$\mathcal {M}\sim 10$| simulations, the rms field strength in physical units is |$\langle |{\boldsymbol {B}}|^{2} \rangle ^{1/2} \sim 4.2\,\mu G \,\sqrt{(\langle n_{\rm gas} \rangle /10\,{\rm cm^{-3}})\,(T/100\,{\rm K})}$|⁠, comparable to observations (see e.g. Elmegreen & Scalo 2004; Brown et al. 2007; Crutcher et al. 2010) in typical clouds.

Even absent Lorentz forces on dust, magnetic fields can alter concentration. For example, fields modify the velocity scalings, imprint local anisotropy, and change the ratio of solenoidal to compressible modes (see e.g. Kowal, Lazarian & Beresnyak 2007; Burkhart et al. 2009; Lemaster & Stone 2009; Kritsuk et al. 2011; Molina et al. 2012; Downes 2012; Federrath & Klessen 2013). All of these can change the ‘response’ of a parcel of grains to an encounter with a turbulent velocity structure (Lazarian & Cho 2004; Yan et al. 2004; Hopkins & Christiansen 2013; Hopkins 2014a, 2016). But these effects depend not just on the rms field strength but also on the mean-field strength, especially as the turbulence goes from super-Alfvénic to sub-Alfvénic (Collins et al. 2012, and references therein).

We therefore consider simulations with α = 0.03 and |$\mathcal {M}\sim 10$| fixed⁵ and varying mean-field |$|\langle {\boldsymbol {B}} \rangle |$|⁠. We consider our ‘default’ case above (⁠|$|\langle {\boldsymbol {B}} \rangle |\approx 0.5$|⁠, with mean Alfvén |$\mathcal {M}_{A}\sim 6$|⁠), a pure hydro case (⁠|$|\langle {\boldsymbol {B}} \rangle |=0$|⁠, |$\mathcal {M}_{A}=\infty$|⁠), and two strong-field cases with |$|\langle {\boldsymbol {B}} \rangle | = (4.3,\,27)$|⁠, producing saturated |$\mathcal {M}_{A} \approx (2.5,\,0.4)$|⁠, i.e. rms field strengths |$\langle |{\boldsymbol {B}}|^{2} \rangle ^{1/2} \sim (10,\,65)\,\mu G \,\sqrt{(\langle n_{\rm gas} \rangle /10\,{\rm cm^{-3}})\,(T/100\,{\rm K})}$|⁠.

Figs 8 and 9 shows the results. Most obviously, increasing the field strength suppresses gas-density fluctuations, especially at very low densities n_gas ≪ 〈n_gas〉 (where magnetic fields dominate the pressure); this effect is well known (see Ostriker, Stone & Gammie 2001; Lemaster & Stone 2009; Burkhart et al. 2009; Collins et al. 2012). Because the dust and gas decoupled at low densities, this can weakly reduce the absolute magnitude of dust concentration (e.g. the clumping factor) averaged over all densities in the simulation.

Figure 9.

Top: bivariate distribution of dust and gas densities as Fig. 3, for the |$\mathcal {M}\sim 10$|⁠, α = 0.03 simulations where we vary the magnetic mean-field strength (see Table 1). Bottom: distribution of dust-to-gas ratios δ as Fig. 6, for all gas (left) and just the high-density gas (right). Increased magnetic field strengths suppress the low-density tail of gas density fluctuations, where the dust and gas decouple, so the highest δ values arising from these low-density regions decrease. However, at a fixed density n_gas, increasing field strength (lower |$\mathcal {M}_{A}$|⁠) actually increases the dust-to-gas ratio fluctuations. The effect is weak compared to variations in α, but true at every n_gas we measure here. This owes to the increased solenoidal/vorticity component of the flows with stronger fields (which produces dust density fluctuations without corresponding gas density fluctuations) and to the existence of additional magnetic pressure to create local ‘pressure traps.’

However, at fixed density (n_gas), variations in the dust-to-gas ratio δ are stronger with larger field strengths. Stronger fields direct more energy into solenoidal as opposed to compressive modes (evident in the coherent filamentary structure and smaller voids in Fig. 9 for |$\mathcal {M}_{A}=0.4$|⁠); the solenoidal (vorticity/strain) modes can still induce large changes in dust density (see references in Section 1), but do not alter the gas density, so they directly alter the dust-to-gas ratio. Moreover, magnetic fields provide another source of pressure which the grains do not feel – this can produce phenomena such as zonal flows which create ‘pressure traps’ (local maxima) in which grains with appropriate stopping times collect (see Whipple 1972; Pinilla et al. 2012; Dittrich et al. 2013). The effects are weak compared to changing α, but are significant at every density n_gas.

We caution, however, that we have neglected Lorentz forces on grains, which of course become larger with increasing field strength, and may therefore reverse some of these effects.

3.2.4 Dust and gas correlation functions

The scale-dependence of grain clustering is also reflected in the autocorrelation function ξ(r), defined in the usual fashion as

\begin{eqnarray} \xi (r)_{\rm dust} &\equiv &\frac{1}{\langle n_{\rm dust} \rangle }\,\left\langle \frac{{\rm d} N_{\rm dust}({\boldsymbol {r}})}{{\rm d}^{3}{\boldsymbol {x}}} \right\rangle _{|{\boldsymbol {r}}|=r} - 1 \end{eqnarray}

(11)

where |$\langle {d N_{\rm dust}({\boldsymbol {r}})}/{d^{3}{\boldsymbol {x}}} \rangle$| is the average number density of dust particles at a distance |$r=|{\boldsymbol {r}}|$| from each dust particle.⁶ Replacing dust with gas, we have the autocorrelation of gas (mass), ξ(r)_gas. These are particularly simple to compute given our Lagrangian numerical method. In cases with a weak mean (coherent) magnetic field, there is no preferred direction in our simulations, so we need consider only the isotropic ξ(r) (even in the strong-mean field case we simulate, the anisotropic corrections are small compared to the effects of different grain sizes). Fig. 10 shows the correlation function for dust and gas in our |$\mathcal {M}\sim 10$|⁠, |$\mathcal {M}_{A}\sim 6$| simulations.

Figure 10.

The auto-correlation functions of the dust (ξ_dust) and gas (ξ_gas). Here we measure the three-dimensional, isotropic, radial ξ(r), per equation (11) (at single random instant in time); this relates to the radial distribution function g(r) = 1 + ξ(r). This measures the variance in density fluctuations (〈δ(r)〉²/〈δ〉²) smoothed on different scales. The gas ξ_gas flattens below a scale of order the sonic length |${\sim } R_{\rm sonic} \sim L_{\rm box}/\mathcal {M}^{2}$|⁠, as expected. The dust clustering is reflected by ξ_dust continuing to rise to much smaller scales. In all cases with α ≪ 1, ξ_dust ∝ r^−0.3–0.5 over the clustering dynamic range (down to the free-streaming scale L_stream at the maximum gas densities). For α ≳ 1 there is little clustering, except in the most extreme high-density regions (the ‘bump’ at small scales).

The gas correlation function for all our runs with the same |$\mathcal {M}$| (and |$\mathcal {M}_{A}$|⁠) is identical (since the gas dynamics are not altered by the dust) and behaves as expected. ξ_gas rises towards small scales, until it flattens completely at a scale r a factor of a few below the sonic scale |$R_{\rm sonic}\sim L_{\rm box}/\mathcal {M}^{2}$|⁠, where pressure effects suppress small-scale density fluctuations.

At large scales ≳ (a couple) α L_box ∼ 〈L_stream〉, the dust is well coupled to the gas (at least around the mean n_gas), so ξ_dust ∼ ξ_gas. But below this scale the two decouple. Initially (scales R_sonic ≲ r ≲ 〈L_stream〉), the dust clusters more weakly than the gas (ξ_dust < ξ_gas). This reflects the dust free-streaming. But dust clustering/trapping in dense gas means ξ_dust continues to rise as a power law towards much smaller scales |${\sim } L_{\rm stream}(n_{\rm gas}^{\rm max})$| (where |$n_{\rm gas}^{\rm max}$| represents the highest gas densities reached by a significant volume fraction, ∼100 〈n_gas〉 in these runs). We can approximate ξ_dust reasonably well over most of its dynamic range with a power law

\begin{eqnarray} \xi _{\rm dust} = \left(\frac{r}{r_{0}} \right)^{-\eta } \end{eqnarray}

(12)

with η ≈ 0.3–0.5 (depending on the simulation) and r₀ ∼ (0.3 − 1) L_box (note ξ must drop more rapidly and eventually become negative as r → L_box, but this is not interesting). For a power-law ξ(r), the variance in the density field averaged within a spherical radius r (or equivalently, the mass enclosed in spheres of radius r), is trivially related to ξ by

\begin{eqnarray} \sigma ^{2}_{n_{\rm dust}}(r) &\equiv & \frac{\langle n_{\rm dust}(r)^{2} \rangle }{\langle n_{\rm dust} \rangle ^{2}} = C_{\eta }\,\left( \frac{r}{r_{0}} \right)^{-\eta } \end{eqnarray}

(13)

where C_η ≈ 1.035 ∼ 1 for all η in the range of interest (Peebles 1993). So the rms dispersion in n_dust scales as |$\sigma _{n_{\rm dust}} \approx (r/L_{\rm box})^{-\eta /2} \sim (r/L_{\rm box})^{-0.2}$|⁠. The small value of η/2 ≈ 0.2 here means that the scale-dependence of the fluctuations is weak; this is why we saw relatively mild changes in the PDF as we increased h_min.⁷

Note, however, that the scaling ξ(r) determines how the volume-weighted, linear density variance scales (not logarithmic variance); moreover it measures this for n_dust not the dust-to-gas ratio (so some of the power comes from gas-density fluctuations). None the less since the power in gas-density fluctuations is small on small scales, this should translate there to dust-to-gas fluctuations. If we also assume a lognormal distribution, the linear and logarithmic variances are related by

\begin{eqnarray} \sigma ^{2}_{\ln (n_{\rm dust})}(r) &=& \ln {\left[ 1 + \sigma ^{2}_{n_{\rm dust}}(r) \right]} \sim \ln {\left[ 1 + \left( \frac{r}{r_{0}} \right)^{-\eta } \right]} \end{eqnarray}

(14)

\begin{eqnarray} \sigma _{\log _{10}(n_{\rm dust})} &=& \frac{\sigma _{\ln {(n_{\rm dust})}}}{\ln {10}}\sim 0.43\,\sqrt{\ln {\left[ 1 + \left( \frac{r}{r_{0}} \right)^{-\eta } \right]}} \end{eqnarray}

(15)

This provides a good approximation to how the PDF of δ scales versus h_min; however, from ξ alone we do not have the dependence on n_gas.

4 DISCUSSION

We have shown that aerodynamic particles (e.g. dust grains in neutral gas) exhibit large dust-to-gas variations, as well as structure and dynamics qualitatively different from the gas, in supersonic, MHD turbulence.

In some respects, this is similar to the well-studied sub-sonic case in protoplanetary discs. However, a key difference is that in the supersonic case, the gas density exhibits large fluctuations (and the gas-dust velocity contributes to the stopping time). This means that the ‘stopping time’ of the dust is no longer constant across the flow, even for a single dust species. We find that this actually enhances the dynamic range of scales which exhibit dust clustering, in contrast to the case of sub-sonic turbulence, where the dust-to-gas ratio fluctuations tend to be concentrated in a narrow range of scales around the ‘resonant’ scale where the eddy turnover time is about equal to the constant stopping time (see references in Section 1).

We show that fluctuations in the dust-to-gas ratio are approximately lognormal, with two regimes. (1) At low densities ρ_gas < α 〈ρ_gas〉, the grains decouple from the gas, so the dust scatters about its mean volume density independent of gas density changes. The nominal dust-to-gas ratios ρ_dust/ρ_gas in this limit can reach extremely large values, with a power law tail towards high ρ_dust/ρ_gas. However, the low-density regime is also the limit in which we expect Lorentz forces to begin dominating over drag forces, so the fluctuations may be suppressed. (2) At high densities, the dust and gas are partially coupled. The mean dust density follows the mean gas density; however, there are approximately lognormal fluctuations owing to non-linear grain clustering. Some of this resembles well-studied grain clustering in the sub-sonic limit, since the clustering scales of the dust can be below the sonic scale of the turbulence. But there are additional effects as well, for example, grains can sediment into very thin filaments within shock fronts, similar dynamically to sedimentation under gravitational forces but here the effective acceleration owes to the pressure forces felt by gas and not dust. The magnitude of these fluctuations is large, ∼0.3–0.5 dex 1σ dispersion, for grains with size parameter over a wide range α ∼ 0.001–0.3, with a maximum around α ∼ 0.01–0.1. Much larger grains (α ≳ 1) are never tightly coupled; much smaller grains (α ≲ 0.001) are too well coupled to gas. The characteristic spatial scales of the grain structures/clustering increase with the grain size (see equation 7).

These clustering effects can have many important consequences, which we will explore in future work. For example:

• Dust formation and growth: because dust is highly clustered, its growth and evolution, particularly via dust–dust collisions (coagulation or shattering) can be dramatically altered. To lowest order, these effects are manifest in the clumping factor |$\langle n_{\rm dust}^{2} \rangle / \langle n_{\rm dust} \rangle ^{2}$|⁠, which governs the dust–dust interaction rate and reaches values ∼50–100. The effects of grain clustering on growth have been extensively studied in protoplanetary discs; however they are not well understood in molecular clouds. At the very least, the large clumping factors imply order-of-magnitude faster evolution of large grains in neutral clouds compared to what is usually assumed. Other effects, for example the non-uniform and size-dependent velocity dispersions of grains, may substantially alter both collision rates and the outcomes of those collisions (sticking versus shattering, for example). These effects, in turn, may dramatically influence the size distribution of dust.

• Extinction mapping and dust emission: visually, it is obvious that the dust and gas are not necessarily co-located. In probes of extinction and dust emission, this may be directly visible; however, we caution that the predictions here correspond to large dust grains. These do not dominate extinction. Rather, one would have to use diagnostics specifically sensitive to large grains (for example, sub-mm observations). Moreover, there is a finite scale which must be resolved in order to see the dust-to-gas fluctuations. On larger scales compared to the critical scale for dust clustering, they will be smoothed out, and one will simply trace the mean dust-to-gas ratio. But such fluctuations, on scales similar to the critical scale predicted, have been observed in many nearby clouds and some centres of nearby galaxies including e.g. Taurus (Padoan et al. 2006; Flagey et al. 2009; Pineda et al. 2010), NGC 1266 (Nyland et al. 2013; Pellegrini et al. 2013), Orion (Abergel et al. 2002), the Ursa Major cirrus (Miville-Deschênes et al. 2002), IC 5146, CGCG525-46, IR04139+2737, and G0858+723 (Thoraval, Boisse & Duvert 1997; Thoraval, Boissé & Duvert 1999). The absolute scales where the fluctuations are observed range from ∼0.003to10 pc, but in each case the critical scale and magnitude of fluctuations appears to agree with the simple scalings expected for turbulent concentration, given the different cloud densities and grain sizes probed (for detailed comparisons, see Padoan et al. 2006 and Hopkins 2014b). Thus, great care is needed, especially as observations push to higher resolution. Both dust and gas have a filamentary morphology, but dust filaments may, in fact, be much narrower than gas filaments (which are characteristically of order the sonic length); in rare cases, dust filaments can exist where there is no gas filament at all (owing to dust concentration by gas vorticity). Very large dust grains ( ≫ 1 μm), on the other hand, may be more uniformly distributed than gas throughout clouds. This may resolve several long-standing puzzles regarding apparently different extinction measurements that have alternatively been attributed to different dust chemistry in different regions.

• Cooling physics and star formation: in dense regions of clouds or galactic nuclei, gas cooling or heating can be regulated by collisions with dust, with the relevant rate proportional to the local dust-to-gas ratio. When this dominates, we therefore predict that there may be order-of-magnitude variations in the cooling physics of some regions. In metal-poor galaxies, regions which are relatively overabundant in dust may be preferentially able to form stars, since low-mass star formation may be difficult without sufficient dust present to act as a coolant.

• Stellar abundances: as proposed in Hopkins (2014b), large fluctuations in the local dust-to-gas ratio should translate to interesting variations in stellar abundances, even for stars formed in the same cluster. Large dust grains contain most of the dust mass (about ∼1/2 the total metal mass), and they are the ones for which these fluctuations are important. Even smoothing on relatively large scales, corresponding to ≳ 0.1 pc (the size of large protostellar cores), we predict significant fluctuations if grains have the appropriate sizes. Specifically, assuming L_box ∼ 10 pc and Σ_GMC ∼ 300 M_⊙ pc⁻², and that 1/3 the metals are in grains with sizes ∼0.1 μ m, we predict an approximately ≈0.05–0.1 dex 1σ dispersion in the total metallicity of the dense regions (owing to dust-to-gas fluctuations); this is small and well within the dispersion observed for nearby clusters (Casagrande et al. 2011; Duran et al. 2013). More interestingly, though, because it is lognormal, the distribution has a long tail, and one dense star-forming region per million could have a total metallicity enhancement of a factor ∼20–50!

Studying these in more detail requires additional simulations with the relevant physics included, which makes our calculations no longer scale-free. However, it is straightforward to extend our models and follow these additional processes. Moreover, applying these simulations to a specific scale and situation allows for the inclusion of additional, non-scale-free physics which we have ignored in this first study (such as Lorentz forces and grain collisions).

We thank Jim Stone, Paolo Padoan, Jessie Christiansen, Charlie Conroy, Evan Kirby, and Paul Torrey for helpful discussions and contributions motivating this work. We also thank the anonymous referee for valuable suggestions. Support for PFH was provided by the Gordon and Betty Moore Foundation through grant no. 776 to the Caltech Moore Center for Theoretical Cosmology and Physics, an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant no. 1411920. Numerical calculations were run on the Caltech compute cluster ‘Zwicky’ (NSF MRI award no. PHY-0960291) and allocation TG-AST130039 granted by the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the NSF.

1

A public version of this code is available at http://www.tapir.caltech.edu/∼phopkins/Site/GIZMO.html.

2

Equation (1) assumes grains are in the Stokes, rather than the Epstein, limit for drag, i.e. a_d < (9/4) λ, where λ = (n_gas σ_gas)⁻¹ is the gas mean free path. For molecular hydrogen cross-sections at ∼30 K, this requires a_d < 10¹³ cm (n_gas/10 cm⁻³), so is obviously satisfied.

3

We use a standard kernel-density estimator to calculate the density of dust particles and gas particles around each point |${\boldsymbol {x}}_{i}$| based on the distribution of dust/gas particle neighbours, e.g. |$n_{{\rm dust},\,i} = \sum _{j} W({\boldsymbol {x}}_{j}-{\boldsymbol {x}}_{i},\,h_{i})$|⁠, where h = MAX(h_min, h_N) with h_N the radius that encloses a finite neighbour number of particles j, and W the kernel function chosen here for consistency to match the same used in our mesh-free hydrodynamics methods. This is much more accurate, given the Lagrangian nature of our code, than a simpler particle-in-cell estimate. We have confirmed that our results are insensitive to the (arbitrary) neighbour number in the estimation kernel and the specific choice of kernel function. We have also confirmed that a direct reconstruction output by our hydrodynamic solver gives indistinguishable results to this estimator applied in post-processing.

4

For dP_vol, we have equation (10) with 〈ln (δ)〉_vol = 〈ln δ(n_gas)〉_dust − S_dust. For dP_gas, we apply the same shift in 〈ln (δ)〉 and also shift the mean 〈ln (n_gas/〈n_gas〉)〉_gas = +S_gas/2. See Appendix A for details.

5

Non-linear interactions make it difficult to maintain exactly |$\mathcal {M}=10$| as we vary |$|\langle {\boldsymbol {B}} \rangle |$|⁠; however the runs in Table 1 all saturate with |$\mathcal {M}\approx 8{\rm -}12$|⁠.

6

Note in the terrestrial literature, it is more common to define the radial distribution function g(r), but this is trivially related to ξ(r) by g(r) = 1 + ξ(r).

7

Note that ξ(r) is three dimensional. The projected/angular correlation function ω(θ) is simply related by integration. Because ξ(r) ∝ r^−η is sufficiently shallow (η < 1), we obtain a nearly flat projected ω(θ) ∼ constant.

8

We note that the errors in the ‘tracer particle limit’ α = 10⁻¹⁰ here are significantly smaller than those shown in Genel et al. (2013). This owes to several effects: we use a smooth (as opposed to discontinuous) interpolation for the velocity field, our method is fully Lagrangian (there are no interparticle mass fluxes to enhance discrepancies in advection), we use the exact solution over the timestep to update particle velocities (as opposed to only the instantaneous acceleration), and we synchronize the time updates between gas and dust (and use a stricter dust timestep criterion). As discussed therein, these all reduce (although do not completely eliminate) the relevant errors.

APPENDIX A: DISTRIBUTIONS PER UNIT VOLUME AND PER UNIT GAS MASS

In the text, we have (unless otherwise specified) shown the PDF of field properties such as n_dust and n_gasaround a random dust grain. This is the dust-mass weighted PDF dP_dust:

\begin{eqnarray} {\rm d}P_{\rm dust} = {\rm Probability(around\ random\ grain)} \end{eqnarray}

(A1)

Because, within an infinitesimally small differential volume |$d^{3}{\boldsymbol {x}} = {\rm d}V$|⁠, the dust properties are sampled N = n_dust dV times, it is trivial to show that this is related to the PDF around a random point in space|${\boldsymbol {x}}$|⁠, i.e. the volume-weighted PDF dP_vol, by

\begin{eqnarray} {\rm d}P_{\rm vol} \propto n_{\rm dust}^{-1}\,{\rm d}P_{\rm dust} \end{eqnarray}

(A2)

Note that this trivial conversion is only strictly true if we measure the point-like density n_dust – if we instead smooth the properties on some finite averaging scale h_min > 0, this is only approximate. However, we can, of course, rigorously calculate dP_Vol by following the same measurement procedure described in the text, beginning from randomly selected, uniformly distributed points in space, as opposed to the locations of dust grains. Finally, it is likewise trivial to show that the probability around a random gas element, i.e. the Lagrangian or gas-mass weighted PDF dP_gas, is given by

\begin{eqnarray} {\rm d}P_{\rm gas} \propto n_{\rm gas}\,{\rm d}P_{\rm vol} \end{eqnarray}

(A3)

Figs A1 and A2 show the full bivariate distribution of n_dust and n_gas as Fig. 3 in the text. However, Fig. A1 shows the volume-weighted probability P_vol and Fig. A2 shows the gas-mass weighted probability P_gas. These are calculated correctly for h_min > 0 but are very close to the approximate values given by the simple relations above.

Figure A1.

Distribution of dust and gas densities in our |$\mathcal {M}\sim 10$| simulations, as Fig. 3 in the text. The only difference is that here, we measure the PDF around random points in space within the box (i.e. measure the volume-weighted PDF dP_vol), instead of the distribution around random dust particles (dP_dust). For h_min = 0, the two are trivially related by |${\rm d}P_{\rm vol} \propto n_{\rm dust}^{-1}\,{\rm d}P_{\rm dust}$|⁠. The contours shift to lower n_dust and n_gas as these have larger volume-filling factors, and the peak of the volume-averaged (h_min > 0) probability density shifts closer to n_gas = 〈n_gas〉, n_dust = 〈n_dust〉, as it must, but the qualitative behaviour and scatter in n_dust at fixed n_gas is similar in all cases.

Figure A2.

Distribution of dust and gas densities in our |$\mathcal {M}\sim 10$| simulations, as Fig. 3 in the text. As Fig. A1, the difference here is that we measure the PDF around random gas elements (i.e. the gas-mass weighted PDF dP_gas) instead of around random dust elements (dP_dust) or random volume elements (dP_vol). For h_min = 0 these are related by dP_gas ∝ n_gas dP_vol ∝ (n_gas/n_dust) dP_dust. Overall, the contours shift noticeably towards the mean dust-to-gas ratio (especially at low densities), i.e. most of the gas sees a ratio closer to the mean, compared to what most of the dust sees, because the dust is more highly clustered than the gas. The scatter about this mean at high n_gas, however, is similar in dP_gas and dP_dust.

First consider dP_vol. As expected, the distribution shifts to lower n_dust and n_gas as these regions have larger volume-filling factors. Similarly the peak of the probability density shifts closer to n_gas = 〈n_gas〉, n_dust = 〈n_dust〉, and rapidly moves on to this point as we increase the volume-averaging scale h_min (since the distribution function for even infinitely clustered dust must converge to a delta function around this point as h_min → L_box). Modulo the mildly reduced scatter towards high n_dust, however, the qualitative behaviour of the distribution functions is similar to dP_dust in all cases.

Next consider dP_gas. Interestingly, in this case, the distribution shifts significantly towards the mean dust-to-gas ratio, especially at low densities. This is partly because the low-n_gas regions (where the dust decouples) are downweighted in the distribution. It cannot entirely be explained by this effect however – even at intermediate/high densities, most of the gas has a local dust-to-gas ratio closer to the mean compared to most of the dust. As noted in the main text, this is expected because the dust is more highly clustered than the gas. But once again, the high-n_gas behaviour and scatter is still similar to dP_dust.

A simple model explains these results. In supersonic turbulence, the gas density is approximately lognormal, with dP_vol ∝ exp [ − (ln n_gas − S_gas/2)²/(2 S_gas)] where |$S_{\rm gas}\approx \ln {[1+b^{2}\,\mathcal {M}^{2}]}$| is the variance (b ∼ 1/2, depending on details of the turbulence). Since the dust does not modify the gas in our runs, the bivariate distribution should reflect this for the gas, with P(n_dust | n_gas) also approximately lognormal, as shown in Fig. 6. But for a lognormal distribution dP(ln x), it is trivial to show that the distribution dP_new ∝ x dP is also lognormal, with the same variance but a shifted mean. Therefore, in this case, we expect the bivariate distributions to be approximately given by

\begin{eqnarray} \frac{{\rm d}P_{\rm vol}}{{\rm d}\ln {n_{\rm gas}}\,{\rm d}\ln {n_{\rm dust}}} &=& \frac{1}{2\pi \sqrt{S_{\rm gas}\,S_{\rm dust}}}\, \exp {\left[ -\frac{\Delta _{\rm gas}^{2}}{2\,S_{\rm gas}} - \frac{\tilde{\Delta }_{\rm dust}^{2}}{2\,S_{\rm dust}} \right]} \end{eqnarray}

(A4)

\begin{eqnarray} \frac{{\rm d}P_{\rm gas}}{{\rm d}\ln {n_{\rm gas}}\,{\rm d}\ln {n_{\rm dust}}} &=& \frac{1}{2\pi \sqrt{S_{\rm gas}\,S_{\rm dust}}}\, \exp {\left[ -\frac{\tilde{\Delta }_{\rm gas}^{2}}{2\,S_{\rm gas}} - \frac{\tilde{\Delta }_{\rm dust}^{2}}{2\,S_{\rm dust}} \right]} \end{eqnarray}

(A5)

\begin{eqnarray} \frac{{\rm d}P_{\rm dust}}{{\rm d}\ln {n_{\rm gas}}\,{\rm d}\ln {n_{\rm dust}}} &=& \frac{1}{2\pi \sqrt{S_{\rm gas}\,S_{\rm dust}}}\, \exp {\left[ -\frac{\Delta _{\rm gas}^{2}}{2\,S_{\rm gas}} - \frac{{\Delta }_{\rm dust}^{2}}{2\,S_{\rm dust}} \right]} \end{eqnarray}

(A6)

\begin{eqnarray} \Delta _{\rm gas} &\equiv & \ln {\left(\frac{n_{\rm gas}}{\langle n_{\rm gas}\rangle }\right)} + \frac{S_{\rm gas}}{2} \end{eqnarray}

(A7)

\begin{eqnarray} \tilde{\Delta }_{\rm gas} &\equiv & {\Delta }_{\rm gas}-S_{\rm gas} = \ln {\left(\frac{n_{\rm gas}}{\langle n_{\rm gas}\rangle }\right)} - \frac{S_{\rm gas}}{2} \end{eqnarray}

(A8)

\begin{eqnarray} {\Delta }_{\rm dust} & \equiv & \ln {(\delta )} - \langle \ln {\delta (n_{\rm gas})} \rangle _{\rm dust} \end{eqnarray}

(A9)

\begin{eqnarray} \tilde{\Delta }_{\rm dust} &\equiv \Delta _{\rm dust} + S_{\rm dust} = \ln {(\delta )} - \delta _{0} \end{eqnarray}

(A10)

where S_gas is constant but S_dust = S_dust(n_gas) and δ₀ = δ₀(n_gas) are functions of n_gas, and

\begin{eqnarray} \delta &\equiv &\frac{n_{\rm dust}/n_{\rm gas}}{\langle n_{\rm dust} \rangle / \langle n_{\rm gas} \rangle } \end{eqnarray}

(A11)

\begin{eqnarray} \delta _{0} &\equiv \langle \ln \delta (n_{\rm gas})\rangle _{\rm vol} = \langle \ln \delta (n_{\rm gas})\rangle _{\rm dust}-S_{\rm dust} \end{eqnarray}

(A12)

It is easy to verify these obey dP_gas ∝ n_gas dP_vol ∝ (n_gas/n_dust) dP_dust. Trivially, therefore, if the distributions are lognormal in δ and n_gas, the logarithmic scatter is identical regardless of how we weight the distributions, and the mean values simply shift.

The values of S_dust and δ₀ can be read off of Fig. 7, noting S_dust = [σ_log 10(δ) ln 10]² ∼ 0.5–2, and that the plotted 〈log₁₀δ(n_gas)〉_dust = (δ₀ + S_dust)/ln (10). If δ₀ and S_dust are constant (approximately true in the high-density limit), then the constraint that the PDF integrates correctly gives δ₀ = −S_dust/2; if we instead assume S_dust is constant but δ₀ = A − ln (n_gas/〈n_gas〉) (i.e. n_dust ∼ constant, approximately true in the low-density limit), we have A = −S_dust/2. These give good approximations in both limits to the results in Fig. 7.

APPENDIX B: RESOLUTION STUDY

Obviously it is important to test that our simulations are numerically converged. Because the resolution we can achieve is more limited in 3D, we consider a resolution study first in 2D reaching much higher resolution than our standard runs in the text.

Fig. B1 plots the time-averaged PDF of the dust-to-gas ratio, in the style of Fig. 5 in the text, for 2D simulations with Mach number |$\mathcal {M}_{2D}\sim 5$| and α = 0.01. We consider resolutions 64²–1024². As expected, the tails of the PDF are better sampled at higher resolution – this follows from simple counting statistics. Remarkably, the core of the PDF appears reasonably well-converged at just ∼64² resolution; by ∼256² the ‘wings’ agree well down to part-per-million amplitudes (there is a small deviation in the 512² run, such that the 1024² run actually agrees slightly better with the 256² run; this appears to depend on how the turbulent driving routines depend on resolution). This justifies our choice in the text of 256³ resolution.

Figure B1.

Top: distribution of dust-to-gas ratios (as Fig. 5), in a 2D resolution study with |$\mathcal {M}_{2{\rm D}}\sim 5$| and α = 0.01. Bottom: Same, in a 3D study with |$\mathcal {M}\sim 10$| and α = 0.03. Owing to the Lagrangian nature of our code, and to the fact that the turbulence is supersonic (and so structures are driven by relatively easily captured shocks and rarefactions), the convergence is remarkably good. Even 64^D runs appear well converged in the core of the PDF; by 256^D the results agree well with our 1024^D simulations (in 2D). We expect our conclusions in the text are robust to resolution effects.

Interestingly, the convergence here is much faster than seen in some previous studies (compare e.g. Bai & Stone 2010a). This owes in part to the Lagrangian nature of our method, which is able, in principle, to capture arbitrarily large fluctuations in density (so long as they involve equal to or larger than some fixed mass scale) at low ‘resolution’ (i.e. there is no fixed spatial resolution). It also owes to the supersonic nature of the turbulence here, where much of the dynamics is driven by shocks and rarefactions (relatively ‘easily’ captured in these methods), as opposed to the streaming instability or details of the vorticity field of small turbulent eddies (the dominant effects in the highly subsonic limit).

We have also considered a limited study in 3D, taking one of our standard runs and re-running at lower resolution. Even at 64³, our qualitative conclusions are essentially identical (although the extremes of the distribution functions are sampled relatively poorly).

APPENDIX C: NUMERICAL EFFECTS OF THE DUST DRAG ALGORITHM AND POISSON ERRORS

Although the scheme used here to integrate the trajectories of dust particles is standard and relatively well tested (similar to Carballido et al. 2008; Hogan et al. 1999; Johansen & Youdin 2007; Johansen et al. 2009; Bai & Stone 2010a; Pan et al. 2011), there are known sources of numerical error.

The advantage of a Lagrangian ‘superparticle’ approach is that, in the limit where the grains are decoupled from the gas (α → ∞), their dynamics (free-streaming) are perfectly recovered, and the only source of error in the density field is Poisson noise (from our finite particle number). This is not true in ‘two-fluid’ approximations, for example, which cannot account for the full velocity distribution function of grains at a single location.

In the opposite limit of perfect coupling (α → 0), the grains should perfectly trace the gas (as tracer particles), up to Poisson noise in the initial tracer field. However, our methods introduce an additional error: when α → 0, the algorithm used to update the particle velocities and positions (interpolating to the particle position) does not, numerically, perfectly match the Godunov-type update to the gas particle velocities (involving the solution of a Riemann problem). In a sufficiently smooth flow, these should be identical, but given numerical noise or physical discontinuities, they can differ (for detailed analysis of these errors, see Genel et al. 2013).

We therefore test both limits here. We take our standard |$\mathcal {M}=10$| simulation and re-run with α = 10¹⁰ (effectively infinite) and α = 10⁻¹⁰ (effectively zero). In Fig. C1, we plot the resulting images, bivariate density distributions, and time-averaged PDF of the dust density (for α = 10¹⁰) and dust-to-gas ratio (for α = 10⁻¹⁰). We take h_min = 0, since the errors of interest rapidly become smaller as the averaging scale becomes larger. For α = 10¹⁰, we confirm that the scatter in dust density is what we expect from Poisson statistics (with smaller residual errors owing to our post-processing kernel density estimator). For α = 10⁻¹⁰, we find the dust traces gas at all densities, with a comparable scatter to the Poisson case.

Figure C1.

Top: images of the dust and gas density as Fig. 2, in a |$\mathcal {M}\sim 10$| simulation with almost perfectly coupled α = 10⁻¹⁰ ( left) and almost perfectly uncoupled α = 10¹⁰ (right). Middle: bivariate dust and gas distribution as Fig. 3 for both cases. We show h_min = 0; the scatter decreases for larger h_min. Bottom: histogram of the dust-to-gas ratio δ (left) and dust density n_dust (right), as Fig. 5. In the perfectly coupled case, dust should track gas exactly (δ = 1), in the un-coupled case, dust should remain at the mean density (n_dust = 〈n_dust〉). Poisson sampling from our finite particle number (laid down randomly in the initial conditions) leads to some scatter. For the strongly coupled case, these errors are enhanced by small numerical differences between the algorithms used to update the gas and dust velocities. However, the distributions do not show any systematic deviation from the expected behaviour. Their widths (σ ∼ 0.05–0.07 dex) are much smaller than any α ∼ 0.001–1 case we consider in the text, so the errors are not significant in our study.

In both cases, the scatter in the core of the distribution is <0.1 dex; much smaller than we see in any of our simulations with 0.001 ≲ α ≲ 1. Moreover, the tails of the distribution are dramatically suppressed – these are many orders of magnitude smaller than we see in the text. And in both cases, the mean dust density behaves as it should and we see no unphysical features (only noise).⁸

We conclude that these sources of error are not significant for the α values in the text. Based crudely on the scaling of the variance in Fig. 7, we estimate that Poisson noise and/or integration errors would, at our current resolution, become significant compared to physical effects at α ≪ 10⁻⁴ or α ≫ 100, necessitating higher resolution studies.