© The Authors 2023

1 Introduction

The properties of dust grains regulate the opacity, ionisation, initial conditions of planet formation in star forming regions, but their evolution remains poorly constrained. Contemporary observing facilities on the ground and in space have shed light on dust-grain properties at various scales from the diffuse interstellar medium (ISM) using Planck (Planck Collaboration XXIV 2011) down to protoplanetary discs using ALMA and SPHERE (e.g. Muro-Arena et al. 2018). It is now understood that the dust-grain population evolves at all scales via dynamical interaction with the gas (drag force) and the magnetic fields (Lorentz force), as well as via interaction within the dust-grains population that leads to growth and fragmentation processes (e.g. Lesur et al. 2022). The study of the evolution of dust grains during star and planet formation has therefore attracted significant interest in recent years. There is growing evidence from observations and theoretical works that dust-grain properties evolve dramatically within dense cores during the very early stages of star and planet formation, suggesting rapid dust growth and dynamical segregation (e.g. Steinacker et al. 2010; Bate & Lorén-Aguilar 2017; Sadavoy et al. 2018; Galametz et al. 2019; Lebreuilly et al. 2020; Guillet et al. 2020; Tsukamoto et al. 2021). It remains unclear as to whether or not dust grains also evolve at larger scales within molecular clouds, prior to protostellar collapse.

Thanks to major advances achieved in computational astrophysics for modelling dust and gas mixtures as well as in hardware development, it is now possible to study the dynamical evolution of dust grains at various scales in the ISM with a number of astrophysical codes (see Teyssier & Commerçon 2019, for a review). In particular, recent works focused on the dynamics of dust grains within molecular clouds and report variation of the dust concentration (Hopkins & Lee 2016; Tricco et al. 2017; Mattsson et al. 2019a). This result is of significant importance because it can have different astrophysical implications. First, the dust-to-gas ratio could vary in molecular clouds beyond the canonical value of 1 percent. Second, dynamical decoupling of grains of different sizes could favour dust growth and fragmentation processes. Variations in dust-grain concentration and size distribution will then have potential implications on both heating and cooling, as well as on mass estimates from observations through opacity variations (Ormel et al. 2011). Furthermore, the non-ideal magnetohydrodynamics resistivities within collapsing dense cores, which regulate the formation of protoplanetary discs, are highly sensitive to dust-grain size distribution (Guillet et al. 2020; Marchand et al. 2021).

However, the critical dust-grain size for measurements of dynamical decoupling varies between these latter studies. Hopkins & Lee (2016) report significant dust decoupling (variation in dust concentration by more than a factor 1000) for dust grains of s_grain ≥ 0.01 μm in size. Similarly, Mattsson et al. (2019a) report dust-ratio variations for sizes s_grain > 1 μm. On the contrary, Tricco et al. (2017) report significant decoupling only for dust grains larger than 10 μm. In addition, Mattsson et al. (2019b) find that very small grains (typically nanometer size) are also clustering, which increases the local grain density by at least a factor of a few in their models.

At first glance, the difference in the critical size for dynamical decoupling might seem insignificant, but it can have dramatic effects on the evolution of dust-size distribution in the ISM. The dust-grain size distribution in the diffuse ISM is relatively well constrained, with maximum sizes of <1 μm (Guillet et al. 2018) and a peak of the distribution at around 0.1−0.3 μm in the Milky Way (Weingartner & Draine 2001; Draine 2003; Guillet et al. 2018). This peak lies just in the range where Hopkins & Lee (2016) and Mattsson et al. (2019a) report decoupling while Tricco et al. (2017) do not.

The cause of the discrepancy in the literature for the critical size for decoupling remains unknown. Both the physical and numerical setup are different from one study to another. In the following, we investigate whether these discrepancies may arise from differences in numerical implementation. We briefly present the current state of the art in the numerical methods used, as well as their main limitations, which we will test in this study.

Hopkins & Lee (2016) and Mattsson et al. (2019a) use a two-fluid model, where the dust and the gas are considered as two separate fluids interacting via a drag force. These authors also use two different fluid descriptions. In Mattsson et al. (2019a), the gas is handled on a uniform Cartesian grid (Eulerian approach), while the dust fluid is handled using inertial particles (Lagrangian approach). Hopkins & Lee (2016) use a mesh-free method where the gas and the dust fluids are handled using two separate particle distributions (Hopkins 2015). In both implementations, the back-reaction of the dust onto the gas is not considered.

In addition, very small grains, which are well coupled to the gas, can be assimilated to gas tracer particles; that is, they mostly move with the gas velocity (with a tiny velocity shift). Price & Federrath (2010) and Cadiou et al. (2019) showed that classical particle-integration schemes based on velocity increments interpolated from the grid (using for instance a Cloud-In-Cell (CIC) algorithm) cannot properly trace the gas dynamics. These authors show that the velocity tracers aggregate in converging flow (high-density region) leading to artificial clustering. The density in the vicinity of converging regions can therefore be overestimated by one order of magnitude, while it is largely underestimated around filaments. We test whether this result affects the dynamics of very small dust grains treated as Lagrangian particles by comparing their distribution to the distribution of dust velocity tracer particles.

On the other hand, Laibe & Price (2014a,b) propose a full monofluid approach, which is well adapted for dust grain dynamics as long as the fluid approximation remains valid (Stokes number (St) of the order unity). The full monofluid system of equations for the gas and dust dynamics and the system of equations of the two-fluid formalism are mathematically equivalent. The monofluid approach consists of a change of variables, which allows us to follow the evolution of the barycentre of the mixture and the relative velocity difference and concentration of the dust. We have further simplified this formalism using the so-called diffusion approximation (Price & Laibe 2015), which is based on the terminal velocity approximation (TVA, Youdin & Goodman 2005). In this approximated monofluid formalism, the drift velocity is set directly by the force budget on the gas and the dust, and only one equation (the mass conservation) needs to be solved per dust size. In the context of molecular clouds, this diffusion approximation has been used in Tricco et al. (2017) and is well suited for well-coupled grains with St < 1; it is not valid for larger grains, because the velocity difference between the dust and the gas is underestimated.

In this study, we perform a comprehensive test of well-controlled numerical experiments which allows the use of both the monofluid formalism in the diffusion approximation and a two-fluid method based on inertial Lagrangian particles for the dust. We focus on the effect of the numerical implementations on the dust dynamics in driven dusty turbulence experiments. We will explore the effect of the physical parameters (amplitude and properties of turbulence and the effect of magnetic fields) in a forthcoming study. The paper is organised as follows. In Sect. 2, we estimate the expected critical dust grain size for decoupling within typical turbulent molecular clouds from simple analytical arguments. The numerical methods and the physical setup are described in Sect. 3. Section 4 is devoted to the results obtained with the monofluid implementation. We then compare with the two-fluid results in Sect. 5. In Sect. 6, we discuss the limitation of this work and the comparison with previous works.

2 Condition for decoupling in molecular clouds

2.1 Turbulence in molecular clouds

In a turbulent cloud of size L and internal velocity dispersion υ_rms, the dynamical time is defined as the time at which a turbulent fluctuation travels through the blob:

$$t_{\text{turb}}\equiv \frac{L}{{\upsilon }_{\text{rms}}},$$(1)

where t_turb can been seen as the time needed for the turbulence to completely renew all the hydrodynamic fields of the flow. In molecular clouds, observations show that the amplitude of the turbulence υ_rms and the gas density n scale with the size L as (Larson 1981; Heyer & Brunt 2004; Roman-Duval et al. 2011; Hennebelle & Falgarone 2012)

$$\begin{array}{c}{\upsilon }_{\text{rms}}\simeq 1{\text{kms}}^{-1}{\left(\frac{L}{1\text{Pc}}\right)}^{0.5},\\ n\simeq 3000{\text{cm}}^{-3}\left(\frac{L}{1\text{Pc}}\right).\end{array}$$(2)

We aim to study dust dynamics in typical conditions that satisfy the two aforementioned scaling relations.

2.2 Dust grain and gas dynamical (de)coupling

Dust grains experience a drag force from the collisions with gas particles. The drag force experienced by a dust grain can be written as

$$F_g=\frac{m_{\text{grain}}}{t_{\text{s,grain}}}\Delta v,$$(3)

where Δν = v_g − v_d is the differential velocity between the dust grain and the gas (locally described by a Maxwellian velocity distribution), and t_s,_grain is the grain stopping time, which characterises the time needed for the dust grain to adjust its velocity to a change of gas velocity. For typical physical conditions of molecular clouds, the mean free path of the gas is much larger than the size of spherical dust grains; this is the Epstein regime (Epstein 1924). In this regime, the stopping time is defined as

$$t_{\text{s,grain,}0}=\sqrt{\frac{\pi \gamma }{8}}\frac{{\rho }_{\text{grain}}}{{\rho }_g}\frac{s_{\text{grain}}}{c_s},$$(4)

where ρ_grain is the grain intrinsic density (typically 1–3 g cm⁻³, Love et al. 1994), s_grain is the grain size, γ the heat specific ratio of the gas, ρ_g the gas density, and c_s the gas sound speed. In molecular clouds, the flow is supersonic, as described above. As a consequence, the velocity drift between the dust and the gas can exceed the sound speed locally. In that case, we apply a correction for supersonic flow to handle this regime and the stopping time is defined as (Kwok 1975; Draine & Salpeter 1979)

$$t_{\text{s,grain}}=t_{\text{s,grain,0}}{\left(1+\frac{9\pi }{128}\frac{\Delta v^2}{c_s^2}\right)}^{-1/2}.$$(5)

Assuming that all grains have the same intrinsic density, the larger the grain, the longer the stopping time, and consequently, the smaller grains are the most coupled to the gas. The degree of coupling between the dust grain and the gas is characterised by the Stokes number, which is defined as

$$S\text{t=}\frac{t_s}{t_{\text{dyn}}},$$(6)

where t_dyn is the typical dynamical time of the system and $t_s\equiv \frac{{\rho }_{\text{g}}}{\rho }{t}_{\text{s,grain}}$ is the fluid stopping time of the dust species (hereafter stopping time), ρ being the total density of the gas and the dust. For St ≪ 1, the gas and the dust are dynamically coupled.

2.3 Critical grain size for dynamical decoupling

