© The Authors 2025

1 Introduction

Interstellar dust is of great importance in many astrophysical environments in regard to the thermodynamics and chemistry of its surrounding medium. In particular, dust grains are essential ingredients for star and planet formation. The wide spectrum of dust grain sizes affects the optical properties of the medium and, consequently, the heating and cooling processes at work (McKee & Ostriker 2007). In addition, dust is a very useful observational probe in the interstellar medium (ISM) and more generally in the context of star and planet formation (Draine 2004). Dust grains affect the ionization degree of the medium and can become the main charge carriers during protostellar collapse (Zhao et al. 2016; Tsukamoto & Okuzumi 2022) as they allow for electron and ion recombination and electron capture on their surface (Ivlev et al. 2015; Marchand et al. 2016). As such, they control the degree of coupling between the gas and the magnetic field (Nakano et al. 2002). All of the above processes are sensitively dependent on local dust concentration levels and on dust size, which is modeled as a power-law Mathis-Rumple-Nordsiek (MRN) distribution in the diffuse ISM (Mathis et al. 1977). However, the dust size distribution is expected to evolve as dust grains undergo coagulation and fragmentation processes (Dominik & Tielens 1997; Ormel et al. 2009). Following such an evolution of the dust distribution is very challenging, but essential. Dust growth is also inherently related to the formation of larger solid bodies such as planetesimal and later on planets (presumably in less than 1 Myr, see Testi et al. 2014).

Dust growth is very sensitive to the amount of dusty material available, that is to dust density. Although it is widely accepted that dust represents only 1% of the gas mass in the ISM and in protostellar dense cores, this average value holds on the large scale but could vary on smaller scales due to various mechanisms such as gas-dust differential dynamics (Lebreuilly et al. 2020), turbulence and (magneto)-hydrodynamical (MHD) instabilities (see Hendrix & Keppens 2014; Squire & Hopkins 2018, respectively for Kelvin-Helmholtz instability and resonant drag instabilities). A number of works have recently studied such concentration of solid particles, unveiling up to orders of magnitude in dust-to-gas ratio fluctuations resulting from aerodynamical drag. Hydrodynamical dust clumping by supersonic turbulence has been studied numerically in an environment reproducing molecular cloud physical conditions (Hopkins & Lee 2016; Tricco et al. 2017; Commerçon et al. 2023). In the latter work, they compare two numerical implementations, namely a monofluid Eulerian approach with terminal velocity approximation for the dust grains, and a bifluid approach treating the dust grains as Lagrangian superparticles. They show that the terminal velocity approximation is well suited for dense regions but dust decoupling and enrichment is not well described for grains of a Stokes number above unity (in regions of low density). In the same Stokes number regime, the bifluid approach capacity to describe dust enrichment is limited by the grid resolution used for the gas. For very well coupled dust grains (St << 1), the Lagrangian treatment of the dust leads to artificial trapping in the densest regions (Price & Federrath 2010; Cadiou et al. 2019). Large grains tend to clump more efficiently as they are less coupled to gas molecules (Mattsson et al. 2019). With both an analytical and numerical approach, Hopkins & Squire (2018) expanded their work on the resonant drag instability (RDI) by including the role of magnetic fields and charged dust. As dust and gas drift relative to one another (via a different force balance including, for example, gravity for both and pressure gradients for the gas only), instabilities develop for all mode wavelengths, leading to strong fluctuations in the dynamics of both fluids. In later phases of the star formation cycle, such as protoplanetary disks, a number of mechanisms and instabilities have been suggested and studied to significantly concentrate the dust locally (Xu & Bai 2022; Lehmann & Lin 2023; Birnstiel 2024), including the so-called streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Drazkowska 2017; Krapp et al. 2019; Li & Youdin 2021; Carrera et al. 2025, and references therein) and the dust settling instability (part of the RDI family, see Squire & Hopkins 2018; Krapp et al. 2020).

Including magnetic fields, Lee et al. (2017); Moseley et al. (2023) and Moseley & Teyssier (2025) have investigated dust concentrations in molecular clouds considering charged dust as Lagragian particles feeling Lorentz forces. In such diffuse environments, most of the charges are mobile ions and electrons and the authors assumed them to be moving as one with gas molecules, therefore solving the induction equation in the ideal regime to compute the magnetic field evolution. Despite of considering charged dust grains, magnetic field evolution is dictated by gas motions, i.e., field lines are attached to the gas fluid in the induction equation. They initially found very large dust-to-gas ratio fluctuations being curtailed to some extent for small particles by Lorentz forces, which offer an extra coupling between gas and dust, while large grains are less affected due to their lower magnetization.

However, to our knowledge, studies of dust clumping with charged dust and beyond-standard MHD models in denser environments of the ISM such as protostellar envelopes and dense cores have never been undertaken. At such densities (>10⁴ cm⁻³), small dust grains are to be considered as main charge carriers (Nishi et al. 1991; Nakano et al. 2002; Zhao et al. 2016) as they collect electrons and recombine positive and negative charges on their surface. More sophisticated MHD models and chemical networks have to be used to model dust and gas dynamics correctly. For this work, we perform 1D simulations of a two-fluid (magnetized dust and neutral gas) mixture in a turbulent box, using a lower computational cost due to our lack of multidimensionality to focus on microphysical processes and clumping on a very small scale. In Sect. 2, we present the methods and MHD equations solved. Numerical setups are described in Sect. 3. In Sect. 4, we first consider a single nonlinear torsional Alfvén wave, emphasizing the role of a parametric-like instability to produce compressible modes (note that we let numerical noise serve as perturbations). Then, we tested the robustness of the presented clumping mechanism in a more realistic environment including 1D turbulence driven on large scales. We discuss those results in Sect. 5. Finally, we summarize our findings in Sect. 6 and provide additional material in the appendix.

2 Methods: Multifluid implementation

We use the code shark (Lebreuilly et al. 2023) to run our simulations. We consider a 1D box with a uniform grid and periodic boundary conditions which extends along the x-axis. We make a full treatment of the gas and dust dynamics which interact through drag terms. The gas evolves isothermally and the dust fluid is considered pressureless. The box is threaded by a uniform magnetic field. We allow transverse (y and z-wise) components to exist.

2.1 Dusty non-ideal MHD

Considering neutral gas and a single charged dust fluid, the magnetohydrodynamical equations are $$\begin{array}{rl}\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [\rho \mathit{v}]=& 0,\\ \frac{\mathrm{\partial }{\rho }_{\mathrm{d}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left[{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}\right]=& 0,\\ \frac{\mathrm{\partial }\rho \mathit{v}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [\rho \mathit{v}\mathit{v}+P\mathbb{I}]=& \frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v})+{\mathit{F}}_{\mathrm{i}\to \mathrm{N}},\\ \frac{\mathrm{\partial }{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left[{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}\right]=& -\frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v})+{\mathit{F}}_{\mathbf{L}\mathbf{o}\mathbf{r},\mathrm{d}},\end{array}$$(1)

where the first term on the right hand side of the momentum equations is the hydrodynamical drag between gas and dust. ρ and ρ_d refer respectively to the gas and dust density. ${\mathit{F}}_{\text{Lor,}\mathrm{d}}=n_{\mathrm{d}}{Z}_{\mathrm{d}}e(\mathit{E}+\frac{{\mathit{v}}_{\mathrm{d}}}{c}\times \mathit{B})$ is the Lorentz force felt by the dust fluid. n_d and Z_d < 0 are the dust numerical density and the number of electrical charges carried by dust grains. e is the fundamental electrical charge, c the speed of light and E the electric field. We also include ions, whose inertia is neglected but we ignore the presence of electrons for the sake of simplicity. We consider the friction exerted by ions on neutrals to be ${\mathit{F}}_{\mathbf{i}\to \mathbf{N}}={\rho }_{\mathrm{i}}\rho {\gamma }_{\mathrm{i},\mathrm{g}}({\mathit{v}}_{\mathbf{i}}-\mathit{v})$. The ions momentum equation is $$s_{\mathrm{d}}\ge 1\mu \mathrm{m}$$(2)

where n_i is the ion numerical density.

We use the momentum equations of the ions along with electroneutrality condition to infer a generalized Ohm's law, which gives us the electric field (the complete derivation is provided in the appendix, Appendix A). We then inject the electric field E in the induction equation to get $$\begin{array}{rl}\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }t}-\mathrm{\nabla }\times [{\mathit{v}}_{\mathrm{d}}\times \mathit{B}-\frac{\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}}({\mathit{v}}_{\mathbf{d}}-\mathit{v})]=& \mathrm{\nabla }\times [\frac{c}{e{n}_{\mathrm{i}}4\pi }(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}]\\ & -\mathrm{\nabla }\times [\frac{c\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}e{n}_{\mathrm{i}}4\pi }(\mathrm{\nabla }\times \mathit{B})],\end{array}$$(3)

where the first term on the right hand side takes the form of a dispersive Hall term and the second of a dissipative Ohm term. The former is treated as a conservative term while the latter is treated as a diffusion equation. More information along with a test of our solver for the treatment of diffusion equations are provided in Appendix D. Note also that the Hall effect makes almost no difference in our simulations. We define the ion Hall factor ${\mathrm{\Gamma }}_{\mathrm{i}}={\omega }_i{t}_{\mathrm{s},\mathrm{i}}$ to be the product of the ion stopping time with the gyromagnetic frequency ${\omega }_i=\frac{e\parallel \mathit{B}\parallel }{c{m}_i}$. The ion stopping time $t_{\mathrm{s},\mathrm{i}}={\left(\rho {\gamma }_{\mathrm{i},\mathrm{g}}\right)}^{-1}$ is given in Pinto & Galli (2008) and Marchand et al. (2016) for HCO⁺ ions. As mentioned earlier, our box is threaded by a uniform background magnetic field B_x = constant. Within this 1D configuration, the divergencefree condition $\mathrm{\nabla }\cdot \mathit{B}={\mathrm{\partial }}_x{B}_x=0$ is naturally complied and no constrained transport scheme is needed.

We also inject the electric field E in the dust Lorentz force and obtain after a few calculation steps (Appendix A) $${\mathit{F}}_{\mathbf{L}\mathbf{o}\mathbf{r},\mathbf{d}}=\frac{{\mathit{J}}_{\mathbf{t}\mathbf{o}\mathbf{t}}}{c}\times \mathit{B}-\frac{\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}{\mathit{J}}_{\mathbf{t}\mathbf{o}\mathbf{t}}-\frac{e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}(\mathit{v}-{\mathit{v}}_{\mathrm{d}}),$$(4)

where ${\mathit{J}}_{\text{tot}}={\mathit{J}}_{\mathbf{d}}+{\mathit{J}}_{\mathbf{i}}$ is the total electric current (contributions from dust and ions) related to the magnetic field as per MaxwellAmpère equation $${\mathit{J}}_{\text{tot}}=\frac{c}{4\pi }\mathrm{\nabla }\times \mathit{B}.$$(5)

The first term of Eq. (4) has to be rewritten in conservative form and moved to the left-hand part of the dust momentum equation to ensure stability of the dust Riemann solver. The third term behaves as an additional friction term between gas and dust coming from the presence of ions. This yields $$\begin{array}{rl}\frac{\mathrm{\partial }{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathit{d}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [& {\rho }_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}+\frac{B^2}{8\pi }-\frac{\mathit{B}\mathit{B}}{4\pi }]\\ & =-\frac{\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}{\mathit{J}}_{\mathbf{t}\mathbf{o}\mathbf{t}}-\frac{e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}(\mathit{v}-{\mathit{v}}_{\mathrm{d}})+\frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}(\mathit{v}-{\mathit{v}}_{\mathrm{d}}).\end{array}$$(6)

Similarly for the gas, we use Eq. (2) to substitute ${\mathit{F}}_{\mathbf{i}\to \mathbf{N}}$ with $e{n}_{\mathrm{i}}(\mathit{E}+\frac{v_{\mathrm{i}}}{c}\times \mathit{B})$ and get $$\begin{array}{rl}& \frac{\mathrm{\partial }\rho \mathit{v}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [\rho \mathit{v}\mathit{v}+P\mathbb{I}]\\ & =\frac{\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}{\mathit{J}}_{\mathbf{t}\mathbf{o}\mathbf{t}}-\frac{e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}c}({\mathit{v}}_{\mathrm{d}}-\mathit{v})+\frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v}).\end{array}$$(7)