In turbulent molecular clouds, the dynamical time is the turbulent time t_turb. One can consider that a dust grain decouples from the gas dynamics if St > 1, which yields the critical dust grain size (ignoring the correction for supersonic flow in the stopping time expression):

$$s_{\text{crit}}=\sqrt{\frac{8}{\pi \gamma }}\frac{{\rho }_{g}L}{{\rho }_{\text{grain}}ℳ}.$$(7)

According to the Larson relations (Eq. (2)), the critical dust grain size scales as

$$s_{\text{crit}}\sim \text{μm}\left(\frac{{\rho }_g}{{10}^{-20}\text{g\hspace{0.17em}}{\text{cm}}^{-3}}\right)\left(\frac{\text{L}}{1\text{\hspace{0.17em}pc}}\right){\left(\frac{{\rho }_{\text{grain}}}{1\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}}\right)}^{-1}{\left(\frac{ℳ}{10}\right)}^{-1},$$(8)

where ℳ = υ_rms/c_s is the turbulent Mach number, and we have taken γ = 5/3. Dust grains of a few 10s microns are therefore expected to decouple significantly from the gas dynamics within a crossing time t_turb. Numerical experiments show that the decoupling is already significant for Stokes numbers St ≥ 0.1 (Dipierro et al. 2015; Lebreuilly et al. 2020), which indicates that dust grains >1 μm are expected to decouple on a dynamical timescale as well.

This value of s_crit is consistent with the maximum dust-grain sizes expected in molecular clouds (<1 μm e.g. Köhler et al. 2015). Recent works have shown that grain growth could be at play in the dense and turbulent ISM (Ormel et al. 2009; Guillet et al. 2020). In this study, we therefore consider grains up to sizes of 10 μm.

3 Numerical methods and initial conditions

3.1 The RAMSES code

We use the adaptive-mesh-refinement (AMR) code RAMSES (Teyssier 2002) which integrates the Euler equations in their conservative form using a second-order finite volume Godunov scheme. In this work, we do not take into account magnetic fields or gravity. In addition, we use a uniform grid, which is well suited to the study of turbulent boxes.

Our setup follows the classical one designed for investigating compressible turbulence in molecular clouds using turbulent boxes. The boundary conditions are periodic for all hydrodynam-ical quantities. The gas thermal evolution is isothermal. We use the HLL Riemann solver combined with a minmod slope limiter for the hyperbolic part and upwind for the dust solver (see e.g. Toro 1999, for more details). In our experience, this combination of solver and slope limiter is the best compromise between accuracy and robustness.

We employ the same module to drive the turbulence as that used by Commerçon et al. (2019). This latter is based on a source term, a force, in the momentum equation. The turbulent force field is generated in Fourier space and its mode is given by a stochastic method based on the Ornstein-Uhlenbeck process (Eswaran & Pope 1988; Schmidt et al. 2006). For more details on the turbulence forcing module, we refer readers to the work of Schmidt et al. (2009) and Federrath et al. (2010) on which our implementation is based. We set the forcing to a mixture of compressive and solenoidal modes, with a compressive force power equal to one-third of the total forcing power. The turbulence is driven on the scale of the computational box, with a peak at half the box length.

3.2 Dust and gas dynamics implementations

3.2.1 Monofluid, two-fluid, and multi-fluid: terminology explained

The terms monofluid, two-fluid, Lagrangian dust, and Eulerian dust are often used in this paper, as well as much of the literature in the field. Therefore, it is useful to clarify what is meant by these terms in this context. The two formulations of dust hydrodynamics, 'two-fluid' or 'multi-fluid' and 'monofluid', differ fundamentally in one aspect, namely in how the equations for the dust-gas mixture are written. In a multi-fluid approach, dust species are treated as additional fluids with their own Navier-Stokes equations coupled to the gas by a drag term. The monofluid description is instead designed to describe the mixture as a single fluid, with velocity equal to the barycentre velocity of the mixture, and evolving the drift velocity, which is the vector field that describes the velocity difference between dust and gas. The two descriptions are mathematically equivalent; however, starting from the monofluid description, it is possible to construct approximate descriptions of the mixture by truncating the drift velocity expansion at different orders. The monofluid method in this paper is a TVA of the complete monofluid equations (see Sect. 3.2.2). A further level of complexity is added to the problem when one considers that, when constructing a multi-fluid solver, it is not required that the same approach be used to solve all the fluid equations (provided the drag term is treated correctly). Hybrid schemes solve the gas dynamics on a Eulerian grid, but solve the dust dynamics using Lagrangian super-particles. This approach may provide some advantages, especially in the limit of St → inf where the dust particle dynamics becomes ballistic mechanics. The two-fluid scheme in this paper is one such 'hybrid scheme', which solves the gas dynamics on a grid, but solves the dust dynamics by integrating the motions of Lagrangian dust particles.

3.2.2 Monofluid implementation

The implementation in RAMSES of dust dynamics in the diffusion approximation - which is based on the TVA - and the monofluid formalism was presented in Lebreuilly et al. (2019). We use the multi-species framework, which allows us to follow the coupled dynamical evolution of 𝒩_d different grain species on a single numerical simulation. Below, we describe the main characteristics of our monofluid implementation. We refer readers to Youdin & Goodman (2005), Laibe & Price (2014b,a), Price & Laibe (2015), Hutchison et al. (2018), and Lebreuilly et al. (2019) for more details on the derivation of the diffusion approximation in the monofluid formalism.

We first define the monofluid hydrodynamical quantities. The total mixture density ρ is

$$\rho \equiv {\rho }_g+{\displaystyle \sum _{k=1}^{{\mathcal{N}}_{\text{d}}}{\rho }_{d,k},}$$(9)

where ρ_d,k is the dust fluid density of the kth size bin. The total dust density is simply ${\rho }_{\text{d}}={\displaystyle {\sum }_{k=1}^{{𝒩}_{\text{d}}}{\rho }_{\text{d,}k}}$. The mixture barycentric velocity ν is the defined as

$$v\equiv \frac{{\rho }_{\text{g}}{v}_{\text{g}}+{\displaystyle {\sum }_{k=1}^{{\mathcal{N}}_{\text{d}}}{\rho }_{\text{d,}k}{v}_{\text{d},k}}}{{\rho }_g+{\rho }_{\text{d}}},$$(10)

where v_d,k is the velocity of the k bin dust fluid. The dust ratio of dust species k is then ∈k = ρ_d,k/ρ and the total dust ratio is ɛ = ρ_d/ρ¹.

The dust and gas mixture dynamical coupling is characterised by the stopping time

$$t_{s,k}\equiv \sqrt{\frac{\pi \gamma }{8}}\frac{{\rho }_{\text{grain,}k}}{\rho }\frac{s_{\text{grain,}k}}{c_s}=\frac{{\rho }_g}{\rho }{t}_{\text{s,grain}},$$(11)

which leads to the effective stopping time T_s,k in the case of multiple dust species. This coupling accounts for the interaction between dust species due to their cumulative back-reaction on the gas:

$$T_{s,k}\equiv \frac{t_{s,k}}{1-{ϵ}_k}-{\displaystyle \sum _{l=1}^{{\mathcal{N}}_{\text{d}}}\frac{{ϵ}_l}{1-{ϵ}_{ls,l}}}$$(12)

Finally, we introduce the mean stopping time,

$${\mathcal{T}}_s\equiv \epsilon {\displaystyle \sum _{k=1}^{{\mathcal{N}}_{\text{d}}}{ϵ}_k}{T}_{s,k}.$$(13)

The set of equations of the evolution of the gas and dust mixture is then

$$\begin{array}{l}\frac{\partial \rho }{\partial t}+\nabla \cdot \left[\rho v\right]=0,\hfill \\ \frac{\partial \rho v}{\partial t}+\nabla \cdot \left[\rho v\otimes v+P_g\mathbb{I}\right]\rho f,\hfill \\ \frac{\partial {\rho }_{\text{d},k}}{\partial t}+\nabla \cdot \left[{\rho }_{\text{d},k}\left(v+\frac{{\text{T}}_{\text{s,k}}\nabla P_{\text{g}}}{\rho }\right)\right]=0,\forall k\in \left[1,{\mathcal{N}}_{\text{d}}\right],\hfill \end{array}$$(14)

where P_g is the gas pressure, and f is the acceleration applied to account for the turbulence driving. The back-reaction of the dust onto the gas is taken into account in this monofluid framework. We note that no energy equation is considered here, because we use the isothermal approximation.

This set of equations is integrated using the classical second-order Godunov scheme of RAMSES, which is complemented by the second-order scheme developed by Lebreuilly et al. (2019) to integrate the dust differential dynamics. We refer readers to Lebreuilly et al. (2019, 2020) for more details on the numerical implementations, as well as tests and applications to prestellar dense core collapse and protostellar disc formation.

3.2.3 Two-fluid implementation

We developed an implementation to account for the dynamics of inertial Lagrangian super-particles in the two-fluid formalism for this study. Our implementation is based on the particle-mesh method used for tracer particles already implemented in RAMSES. The method is well tested, and is used, for example, to run large chemical networks as postprocessing, as it stores the density and temperature of individual parts of the fluid (e.g. Coutens et al. 2020).

The dust fluid here is represented by Lagrangian super-particles with constant mass moving at the dust velocity v_d. Their motion is described by

$$\frac{\text{d}{v}_{\text{d}}}{\text{d}t}=\frac{v_{\text{g}}-v_{\text{d}}}{t_s}.$$(15)

In the case of strong dynamical coupling, the drag term becomes stiff, which leads to prohibitive time-step constrains. We therefore solve the equation of motion analytically, which leads to the integration scheme

$$v_{\text{d}}^{n+1}=v_{\text{g}}^{\text{n+1}}\left(1-e^{\Delta t/t_s^n}\right)+v_{\text{d}}^n{e}^{-\Delta t/t_s^n}.$$(16)

This analytic formulation naturally yields the limits St ≪ 1 and St ≫ 1. The particle velocity is updated using a CIC interpolation of the gas density and velocity from the grid. The particles are then moved using the second-order midpoint scheme of RAMSES (Teyssier 2002). This scheme is very similar to the one used by Mattsson et al. (2019a) using the PENCIL code (Pencil Code Collaboration 2021). As mentioned earlier, with this implementation, the numerical resolution of the dust depends on the number of particles used, but the accuracy of the drag estimate is limited to the grid resolution used for the gas. Finally, this two-fluid implementation does not take into account the back-reaction of the dust onto the gas.

3.3 Initial conditions and numerical setup

The initial physical setup consists of a uniform density box, which represents a portion of a large molecular cloud. We use a uniform grid with resolution ranging from 128³ to 512³, always assuming periodic boundary conditions. The initial density is ρ₀ = 4.4 × 10⁻²¹ g cm⁻³ and the box length is set to 4 pc, which satisfies the Larson relations (Eq. (2)). The gas thermal evolution is forced to remain isothermal with a temperature of T₀ = 10 K. In all models, the grain intrinsic density is ρ_grain = 1 g cm⁻³. This value is rather low if one considers that the dust grains are made of a mixture of graphite (2.2 g cm⁻³) and crystalline olivine with composition MgFeSiO₄ (3.6 g cm⁻³). Nevertheless, we note that at a given Stokes number, the dust size and intrinsic density are degenerate. If we were to chose an intrinsic density of 1 g cm⁻³ as used in other studies such as Tricco et al. (2017), we would have to divide the dust size by a factor of 3.

In this study, we explore only one level of turbulence in order to focus on the numerical methods. We set the RMS velocity of the dust and gas mixture in accordance with the Larson relation (Eq. (4)), i.e. υ_rms ≃ 1.9 km s⁻¹, which corresponds to a Mach number M ≃ 10 . The turbulent crossing time is then t_turb ≃ 2.4 Myr. For the turbulence driving, we use a unique autocorrelation timescale 𝒯 (timescale for a full change of the turbulent field), set to one-quarter of the turbulent crossing time, namely 𝒯 ≃ 0.6 Myr.

In the following, we present two types of numerical simulations. In the first, we identify the critical dust grain size for decoupling using the monofluid solver presented above. In this setup, we have introduced 𝒩_d = 10 bins of dust-grain size, ranging from 1 nm to 12 μm, which corresponds to Stokes numbers ranging from 10⁻⁵ to 0.18, respectively. The fiducial resolution is 256³. In addition, in Appendix B, we investigate the effect of the back-reaction by varying the initial total dust ratio ɛ₀ from 0.01 (fiducial value) to 10⁻⁵ (a negligible amount of dust should not affect the gas dynamics). In Appendix A, we also propose a resolution convergence study where the grid size ranges from 128³ to 512³. In the second setup, we compare the monofluid results with the ones obtained using the dust Lagrangian particles (2F runs). To this end, we use 𝒩_d = 5 bins of dust-grain size ranging from 1 nm to 10 μm. Thanks to our numerical tool, we can compare the two solvers using a single numerical simulation. As our two-fluid formalism does not allow us to account for the back-reaction, we have set the total dust ratio to the very small value of 10⁻⁶ in the monofluid solver. On top of this, we have added N_p dust Lagrangian super-particles in the two-fluid solver. These particles do not impact the hydrodynamical fields. Table 1 summarises the numerical parameters of the various simulations we run in the context of this study.

3.3.1 Expected limits to the monofluid and two-fluid methods in this setup

From our analysis in Sect. 2.2, the expected critical dust-grain size corresponding to St = 1 is $s_{crir}\underset{˜}{\underset{˜}{>}}70\text{ μm}$. As already mentioned, in our experience, grains with Stokes number St ≥ 0.1 already decouple significantly on a dynamical timescale. We therefore expect grains with size $s\underset{˜}{\underset{˜}{>}}7\text{ μm}$ to decouple in our setup.

3.3.2 Maximum dust velocity

In the low-gas-density regions, the Stokes number can exceed unity even for the smallest grains. In these conditions, the TVA approximation breaks down, leading to large dust velocity. We therefore limit the dust velocity estimated from our monofluid module to the mean RMS turbulent velocity, that is, Mcs. We also apply the correction for supersonic flow when the velocity drift between the dust and the gas exceeds the sound speed (Kwok 1975; Draine & Salpeter 1979).

3.4 Postprocessing of the grid versus particle data

In this study, we analyse sets of data that encompass both Eulerian and Lagrangian variables. Our analysis is based mostly on volume-weighted probability density functions (PDFs) and density cuts through the computational domain. We use well-tested suites of postprocessing tools, as follows.

For the grid quantities (monofluid), we simply bin the grid cells according to the value of the density for the PDF and the slices are straightforward. We use the OSYRIS² visualisation package for RAMSES. For the particle quantities (Lagrangian dust and tracers), we follow the procedure described in Price & Federrath (2010) and Price et al. (2011). First, all particles have the same mass, m_p, which is equal to the total mass of the dust species divided by the number of particles). To produce density maps, we project the particle distribution onto a regular grid, which has the same dimensions as the grid used for the monofluid. We then use the Smoothed-particle hydrodynamics (SPH) density calculation routine from the PHANTOM code (Price et al. 2018), where the density and smoothing length are iterated self-consistently using classical SPH procedures. We produce cross-section slices of the density field using the SPLASH visualisation software (Price 2007). To obtain the volume-weighted PDF from the particles, we follow Price et al. (2011) and again interpolate the density on an adaptive grid.

4 Monofluid results

4.1 Time evolution

Figure 1 shows the time evolution of the RMS velocity $\langle \sqrt{{\upsilon }^2}\rangle$ and of the density clumping factor of the various fluids (dust, gas, barycentre). The clumping factor C_X of a quantity X is defined as C_X = 〈X²〉/〈X〉². A clumping factor C_X > 1 indicates that the field X shows significant variation. The RMS velocity and the density clumping factor of the barycentre and of the ten different dust species are shown. In addition, we show the evolution of the total dust density and gas density clumping in the clumping factor plot. After roughly one turbulent crossing time, the various quantities oscillate around a mean value and do not show any trend as a function of time, which indicates that we have reached a steady state. While dust species with size s < 1 μm do not show any significantly different evolution from the barycentre, the largest dust grains tend to show a decoupling with amplitude in the measured quantities, which remains constant with time. In addition, the fluids consisting of dust grains with sizes s ≥ 4 μm have RMS velocities and clumping factors that increase with size. The larger the grain, the larger the turbulence amplitude and the stronger the clumping. Interestingly, the total dust density also shows a stronger clumping than the barycentre and the gas. This is mostly due to the fact that under the MRN size distribution, most of the dust mass is contained in the largest grain species. In summary, this first qualitative analysis shows that grains with St > 0.1 decouple from the gas dynamics, in good agreement with our theoretical estimate. In addition, the gas and barycentre quantities exhibit almost identical time evolution, meaning that the gas-phase evolution corresponds to the evolution of the barycentre.

In addition, we verified that the evolution of the total dust mass contained above given the total density threshold remains roughly constant with time. This indicates that there is no cumulative trapping of the dust grains at any density as a function of time, and that the degree of decoupling also remains constant with time. The quantities we derive in the following are then independent of time, and we can safely analyse our models at a given time snapshot, even though the absolute amplitude of the different quantities might fluctuate with time.

Fig. 1 Time evolution of the barycentre and dust-species RMS velocity (top) and of the clumping factor of the density of the different species (bottom) in the fiducial model 256MRN. The time evolution has been normalised according to the turbulent crossing time tcross. In the bottom panel, the dotted line shows the clumping factor evolution of the total dust density ρd and the dashed-dotted line shows the clumping factor evolution of each dust species concentration ϵi. The evolution of the clumping factor of the gas density (dashed) is indistinguishable from that of the barycentre (solid) density.

4.2 Characteristic features of the decoupled regions

Figure 2 shows density maps at time 3t_cross of the variations of the gas density and dust ratio ϵ_i of each dust species relative to their initial value. The gas density varies over more than four orders of magnitude with large velocities in the low-density regions. The amplitude of the density variation is typical of isothermal compressible turbulence, with density jumps at shocks varying as ∝ ℳ². The amplitude of the dust ratio variations increases with dust-grain size, up to very large values (more than eight orders of magnitude) for the largest grains. The region that show a large variation in dust ratio, for both depletion and enrichment³, are associated to densities ρ_g < ρ_g,0 where the largest density gradients develop. As we use an isothermal equation of state, the pressure gradients are proportional to the density gradients. The dust velocity shift is also inversely proportional to ρ², and so the smaller the density, the larger the velocity drift, which is consistent with our findings. We also observe strong depletion of dust density for dust grains >0.1 μm These depleted regions correspond to low gas densities and therefore to low dust densities. As a consequence, the dust enrichment is mostly at higher densities, because the total dust mass is conserved. Interestingly, we observe that the strongly depleted regions (shown in dark blue in Fig. 2, corresponding to log(ϵ_i;/ϵ_i,0) < −2) always sit in low gas densities where log(ρ_g /ρ_g,0) < −1. These regions also exhibit the largest dust velocity. On the other hand, the dust-enriched regions do not correspond to higher densities, but rather are found where −1 < log(ρ_g/ρ_g,0) < 1. The high-gas-density material shows a nearly uniform dust ratio, which corresponds to the initial value for all dust species, except the largest dust species where variations of a factor of a few can be observed.

In order to better characterise the regions of decoupling, Fig. 3 shows the dust-to-gas ratio variations as a function of the gas density for the total amount of dust, as well as the dust ratio for the various dust-grain species we modelled. All dust species show the largest variations in the dust ratio for gas density ρ_g < ρ_g,0. The amplitude of the variations, in particular depletion, increases with dust size. Interestingly, we see that the smallest dust grains (s < 0.1 μm) are only strongly depleted at the smallest gas density, that is, ϵ_i;/ϵ_i,0 < 10⁻², while the enrichment is effective for gas densities −4 < ktg(ρ_g/ρ_g,0) < 1. The bulk of the mass remains at the initial dust ratio for the smallest grains. These grain species also exhibit Stokes numbers St < 1 for the bulk of the mass, which indicates that our monofluid approximation is well suited for these species. The grains with sizes >0.1 μm show stronger depletion even at gas densities corresponding to ρ_g,0. The amplitude of the dust enrichment is lower compared to the smaller grains, typically ϵ_i;/ϵ_i,0 < 100, but corresponds to a larger bulk mass. In terms of total mass budget, as the largest grains carry the mass, the total dust-to-gas ratio exhibits some variation, but this is only moderate (less than one order of magnitude for the bulk). For grains with sizes s > 1 μm the maximum dust variation occurs at St > 1, which is outside the range of validity of the diffusion approximation. In Sect. 4.3, we show that this is not a problem for the total mass budget because the bulk of the mass is found in the St < 1 regions. However, the amplitude of the dust variations might be inaccurately estimated.