Both drag terms between gas and dust are treated implicitly using the solver presented in Krapp & Benítez-Llambay (2020). We assume the mean free path of gas molecules to be large compared with the dust grain size. This is known as the Epstein regime (Epstein 1924) for which the stopping time of the grain is expressed as $$[t_{\mathrm{s},\mathrm{d}}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{\mathrm{g}\mathrm{r}}{s}_{\mathrm{d}}}{\rho c_s},$$(8)

where ρ_gr=2.3 g cm⁻³ is the internal density of a (spherical) dust grain and s_d its size (diameter). c_s is the isothermal sound speed. t_s, d is the characteristic time needed for a dust grain to adjust its velocity to a change in gas velocity. We define the dimensionless Stokes number St as $$St=t_{\mathrm{s},\mathrm{d}}/t_{\mathrm{d}\mathrm{y}\mathrm{n}}=t_{\mathrm{s},\mathrm{d}}{c}_s/L.$$(9)

The dynamical timescale of interest t_dyn is taken to be the time required for a sound wave to cross our numerical box of length L. St=St₀ ρ₀ / ρ where St₀ is the initial Stokes number (associated with the background state of initial uniform gas density ρ₀ and scale L).

2.2 Dusty ideal MHD: Perfect coupling between dust and magnetic field

If we neglect ion-neutral friction and have the ion Hall factor ${\mathrm{\Gamma }}_{\text{ion}}⟶\mathrm{\infty }$, as well as drop the Hall term in the induction equation, the latter takes the classic form encountered in ideal MHD, however with the dust velocity in place of the gas velocity $$\frac{\mathrm{\partial }\mathit{B}}{\mathrm{\partial }t}-\mathrm{\nabla }\times [{\mathit{v}}_{\mathit{d}}\times \mathit{B}]=0.$$(10)

The dust fluid is then considered as the sole charge carrier and is frozen to the magnetic field lines. The gas and dust momentum equations become $$\begin{array}{rlr}\frac{\mathrm{\partial }\rho \mathit{v}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [\rho \mathit{v}\mathit{v}+P\mathbb{I}]& =& \frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v})\\ \frac{\mathrm{\partial }{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathit{d}}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot [{\rho }_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}{\mathit{v}}_{\mathrm{d}}+\frac{B^2}{8\pi }-\frac{\mathit{B}\mathit{B}}{4\pi }]& =& -\frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v})\end{array}$$(11)

This approach is used to construct the dispersion relation depicted in Appendix B.

2.3 Ionization

In order to compute the charge abundances and stay consistent with the assumptions made within our MHD model, we need to neglect the presence of electrons in our chemical network (in the electroneutrality condition). To avoid further complicating the physics, we use a very simple approach where the abundance of ions (HCO⁺) is given by a balance between cosmic-ray ionization (assuming an ionization rate of ζ = 1 × 10⁻¹⁷ s⁻¹) and two-body recombination of charged particles (Elmegreen 1987; Shu et al. 1987): $$x_i=\frac{n_{\mathrm{i}}}{n_{\mathrm{H}}}={10}^{-7}{\left(\frac{n_{\mathrm{H}}}{{10}^3{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-1/2},$$(12)

and the electroneutrality condition is $$n_{\mathrm{i}}+n_{\mathrm{d}}{Z}_{\mathrm{d}}=0.$$(13)

3 Numerical setup

3.1 Initial conditions

Our setup is intended for the exploration of dust clumping in a magnetized turbulent environment, such as protostellar dense cores and envelopes. Our 1D box extends only in the x direction, but we allow transverse components (y and z-wise) of magnetic field and velocities to exist as well. The simulations are initiated with uniform gas density ρ₀ and dust density such that ρ_d, 0=0.01 ρ₀ (i.e., a dust-to-gas ratio ε=0.01), a uniform x-wise background magnetic field inferred from the plasma parameter $\beta =c_{\mathrm{s}}^2/c_{\mathrm{a}}^2$ (c_a being the Alfvén velocity), which gives $$B_x=B_{\parallel }=c_{\mathrm{s}}\sqrt{4\pi {\rho }_0/\beta },$$(14)

and a constant isothermal sound-speed $c_{\mathrm{s}}=\sqrt{\frac{k_{\mathrm{B}}T}{\mu m_{\mathrm{H}}}}$ with a temperature T=10 K (because gas is rapidly cooling). μ=2.3 is the mean molecular weight of the gas molecules in terms of units of hydrogen mass m_H. Introducing the gas numerical density n_H=ρ / μm_H, the box length L is taken equal to the Jeans length which is defined as the distance crossed by a sound wave in a free-fall time $$L=c_s{t}_{\mathrm{f}\mathrm{f}}=\frac{c_s}{\sqrt{\mathcal{G}\mu m_{\mathrm{H}}{n}_{\mathrm{H}}}}.$$(15)

Fig. 1 Maximum dust density fluctuations and maximum transverse magnetic field as a function of time upon propagation of a circularly polarized Alfvén wave, for different values of wavelength (as a fraction of the box length L) and wave amplitudes δB within the dusty ideal MHD regime and dusty non-ideal MHD regime (see Sects. 2.2 and 2.1). The early sharp rise in dust density is due to a mechanism similar to the parametric instability (Del Zanna et al. 2001). Note that growth rates measured here compare well with analytical predictions of the standard parametric instability (see Appendix C). Afterwards, the dust density decreases as a consequence of gas dust friction. The gas density is not displayed because subject to negligible variations. Parameters: nH=106 cm−3, sd=10 μm(St = 0.01) and β=0.1.

3.2 Driven turbulence

We drive longitudinal (x-wise) and transverse (y and z-wise) turbulent motions in the gas by adding an acceleration kick at each timestep, composed of 10 sinusoidal modes of wavenumber lying in the range $L/3\le k_{\text{mode,}\mathrm{i}}\le L/2$ (i.e., the turbulence is injected on a scale which is a fraction of the Jeans length) $$F_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k},\mathrm{x}\mathrm{y}\mathrm{z}}=\rho \sum _{\mathrm{i}=1}^{10}{A}_{\mathrm{x}\mathrm{y}\mathrm{z},\mathrm{i}}\sin (k_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e},\mathrm{i}}x+{\varphi }_{\mathrm{x}\mathrm{y}\mathrm{z},\mathrm{i}}).$$(16)

The amplitude A_i of each mode is generated randomly from a Gaussian distribution and we allow the phase to randomly (from uniform distribution) drift over a turnover time taken to be equal to the sound crossing time, that is t_turnover=t_cross=L/c_s. We measure the Mach number through the volume-averaged root mean squared (RMS) velocity $v_{\text{RMS}}^2={\int }_0^L(\mathit{v}-\overline{\mathit{v}}{)}^2\mathrm{d}x/{\int }_0^L\mathrm{d}x$, (with $\overline{\mathit{v}}={\int }_0^L\mathit{v}\mathrm{d}x/{\int }_0^L\mathrm{d}x$) and calibrate the modes amplitude with a correcting factor to approach equipartition and to vary the Mach number when needed. We do this for each component, which leads to the total Mach number $\mathcal{M}=v_{\mathrm{R}\mathrm{M}\mathrm{S}}/c_{\mathrm{s}}$ where $v_{\text{RMS}}^2={\int }_0^L{(v_x-\overline{v_x})}^2\mathrm{d}x/{\int }_0^L\mathrm{d}x+{\int }_0^L{(v_y-\overline{v_y})}^2\mathrm{d}x/{\int }_0^L\mathrm{d}x+{\int }_0^L{(v_z-\overline{v_z})}^2\mathrm{d}x/{\int }_0^L\mathrm{d}x=v_{\mathrm{x},\text{RMS}}^2+v_{\mathrm{y},\text{RMS}}^2+v_{\mathrm{z},\text{RMS}}^2$. It follows that ${\mathcal{M}}^2=(v_{\mathrm{x},\mathrm{R}\mathrm{M}\mathrm{S}}^2+v_{\mathrm{y},\mathrm{R}\mathrm{M}\mathrm{S}}^2+v_{\mathrm{z},\mathrm{R}\mathrm{M}\mathrm{S}}^2)/c_{\mathrm{s}}={\mathcal{M}}_{\parallel }{}^2+{\mathcal{M}}_{\perp }{}^2$ where we denote ${\mathcal{M}}_{\parallel }={\mathcal{M}}_x$ the Mach number in the longitudinal direction (x-wise) and ${\mathcal{M}}_{\perp }=\sqrt{{\mathcal{M}}_y^2+{\mathcal{M}}_z^2}$ the transverse one, perpendicular to the background uniform magnetic field B_x. Both the dust and gas are initially at rest.

4 Results

In this section, we investigate dust and gas clumping. We first provide a known and documented framework without turbulence in order to give a simple description of the clumping mechanism observed in our simulations. The remainder of this section is then dedicated to the influence of previously mentioned parameters in simulations with turbulence and to the question whether dust clumping endures in such a turbulent environment.

4.1 Dust clumping due to a parametric-like instability

For the sake of simplicity, we initiate a nonlinear circularly polarized Alfvén wave of the form $$B_y={\delta }_B{B}_x\sin \left(\frac{2\pi }{\lambda }x\right)B_z={\delta }_B{B}_x\cos \left(\frac{2\pi }{\lambda }x\right),$$(17)

where the amplitude of the wave is taken to be a fraction δ_B of the background magnetic field B_x. We find the so-called parametric instability (Del Zanna et al. 2001; Hennebelle & Passot 2006) to likely be the underlying mechanism responsible for the magnetically induced clumping of the dust. Figure 1 shows dust maximum density fluctuations as a function of time upon propagation of the Alfvén wave, for different values of wavelength (taken to be a fraction of the box length, either L / 10, L / 50 or L / 100) within both the dusty ideal MHD regime (presented in Sect. 2.2) and the dusty non-ideal MHD regime (presented in Sect. 2.1). Note that perturbations of the background state come from numerical noise.

The stability of Alfvén waves was first studied long ago by Galeev & Oraevskii (1963) and then by many authors. Later, the parametric instability was investigated numerically in the context of a single fluid in the ideal MHD regime with 1D to 3D simulations by Del Zanna et al. (2001) who showed that a circularly polarized mother Alfvén wave (which is an incompressible mode since its magnetic pressure is constant) is unstable and couples to density and compressible fluctuations in the plasma provided that its amplitude is large enough (nonlinear regime). Under such circumstances, the wave decays and gives rise to a whole spectrum of “daughter” waves, namely Alfvén modes and compressible modes (prone to nonlinear steepening and shock dissipation) propagating in both directions. Such compressible modes feed on the energy of the mother wave and lead to remarkable density fluctuations, sensitive to the frequency and amplitude of the mother wave.

In the context of our two-fluid work, triggering the instability is much more challenging since the magnetized species is the dust, which is altogether 100 times less massive than neutral gas (ε=1%). Due to hydrodynamical drag, dust gyro-motions and propagation of Alfvén waves are altered, and are even completely curtailed for certain frequencies (this very dissipative regime is described in the Appendix Appendix B within the dusty ideal MHD setup). When considering dust inertia, propagation is recovered at sufficiently high frequencies (Hennebelle & Lebreuilly 2023). Looking at Fig. 1 in the dusty ideal regime, we see an early sharp increase in dust density followed by a smoother decrease due to hydrodynamical drag (gas is almost not affected for it is not magnetized). In addition, there is a clear dependency of dust density fluctuations on wavelength. Dust size was chosen to give a Stokes number of St=0.01, meaning that dust should strongly decouple hydrodynamically from the gas on scales smaller than a hundredth of the box length ($\lambda \le L/100$). Indeed, while dust density fluctuations drop below 10^0.75=5.62 over a few Alfvén crossing times for $\lambda \ge L/10$, they remain much higher for shorter wavelengths, showing that in this case dust grains are free to concentrate locally. We point out that in the dusty ideal regime, strong dust density enhancements are observed (${\rho }_{\mathrm{d},max}/{\rho }_{\mathrm{d},0}\ge {10}^{1.5}=31.6$) in spite of a modest wave amplitude of δ_B=0.2. Once the decoupling scale has been reached, further reducing the wavelength only leads to a modest increase in maximum dust density as can be seen when comparing λ=L / 50 to λ=L / 100. In addition, the decrease in maximum transverse magnetic field indicates that the Alfvén wave is dissipated due to dust-gas collisions, with a more severe dissipation for large wavelengths.

In addition to hydromagnetic drag, dissipative non-ideal MHD effects (see Sect. 2.1) further complicate the development of such an instability. Considering now the non-ideal MHD regime in Fig. 1, it is clear again that reducing the wavelength leads to a higher saturation density as hydrodynamical friction is reduced on smaller scales (when wave frequency starts exceeding the collision frequency). However, magnetic drag will persist even on small scales since the ion stopping time t_{s, i} is very low compared to that of dust (see Sect. 5.1), leading to overall dust density enhancements lower than in the dusty ideal MHD regime (comparing both models with δ_B=0.2). In addition, increasing the wave amplitude δ_B does yield a higher level of dust concentration, but essentially only at early times (t < 1) before magnetic drag disrupts dust clumps. A key ingredient to the development of the parametric instability is the transverse-to-longitudinal magnetic field ratio $B_{\perp }/B_{\parallel }$ which is here directly controlled by δ_B (in simulations with turbulence, the parameters controlling this ratio are depicted in Fig. 3). Similarly, we see that again the transverse magnetic field intensity drops more steeply when considering larger wavelengths. Within the dusty non-ideal regime, the transverse magnetic field approaches a final value that is lower than that within the ideal regime even for a wavelength as small as λ=L / 100 confirming that dissipative effects persist. Noteworthy, we show in Appendix C that growth rates measured in simulations compare well with analytical predictions of the standard parametric instability. To summarize, although partially curtailed due to non-ideal MHD effects, dust clumping due to a parametric-like instability persists to some extent under specific conditions, namely high wave amplitude δ_B and low wavelength.

In that regard, turbulence plays an essential role as it allows for the continuous and sustained production of Alfvén waves that span a large spectrum of wavelengths. Unstable Alfvén waves are rejuvenated over time and small wavelengths (inducing strong clumping) are reached, although the corresponding amplitudes of such modes cannot be arbitrary large due to energy decreasing with wavenumber (if Kolmogorov spectrum, energy cascades as $E(k)\propto k^{-\frac{5}{3}}$, see Kim & Ryu 2005). Note however that because of the intermittent nature of the turbulence, energetic but localized events occur over very small scales in shocks, which is exactly where dust clumps form, as we shall see. Note however that our 1D simulations cannot produce all the important features of a 3D turbulent cascade.

In this work (when turbulent driving is on), the transverse-to-longitudinal ratio $B_{\perp }/B_{\parallel }$ is controlled by two parameters: the plasma parameter β and the transverse Mach number ${\mathcal{M}}_{\perp }$. On one hand, increasing the former to a certain extent (in such a way as to keep a decent level of magnetization) reduces the background magnetic field $B_{\parallel }$ and thus increases the said ratio. On the other hand, a larger transverse Mach number implies higher transverse dust velocities which in turn (dust being the charged species) produce a stronger transverse magnetic field by bending and distorting $B_{\parallel }$. Another point has to be made regarding the impact of the plasma parameter. From theory and simulations of a single fluid found in the literature, β is known to suppress the parametric instability for values $\beta \ge 1$ as thermal pressure begins to prevail over magnetic pressure, enabling sound waves to smooth out density fluctuations. However in this work, the magnetized fluid (dust) is pressureless and thus not subject to this process. We again stress that here, low values of β are stabilizing while high values favor larger $B_{\perp }/B_{\parallel }$ and thus strong clumping of the dust. An informative illustration is provided in Fig. 3.