Figure 4 shows the cumulative volume-weighted PDFs of the dust-ratio variations for various dust species, including the three largest ones. For grains with sizes s < 0.1 μm, the bulk of the volume does not show variation in dust ratio. The dust-enriched region corresponds to less that 5% of the total volume for the 0.03 μm dust grains. As explained in the previous paragraph, this material corresponds to low gas density. For larger grains, s > 1 μm, more than half of the total volume shows dust-enriched regions. About 4% of the 12 μm dust-grain material shows dust-ratio variations of greater than a factor of two, and about 1.3% shows variations of greater than a factor three.

Figure 5 shows the variation of the total density along a line of sight, as well as the drift velocity and dust ratio for the 12 μm dust grains in the fiducial 256MRN model. As our model is isothermal, the density directly traces the pressure, meaning that the density variations correspond to the pressure ones. In the TVA approximation, the drift velocity is proportional to the pressure gradient. The pressure and drift velocity evolutions are anti-correlated. The drift-velocity minima correspond to pressure maxima. The dust-enriched regions are found in the pressure maxima. As already reported, the highest concentrations are measured in the regions of intermediate density (ρ < ρ₀). This behaviour confirms that our implementation of the TVA correctly captures the predicted physics.

Fig. 2 Gas density (top-left) and dust ratio variation maps in the xz-plane for the ten dust bins in the 256MRN run. The variations are given relative to the initial value of the quantity and are shown in logarithmic scale. In the dust-ratio panels, the red colour shows dust-ratio enhancement, while blue indicates a dust-ratio decrease. The isocontours correspond to gas-density variations of −1 (orange), 0 (black), and +1 (green) in logarithmic scale. The arrows represent the barycentric velocity (top-left), and the dust-grain drift velocity vectors in the plane for all the other plots. All plots are made at a time corresponding to 3tcross.

Fig. 3 2D histograms of the dust-ratio variations as a function of gas density for the fiducial model 256MRN. The top-left panel shows the variation in total dust-to-gas ratio ɛ, while the other panels show the dust-ratio variations of each dust species. The colour coding indicates the mass of the gas in log scale (in solar units). The horizontal grey dotted line separates the dust-enriched and dust-depleted regions. The vertical line indicates St = 1 (solid) for each dust species. All quantities are averaged over more than one dynamical time tcross and the variations are normalised to the initial value of each dust ratio.

Fig. 4 Volume-weighted cumulative PDF of the dust ratio variation for the 1.6 nm, 0.03, 1.7, 4.5, and 12 μm dust grains in the fiducial model 256MRN. All quantities are averaged over more than one dynamical time tcross.

Fig. 5 Profile of the total density ρ (blue) and the drift velocity (orange) and dust ratio (dotted green) for the 12 μm dust grains in the fiducial model 256MRN. The density is normalised by its initial value. The norm of the drift velocity with respect to the barycentre is normalised to 1 km s−1. The variations in dust ratio are computed with respect to its initial value. The quantities are plotted against the x-axis at a random time (>2tcross) and at a random location in the yz-plane. Regions with positive relative dust ratio correspond to dust-enriched material.

4.3 Validity of the diffusion approximation

Figure 6 shows the volume- and density-weighted PDFs of the Stokes number for various dust species, including the three largest ones. The bulk of the mass of all dust species sits in St < 0.1 regions. However, the volume-weighted PDFs of the largest grains show large volume fractions with St > 1. This corresponds to ≃60% of the total computational volume for the 12 μm grains, ≃25% for the 4.5 μm grains, and ≃15% for the 1.7 μm grains. In terms of mass, the fraction of the total dust mass of the largest bin (12 μm) where the Stokes number is St > 1 is ≃1% (4 ≃22% for St > 0.1).

We also checked the fraction of the total volume and mass of each dust species where the limitation of the dust velocity is negligible. For the most decoupled dust grains (12 μm), the fraction of the mass where the velocity floor is active is always $\underset{˜}{\underset{˜}{<}}{10}^{-6}$, which corresponds to a fraction on the order of a few times 10⁻⁵ of the total volume. This confirms that the dust diffusion approximation and our numerical implementation are well suited for our purposes.

4.4 Summary of the monofluid results

The main results of our analysis of the monofluid results are as follows:

• Dust grains of size <0.1 μm remain well coupled to the gas.

• Dust grains of size >1 μm show strong variation in concentration. The decoupling regions correspond mainly to low-density material.

• Dust concentrates in the pressure maxima.

• The total dust fraction does not show variation in the high-density regions.

• Our monofluid and TVA implementation is well-suited for the bulk of the mass (99% for the largest grains). However, the maximum dust variation for large grain sizes occurs for physical conditions that are outside the range of validity of the diffusion approximation.

Fig. 6 Cumulative PDF of the Stokes number for the 1.6 nm dust grains and the 0.03, 1.7, 4.5, and 12 μm dust grains in the fiducial model 256MRN. The coloured dotted lines show the mass-weighted PDF and the coloured solid lines show the volume-weighted PDFs. The Stokes number is computed using the local gas density and the global turbulent crossing time turb. All quantities are averaged over more than one dynamical time tcross. The vertical lines indicate the St = 1 (solid) and St = 0.1 (dashed).

5 Comparison between monofluid and two-fluid formalisms

In this section, we focus on the comparison of the above results based on the monofluid formalism with those produced using our two-fluid implementation based on Lagrangian superparticles. We also present a convergence study where we vary both the grid resolution N and the number of particles N_p.

Our comparison focuses on the evolution of five dust-grain size populations, which are treated using both the monofluid and the two-fluid solvers in single simulation runs. The dust-grain size varies from 1 nm (well-coupled grains) to 10 μm (least coupled), and the dust-to-gas-mass ratio of each dust species is set to 10⁻⁶. Indeed, as already mentioned, as our two-fluid implementation cannot account for the back-reaction of the dust-grain dynamics on the gas dynamics whereas the monofluid implementation can, choosing a low value of the dust-to-gas ratio means we are mostly neglecting the back-reaction in the monofluid solver. All other aspects of the monofluid setup are identical to the ones presented in the previous section (except for grain size and total dust mass). In addition to the two solvers, we have included classical tracer particles designed to trace the dust fluids handled in the monofluid runs, that is, with the velocity computed in the terminal-velocity approximation. The tracer-particle populations are used for comparison with both the monofluid and the two fluid results.

Fig. 7 Density slice in the z = 0 plane obtained from the 256_2F run for the 1 nm (top) and 10 μm (bottom) dust grains at a time corresponding to ≃ 4.5tcross. The left column represents the density field of the monofluid dust, while the right column represents the density field interpolated from the Lagrangian dust particles. The density is normalised by the initial density ρ0. The black lines represent the contour for which the gas density is equal to its initial value ρ0.

5.1 Density fields

Here, we investigate the results of the 256_2F run. In this run, each dust-grain size is represented by 10⁶ 'two fluid' Lagrangian particles and another 10⁶ monofluid 'tracer' particles. Figure 7 shows slices of dust-density variations for the 1 nm and 10 μm dust-grain fluids for both the monofluid and the two-fluid approximations. All quantities are normalised to the initial dust density of each species. We also show the isocontour, which indicates the threshold where the gas density is equal to its initial value. We can therefore compare the dust variations to those of the gas. For the 1 nm dust grains, the dust-density variations from the monofluid closely match the gas-density variations, which is in very good agreement with what we observe in the previous section for the most coupled dust species. The 1 nm dust-grain density obtained with the dust as Lagrangian particles globally matches the density of the gas. The gas-depleted regions correspond to dust-depleted regions. However, one can distinguish small islands of dust-depleted regions within gas-enriched regions and dust-enriched filaments in some gas-depleted regions, which is indicative of decoupling. For the 10 μm dust grains, the dust-density variations roughly match the gas variations for the monofluid, with large dust-depleted regions of more than a factor of two. Conversely, the variations in the dust-density field from the Lagrangian particles are less important than those found in the monofluid, with smaller enriched and depleted regions. Large dust-depleted islands are also observed within the gas-enriched regions.

5.1.1 Tightly coupled dust grains (1 nm)

Figure 8 shows the averaged PDFs of the density obtained from the gas and the monofluid, the Lagrangian dust particle, and the monofluid tracer-particle populations for the 1 nm and 10 μm dust grains. We first focus on the PDFs of the 1 nm dust grains.

We observed two clearly distinct behaviours. On the one hand, the gas and monofluid PDFs match relatively well, in particular for the largest densities, while on the other, the tracer and Lagrangian-dust-particle PDFs are identical. The dust-tracer-particle PDF does not match the monofluid PDF at all, while they should match.

First, comparing the monofluid and the gas PDFs, there are more dust-depleted regions as well as gas-depleted ones. As a consequence, there is a slight shift of the peak to high density for the dust population. This is seen in the previous section for tightly coupled dust grains. For the 1 nm dust grains, the expected behaviour would be for the dust density PDFs to be similar to the gas ones, as the gas drag is strongest. However, at low gas density, the Stokes number becomes large enough for the dust to slightly decouple.

Second, the monofluid tracer particles should trace the monofluid dust species, but as the two PDFs are very different, they do not. We observe an excess in both the low-and high-density tails of the tracer particles distribution, meaning that some particles are trapped in the high-gas-density regions, resulting in a larger depletion in the regions of low gas density. This result suggests that particles may be subject to artificial clustering, which we investigate in the following subsection. From the identical tracer- and Lagrangian-dust-particle population PDFs, we see that both particle populations experience exactly the same velocity field, computed using either the terminal velocity approximation (monofluid) or our analytical integrator scheme (two-fluid Lagrangian particles). Given the small Stokes numbers of the 1 nm dust grains, the computed velocity is equal to the gas velocity, which further verifies that the dust velocity computed by the monofluid approximation is correct in this regime. These two PDFs depart significantly from the gas-density PDF as well.

Fig. 8 PDFs of the logarithm of density ln ρ in the 256_2F run for the 1 nm (top) and 10 μm (bottom) dust grains. The density variations are normalised to the initial density values for each field.

5.1.2 Least-coupled dust grains