We now turn to turbulent simulations and look at Fig. 2 to follow the formation of dust clumps. This figure depicts various fields (dust and gas densities and velocities) at three different times; at t=0.23, t=0.44 and t=0.82. At first (t=0.23), dust and gas are tightly coupled, as seen in the nearly perfectly overlain profiles (in density and velocity). Shocks form in the gas and subsequently in the dust, giving rise to over-densities. In these high density regions, assuming the relevant Alfvén velocity to be that of the gas (valid in the strong coupling regime, see Eq. (B.6), the said velocity is lowered (scales as $c_a\propto {\rho }^{-\frac{1}{2}}$) making it easier for the dust to be super-Alfvénic and consequently to distort magnetic field lines. Indeed, it is clear that the dust Alfvénic Mach numbers are significantly over unity there. In addition, we see that high dust Alfvénic Mach numbers spatially coincide with sharp magnetic pressure gradients, suggesting that field lines have been bent. Later, at t=0.44 and t=0.82, we see that prominent dust clumps have indeed formed in these regions. To summarize, turbulence in the gas aids in triggering what looks like a dust parametric instability in shocks (high density regions) where Alfvén velocity is lower and high transverse-to-longitudinal magnetic field ratio can arise as charged dust grains efficiently distort the background magnetic field. Long-lived high density dust clumps emerge owing to the absence of dust internal pressure (see Appendix E). Driven turbulent motions sustain the production of such clumps whose survival is imperilled by drag effects, especially magnetic drag which recouples both fluids. There is a competition between this instability and magnetic drag.

4.2 Impact of parameters on dust and gas density fluctuations in turbulent simulations

In this section, we drive turbulent motions and choose to vary four relevant parameters, namely the plasma parameter β, the gas initial numerical density n_H=ρ / μ m_H, the dust grain size s_d (corresponding to a given St) as well as the parallel and perpendicular Mach numbers ${\mathcal{M}}_{\parallel }$ and ${\mathcal{M}}_{\perp }$ (which we vary jointly). Table 1 shows the values taken for each parameter and in bold text the fiducial ones.

We point out that those are common astrophysical parameters that facilitate discussion and applications of our results to astrophysical environments. Nevertheless, based on Sect. 4.1, it is worth keeping in mind that their influence on dust clumping can be explained through their impact on two essential ingredients: the Stokes number of dust grains St and the transverse-to-longitudinal magnetic ratio $B_{\perp }/B_{\parallel }$. The latter controls the intensity of density fluctuations via the parametric-like instability, while the former encodes the degree of (hydrodynamical) coupling between gas and dust, crucial for the formation and survival of dust clumps.

For clarity, we provide at the end of this section in Table 2 the equivalent St and maximum (time-averaged) B_⊥ / B_∥ for each model (including those investigating the impact of a lower ion abundance; see Fig. 7). Other useful data are gathered, including the time-averaged dust clumping fraction C_{f, 10} (defined as the fraction of dust mass whose density has increased by a factor of more than 10), along with the maximum and mean dust density fluctuation (spatially averaged first as in Eq. (18) and then time averaged over steady-state regime) and the dust clumping factor $⟨{\rho }_{\mathrm{d}}^2⟩/{⟨{\rho }_{\mathrm{d}}⟩}^2$.

1. The plasma parameter β takes values from 0.1 to 1, corresponding to the Alfvén velocity being 10 times the sound speed for the former (i.e., strong magnetization), and both being equal for the latter. For a gas density n_H=10⁴ cm⁻³, values in the range β=[0.1, 0.3, 0.7, 1] translate into the following range of magnetic field intensity B_x=[41.6, 24, 15.7, 13.2] μ G. For n_H=10⁶ cm⁻³, we get B_x=[416, 240, 157, 131] μ G, which are typical values observed in dense cores (see Crutcher 2012; Pattle et al. 2023; Whitworth et al. 2025).

2. The initial numerical gas density ranges between n_H= 10⁴ cm⁻³ and n_H=10⁸ cm⁻³. We keep β fixed, which implies that B_x increases as the square root of n_H when varying it, as seen in Eq. (14). This is a reasonable assumption since the magnetic field is expected to behave roughly this way as collapse proceeds (Li et al. 2011; Pattle et al. 2023). Therefore, when increasing n_H, we probe different stages of gravitational collapse of dense cores.

3. Dust grains have sizes ranging from 0.1 to 100 μm. The lowest value lies within the size range of the Mathis-Rumple-Nordsiek distribution (Mathis et al. 1977) observed in the ISM. 1 μm corresponds to a few times the upper value of the MRN distribution, i.e., to dust grains that have undergone some growth beforehand with respect to ISM dust. We also explore larger sizes to gain insight. The size range translates into a Stokes number range: St=[10⁻⁴, 10⁻³, 10⁻², 10⁻¹].

4. The parallel Mach number values we explore lie between ${\mathcal{M}}_{\parallel }=0.12$ and ${\mathcal{M}}_{\parallel }=0.92$ while the perpendicular Mach number covers larger values, ranging from ${\mathcal{M}}_{\perp }=0.17$ to ${\mathcal{M}}_{\perp }=1.91$. This choice of range is justified by turbulence in dense cores being observed as subsonic (Barranco & Goodman 1998; Choudhury et al. 2021). Again, it should be kept in mind that we drive 1D turbulence, which is inherently distinct from 3D turbulence.

We find that the plasma parameter β is the most influential. For a value as low as β=0.1, the other parameters are very little impactful, and the associated plots are supplemented in the appendix, Appendix J, while we provide in the present section the results for β=0.7 for which dust clumping is significant. Every simulation has been run with a resolution of 4096 cells and we note that convergence is not reached regarding dust density. The discussion of this point is delayed to Sect. 5.3.

Figure 4 depicts the maximum dust density (solid lines) and gas density (dotted lines) fluctuations as a function of time for different values of the dust grain size, longitudinal Mach number (perpendicular also varied), gas initial density and plasma parameter. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines. Figure 5 displays the fluctuations of the same quantities, but space-averaged (mass weighted) as follows $$⟨a{⟩}_a={\int }_0^{L}a\cdot a\mathrm{d}x/{\int }_0^{L}a\mathrm{d}x={\int }_0^L{a}^2\mathrm{d}x/{\int }_0^{L}a\mathrm{d}x,$$(18)

where a is either the dust density or gas density. Both are self-weighted, allowing us (for instance if a=ρ_d) to capture regions of dust enrichment. The dust-to-gas ratio average is however weighted by the dust density. Besides, in order to discard high dust-to-gas ratios induced by severe gas depletion (where the dust could be depleted as well but to a lesser extent), we apply a filter and retrieve the dust-to-gas ratio only in cells where ρ_d > ρ_d, 0. Indeed, this corresponds to regions of dust concentration, which is what is of interest when it comes to dust growth. Time is measured in Alfvén crossing time, defined as the time needed for an Alfvén wave to cross the box of length L, i.e., as t_a=L / c_a.

Figure 6 shows the time-averaged (over driven turbulence steady-sate regime) probability distribution functions (PDF) of dust-to-gas ratio fluctuations (first row), gas density fluctuations (second row) and dust density fluctuations (third row).

Fig. 2 Different fields at times t=0.23, t=0.44, and t=0.82 in a simulation with driven turbulence. First panel: dust (solid lines) and gas (dotted lines) density fluctuations. Second panel: x-wise dust (solid) and gas (dotted) Mach number. Third panel: magnetic pressure gradient. Fourth panel: total (x, y, z-wise) dust Alfvénic Mach number. Parameters: nH=106 cm−3, sd=10 μm(St=0.01) and β=1.

Fig. 3 Schematic illustration of the parameters controlling the transverse-to-longitudinal magnetic field ratio $B_{\perp }/B_{\parallel }$ and thus the intensity of dust clumping via the parametric instability in the simulations with turbulence.

Fig. 4 Maximum dust density (solid lines) and gas density (dotted lines) fluctuations as a function of time for different values of the dust grain size sd, longitudinal Mach number (perpendicular one being varied too) ${\mathcal{M}}_{\parallel }$, gas initial density nH, and plasma parameter β. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines. When not varied, β=0.7.

Fig. 5 Same as Fig. 4 but depicting space averages (mass weighted) instead of maximum values. On the first row is displayed the dust-to-gas ratio (dust density weighted), on the second row the dust density (dust density weighted) and on the third row the gas density average (gas density weighted).

Fig. 6 Time-averaged (over driven turbulence steady-sate regime) probability distribution functions (PDFs) of dust-to-gas ratio fluctuations (first row), gas density fluctuations (second row) and dust density fluctuations (third row).

4.2.1 Impact of plasma parameter

We first explore the impact of the plasma parameter β which, as mentioned above, is found to be the most influential one. Looking at Fig. 4, we notice a high sensitivity of the maximum dust density to this parameter. Indeed, at β=0.1, the dust maximum density is enhanced by a factor of 10^0.8=6.31 in the steady-state regime relative to its initial value. For larger values of β, we see a stronger clumping of dust whose density is enhanced by a factor of over 10^1.2=15.85(β=1). In contrast, gas density fluctuations do not vary as we increase β and remain near a value slightly below 10^0.8. This is also the average value of both gas and dust density fluctuations in pure hydrodynamics simulations (pictured as a dashed line) suggesting that neutral gas concentration is controlled solely by hydrodynamical turbulent motions while dust concentration is driven by magnetic effects. This behavior is confirmed in Fig. 5.

As explained in detail in Sect. 4.1, the process responsible for dust concentration is most likely the parametric instability, which develops more efficiently for large transverse-to-longitudinal magnetic field ratios B_⊥ / B_∥. As we increase β, we reduce the longitudinal magnetic field component B_∥=B_x and thus increase the said ratio. When a low value of β=0.1 is considered, magnetic effects do not contribute to dust concentration. Instead, dust and gas remain tightly coupled as ions force neutral gas molecules to move along magnetic field lines (magnetic drag, see Eq. (6)) and only modest density fluctuations are observed (see Fig. J.1). When β increases and magnetic effects prevail, sharper dust clumps are allowed to form and decouple from the gas, resulting in a broader spectrum of dust-to-gas ratios, as seen in Fig. 6.

We show in the Appendix (Fig. G.1) the evolution of the maximum B_⊥ / B_∥ as a function of time in simulations with turbulence for different values of β. We see clearly that stronger B_⊥ / B_∥ are achieved for larger values of β.

4.2.2 Impact of grain size

We investigate here the impact of grain size. As a reminder, when not varied, the plasma parameter takes its fiducial value β=0.7 allowing dust density fluctuations to be enhanced. As readily seen in Fig. 4, larger grains are prone to stronger density fluctuations. For small grains of s_d=0.1 μm in size, a modest level of dust concentration is reached which is equal to that of gas and comparable to what we get from pure hydrodynamics simulations. Both fluid are very well coupled and magnetic forces cannot produce sharp dust clumps. However, for marginally coupled dust grains (s_d ≥ 10 μm, that is St ≥ 10⁻²), an increase of more than a factor of 10 in dust density is easily attained, which continues to increase even after 10 Alfvén crossing times. Even relatively small particles of size s_d ≥ 1 μm experience remarkable density enhancements, close to a factor of 10. Interestingly, the gain in clumping efficiency with size seems to approach saturation as we go from s_d=10 μm to s_d=100 μm, meaning that the largest grains reach a state of full decoupling from the gas, which thus ceases to influence their dynamics. Noteworthy, dust-to-gas ratio fluctuations also increase with size, reaching values above 100 and as high as 1000 for s_d ≥ 100 μm (see Fig. 5).

This dependency of dust concentration on grain size has interesting consequences in terms of dust coagulation. As discussed in Sect. 5.2.2, this could lead to a runaway dust growth.

4.2.3 Impact of Mach number

We now turn our attention to the influence of the Mach number. Note that computational cost of simulations sharply increases with this parameter, therefore we chose not not to integrate the run for ${\mathcal{M}}_{\parallel }=0.92$ beyond 4 Alfvén crossing times.

First of all, we notice that unlike the previous parameters, varying the Mach number does affect gas density fluctuations. Indeed, increasing the turbulent Mach number leads to stronger fluctuations in the gas density, reaching higher maximum values and strikingly lower minimum values, i.e., a much more inhomogeneous distribution of matter. Figure 6 clearly shows that gas density fluctuations as low as 10⁻⁶ are not uncommon in regions of gas depletion. The maximum gas density increases by a factor of 10^0.4=2.5 for ${\mathcal{M}}_{\parallel }=0.12$ to up to 10 for ${\mathcal{M}}_{\parallel }=0.92$ (third column, third row). Again, those values match those from hydrodynamics simulations, showing that gas clumping is due to non magnetic turbulent effects only.

As for the dust, we observe a level of clumping similar to that of gas in the lowest Mach number run (${\mathcal{M}}_{\parallel }=0.12$), implying that both fluids are well coupled and concentrate jointly. However, as the Mach number increases, a substantial departure arises and dust is allowed to exhibit stronger density enhancements, by a factor of about 10 for ${\mathcal{M}}_{||}=0.56$ and over 10^1.4= 25.1 for ${\mathcal{M}}_{\parallel }=0.92$ (third column, second row). Here again, this behavior should be attributed to magnetic effects, keeping in mind that the transverse Mach number ${\mathcal{M}}_{\perp }$ is varied along with its longitudinal counterpart. As explained in Sect. 4.1 and discussed in Sect. 5.2.1, the transverse-to-longitudinal magnetic ratio B_⊥ / B_∥ (and thus the parametric instability) is controlled by ${\mathcal{M}}_{\perp }$ (see Table 2).

Finally, while it is striking that the maximum dust-to-gas ratio significantly increases with Mach number, the minimum value does not drop as low for ${\mathcal{M}}_{\parallel }=0.92$ as for ${\mathcal{M}}_{\parallel }=0.56$. This means that in supersonic regime, the gas is less inclined to concentrate in regions devoid of dust.

Fig. 7 Average (dust density weighted) dust density (solid lines) and average (gas density weighted) gas density (dotted lines) fluctuations as a function of time for different values of the dust grain size sd. Upper left panel: β=0.1. Lower left panel: β=1. Second column: same but with ten times less ions. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines.

4.2.4 Impact of background gas density

Finally, we explore the impact of the initial background gas density n_H. Although local gas density fluctuations remain unaffected, dust clumping appears to be reduced in denser environments. Indeed, while maximum fluctuations by a factor of 10 are easily reached at $n_{\mathrm{H}}={10}^4{\mathrm{c}\mathrm{m}}^{-3}$ and n_H=10⁶ cm⁻³, they remain below this threshold at n_H=10⁸ cm⁻³ (see Fig. 4, fourth panel). This is caused by the competition of two different effects, one inhibiting dust concentration and the other promoting it. First, a higher background gas density implies a lower Stokes number St (a lower stopping time t_s, d, see Eq. (8)) for a given dust grain size, that is, a better coupling between dust and gas, preventing efficient dust clumping. Second, we expect on average that the ion numerical density relative to dust density is lower in denser environments since ${\rho }_i\propto \sqrt{n_{\mathrm{H}}}$ (see Eq. (12)) and ρ_d, 0=ε ρ₀ ∝ n_H. As discussed in Sect. 5.1, the impact of magnetic drag due to ions should be reduced compared to that of hydrodynamical drag, allowing for an easier dust clumping. However, this effect is negligible with respect to the influence of a lower St, and the overall tendency is a reduction in dust density fluctuations.

Note that, as per Eq. (14), a higher density with a fixed β yields a stronger longitudinal magnetic field B_∥=B_x, and one could think that for a given level of turbulence (${\mathcal{M}}_{\perp }$), higher transverse-to-longitudinal magnetic ratio B_⊥ / B_∥ could be reached. However, it is not the case, because a fixed β implies that the initial gas density n_H and the background B_x are varied jointly in order to keep the same Alfvén velocity c_a. As a consequence, similar dust Alfvénic Mach numbers are achieved, giving rise to very similar B_⊥ / B_∥ (see Table 2). Thus, it is only the decrease in St that leads to a less efficient clumping of the dust at higher n_H. In Fig. 6, we see that the dust-to-gas ratio PDF is slightly narrower at n_H=10⁸ cm⁻³ due to tighter coupling between gas and dust.

5 Discussion

5.1 non-ideal effects and chemical network

As seen in Fig. F.1 in Appendix F (see also Fig. 1), dust clumping is curtailed to a certain extent when considering non-ideal MHD effects. We thus see that there is a need to consider microphysical processes along with an accurate description of MHD equations to predict relevant physical quantities such as density fluctuations. First, in simulations, ions cannot be considered as main charge carriers in dense cores, and therefore gas should be considered neutral (contrary to more diffuse environments such as molecular clouds, see Lee et al. 2017; Moseley et al. 2023). Second, the present work demonstrates that assuming a perfect coupling between dust particles and magnetic field lines (Sect. 2.2) is inaccurate. non-ideal MHD effects do make a difference.

We found the discrepancies between dusty ideal MHD and dusty non-ideal MHD models to be primarily attributed to the additional magnetic drag term induced by the presence of ions in the dust momentum equation (Eq. (6)) that we rewrite $\frac{e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}B}{{\mathrm{\Gamma }}_{\mathrm{i}}c}(\mathit{v}-{\mathit{v}}_{\mathrm{d}})=\frac{{\rho }_i}{t_{\mathrm{s},\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{d}})$. It is instructive to compare this magnetic drag term with the hydrodynamical one, i.e., to compare $\frac{{\rho }_i}{t_{\mathrm{s},\mathrm{i}}}$ with $\frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}$. We first compute stopping times. Using the fiducial background gas density ρ=3.9 × 10⁻¹⁸ g cm⁻³ (corresponding to n_H=10⁶ cm⁻³) and an ion (HCO⁺) gas collision rate γ_{i, g} = 3.5 × 10¹³ cm⁻³ g⁻¹ s⁻¹ given in Pinto & Galli (2008), we find $t_{\mathrm{s},\mathrm{i}}={\left(\rho {\gamma }_{\mathrm{i},\mathrm{g}}\right)}^{-1}\simeq {10}^4\mathrm{s}$. From Eq. (8) for s_d=100 μm (St = 0.1) we compute the dust stopping time t_{s, d} ≃ 2 × 10¹¹ s. In terms of orders of magnitude, it yields $$\begin{array}{rl}& \frac{{\rho }_{\mathrm{d}}}{t_{\mathrm{s},\mathrm{d}}}\simeq \frac{{10}^{-20}}{{10}^{11}}={10}^{-31}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}{\mathrm{s}}^{-1}\\ & \frac{{\rho }_i}{t_{\mathrm{s},\mathrm{i}}}\simeq \frac{{10}^{-25}}{{10}^4}={10}^{-29}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}{\mathrm{s}}^{-1}\end{array}$$(19)

where ρ_i was computed from Eq. (12) (see also Fig. I.1). The higher mass density of dust grains does not compensate for the enormous gap in stopping times. Magnetic drag dominates by a factor of ≃ 100. Since t_{s, d} ∝ s_d, the two drag terms become comparable if we decrease the dust size by two orders of magnitude, that is, for s_d ≤ 1 μm.

Based on the previous calculations, we see that magnetic drag prevails under typical conditions. As long as it does, dust grains remain very well coupled to the gas and cannot concentrate independently. Unlike the dust stopping time, the ion stopping time is not a function of grain size. Consequently, in this regime, larger grains do not decouple from the gas as we usually expect from hydrodynamical drag. This explains the modest dust density fluctuations observed in Fig. J.1 (where β=0.1) and in particular why we do not see a dependency of the dust density fluctuations on grain size. This behavior is supported by the impact of Mach number, where we see that regardless of its value, both dust (solid lines) and gas (dotted lines) clump by the same amount (they are tightly coupled by magnetic drag) and reach the same values as those obtained in pure hydrodynamical simulations (dashed lines).