For the 10 μm dust grains, we observe four different PDFs in Fig. 8. First, the tracer particles and monofluid PDFs are different, indicating that the tracer particles do not trace the dust fluid as they should for all Stokes regimes. This brings into question the reliability of the results obtained with the classical Lagrangian tracer particles. Further, the monofluid PDF departs from the gas PDF for low density. At high density, where the Sokes number is St ≪ 0.1, the two PDFs (gas and monofluid) are a close match, as indicated by the small dust-ratio variations observed previously. The Lagrangian dust-particle PDF falls in the middle of the gas and monofluid PDFs, close to the PDF of the tracer particles at high density and to the gas PDF at low density. This last observation seems to indicate that the two-fluid solver based on Lagrangian particles also has problems handling the large gas density corresponding to small Stokes numbers (clustering). The relatively good match with the gas PDF at low density is also questionable, because the dust-grain dynamics should be decoupled from the gas dynamics given the high Stokes numbers. We see in the following sections that this match is not robust, and that the low-density slope of the PDF depends on the number of particles N_p.

5.1.3 Intermediate summary

The main results of the analysis of the density PDFs are the following:

• Two-fluid results exhibit stronger decoupling than the monofluid ones, for all dust sizes.

• For small dust-grain sizes, the tracer particles and very small dust Lagrangian super-particles have identical density PDFs.

• For large dust-grain sizes, the density PDFs obtained with the monofluid and the two-fluid formalisms do not agree.

5.2 Dust ratio

In this section, we focus on the dust variations obtained in the 256_2F run. We computed the dust enrichment in the two-fluid case with dust as Lagrangian particles by projecting the dust-particle distribution onto a uniform grid of the same resolution as that used for the gas, and divided by the gas density. For the dust species handled with the two-fluid solver, we can determine a maximum dust enrichment by assuming that the minimum length scale of the particle distribution cannot go beyond the gas resolution set by the grid for the gas. We define the maximum dust enrichment ϵ_max by setting the minimum adaptive resolution length h_min (see definition in the following section) equal to the grid size Δx:

$${ϵ}_{\max }\equiv \frac{{\rho }_{\text{d,max}}}{\rho }=\frac{m_{\text{p}}}{\Delta x^3\rho },$$(17)

where m_p is the mass of a particle.

Figure 9 shows 2D histograms of the dust enrichment obtained in the 256_2F run for the monofluid and Lagrangian-particle dust species of 1 nm and 10 μm in size. First, the monofluid results are different from those observed in Fig. 3 for the 256MRN run. In particular, there is almost no dust enrichment for the 1 nm dust grains. In Appendix A, we show that this is due to the neglect of the dust back-reaction onto the gas, which mainly affects the dust-enriched tail of the PDFs. The monofluid results for the 10 μm dust grains are similar to the ones observed in the 256MRN run.

We now compare with the results obtained in the two-fluid case with dust as Lagrangian particles. For both grain sizes, monofluid and two-fluid solvers give very different results. For the well-coupled dust grains, there is a huge spread in dust enrichment (up to at least two orders of magnitude) in the full density range. As a consequence, the depletion is also large at all densities. For the 10 μm dust grains, the enrichment is higher at low gas density compared to the monofluid results. In particular, the regions of lowest gas density are only dust enriched. This corresponds to regions where St > 1 . At high gas density, the enrichment is similar to that observed for the smallest grains, which indicates that they undergo the same dynamical (de-)coupling with the high-density gas at low Stokes number. The red line in Fig. 9 indicates the maximum dust enrichment, which corresponds to the maximum dust density allowed by the gas grid resolution. The entire region above ϵ_max corresponds to dust enrichment with resolution length smaller than the gas resolution, that is, unresolved dynamical coupling. For 1 nm dust grains, only the regions of high gas density are affected, while for the 10 μm dust grains, the upper part of the distribution follows a trend inversely proportional to the gas density. This suggests that the maximum dust enrichment is regulated by the grid resolution used for the gas dynamics. We observe a maximum enrichment of about four orders of magnitude in the low-gas-density material, which corresponds to the maximum enrichment permitted by the gridresolution usedforthe gas dynamics. Even if the two-fluid implementation with dust as Lagrangian particles is better suited for large Stokes numbers than our monofluid implementation, it remains severely limited by the grid resolution, leading to numerical artefacts in of high dust density.

Fig. 9 Same as Fig. 3 but for the 256_2F run. The left column shows the monofluid results and the right column shows the results when dust is modelled as Lagrangian particles. The top row shows the 1 nm dust grains and the bottom row the 10 μm dust grains. The red dashed line indicates the maximum dust enrichment as defined in Eq. (17). All the material that is above this line has an excessive dust concentration.

Fig. 10 PDF of the effective adaptive resolution length h of the 1 nm (dotted lines) and 10 μm (solid lines) dust-grain fluid particle distributions for the 256_2F run. The colour coding indicates the particle distribution (two-fluid dust as Lagrangian particles in blue, tracer in red). The vertical line indicates the numerical resolution of the grid (2563).

Intermediate summary

The main results of our analysis of the dust-ratio variations are as follows:

• Our two-fluid implementation with dust as Lagrangian particles leads to dust enrichment that exceeds the maximum value set by the grid resolution used to compute the drag from the gas (see Fig. 9).

• The high-density regions show dust-ratio variations at all sizes with the two-fluid solver that were not observed in the monofluid results.

5.3 Adaptive resolution length PDFs

We define the adaptive resolution length of the dust-particle distribution as

$$h={\left(\frac{m_{\text{p}}}{{\rho }_{\text{d,i}}}\right)}^{1/3},$$(18)

where ρ_d,i is the local density associated with each particle, which is computed following the procedure described in Sect. 3.4. Therefore, h is a measure of the highest resolution reached in the particle distribution, to be compared with the grid scale Δx.

Figure 10 shows the PDFs of the adaptive resolution length h obtained for the particle distributions of run 256_2F. The scales covered range over more than six orders of magnitude. At large values (high densities), the PDFs of the 10 μm behave very differently. The two-fluid dust Lagrangian-particle distribution is more compact, while the particle distribution of the tracer exhibits larger voids. Conversely, all the PDFs show the same trend at small scales, namely high densities. Interestingly, when we compare the adaptive resolution length to the resolution of the grid used for the gas dynamics, we observe that for h < 2Δx, the PDFs are almost identical, independently of the dust-grain size and integration method (two-fluid dust Lagrangian or tracer particles). This indicates that the dust-grain dynamics at scales smaller than the grid size is dominated by numerical artefacts and the clustering is therefore artificial. Indeed, the gas drag is the only mechanism that can move dust particles, and the gas drag cannot be accurately computed on scales smaller than Δx. As a consequence, the clusters of high dust density are characteristic of regions dominated by numerical diffusion, which leads to the conclusion that tiny grain clustering observed in the particle distributions has a numerical origin. It remains unclear as to what extent this clustering affects the rest of the simulation volume, in particular the low-density tail of PDF.

Intermediate summary

The main results of the analysis of the Lagrangian-particle clumping scale are as follows:

• The two-fluid Lagrangian-particle dust results show clustering on scales smaller than the grid size, where the physics of the dust dynamics is not resolved. This behaviour is not physical, but numerical in nature.

• All dust sizes show the same trend of artificial clustering below the grid scale.

5.4 Effects of the grid resolution and the number of particles

In this section, we investigate the effect of varying the grid resolution N and the number particles N_p in order to test the effect of the relative resolution used for dust particles in comparison with the relative resolution of the gas. First, we compare the 256_2F and 128_2F runs, in which the only difference is the resolution used for the grid. Figure 11 shows the PDFs of the 1 nm and 10 μm density and the adaptive resolution length of the dust-particle distributions in these two models. For the 1 nm dust grains, the PDFs of the dust density are identical regardless of the treatment (two fluid vs. tracer) at a given resolution. The PDFs of the 128_2F are wider compared to the 256_2F ones, which indicates a more important decoupling from the gas dynamics as the gas resolution decreases. For the 10 μm dust grains, we observe the same trend, but with different PDFs for the two-fluid Lagrangian-particle distribution and tracer particle distribution, as already reported in Sect. 5.1. From the adaptive length PDFs, we observe that the 128_2F exhibits the largest clustering effect. One would have naively expected the 128_2F model to give better results than the 256_2F model for the dust-particle distribution, as the particle resolution per grid cell is eight higher in the 128_2F model. However, as demonstrated in the previous section, the dust-particle dynamics is strongly affected by the numerical resolution used for the gas on the grid, and the amount of artificial clustering is higher.

Figure 12 shows the PDFs of the 1 nm and 10 μm dust-grain density and adaptive resolution length of the dust-particle distributions in the 128_2F (10⁶ particles per grain size bin), 128_2F_LR (125 × 10³ particles), and 128_2F_VLR (10³ particles) models. The aim is to test the effect of the number of particles used to describe the dust fluid at a given resolution for the gas. We see that the dust-density PDFs get tighter as the number of particles decreases. The best match between the two-fluid and the monofluid results is found for the 128_2F_LR run for the 1 nm dust grains. For the 128_2F_VLR model, the results for the 1 nm and 10 μm dust grains are qualitatively similar: the resolution is insufficient to correctly capture the dust-grain dynamics. Finally, the minimum adaptive resolution length increases when the number of particles decreases.

Our short resolution study does not indicate a convergence on the particle-distribution properties. However, it is clear that the properties of the two-fluid Lagrangian-particles and tracer-particle evolution strongly depend on the grid resolution, as well as on the total number of particles used to sample the dust fluid. We discuss the limits of our different methods used to follow the dust dynamics in the following section.

Fig. 11 PDF of the dust density of the 1 nm (left) and 10 μm (middle). The colour coding indicates the dust fluid treatment (two fluid Lagrangian particles in blue, tracer particles in red, monofluid in black). For clarity, only the monofluid quantities of the 256_2F are shown. The solid lines show the results of the 256_2F run, while the dotted lines show the results of the 128_2F. Right: PDF of the adaptive length of the two-fluid dust Lagrangian-particle distributions (solid lines for the 10 μm dust grains, dotted line for the 1 nm) for the 256_2F (blue) and 128_2F (red) models. The vertical dashed lines indicate the grid numerical resolution (black for 256_2F, grey fro 128_2F).

Fig. 12 Same as Fig. 11 for the 128_2F, 128_2F_LR, and 128_2F_VLR models. In the density PDFs, the blue lines indicate the results for the two fluid dust Lagrangian-particle distribution while the red lines indicate the tracer particle ones.