There is a competition between the parametric-like instability identified in Sect. 4.1 and magnetic drag. The former favors dust concentration, while the latter inhibits it. As already mentioned, we need to increase the transverse-to-longitudinal magnetic field ratio B_⊥ / B_∥ to enhance instability and hence dust clumping. To do so, we can either increase the plasma parame ter β (i.e., decrease B_x) or increase the transverse Mach number ${\mathcal{M}}_{\perp }$ (see Fig. 3). Indeed, the induction equation tells us that the background longitudinal magnetic field B_x is distorted, producing transverse components B_y and B_z if the main charge carrier (here the dust) is allowed to move in transverse directions. On the other hand, to reduce magnetic drag, one could play with the ion numerical density n₁. Decreasing it would allow dust grains to decouple more easily from the gas and clump. This point leads us to discuss our choice of chemical network.

Although our approach (see Sect. 2.3) offers simplicity, it is not fully realistic. First of all, there is no dependency of the total numerical density of charge n_d Z_d of the dust fluid on grain size. Although it is expected to increase with dust size (Draine & Sutin 1987), the exact charge scaling with size is not straightforward when considering various reactions such as dust-electron capture, dust-ion capture, ion-electron recombination on grain surface or even grain-grain reactions (e.g., charge transfer, see Waitukaitis et al. 2014). Second, a more sophisticated chemical network including a dust size distribution would lead to a lower ion abundance owing to very small grains offering a large total cross section for ion-electron recombination on their surface. Upon comparison with the chemical network of Marchand et al. (2021) (see Fig. I.1) where a dust size distribution (MRN like) is considered, ions should be a few times less than predicted by our chemical network at about n_H=10⁶ cm⁻³, leading to a weaker magnetic drag. The departure from one model to the other is reduced at lower density but is more severe at higher density. Note that abundances also depend on the cosmic ray ionization rate ζ.

5.2 Dust growth in protostellar envelopes

Recent observational studies have revealed low dust emissivity index (β < 1)¹ in protostellar envelopes (at radii from 100 to 500 AU), suggesting the presence of dust grains ranging between 100 μm and 1 mm in size (Galametz et al. 2019; Cacciapuoti et al. 2023, 2025), supported by polarization measurements (Valdivia et al. 2019; Sadavoy et al. 2019). Indeed, the emissivity index β is anticorrelated with the maximum grain size of the distribution where values β ≥ 1.5 are consistent with dust sizes s_d ≤ 10 μm, while β ≤ 1 indicate the presence of s_d ∼ 100 μm dust grains. Although other considerations could explain this particularly low value (effects of temperature and optical depth, contamination by synchrotron radiation, free-free emission and spinning dust emission, and in particular grain structure and composition effects, see Carpine et al. 2025), the actual presence of significantly grown dust grains is more than plausible (Demyk et al. 2017; Ysard et al. 2019).

Two (not mutually exclusive) scenarios have been proposed to explain these observations, namely the possibility for the grains to grow in situ by means of coagulation, and the ejection of already grown grains via winds and outflows from the inner regions of the protoplanetary disk into the inner and intermediate regions of the envelope (see Wong et al. 2016; Tsukamoto et al. 2021, for the latter). This mechanism depends on some criteria including the mass of the central stellar embryo, the gas velocity, the opening angle of the outflow and the mass loss rate associated. The mass of the central nascent star needs to be sufficiently low, and the outflow mass ejection rate high enough to ensure an efficient transportation of the grains.

When it comes to in situ grain growth, current coagulation models face difficulties to form such large grains in a reasonable amount of time under standard assumptions (see Silsbee et al. 2022). However, we stress here that this scenario cannot be ruled out yet as it seems possible to form millimeter grains in situ provided that the dust is concentrated locally. In this section, we discuss the clumping measurements made in our simulations and the likelihood of this scenario.

5.2.1 Turbulent clumping is a good candidate

We saw that the two main ingredients when it comes to dust clumping are the transverse-to-longitudinal magnetic field ratio B_⊥ / B_∥ and the ion numerical density n_i. One way to increase the former, is by increasing the transverse Mach number ${\mathcal{M}}_{\perp }$. The corresponding effect is readily seen in Fig. 4. Although dense core environments are expected to be sub-/transonic (Barranco & Goodman 1998; Choudhury et al. 2021), certain mechanisms could generate supersonic turbulence (Hennebelle 2021; Higashi et al. 2021, explore the generation and amplification of turbulence in collapsing environments (gravo-turbulence)). Studies show that a significant fraction of the velocity dispersion appears as a systematic drift that cannot be interpreted as turbulence. However, systematic drift would lead to instabilities (resonant drag instabilities, Hendrix & Keppens 2014; Squire & Hopkins 2018) favoring dust concentration. The rest is in the form of random velocity fluctuations which can be seen as supersonic turbulence on smaller scales. In massive dense cores, turbulence is expected to reach supersonic regime, although this claim has recently been questioned (see Wang et al. 2024). Figure 4 shows that for ${\mathcal{M}}_{\parallel }=0.56({\mathcal{M}}_{\perp }=0.99)$ and β=0.7, an increase by a factor of 0.8-1.0 in the average dust density is expected, which would allow for rapid grain growth. Going supersonic leads to even stronger dust clumping.

Another way to increase the transverse-to-longitudinal magnetic field ratio for a given Mach number, is by increasing the plasma parameter β. Although too weak a magnetic field would suppress the clumping observed in our simulations $(\beta ⟶\mathrm{\infty }$ implies $\mathit{B}⟶\mathbf{0}$), reducing the magnetization of the background magnetic field to a certain extent does favor dust clumping via parametric instability (Sect. 4.1). Indeed, for a given turbulent driving (i.e., a given transversal Mach number ${\mathcal{M}}_{\perp }$), the level of dust concentration is very sensitive to β as it controls the transverse-to-longitudinal magnetic ratio B_⊥ / B_∥(see Table 2 and Fig. G.1). Figure J.1 indicates that the average dust density increases only slightly (by a factor of $\sqrt{10}$) relative to initial value for β=0.1 compared with a magnetization seven times lower (β=0.7), which yields a dust density enhancement of more than a factor of ten (Fig. 4). Although β=0.1 is a typical value for protostellar cores, a distribution of values is measured observationally (Crutcher 2012; Pattle et al. 2023) and found in simulations (Hennebelle et al. 2011; Tu et al. 2022). Cores threaded by weaker background magnetic fields would be more sensitive to this mechanism and prone to significant dust clumping. Nevertheless, one should keep in mind the timescales of interest. Our results suggest that significant clumping is achieved after about 1 to 2 Alfvén crossing times. When β=1, by definition, the Alfvén crossing time corresponds to the sonic crossing time, and since we defined our box length L as the Jeans length, such a time is equal to the free-fall time. Fortunately, protostellar envelopes should persist over a few free-fall times owing to other support mechanisms, thus offering enough time for dust grains to grow significantly. However, this conclusion should not hold for larger values of β.

We explained in Sect. 5.1 that fewer ions would allow dust grains to decouple from the gas and thus accumulate more efficiently. Based on the associated discussion, we explore here the impact of considering a 10 times lower ion abundance: n_i → n_i / 10. This is readily seen in Fig. 7 which depicts a striking soaring in the dust density when combining a lower magnetization (β=1) and ten times less ions. Micron-sized (s_d=1 μm) grains display a increase by a factor of 10 in density, while a factor of about 30 is reached for grains ten times bigger (S_d=10 μm). Dust density fluctuations climb as high as 10^2.5=300 for grains of size S_d=100 μm. Considering ion numerical densities ten times lower than those predicted by our chemical network is not unreasonable since it matches the values obtained with a much more sophisticated chemical network including a distribution of grain size (see Appendix I), where very small grains would efficiently recombine ions and electrons on their surface. Despite a swift coagulation owing to ambipolar drift, the small grain population is expected to persist at low densities in dense cores due to electrostatic repulsion (investigated in protoplanetary disks in Akimkin et al. 2023) and grain-grain erosion (Vallucci-Goy et al. 2024). Note however, that if sufficiently well coupled to the magnetic field, very small grains could behave as ions to some extent and maintain a certain level of magnetic drag between gas and larger grains, making efficient clumping more difficult to achieve. On the other hand, we should keep in mind that very small dust grains can carry only a single electric charge (just as the ions we consider here) but are substantially more massive than ions. This means that their (small dust grains) Hall factor would be higher and thus the subsequent induced magnetic drag lower (according to Eq. (6)) than that induced by ions in the present work. Simulations with a distribution of dust sizes and a more realistic chemical network are needed to address this point, but require substantial developments that are beyond the scope of this work.

In addition, a different ionization rate induced by cosmic rays ζ would change the abundance of ions as well. In Eq. (12), we assumed ζ=1 × 10⁻¹⁷ s⁻¹. However, this parameter is not tightly constrained and values ranging from 5 × 10⁻¹⁸ to 5 × 10⁻¹⁶ s⁻¹ are measured in dense cores (Padovani et al. 2022). A lower ionization rate would imply fewer ions and thus possibly a more efficient dust concentration, although less charges would be available for dust grains to carry.

In summary, specific conditions are needed for strong clumping that are not systemically met. However, the diversity of physical conditions (level of turbulence and magnetization, cosmic-ray ionization rate) encountered in dense core environments and protostellar envelopes make this mechanism a promising pathway to the formation of the large dust grains observed, at least in certain cases.

5.2.2 Toward a coagulation instability

In this section, we mention the potential development of a “coagulation instability” as defined in Tominaga et al. (2021, 2022) except that the said instability was studied in protoplanetary disks and resulted from a combination of dust growth and radial drift. In the present work, the instability would be driven by the interplay between magnetic clumping and dust growth.

On one hand, an increase in dust density as a result of magnetic clumping leads to a lower coagulation timescale. On the other hand, as dust grains grow in size (and in Stokes number), they decouple more efficiently from the gas. As showed in Sect. 4.2.2 and discussed in Sect. 5.2.1, larger grains concentrate more strongly and reach higher density fluctuations (in the regime where the parametric instability can develop). In turn, a higher density further reduces the coagulation timescale, allowing dust grains to grow faster and clump even more efficiently, etc. The positive feedback of dust growth on dust concentration is what triggers the instability. A very simple illustration of this concept is displayed in Fig. 8. This leads to a runaway process that could allow for a very rapid growth and the formation of large grains over reasonable timescales. However, a point that remains to be elucidated is the time that would be needed to trigger this instability.

Again, if the parametric-like instability does not overcome magnetic drag due to transverse-to-longitudinal magnetic field ratio being too low (case β=0.1 in Fig. 7), larger grains do not exhibit a higher level of clumping and the coagulation instability would thus not operate. Interestingly, local magnetic flux expulsion observed in 3D simulations of collapsing dense cores would induce inhomogeneous dust clumping and growth. Indeed, in regions of lower magnetization, the coagulation instability would be at play and dust growth would proceed very rapidly, resulting in fine in a spatial segregation/sorting of dust by size.

Fig. 8 Schematic illustration of a possible coagulation instability in turbulent magnetized dense cores. Dust concentration leads to a more rapid dust growth, which in turn results in larger grains prone to stronger concentration.

5.3 Absence of convergence in the dust density

We show in Fig. E.1a convergence test by plotting the dust and gas average density fluctuations with increasing resolution. While resolution has no effect on the gas density, convergence is not reached in the dust even for a number of cells as high as NX =8192. This results from the absence of any pressure in the dust momentum equation (Eq. (6)), which enables the formation of smaller and sharper clumps. This further supports the capacity of magnetic effects to induce substantial dust clumping, which could be significantly higher than that presented in this work. Note that similar trends are observed in simulations of instabilities involving significant levels of dust concentration when the dust fluid is considered as pressureless (for instance with the streaming instability, see Johansen & Youdin 2007; Benítez-Llambay et al. 2019).

However, on sufficiently small scale, we would expect dust velocity dispersion to act as a pressure and halt further increase in dust density, even in absence of grain-grain collisions. See the work of Garaud et al. (2004); Hersant (2009) and Lynch & Laibe (2024) for more information on the relation between velocity dispersion and pressure for a dust fluid. External forces such as gravity or magnetic forces, turbulent converging flows or even residual friction with gas (if some of it is being dragged along) could provide such velocity dispersion. The scale on which a pressure would emerge remains unclear. Marginally coupled dust grains could even concentrate to the point of gravitational instability if not opposed by pressure before, leading potentially to early planetesimal formation.

In that regard, we test the robustness of our results by showing in Fig. H.1 the impact of a dust pressure added a posteriori in the corresponding momentum equation. We parametrize the dust sound speed as a fraction of that of the gas. We see no effects on dust density fluctuations for values c_{s, d} < 0.5 c_s for the chosen spatial resolution (NX=4096). For the highest value c_{s, d}=c_s, the magnetic clumping is suppressed to a certain extent. However, it is reasonable to think that such a high incompressibility in the dust is unrealistic within our setup. To show that, we use the subgrid model of Ormel & Cuzzi (2007) which provides velocity dispersions for dust grains coupled to turbulent vortices. For a dust grain of Stokes number St, it gives $$\mathrm{\Delta }{v}_{\mathrm{d}}\simeq {\sigma }_{\mathrm{v},\mathrm{d}}=v_{\mathrm{g}}\sqrt{1.97\mathrm{S}\mathrm{t}}=\sqrt{\alpha }{c}_{\mathrm{s}}\sqrt{1.97\mathrm{S}\mathrm{t}},$$(20)

where α ≃ 1 is a suitable measure of the (subsonic) turbulence intensity in dense cores. Velocity dispersions equal to σ_{v, d} ≃ c_{s, d}=[0.5 c_s, c_s] are associated with St=[0.13, 0.51], which are values above those explored in this work (the fiducial value being St=0.01, see Table 2). Those simple estimates suggest that the dust concentration measured in the simulations should hold even in presence of a dust pressure in the equations. In addition, one has to keep in mind that this subgrid model ignores magnetic effects and dust backreaction. Accounting for mass loading of dust particles, we expect gas turbulence to be attenuated and thus dispersion velocities to be further reduced (Tominaga & Tanaka 2025; Carrera et al. 2025).

5.4 Caveats

A major caveat is the one dimensional nature of our simulations because 1D and 3D turbulence are known to behave differently. On one hand, geometrical effects in 3D simulations could reduce dust concentration. On the other hand, turbulent solenoidal modes would produce vortices known to be capable of trapping solid particles. Indeed particles are centrifuged out of turbulent eddies and concentrate in regions of low vorticity. This happens when the local Stokes number reaches unity (Squires & Eaton 1991).

As already mentioned throughout the text, there is a need for a more realistic chemical network to properly compute grain charges and ion abundances. To do so, a distribution of dust sizes will have to be implemented along with inertialess electrons in the equations, which we will do in a future work. In addition, a numerical treatment of dust growth (coagulation and fragmentation) will have to be included to follow self-consistently the formation of large grains.

6 Conclusion

In this work, we have investigated dust dynamics and clumping in the dense ISM (n_H=10⁴−10⁸ cm⁻³) performing multifluid 1D simulations of magnetized dust and neutral gas within two different MHD regimes:

1. Dusty ideal MHD: dust grains are assumed to be perfectly coupled to magnetic field lines.

2. Dusty non-ideal MHD: ions (whose inertia is neglected) are added in the equations along with ion-neutral friction, introducing an Ohmic dissipation and a dispersive Hall effect in the induction equation as well as an additional “magnetic drag” in the momentum equations that recouples dust and gas.

We first propagated nonlinear circularly polarized Alfvén waves and identified a mechanism similar to the parametric instability to be an efficient pathway to dust concentration. Strong dust clumping is obtained in the dusty ideal regime on scales where hydrodynamical dust-gas friction is negligible. Lower but remarkable levels of dust concentration are maintained in the dusty non-ideal regime despite non-ideal MHD effects (magnetic drag) imperilling the survival of dust clumps. Short wavelengths and sufficiently high wave amplitudes (i.e., high transverse-to-longitudinal magnetic ratio B_⊥ / B_∥) are needed to achieve significant clumping, conditions that are met and sustained when considering a turbulent cascade driven on large scale.

Turning then to simulations of 1D driven turbulence, we first described how dust clumps form. They develop in shocks, regions of high gas density where the dust can easily reach super-Alfvénic regime and distort the background longitudinal magnetic field lines to produce transverse components on small scales. We then investigated the influence of four different parameters on dust and gas density fluctuations, namely dust size S_d, longitudinal ${\mathcal{M}}_{\parallel }$ and transverse ${\mathcal{M}}_{\perp }$ Mach numbers, background gas density n_H and plasma parameter β. Significant dust concentration, larger than that obtained with pure hydrodynamical simulations is achieved because of small scale compressive magnetic effects. There is a competition between the parametric-like instability which promotes dust concentration, and magnetic drag which hinders it. The former is sensitive to the transverse-to-longitudinal magnetic field ratio B_⊥ / B_∥ while the latter can be weakened by reducing the ion abundance n_i. As a consequence, stronger dust clumping is expected in regions of intense turbulence (${\mathcal{M}}_{\perp }$ controls B_⊥) and lower (but not too low) magnetization (β controls B_∥). Indeed, for a high magnetization, field lines are more difficult to bend and lower B_⊥ / B_∥ are produced for a given level of turbulence.

We found that a substantial fraction of the dust mass is easily and rapidly concentrated to a density more than 10 times its initial value, provided that β ≥ 0.7 with subsonic turbulence and dust size S_d ≥ 1 μm, making this novel mechanism a promising candidate for the formation of large grains in situ in protostellar envelopes.

Acknowledgements

We thank the referee who helped us greatly in improving the clarity and relevancy of this paper. We thank the consortium and ERC (European Research Council) synergy grant ECOGAL (grant 855130)) for their financial support, and members for their contribution through insightful ideas and remarks.

References

Appendix A Derivation of generalized Ohm's law

We give here the few steps leading to the expression of the electric field in the dusty non-ideal MHD regime. Starting from the ion momentum equation (2), one gets $$\mathit{E}=-\frac{{\mathit{v}}_{\mathit{i}}}{c}\times \mathit{B}-\frac{{\rho }_{\mathrm{i}}\rho {\gamma }_{\mathrm{i},\mathrm{g}}}{e{n}_{\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{i}})=-\frac{{\mathit{v}}_{\mathit{i}}}{c}\times \mathit{B}-\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{i}}).$$(A.1)

The electroneutrality condition for ions and dust is $$n_{\mathrm{d}}{Z}_{\mathrm{d}}+n_{\mathrm{i}}=0.$$(A.2)

From the total electric current ${\mathit{J}}_{\text{tot}}=e{n}_{\mathrm{i}}{\mathit{v}}_{\mathit{i}}+e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}{\mathit{v}}_{\mathbf{d}}$, using Maxwell-Ampère's law (Eq. 5) along with Eq. (A.2), the ion velocity ensues $${\mathit{v}}_{\mathit{i}}=\frac{c}{e4\pi n_{\mathrm{i}}}\mathrm{\nabla }\times \mathit{B}+{\mathit{v}}_{\mathrm{d}},$$(A.3)

which we inject into Eq. (A.1) $$\mathit{E}=-\frac{{\mathit{v}}_{\mathbf{d}}}{c}\times \mathit{B}-\frac{1}{e4\pi n_{\mathrm{i}}}(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}+\frac{\parallel \mathit{B}\parallel }{{\mathrm{\Gamma }}_{\mathrm{i}}e4\pi n_{\mathrm{i}}}(\mathrm{\nabla }\times \mathit{B})+\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}({\mathit{v}}_{\mathbf{d}}-\mathit{v}).$$(A.4)

In addition, starting from $F_{\text{Lor,}\mathrm{d}}=n_{\mathrm{d}}{Z}_{\mathrm{d}}e(\mathit{E}+\frac{{\mathit{v}}_{\mathbf{d}}}{c}\times \mathit{B})$ and using Eq. (A.1), we can rewrite the dust Lorentz force the following way $$\begin{array}{rl}{\mathit{F}}_{\text{Lor,}\mathrm{d}}& =-e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\frac{{\mathit{v}}_{\mathrm{i}}}{c}\times \mathit{B}-e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{i}})+\frac{{\mathit{J}}_{\mathrm{d}}}{c}\times \mathit{B}\\ & =\frac{{\mathit{J}}_{\text{tot}}}{c}\times \mathit{B}-\frac{{\mathit{J}}_{\mathrm{d}}}{c}\times \mathit{B}+\frac{{\mathit{J}}_{\mathrm{d}}}{c}\times \mathit{B}-e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-\frac{{\mathit{J}}_{\text{tot}}}{e{n}_{\mathrm{i}}}-{\mathit{v}}_{\mathrm{d}})\\ & =\frac{{\mathit{J}}_{\text{tot}}}{c}\times \mathit{B}-\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}{\mathit{J}}_{\text{tot}}-e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{d}}).\end{array}$$(A.5)

Similarly for the drag felt by the gas due to ion-neutral friction, starting from Eq. (2) and using the ion velocity (Eq. A.3) $$\begin{array}{rl}{\mathit{F}}_{\mathrm{i}\to \mathbf{N}}& =n_{\mathrm{i}}e\mathit{E}+\frac{{\mathit{J}}_{\mathrm{i}}}{c}\times \mathit{B}\\ & =-n_{\mathrm{i}}e{\mathit{v}}_{\mathrm{i}}\times \mathit{B}-\frac{n_{\mathrm{i}}e\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-{\mathit{v}}_{\mathrm{i}})+\frac{{\mathit{J}}_{\mathrm{i}}}{c}\times \mathit{B}\\ & =\frac{{\mathit{J}}_{\mathrm{i}}}{c}\times \mathit{B}-\frac{{\mathit{J}}_{\mathrm{i}}}{c}\times \mathit{B}-\frac{n_{\mathrm{i}}e\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}(\mathit{v}-\frac{{\mathit{J}}_{\mathrm{t}\mathrm{o}\mathrm{t}}}{e{n}_{\mathrm{i}}}-{\mathit{v}}_{\mathrm{d}})\\ & =\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}{\mathit{J}}_{\mathrm{t}\mathrm{o}\mathrm{t}}+\frac{n_{\mathrm{i}}e\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v})\\ & =\frac{\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}{\mathit{J}}_{\mathrm{t}\mathrm{o}\mathrm{t}}-\frac{e{n}_{\mathrm{d}}{Z}_{\mathrm{d}}\parallel \mathit{B}\parallel }{c{\mathrm{\Gamma }}_{\mathrm{i}}}({\mathit{v}}_{\mathrm{d}}-\mathit{v}).\end{array}$$(A.6)

Appendix B Dispersion relation of the dusty ideal Alfvén wave

We study the propagation of a linear Alfvén wave within the dusty ideal MHD model (Sect. 2.2) and compare in Fig. B.1 the analytical dispersion relation with the one measured in our simulations. This analysis serves as a test for the code and as a source of physical insight. The circularly polarized Alfvén wave is initialized as $$B_y={\delta }_B{B}_x\sin \left(\frac{2\pi }{\lambda }x\right)B_z={\delta }_B{B}_x\cos \left(\frac{2\pi }{\lambda }x\right),$$(B.1)

with the wavelength λ being a fraction of the box length L. The wave amplitude is taken to be δ_B B_x=10⁻⁴ B_x so that we work in the linear regime. The plasma parameter is taken to be β=0.1 and we set the initial Stokes number to be St=0.001. Both gas and dust are initially moving in the transverse directions with velocity components imposed by the Alfvén mode targeted. Note that we work with code units here: we consider a uniform gas initial density ρ₀=1 and dust initial density ρ_d, 0=0.01 ρ₀ (i.e., a dust-to-gas ratio ε=0.01), a uniform x-wise background magnetic field inferred from the plasma parameter β as $B_x=c_{\mathrm{s}}\sqrt{4\pi {\rho }_0/\beta }$, and a constant soundspeed c_s=0.316.

Fig. B.1 Dispersion relation for an Alfvén wave in a dusty magnetized medium within the dusty ideal MHD model (Sect. 2.2). Comparison between the analytical solution and numerical measurement from our simulation. The real part of the frequency appears as solid line and the imaginary part as dashed line. St =0.001, β=0.1 and δB=10−4.

Appendix B.1 Analytical solution