Intermediate summary

The main results of our analysis of the two-fluid Lagrangian-particle dust method convergence are as follows:

• The two-fluid results are highly sensitive to both the grid resolution and the total number of particles.

• Our resolution study of the grid and particle resolutions does not show convergence.

6 Discussion

6.1 Limits of the monofluid and two-fluid formalisms

Our results presented in Sect. 4 show that the monofluid approximation is well suited for our setup, even for the largest grains, which exhibit large Stokes numbers St > 1 in significant parts of the computational volume. For all dust sizes, the regions where St > 0.1 show dynamical size sorting and dust-ratio variations of orders of magnitude. However, the amplitude of the largest dust-ratio variations we report for grains with sizes s > 1 μm is questionable, because the diffusion approximation we use is not well suited for large Stokes numbers St > 1 . However, the mass enclosed in these regions remains negligible compared to the total mass (≃10⁻²). In addition, we find no evidence of systematic effects on the time evolution of dust enrichment, because the turbulence resets the flow within a crossing time. More importantly, we show that dust grains with sizes s ≥ 4 μm significantly decouple from the gas, which is consistent with the critical size of $s\underset{˜}{\underset{˜}{>}}7$ that we derive in Sect. 3.3.1 for our setup and St > 0.1. Finally, our monofluid solver does account for the back-reaction of the dust on the gas dynamics. We show in Appendix B that the back-reaction crucially affects the dust enrichment of the smallest grains (nanometers), which is concentrated in the regions of low gas density. Indeed, as large dust grains accumulate, the gas gets pushed away from the regions of high dust concentration. The smallest dust grains, which are well coupled to the gas, are therefore expelled. We do not observe this high concentration at low density in our models that neglect the back-reaction. Our numerical experiments and implementation of the monofluid equations using the terminal velocity approximation do not allow us to investigate the dynamics of grains larger than 10 μm. An alternative implementation, such as full monofluid (Laibe & Price 2014a,b), a full two-fluid Eulerian treatment (Benitez-Llambay et al. 2019), or the two-fluid dust as Lagrangian superparticles we employed (with the addition of dust back-reaction), should better handle high Stokes numbers. We note that for Stokes numbers a few times above unity, the fluid formalism itself becomes questionable.

Using the two-fluid dust as Lagrangian particles, our results are at odds with the monofluid ones for all dust sizes. Firstly, for the smallest dust-grain sizes, the bulk and the peak of the monofluid and two-fluid density PDFs are quantitatively comparable. However, the two-fluid PDF tails extend to much higher and lower densities (more than one order of magnitude). The two-fluid particles indeed encounter the same limitation as classical tracer particles, which experience artificial trapping in the high-density regions (Cadiou et al. 2019). This leads to an apparent decoupling of the dynamics of the dust particles and the gas, which then appears as large variations in the dust-to-gas ratio. However, this decoupling is of numerical origin, and is mostly due to the finite resolution used for the gas, that is, the grid cell size is larger than the adaptive smoothing length computed from the Lagrangian particle distributions. In addition, we neglect the dust back-reaction in our two-fluid Lagrangian-particle implementation, and our monofluid results indicate that the smallest grains should not show evidence of large dust enrichment at low gas density. Secondly, for the large grains, the PDF of the density distribution does not converge between the monofluid and the two-fluid dust as Lagrangian particles. There is an important discrepancy between the low-density PDF of the two methods. Considering only the two-fluid Lagrangian-particle distribution, at high density, the PDF of the largest grains matches that of the smallest grains, with the same issue of h < Δx. We also defined a maximum dust enrichment as a function of the gas density, which is a function of the grid size. From Fig. 9, it is clear that the regions of maximum dust enrichment, that is, the regions with the highest dust density, are regulated by the grid size Δx, in particular for the largest grains.

Laibe & Price (2012, 2014a) showed that two-fluid methods are limited by both temporal and spatial resolution in the strong drag regime. Firstly, the spatial resolution should verify the criterion Δx < c_st_s in order to accurately capture the dust and gas dynamical phase separation, and to avoid artificial velocity difference damping. It is straightforward to estimate the dust-grain size s_bifluid which satisfies Δx < c_st_s using Eq. (4):

$$s_{\text{bifluid}}>\Delta x\frac{{\rho }_{\text{g}}}{{\rho }_{\text{grain}}}\sqrt{\frac{8}{\pi \gamma }},$$(19)

or, in typical conditions of molecular clouds,

$$s_{\text{bifluid}}>380\text{μm}\left(\frac{\Delta x}{1\text{\hspace{0.17em}\hspace{0.17em}pc}}\right)\left(\frac{\rho }{{10}^{-20}g{\text{cm}}^{-3}}\right){\left(\frac{{\rho }_{\text{grain}}}{1g{\text{cm}}^{-3}}\right)}^{-1},$$(20)

where Δx is the characteristic resolution length of the gas fluid. This criterion might have been overlooked in previous studies. In our experiments, we neglected the dust back-reaction, meaning that the over-damping is not problematic (Lorén-Aguilar & Bate 2015). We verified this by running some DUSTYWAVE experiments (Laibe & Price 2012) with our monofluid and two-fluid solvers with a tiny (ϵ = 10⁻⁵) dust fraction. In the strong drag regime, we find that both methods give accurate results. In addition to the spatial resolution criterion, Laibe & Price (2012, Laibe & Price 2014a) showed that the integration time-step constraint due the very short stopping time conditions can be very restrictive. In this work, we solve this issue by using an analytical update of the dust particle velocity. In addition, the characteristic velocity of the gas in supersonic turbulence experiments is larger than the sound speed, meaning that the condition in Eq. (20) relaxes towards smaller values, because Δx < υ_rmst_s. Assuming Larson scaling relations, we get

$$s_{\text{bifluid}}>84\text{μm}\left(\frac{\Delta x}{1\text{\hspace{0.17em}\hspace{0.17em}pc}}\right){\left(\frac{\rho }{{10}^{-20}g{\text{cm}}^{-3}}\right)}^{12/7}{\left(\frac{{\rho }_{\text{grain}}}{1g{\text{cm}}^{-3}}\right)}^{-1},$$(21)

For the two-fluid results, it is important to note that we have found no systematic accumulation of the particles in the regions of high gas density in the two-fluid Lagrangian-particle population. Indeed, the gas flow, being highly turbulent, changes drastically over one crossing time, such that the artificially trapped particles can be released in the low-density material. A hybrid treatment of the particle motion in the high-density regions could help to alleviate this problem by superimposing the dust-mass fluxes from the Eulerian description at small Stokes numbers onto the Lagrangian model (see e.g. Cadiou et al. 2019).

From the study we conduct here, with its limited resolution, we observe a strong dependence of the dust density PDFs and clustering degree on the number of particles and grid resolution. Indeed, the number of particles as well as the number of grid cells should increase, as should the ratio between the two, similarly to what is expected in SPH codes with the number of particles and the number of neighbours (e.g. Commerçon et al. 2008).

6.2 Comparison to previous works

Using SPH, Tricco et al. (2017) employed the same approximations as ours, that is, monofluid and diffusion approximation, and used very similar initial conditions and numerical resolution. In the context of compressible turbulence, it is worth remembering that Price & Federrath (2010) report that SPH and grid codes give roughly comparable results when the number of SPH particles is approximately equal to the number of grid cells. We can therefore directly compare our results with those of Tricco et al. (2017). The models of these latter authors use a grain density of 3 g cm⁻³, while we use 1 g cm⁻³. As a consequence, at equal Stokes number, our dust size has to be multiplied by a factor of three in order to compare with their results. Our results qualitatively agree with those of Tricco et al. (2017), who report that grains with size >10 μm decouple significantly from the gas dynamics. We find that the dust grains with size >4 μm decouple significantly (which thus corresponds to 12 μm in the experiments of Tricco et al. 2017), with large variations in the dust concentration (more than one order of magnitude). In terms of initial Stokes number, this corresponds to St ≳ 0.1. This behaviour is referred to as size-sorting in Tricco et al. (2017). We find that the size-sorting most efficiently occurs in regions of intermediate gas density (−4 < log(ρ_g/ρ_g,0) < 1), whereas the densest parts of the molecular clouds exhibit uniform dust concentration, equal to the initial dust concentration. For the small grains, with sizes <4 μm we report some variation in the local dust ratio, but these regions represent a negligible amount of mass and correspond to gas densities where ρ_g/ρ_g,0 < 1. We can therefore expect apparent variation of the dust-size distribution in these intermediate densities. However, when observed in real conditions with an integration along the line of sight, this inho-mogeneity will be smoothed out by the bulk of the mass, which has a uniform dust concentration.

In order to compare our results with those of Hopkins & Lee (2016) and Mattsson et al. (2019a,b), we introduce the dimensionless parameter α as defined in Hopkins & Lee (2016):

$$\alpha =\frac{{\rho }_{\text{grain}}s}{\langle {\rho }_{\text{gas}}\rangle L},$$(22)

where 〈ρ_gas〉 ≃ ρ₀ is the mean density of the gas. In our numerical setup, this translates to

$$\alpha \simeq 0.002\left(\frac{s}{1\text{\hspace{0.17em}μm}}\right)\simeq 0.12St.$$(23)

Hopkins & Lee (2016) and Mattsson et al. (2019a,b) report a small-scale clustering reaching its maximum for grains in the nanometer range (a few tens of nanometers), which corresponds to α ≃ 0.01−0.1 for their parameter choice. In our setup, we investigate a parameter space for which α < 0.02, corresponding to the very small grains in Hopkins & Lee (2016) and Mattsson et al. (2019a,b).

Hopkins & Lee (2016) designed their numerical experiment to investigate gas and dust dynamics in giant molecular clouds. These authors used a 256³ particle resolution for both the dust and the gas fluids, using the mesh-free GIZMO framework (Hopkins 2015). The numerical setup used by Mattsson et al. (2019a,b) is more similar to ours (gas dynamics on a fixed grid, dust treated as Lagragian superparticles), but they investigate dust grains with α > 0.1. Our results for the two-fluid dust as Lagrangian particles are globally consistent with those of Hopkins & Lee (2016) at very low α, with the largest variation of the dust ratio observed in the regions of low gas density. In our two-fluid models, the maximum dust ratio is set by the grid resolution used for the gas dynamics. In the regions of low gas density, we report a dust ratio increase of up to four orders of magnitude for the α ≃ 0.02 dust grains, while Hopkins & Lee (2016) found variations of up to six orders of magnitude for α = 0.03. Hopkins & Lee (2016) used a completely different numerical framework from ours, with adaptive resolutions of both the gas and the dust particles, meaning that they can better manage different resolutions for the gas and the dust (despite still being limited by the gas resolution for the drag). Further work at higher resolution is needed in order to investigate the behaviour of the α ≃ 0.02 dust grains in the material of low gas density with our setup. At high gas density, we observe dust-ratio variations over more that two orders of magnitude, also similar to Hopkins & Lee (2016) results. We attribute these variations to spurious effects due to trapping of Lagrangian superparticles beyond the grid scale, which gives rise to an artificial increase in the dust-fluid density. Our results therefore question these fluctuations at high density, which are typically regions where St < 1. These correspond to adaptive resolution lengths of the Lagrangian-particle distributions smaller that the grid size and where the grid size is larger than the coupling length. The effects of these spurious variations at small Stokes number on the dust ratio at low density remains unclear from our comparison. Hopkins (2014) and Hopkins & Lee (2016) suggest that the variations of stellar abundances within clusters could be a consequence of the dust dynamical coupling at the scales of giant molecular cloud. In particular, Hopkins (2014) proposes that decoupling due to dust dynamics could lead to the formation of totally metal stars (1 for 10⁴ in clusters). Our two-fluid results show that the maximum dust enrichment for the largest dust grains is regulated by the size of the grid used to interpolate the friction with the gas. On the other hand, our monofluid results show that the maximum dust enrichment for the largest grain corresponds to regions where St > 1, that is, regions where the TVA approximation is out of range. In summary, our results do not allow us to make a firm estimation of the amplitude of dust-ratio variations for dust grains with St > 1. Greater resolution is still clearly needed in order to quantify the occurrence of metal stars as suggested by Hopkins (2014). Two-fluid Eulerian dust simulations should also be used as the method does not cause artificial clustering for particles with St ≳ 1.

More recently, Mattsson et al. (2019b) performed similar experiments to ours, using a 1024³ grid resolution for the gas and 10⁶−10⁷ particles for each dust-size bin. These authors compare the properties of the two-fluid Lagrangian superparticles with the those of the tracer particles for the smallest grains they consider (α ≃ 0.001) and find that the inertia still makes a difference in the grain dynamics and small-scale clustering. Our comparison between the tracer particles and two-fluid Lagrangian superparticles shows that this is not the case for very small Stokes number, as both particle distributions are identical. More importantly, we show that independently of the dust size, the regions where St < 1 show evidence of artificial dust-particle trapping (or clustering) similar to that of the tracer particles. Our findings therefore bring into question the small-scale clustering properties of nano-dust grains reported by Mattsson et al. (2019b), in particular with respect to the minimum resolution set by the grid.

This short comparison points towards the need for a hybrid numerical method for dust-particle mesh codes in order to handle the large variety of Stokes numbers encountered locally in turbulent molecular clouds. For dust-grain species with both St > 1 in the material with low gas density and St < 1 in that with high gas density, one could use our two-fluid implementation for St > 1 to get the particle acceleration, but a Monte Carlo scheme based on the Godunov flux at cell interface for St < 1 as proposed by Cadiou et al. (2019).

6.3 Neglected physics

In this work, we studied the dynamical interplay between dust grains and gas in conditions typical of molecular clouds. We made a number of simplifications in order to focus on numerical implementation effects rather than on physical effects. The first big assumption we made concerns magnetic fields. Magnetic fields are ubiquitously observed in molecular clouds with micro-Gauss amplitudes (Hennebelle & Falgarone 2012; Pattle & Fissel 2019) and can interplay with the turbulent motions of the gas (Federrath 2016). In our setup, accounting for magnetic fields would have two implications. First, the properties of the turbulent flows change in the presence of magnetic fields, in particular in the strong field case, with anisotropies arising from the strong magnetic tension. However, this should not alter our findings as to the reliability of the monofluid and two-fluid formalisms, as long as dust grains are considered to be neutral. Hopkins & Lee (2016) report that including magnetic fields has a weak effect on the amplitudes of the gas-to-dust-ratio variations, but that the dust density can vary in regions where the gas density does not vary. Indeed, the neutral dust grains are sensitive to the gas-pressure gradients and Lorentz force of the gas through the drag. In shearing-box experiments of dust dynamics within protoplanetary discs, Lebreuilly et al. (2021) show that dust grains preferentially concentrate in regions where the pressure-gradient force and Lorentz force cancel out. Second, dust grains are expected to be charged in molecular clouds (e.g. Draine & Sutin 1987; Guillet et al. 2007), and depending on the dust-grain size, either the Lorentz force or the gas drag should regulate the dust-grain motions. The balance between the gyration time and the stopping time of the dust grains of different sizes can therefore favour dust-size sorting. Lee et al. (2017) extended the work of Hopkins & Lee (2016) to include the dynamics of charged dust grains and find that the Lorentz force suppresses the dust-ratio variations for very small grains (s ≪ 1 μm), while the large grains (s ≳ 1 μm) are insensitive to magnetic fields. These authors concluded that the effects of Lorentz forces are subdominant, and so the qualitative conclusions from their previous studies remain unchanged. However, they expect that the dynamics of charged dust grains would be more important in physical conditions other than the cold neutral medium, such as in the warm neutral medium. Exploration of the dynamics of charged dust grains with the monofluid formalism still lacks a proper derivation of the monofluid equation accounting for the Lorentz acceleration and the dust back-reaction.

The second strong assumption is on dust coagulation and fragmentation processes. We limit our physical model to dust dynamical evolution, in which only the interactions between the dust and the gas are taken into account. The supersonic motions within molecular clouds associated with the dust size sorting favour strong interactions between dust grains, in particular collisions that can lead to dust growth and fragmentation. However, we show that for dust-grain sizes typically observed in molecular clouds (<1 μm, Köhler et al. 2015), only a small mass fraction of the dust grains decouples from the gas and could therefore be undergoing grain-grain interactions. Further work should account for dust-grain interaction in order to test whether or not dust growth can already occur in molecular clouds thanks to the dynamical size sorting as suggested by previous work by Ormel et al. (2009), who report growth up to 100 μm in clouds with lifetimes larger than a freefall time.

7 Conclusion

We present a unique suite of numerical experiments designed to investigate neutral dust-grain dynamics in turbulent molecular clouds. We considered typical turbulence driving parameters, which satisfy the classical Larson relations (Larson 1981; Hennebelle & Falgarone 2012. We compared two different numerical implementations of the dust dynamics. The first one relies on the monofluid formalism and the diffusion approximation. The second implementation treats the dust fluid using Lagrangian super-particles, while the gas dynamics is handled on a Eulerian grid. For each implementation, we compared the dust-density distributions and the spatial dust-ratio variations as a function of dust-grain size. Our main findings are as follows:

• Our implementation of the monofluid formalism, which relies on the dust diffusion approximation (Price & Laibe 2015; Lebreuilly et al. 2019), is well suited to studying the dynamics of neutral dust grains with sizes 1 nm < s < 10 μm. We report dust dynamics decoupling for Stokes numbers St > 0.1, that is, dust grains of s > 4 μm in size, which matches the theoretical expectations provided in Eq. (8). Our results are in very good agreement with previous work by Tricco et al. (2017). The dust concentrates in the pressure maxima, with a higher concentration at low density (ρ < ρ₀).

• The two-fluid results with dust as Lagrangian particles are affected by numerical artefacts in both the strong and weak drag regimes as follows. For tightly coupled dust grains, with St ≪ 1, we show that the apparent large dust-to-gas ratio variations are spurious, mostly because of artificial trapping in the high-density regions. For weakly coupled dust grains in the St > 1 regime, the maximum dust enrichment we measure is strongly affected by the grid resolution, from which the drag from the gas is interpolated. These results raise questions about the robustness of dust concentration and clustering predictions measured in numerical experiments of dusty turbulence using a hybrid multifluid numerical method (gas on a grid, dust as Lagrangian particles). Our resolution study shows that results from a setup with dust as Lagrangian particles do not converge to the correct results with increasing grid and particle resolutions. Instead of deriving resolution criteria, a more fundamental study of the correctness of the Lagrangian implementation at different drag regimes should be carried out.

• We discuss the comparison with previous works and we show that there is no tension in terms of the critical size for decoupling between the results reported by Tricco et al. (2017) on one side and Hopkins & Lee (2016) and Mattsson et al. (2019a) on the other. These studies did not explore the same parameter space as suggested in footnote 13 of Hopkins et al. (2020). The controversy arose from trying to conclude on a universal dust size for decoupling, independently of the physical conditions. Our results show that these previous studies are in good agreement when the comparison is done at equivalent Stokes number or α parameter.

• If one takes typical dust grains in molecular clouds to have sizes in agreement with the standard MRN size distribution, namely a maximum dust size of <1 μm, then our results indicate that the dust-to-gas ratio should not exhibit large variations of more than a factor of two in regions of high gas density. We stress that these results are valid for (1) neutral dust grains and (2) turbulent properties following the Larson relations at parsec scales.

Acknowledgements

This work was granted access to the HPC resources of CINES (Occigen) under the allocation 2018-047247 made by GENCI. We gratefully acknowledge support from the PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon. B.C. and F.L. have received fundings from the French national agency ANR DISKBUILD ANR-20-CE49-0006. This work was supported with funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska Curie grant agreement No 823823 (RISE DUSTBUSTERS project). B.C. acknowledges financial support from "Programme National de Physique Stellaire" (PNPS) of CNRS/INSU, CEA and CNES, France. U.L. acknowledges financial support from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130).

All the figures were created using the OSYRIS visualisation package for RAMSES.

Appendix A Resolution convergence

In this Appendix, we test the resolution convergence of the monofluid models. We repeat our fiducial setup 256MRN with resolutions of 128³ and 512³. These models include ten dust species with sizes ranging from 1 nm to 20 μm, distributed according to the MRN power law.