To derive the analytical dispersion relation, we start from the gas and dust momentum equations (Eq. 11) as well as the induction equation (Eq. 10) and carry out a linear perturbation analysis. Since Alfvén modes are incompressible, we discard the mass conservation equations. In addition, we consider no background quantities except for B_x and we work with transverse perturbations only. Combining the y and z components, we use the following change of variables and work with the polarization modes “+” and “-” $$V_{\pm }=V_y\pm i{V}_z{B}_{\pm }=B_y\pm i{B}_z$$(B.2)

We consider the following small perturbations $$\begin{array}{rlr}\delta v_{\pm }& =& V^{\pm }\exp (i(\omega t-kx)),\\ \delta v_{\mathrm{d},\pm }& =& V_{\mathrm{d}}^{\pm }\exp (i(\omega t-kx)),\\ \delta b_{\pm }& =& B^{\pm }\exp (i(\omega t-kx)),\end{array}$$(B.3)

where w and k are the perturbation frequency and wavenumber. The momentum and induction equations then become $$\begin{array}{rlr}i\omega V^{\pm }& =& \frac{{\rho }_{\mathrm{d}}}{\rho t_{\mathrm{s}}}(V_{\mathrm{d}}^{\pm }-V^{\pm })\\ i\omega V_{\mathrm{d}}^{\pm }& =& -\frac{(V_{\mathrm{d}}^{\pm }-V^{\pm })}{t_{\mathrm{s}}}-\frac{ik{B}_x{B}^{\pm }}{{\rho }_{\mathrm{d}}}\\ i\omega B^{\pm }& =& -B_{x}ik{V}_{\mathrm{d}}^{\pm }\end{array}$$(B.4)

It is straightforward to show that the system reduces to $$\left(\begin{array}{cc}i\omega +{\rho }_{\mathrm{d}}/\rho t_{\mathrm{s}}& {\rho }_{\mathrm{d}}\omega /\rho t_{\mathrm{s}}{B}_{x}k\\ -1/t_{\mathrm{s}}& ik{B}_x/{\rho }_{\mathrm{d}}-(i\omega +1/t_{\mathrm{s}})\omega /B_{x}k\end{array}\right)(\genfrac{}{}{0ex}{}{V^{\pm }}{B^{\pm }})=0.$$(B.5)

Equating the matrix determinant to zero, the dispersion relation ensues $$\frac{{\omega }^3}{k}-\frac{i(\rho +{\rho }_{\mathrm{d}})}{k\rho t_{\mathrm{s}}}{\omega }^2-\frac{k{B}_x^2}{{\rho }_{\mathrm{d}}}\omega +i\frac{k{B}_x^2}{\rho t_{\mathrm{s}}}=0,$$(B.6)

which is a third degree polynomial. One mode is evanescent and thus does not propagate. The other two are degenerate Alfvén modes which propagate in opposite directions, which is expected since there is no Hall effect breaking the symmetry.

Appendix B.2 Comparison with simulation

The evolution of the frequency (real part in solid line and imaginary part in dashed line) with respect to the wavelength for the progressive Alfvén mode is plotted in Fig. B.1 with St=0.001, β=0.1 and δ_B=10⁻⁴. We compare the analytical solution with the frequencies and damping rates measured in our simulation. To ensure that the right mode is excited, we inject the associated solution of Eq. (B.6) in the set of perturbed equations Eq. (B.4) and use them to infer the initial conditions for the gas and dust velocity perturbations. We notice a striking match between the analytical solution and the numerical measurements. At large wavelengths, the dynamical timescale of the Alfvén wave (t_a= L / c_a) is much larger than the grain stopping time (i.e., t_a ≫ t_s), and we thus recover the strong coupling regime where the Alfvén velocity is determined by the gas density, $c_{\mathrm{a}}=\parallel \mathit{B}\parallel /\sqrt{4\pi \rho }=1$ in code units. As we decrease the wavelength, we enter a very dissipative regime where the imaginary part becomes much higher than the real part, preventing any propagation. This is due to a resonance effect (t_a ≃ t_s) where most of the magnetic energy is lost to friction with the gas. We note that the frequency of the signal could not be measured from the corresponding simulation in this regime, and we did not display the near-zero real part from the analytical solution for the sake of clarity. Finally, at low wavelength, we see that Alfvén waves do propagate with weak dissipation, which is a feature already identified analytically in Hennebelle & Lebreuilly (2023) when accounting for dust inertia. Since t_a ≪ t_s, gas-grain collisions are scarce and Alfvén waves are carried by the grains. The Alfvén waves thus propagate at a velocity ten times higher $c_{\mathrm{a},\mathrm{d}}=\parallel \mathit{B}\parallel /\sqrt{4\pi {\rho }_{\mathrm{d}}}=10$ as we consider a dust-to-gas ratio ε=0.01.

Appendix C Growth rate of the parametric instability

We aim to support the statement that the instability observed and described in Sect. 4.1 is very similar to the so-called parametric instability (Del Zanna et al. 2001) but in a two-fluid configuration. To do so, we provide the analytical growth rate of the fastest growing mode derived in the standard one-fluid configuration, and compare it to the numerical results (see Fig. C.1) in the early growth phase.

For the standard parametric instability, the dispersion relation for the daughter waves (with a generic β) is (Derby 1978; Goldstein 1978; Del Zanna et al. 2001) $$\begin{array}{rl}({\tilde{\omega }}^2-\beta {\tilde{k}}^2)(\tilde{\omega }-\tilde{k}{)}^2& (\tilde{\omega }+\tilde{k}+2)(\tilde{\omega }+\tilde{k}-2)\\ & -{\delta }_B^2{\tilde{k}}^2(\tilde{\omega }-\tilde{k})({\tilde{\omega }}^3+\tilde{k}{\tilde{\omega }}^2-3\tilde{\omega }+\tilde{k})=0,\end{array}$$(C.1)

where the dimensionless frequency and wavenumber are respectively $\tilde{\omega }=\omega /{\omega }_0$ and $\tilde{k}=k/k_0$. We choose the mother Alfvén wave frequency to be ω₀=c_a, d k₀ with $c_{\mathrm{a},\mathrm{d}}=\parallel \mathit{B}\parallel /\sqrt{4\pi {\rho }_{\mathrm{d}}}$ being the dust Alfvén speed. Indeed, we account for the fact that the magnetized fluid sensitive to the instability is the pressureless dust, and set β=0. Although the study of a two-fluid parametric instability is beyond the scope of this work, Eq. (C.1) is relevant when friction between gas and dust is negligible, that is when considering an Alfvén wave frequency large with respect to the collision frequency of the two fluids. This condition is met in Fig. C.1 where the growth rate of the fastest growing (unstable) mode computed from the roots of Eq. (C.1) is displayed as a dashed gray line.

The analytic growth rate is little sensitive to the wavelength of the mother Alfvén wave and we write the dust density perturbation as δ ρ_d ∝ exp (t / τ) where we computed the growth rate τ ≃ 0.015 Alfvén crossing-time. We see that there is a good match between the analytic prediction and the numerical results, reinforcing the idea that a parametric-like instability is at play. The modest deviation comes from the finite friction between the charged dust and the neutral gas and the fact that initial perturbations from the background state are induced by numerical noise (meaning that several modes are probably triggered).

Fig. C.1 Comparison between the analytic growth rate of the parametric instability (dashed gray line) and numerical results of the two-fluid parametric-like instability presented in Sect. 4.1 in the dusty ideal MHD model for different values of the mother Alfvén wave wavelength (same configuration as in Fig. 1).

Appendix D Numerical treatment of non-ideal MHD terms

We show here the treatment made for the diffusive ohmic dissipation and the dispersive Hall effect appearing in the right-hand side of the induction equation (Eq. 3). In our cartesian coordinate system, The Ohm-like term is written $$\left(\begin{array}{l}\mathrm{\partial }{B}_x/\mathrm{\partial }t\\ \mathrm{\partial }{B}_y/\mathrm{\partial }t\\ \mathrm{\partial }{B}_z/\mathrm{\partial }t\end{array}\right)=\left(\begin{array}{c}0\\ {\mathrm{\partial }}_x\left(\frac{c\parallel \mathit{B}\parallel }{4\pi {\mathrm{\Gamma }}_i{e}_i}{\mathrm{\partial }}_x{B}_y\right)\\ {\mathrm{\partial }}_x\left(\frac{c\parallel \mathit{B}}{4\pi {\mathrm{\Gamma }}_{i}e{n}_i}{\mathrm{\partial }}_x{B}_z\right)\end{array}\right),$$(D.1)

where $\parallel \mathit{B}\parallel =\sqrt{B_x^2+B_y^2+B_z^2}$ denotes the norm of the magnetic field. This is a diffusion equation which is nonlinear owing to the diffusion coefficients (magnetic resistivity) being spacedependent. Figure D.1 provides a comparison of the numerical results with the analytical solution at different times for the following simple nonlinear diffusion equation $$\mathrm{\partial }{B}_y/\mathrm{\partial }t={\mathrm{\partial }}_x\left(2B_y{\mathrm{\partial }}_x{B}_y\right).$$(D.2)

A self similar solution to this equation is the Barenblatt-Pattle solution (Barenblatt 1952) that writes $$B_y(t,x,C)=(2t{)}^{-1/3}(C-\frac{1}{6}\frac{x^2}{(2t{)}^{2/3}}),$$(D.3)

with the constant $C={\left(B_{y,0}{x}_c/\sqrt{6}\right)}^{2/3}$. At initial time $t_0=C^3/2B_{y,0}^3$, we consider the following initial profile $$B_y(t_0,x)=B_{y,0}(1-{\left(\frac{x}{x_c}\right)}^2).$$(D.4)

There is a great agreement between the numerical evolution of the profile and the analytical prediction, showing that the solver works properly and is suited to our numerical investigation. Now, we turn to the Hall-like term which is expressed as $$\left(\begin{array}{l}\mathrm{\partial }{B}_x/\mathrm{\partial }t\\ \mathrm{\partial }{B}_y/\mathrm{\partial }t\\ \mathrm{\partial }{B}_z/\mathrm{\partial }t\end{array}\right)=\left(\begin{array}{c}0\\ {\mathrm{\partial }}_x\left(\frac{c{B}_x}{4\pi e_i}{\mathrm{\partial }}_x{B}_z\right)\\ -{\mathrm{\partial }}_x\left(\frac{c{B}_x}{4\pi e{n}_i}{\mathrm{\partial }}_x{B}_y\right)\end{array}\right).$$(D.5)

This term cannot be treated as a diffusion equation since B_y and B_z are coupled. Instead, we adopt the approach presented in Marchand et al. (2018) and Marchand et al. (2019) and rewrite it in conservative form using Maxwell-Ampère equation (Eq. 5) $$\left(\begin{array}{l}\mathrm{\partial }{B}_x/\mathrm{\partial }t\\ \mathrm{\partial }{B}_y/\mathrm{\partial }t\\ \mathrm{\partial }{B}_z/\mathrm{\partial }t\end{array}\right)=\left(\begin{array}{c}0\\ -{\mathrm{\partial }}_x\left(\frac{B_x}{e{n}_1}{J}_y\right)\\ -{\mathrm{\partial }}_x\left(\frac{B_x}{e{n}_1}{J}_z\right)\end{array}\right),$$(D.6)

and include it as a flux in a separate Riemann solver dedicated to the induction equation (Eq. 3). The Hall effect is known to introduce new dispersive waves in the system of hyperbolic equations, namely Whistler waves. The dispersive nature of this term compels us to cut-off the velocity of Whistler waves on a scale that is determined by the resolution, that is the cell size $\mathrm{\Delta }x=\frac{NX}{L}$. The Whistler velocity is thus defined as $$c_w=\frac{{\eta }_{\mathrm{H}}\pi }{2\mathrm{\Delta }x}+\sqrt{{\left(\frac{{\eta }_{\mathrm{H}}\pi }{2\mathrm{\Delta }x}\right)}^2+c_a^2}$$(D.7)

where the Hall resistivity ${\eta }_{\mathrm{H}}=\frac{c{B}_x}{4\pi e{n}_{\mathrm{i}}}$. We incorporate it in the wave fan of the Riemann solver of the magnetic field, but do not use it to update non-magnetic flow variables. As pointed out in Marchand et al. (2019), doing otherwise would introduce significant truncation errors and would thus lead to a much lower effective resolution.

Fig. D.1 Comparison between analytical solution and numerical resolution of the nonlinear diffusion equation given by Eq. (D.2). The analytical solution (Barenblatt 1952) is given by Eq. (D.3). Each curve shows By profile at a given time.

Appendix E Convergence test

In Fig. E.1, we show a convergence test by plotting the dust and gas average density fluctuations with increasing resolution. While resolution has no effect on the gas density, convergence is not reached in the dust even for a number of cells as high as NX=8192. This results from the absence of any pressure in the dust, which is allowed to form smaller and smaller clumps. This further supports the capacity of magnetic effects to induce substantial dust clumping, which could be significantly higher than that shown in the present work. However, on sufficiently small scale, we would expect a remnant dust velocity dispersion to act as a pressure and halt further increase in dust density.

Fig. E.1 Convergence test. Average dust density (solid) and gas density (dotted) (as defined in Eq. 18) as a function of time for different spatial resolutions. β=1 and sd=100 μm(St=0.1) were used. The vertical dashed black line indicates the time at which Fig. E.2 was made. While convergence is reached for the gas, it is not for the dust.

Fig. E.2 Dust (solid lines) and gas (dotted lines) density fluctuation field throughout the box at t=1.6 (see Fig. E.1) for different resolutions.

Appendix F Ideal vs non-ideal dusty MHD

A comparison between dusty ideal MHD model (Sect. 2.2) and dusty non-ideal MHD model (Sect. 2.1) is depicted in Fig. F. When plasma parameter β=0.1, there is no difference between both models. Indeed, the fiducial transverse Mach number ${\mathcal{M}}_{\perp }$ = 0.99 produces too weak transverse-to-longitudinal magnetic field ratio B_⊥ / B_∥, keeping the parametric instability from inducing a substantial dust concentration beyond that of pure hydrodynamical turbulent processes. However, when going to β=0.7, significant deviation is observed between the two models: maximum dust density fluctuations are one order of magnitude apart. Magnetic drag in Eq. (6) due to the presence of ions is to be held responsible for the reduction in dust concentration within the dusty non-ideal regime. Consequently, an accurate treatment of non-ideal MHD effects cannot be ignored and have to be included in numerical simulations.

Fig. F.1 Comparison of time evolution of the maximum dust density within dusty ideal MHD setup (dashed line. See Sect. 2.2) and dusty non-ideal MHD (solid line. See Sect. 2.1). Plasma parameter β=0.1 (blue) and β=0.7 (purple). Dust size: sd= 10 μm.

Appendix G How the plasma parameter β controls the transverse-to-longitudinal magnetic ratio B_⊥ / B_∥

Figure G.1 depicts the influence of the of the plasma parameter β on the maximum transverse-to-longitudinal magnetic ratio B_⊥ / B_∥ in simulations with turbulence. We clearly see that within our setup, for a given level of turbulence (${\mathcal{M}}_{\perp }$), a higher β (i.e., lower longitudinal B_x) leads to higher values of B_⊥ / B_∥ which in turn leads to higher dust density fluctuations (see Sect. 4.1).

Fig. G.1 Maximum transverse-to-longitudinal magnetic ratio as a function of time for different values of plasma parameter β in simulations with turbulence.

Appendix H Exploring the effect of a dust pressure on density fluctuations

Figure H.1 shows the impact of a dust pressure added a posteriori in the corresponding momentum equation. Parametrizing the dust sound speed as a fraction of that of the gas, we see no effects on dust density fluctuations for values c_{s, d} < 0.5 c_s for the spatial resolution chosen (NX=4096). As discussed in Sect. 5.3, such high values seem unrealistic for the grain sizes explored in this work, meaning that our results are most likely robust.

Fig. H.1 Maximum dust density (solid lines) and gas density (dotted lines) fluctuations as a fraction of time for different values of the dust sound speed (defined as a function of the gas sound speed cs) in simulations with turbulence (the red and black curves are perfectly overlain). Note that the spatial resolution is NX =4096.

Appendix I Chemical network comparison

We provide here a comparison of the chemical network used in this paper (see Sect. 2.3) with the more sophisticated one presented in Marchand et al. (2021). The latter includes the effects of a dust size distribution, with small grains offering a large cross section for reactions such as electron capture and ion-electron recombination at their surface. Figure I.1 shows the ion abundance n_i and ionization fraction x_i=n_i / n_H as a function of gas density n_H. Although both networks seem to agree at lower density (n_H < 10⁵ cm⁻³), discrepancies are exacerbated as density increases.

Appendix J Impact of parameters on density fluctuations for a lower plasma parameter β = 0.1

Similarly to what was done in Sect. 4.2, Fig. J.1 and Fig. J.2 depict the influence of dust size, Mach number, background gas density and plasma parameter on maximum and average density fluctuations as a function of time, however for a lower fiducial value of plasma parameter β=0.1. The fiducial values for the other parameters are those presented in Table 1. As explained in Sect. 4.1 and Sect. 5.2.1, a lower β (i.e., a higher background B_x) makes the development of high transverse-to-longitudinal magnetic field ratios B_⊥ / B_∥ more challenging. As a consequence, magnetic effects are incapable of pushing dust concentration levels beyond what is permitted by pure hydrodynamical effects. In addition, there is a loss of dust clumping dependency on dust size because of ion-induced magnetic drag being insensitive to this parameter, as seen in Eq. (6).

Fig. I.1 Ion density (solid lines, left axis) and ionization fraction ni / nH (dotted lines, right axis) as a function of gas density. Comparison between our chemical model (Sect. 2.3, Shu et al. 1987) and that of Marchand et al. (2021). For the latter, an MRN dust size distribution was considered. For both, the ionization rate by cosmic rate is ζ=10−17 s−1.

Fig. J.1 Dust density (solid lines) and gas density (dotted lines) maximum fluctuations as a function of time for different values of the dust grain size sd, longitudinal Mach number (perpendicular one being varied too) ${\mathcal{M}}_{\parallel }$, gas initial density nH and plasma parameter β. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines. When not varied, β=0.1.

Fig. J.2 Same as Fig. J.1 but depicting dust-to-gas ratio average (dust density weighted), dust density average (dust density weighted) and gas density average (gas density weighted).

All Tables

All Figures

Fig. 1 Maximum dust density fluctuations and maximum transverse magnetic field as a function of time upon propagation of a circularly polarized Alfvén wave, for different values of wavelength (as a fraction of the box length L) and wave amplitudes δB within the dusty ideal MHD regime and dusty non-ideal MHD regime (see Sects. 2.2 and 2.1). The early sharp rise in dust density is due to a mechanism similar to the parametric instability (Del Zanna et al. 2001). Note that growth rates measured here compare well with analytical predictions of the standard parametric instability (see Appendix C). Afterwards, the dust density decreases as a consequence of gas dust friction. The gas density is not displayed because subject to negligible variations. Parameters: nH=106 cm−3, sd=10 μm(St = 0.01) and β=0.1.

In the text

	Fig. 2 Different fields at times t=0.23, t=0.44, and t=0.82 in a simulation with driven turbulence. First panel: dust (solid lines) and gas (dotted lines) density fluctuations. Second panel: x-wise dust (solid) and gas (dotted) Mach number. Third panel: magnetic pressure gradient. Fourth panel: total (x, y, z-wise) dust Alfvénic Mach number. Parameters: nH=106 cm−3, sd=10 μm(St=0.01) and β=1.
In the text

	Fig. 3 Schematic illustration of the parameters controlling the transverse-to-longitudinal magnetic field ratio $B_{\perp }/B_{\parallel }$ and thus the intensity of dust clumping via the parametric instability in the simulations with turbulence.
In the text

Fig. 4 Maximum dust density (solid lines) and gas density (dotted lines) fluctuations as a function of time for different values of the dust grain size sd, longitudinal Mach number (perpendicular one being varied too) ${\mathcal{M}}_{\parallel }$, gas initial density nH, and plasma parameter β. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines. When not varied, β=0.7.

In the text

	Fig. 5 Same as Fig. 4 but depicting space averages (mass weighted) instead of maximum values. On the first row is displayed the dust-to-gas ratio (dust density weighted), on the second row the dust density (dust density weighted) and on the third row the gas density average (gas density weighted).
In the text

	Fig. 6 Time-averaged (over driven turbulence steady-sate regime) probability distribution functions (PDFs) of dust-to-gas ratio fluctuations (first row), gas density fluctuations (second row) and dust density fluctuations (third row).
In the text

Fig. 7 Average (dust density weighted) dust density (solid lines) and average (gas density weighted) gas density (dotted lines) fluctuations as a function of time for different values of the dust grain size sd. Upper left panel: β=0.1. Lower left panel: β=1. Second column: same but with ten times less ions. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines.

In the text

	Fig. 8 Schematic illustration of a possible coagulation instability in turbulent magnetized dense cores. Dust concentration leads to a more rapid dust growth, which in turn results in larger grains prone to stronger concentration.
In the text

	Fig. B.1 Dispersion relation for an Alfvén wave in a dusty magnetized medium within the dusty ideal MHD model (Sect. 2.2). Comparison between the analytical solution and numerical measurement from our simulation. The real part of the frequency appears as solid line and the imaginary part as dashed line. St =0.001, β=0.1 and δB=10−4.
In the text

	Fig. C.1 Comparison between the analytic growth rate of the parametric instability (dashed gray line) and numerical results of the two-fluid parametric-like instability presented in Sect. 4.1 in the dusty ideal MHD model for different values of the mother Alfvén wave wavelength (same configuration as in Fig. 1).
In the text

	Fig. D.1 Comparison between analytical solution and numerical resolution of the nonlinear diffusion equation given by Eq. (D.2). The analytical solution (Barenblatt 1952) is given by Eq. (D.3). Each curve shows By profile at a given time.
In the text

	Fig. E.1 Convergence test. Average dust density (solid) and gas density (dotted) (as defined in Eq. 18) as a function of time for different spatial resolutions. β=1 and sd=100 μm(St=0.1) were used. The vertical dashed black line indicates the time at which Fig. E.2 was made. While convergence is reached for the gas, it is not for the dust.
In the text

	Fig. E.2 Dust (solid lines) and gas (dotted lines) density fluctuation field throughout the box at t=1.6 (see Fig. E.1) for different resolutions.
In the text

	Fig. F.1 Comparison of time evolution of the maximum dust density within dusty ideal MHD setup (dashed line. See Sect. 2.2) and dusty non-ideal MHD (solid line. See Sect. 2.1). Plasma parameter β=0.1 (blue) and β=0.7 (purple). Dust size: sd= 10 μm.
In the text

	Fig. G.1 Maximum transverse-to-longitudinal magnetic ratio as a function of time for different values of plasma parameter β in simulations with turbulence.
In the text

	Fig. H.1 Maximum dust density (solid lines) and gas density (dotted lines) fluctuations as a fraction of time for different values of the dust sound speed (defined as a function of the gas sound speed cs) in simulations with turbulence (the red and black curves are perfectly overlain). Note that the spatial resolution is NX =4096.
In the text

	Fig. I.1 Ion density (solid lines, left axis) and ionization fraction ni / nH (dotted lines, right axis) as a function of gas density. Comparison between our chemical model (Sect. 2.3, Shu et al. 1987) and that of Marchand et al. (2021). For the latter, an MRN dust size distribution was considered. For both, the ionization rate by cosmic rate is ζ=10−17 s−1.
In the text

Fig. J.1 Dust density (solid lines) and gas density (dotted lines) maximum fluctuations as a function of time for different values of the dust grain size sd, longitudinal Mach number (perpendicular one being varied too) ${\mathcal{M}}_{\parallel }$, gas initial density nH and plasma parameter β. The dust density mean value for pure hydrodynamics simulations is displayed for reference as dashed lines. When not varied, β=0.1.

In the text

	Fig. J.2 Same as Fig. J.1 but depicting dust-to-gas ratio average (dust density weighted), dust density average (dust density weighted) and gas density average (gas density weighted).
In the text