Figure A.1 shows the gas-density and dust-ratio variations of the two extreme dust-size bins. As expected in isothermal turbulence, the higher the resolution, the stronger the density contrasts. The dust-ratio variations are also sensitive to the resolution, with larger variations in the 512MRN model. Figure A.2 portrays the PDFs of the total dust-to-gas ratio and of the dust ratios of the 1.6 nm, 0.08 μm, 0.6 μm, 4.5 μm, and 12 μm dust grains for the three resolutions. The PDFs get wider as resolution increases, except in the case of the 1.6 nm grains. Indeed, for the largest grains, the pressure gradients increase with resolution and thus the velocity shift between the dust and the gas. As a consequence, the dust grains tends to decouple more. For the smallest grains, the numerical diffusion decreases with resolution. The 1.6 nm grains should trace almost perfectly the gas since St < 1, see Fig. 3). As a consequence, as resolution increases, the coupling gets stronger.

Our results are qualitatively very similar, and are independent of the resolution. The monofluid results we report in the main part of this study are thus robust against resolution effects.

Fig. A.1 Total density (left), 1.6 nm dust-grain-ratio variations (middle), and 12 μm dust-grain-ratio variation maps in the xz-plane for different resolutions: 1283 (top), 2563 (middle) and 5123 (bottom). The dust variation is given relative to the initial dust ratio value and is shown in logarithmic scale. The red colour shows dust-ratio enhancement while blue means a dust-ratio decrease. The arrows represent the barycentric velocity (left), and the 1.6 nm (middle) and 12 μm dust-grain (right) velocity vectors in the plane. All plots are made at a time corresponding to 2tcross.

Fig. A.2 Dust-ratio PDFs for the total dust-to-gas ratio, the 0.6 nm dust grains, the 1.6 nm dust grains, the 4.5 nm dust grains, the 0.08 μm dust grains, and the 12 μm dust grains. The black, red, and blue lines show respectively the 128MRN, the 256MRN, and the 512MRN runs. The dust-ratio PDFs are normalised to the initial dust ratio of each dust species (ϵi/ϵi,0) and the PDFs of the total dust-to-gas ratio is not normalised to the initial dust ratio and thus peaks around 0.01.

Appendix B Back-reaction

Our two-fluid implementation does not include the back-reaction of the dust onto the gas. In order to compare the monofluid and two-fluid results, we neglected the back-reaction in the monofluid models presented in Sect. 5. To do so, we simply set the initial dust-to-gas ratio to 10⁻⁶ in the monofluid runs. In this Appendix, we compare the impact of the back-reaction in the monofluid models in more depth than in the main part of the manuscript. Our comparison is based on the 256MRN and on a simlar one but with an initial dust-to-gas ratio to 10⁻⁶

Figure B.1 shows the PDF of the total (gas + dust), gas, and total dust densities. First, the PDFs of the total and gas densities are almost identical, because the dust mass is always much less that the gas mass for the bulk of the material. The dust-enriched regions we report in the main text do not enclose sufficient mass to affect the global quantities. Second, all PDFs are very similar between the two runs.

Figure B.2 shows the PDFs of the total dust-to-gas ratio and of the same dust ratios as in Fig. A.2 for the two runs. Only the tails of the PDFs are effected by the back-reaction, while the bulk of the mass remains unaffected.

Fig. B.1 PDFs of the total (gas + dust) density (top), total dust density (middle), and gas density (bottom) variations for the 256MRN (fiducial, red) run and for the same model but with an initial dust-to-gas ratio to 10−6 (black). All quantities are averaged over more than one crossing time. The different densities are normalised to their initial values.

Fig. B.2 Same as Fig. A.2 for the comparison between the 256MRN run and the same one but with an initial dust-to-gas ratio of 10−6.

References

All Tables

All Figures

Fig. 1 Time evolution of the barycentre and dust-species RMS velocity (top) and of the clumping factor of the density of the different species (bottom) in the fiducial model 256MRN. The time evolution has been normalised according to the turbulent crossing time tcross. In the bottom panel, the dotted line shows the clumping factor evolution of the total dust density ρd and the dashed-dotted line shows the clumping factor evolution of each dust species concentration ϵi. The evolution of the clumping factor of the gas density (dashed) is indistinguishable from that of the barycentre (solid) density.

In the text

Fig. 2 Gas density (top-left) and dust ratio variation maps in the xz-plane for the ten dust bins in the 256MRN run. The variations are given relative to the initial value of the quantity and are shown in logarithmic scale. In the dust-ratio panels, the red colour shows dust-ratio enhancement, while blue indicates a dust-ratio decrease. The isocontours correspond to gas-density variations of −1 (orange), 0 (black), and +1 (green) in logarithmic scale. The arrows represent the barycentric velocity (top-left), and the dust-grain drift velocity vectors in the plane for all the other plots. All plots are made at a time corresponding to 3tcross.

In the text

Fig. 3 2D histograms of the dust-ratio variations as a function of gas density for the fiducial model 256MRN. The top-left panel shows the variation in total dust-to-gas ratio ɛ, while the other panels show the dust-ratio variations of each dust species. The colour coding indicates the mass of the gas in log scale (in solar units). The horizontal grey dotted line separates the dust-enriched and dust-depleted regions. The vertical line indicates St = 1 (solid) for each dust species. All quantities are averaged over more than one dynamical time tcross and the variations are normalised to the initial value of each dust ratio.

In the text

	Fig. 4 Volume-weighted cumulative PDF of the dust ratio variation for the 1.6 nm, 0.03, 1.7, 4.5, and 12 μm dust grains in the fiducial model 256MRN. All quantities are averaged over more than one dynamical time tcross.
In the text

Fig. 5 Profile of the total density ρ (blue) and the drift velocity (orange) and dust ratio (dotted green) for the 12 μm dust grains in the fiducial model 256MRN. The density is normalised by its initial value. The norm of the drift velocity with respect to the barycentre is normalised to 1 km s−1. The variations in dust ratio are computed with respect to its initial value. The quantities are plotted against the x-axis at a random time (>2tcross) and at a random location in the yz-plane. Regions with positive relative dust ratio correspond to dust-enriched material.

In the text

Fig. 6 Cumulative PDF of the Stokes number for the 1.6 nm dust grains and the 0.03, 1.7, 4.5, and 12 μm dust grains in the fiducial model 256MRN. The coloured dotted lines show the mass-weighted PDF and the coloured solid lines show the volume-weighted PDFs. The Stokes number is computed using the local gas density and the global turbulent crossing time turb. All quantities are averaged over more than one dynamical time tcross. The vertical lines indicate the St = 1 (solid) and St = 0.1 (dashed).

In the text

Fig. 7 Density slice in the z = 0 plane obtained from the 256_2F run for the 1 nm (top) and 10 μm (bottom) dust grains at a time corresponding to ≃ 4.5tcross. The left column represents the density field of the monofluid dust, while the right column represents the density field interpolated from the Lagrangian dust particles. The density is normalised by the initial density ρ0. The black lines represent the contour for which the gas density is equal to its initial value ρ0.

In the text

	Fig. 8 PDFs of the logarithm of density ln ρ in the 256_2F run for the 1 nm (top) and 10 μm (bottom) dust grains. The density variations are normalised to the initial density values for each field.
In the text

Fig. 9 Same as Fig. 3 but for the 256_2F run. The left column shows the monofluid results and the right column shows the results when dust is modelled as Lagrangian particles. The top row shows the 1 nm dust grains and the bottom row the 10 μm dust grains. The red dashed line indicates the maximum dust enrichment as defined in Eq. (17). All the material that is above this line has an excessive dust concentration.

In the text

	Fig. 10 PDF of the effective adaptive resolution length h of the 1 nm (dotted lines) and 10 μm (solid lines) dust-grain fluid particle distributions for the 256_2F run. The colour coding indicates the particle distribution (two-fluid dust as Lagrangian particles in blue, tracer in red). The vertical line indicates the numerical resolution of the grid (2563).
In the text

Fig. 11 PDF of the dust density of the 1 nm (left) and 10 μm (middle). The colour coding indicates the dust fluid treatment (two fluid Lagrangian particles in blue, tracer particles in red, monofluid in black). For clarity, only the monofluid quantities of the 256_2F are shown. The solid lines show the results of the 256_2F run, while the dotted lines show the results of the 128_2F. Right: PDF of the adaptive length of the two-fluid dust Lagrangian-particle distributions (solid lines for the 10 μm dust grains, dotted line for the 1 nm) for the 256_2F (blue) and 128_2F (red) models. The vertical dashed lines indicate the grid numerical resolution (black for 256_2F, grey fro 128_2F).

In the text

	Fig. 12 Same as Fig. 11 for the 128_2F, 128_2F_LR, and 128_2F_VLR models. In the density PDFs, the blue lines indicate the results for the two fluid dust Lagrangian-particle distribution while the red lines indicate the tracer particle ones.
In the text

Fig. A.1 Total density (left), 1.6 nm dust-grain-ratio variations (middle), and 12 μm dust-grain-ratio variation maps in the xz-plane for different resolutions: 1283 (top), 2563 (middle) and 5123 (bottom). The dust variation is given relative to the initial dust ratio value and is shown in logarithmic scale. The red colour shows dust-ratio enhancement while blue means a dust-ratio decrease. The arrows represent the barycentric velocity (left), and the 1.6 nm (middle) and 12 μm dust-grain (right) velocity vectors in the plane. All plots are made at a time corresponding to 2tcross.

In the text

Fig. A.2 Dust-ratio PDFs for the total dust-to-gas ratio, the 0.6 nm dust grains, the 1.6 nm dust grains, the 4.5 nm dust grains, the 0.08 μm dust grains, and the 12 μm dust grains. The black, red, and blue lines show respectively the 128MRN, the 256MRN, and the 512MRN runs. The dust-ratio PDFs are normalised to the initial dust ratio of each dust species (ϵi/ϵi,0) and the PDFs of the total dust-to-gas ratio is not normalised to the initial dust ratio and thus peaks around 0.01.

In the text

	Fig. B.1 PDFs of the total (gas + dust) density (top), total dust density (middle), and gas density (bottom) variations for the 256MRN (fiducial, red) run and for the same model but with an initial dust-to-gas ratio to 10−6 (black). All quantities are averaged over more than one crossing time. The different densities are normalised to their initial values.
In the text

	Fig. B.2 Same as Fig. A.2 for the comparison between the 256MRN run and the same one but with an initial dust-to-gas ratio of 10−6.
In the text