# A new multifluid method for dusty astrophysical flows

G. Verrier, U. Lebreuilly, P. Hennebelle

► To cite this version:

G. Verrier, U. Lebreuilly, P. Hennebelle. A new multifluid method for dusty astrophysical flows. Astronomy & Astrophysics - A&A, 2025, 701, pp.A174. 10.1051/0004-6361/202554662. cea-05262653

HAL Id: cea-05262653

https://cea.hal.science/cea-05262653v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# A new multifluid method for dusty astrophysical flows
## Application to turbulent protostellar collapses

G. Verrier<fnref n="1" />,*, D. U. Lebreuilly<fnref n="2" />, and P. Hennebelle<fnref n="2" />

Received 20 March 2025 / Accepted 11 July 2025

## ABSTRACT

**Context.** Stars and planets form in collapsing clouds of gas and dust. The presence of dust grains and their local distribution play a significant role throughout the protostellar sequence, from the thermodynamics and the chemistry of molecular clouds to the opacity of collapsing protostellar cores and the coupling between the gas and the magnetic field and down to planet formation in young and evolved disks.

**Aims.** We aim to simulate the dynamics of the dust, considering the whole range of grain sizes, from few nanometers to millimeters.

**Methods.** We implemented a neutral pressureless multifluid that samples the dust size distribution in the RAMSES code. This multi-fluid is dynamically coupled to the gas via a drag source term and self-gravity, relying on the Eulerian approach.

**Results.** We designed a Riemann solver for the gas and dust mixture that prevents unphysical dust-to-gas ratio variations for well-coupled grains. We illustrated the capacities of the code by performing simulations of a protostellar collapse down to the formation of a first hydrostatic core, both for small and large dust grains. Grains over 100 microns significantly decouple from the gas. The spatial maps and the probability density functions indicate that dust enrichment within the first hydrostatic core and in some locations of the envelope increases as a function of the grain size and the level of initial turbulence.

**Conclusions.** Thanks to the novel Riemann solver, we recovered the terminal velocity regime, even at low resolution. Moreover, we successfully extended it to regimes where the grain inertia matters. The multifluid module performs the coupling between the dust and the gas self-consistently all through the dynamical scales. The dust enrichment in the first hydrostatic core and the envelope have been revised here, assuming the initial turbulence and grain sizes. This enables us to probe new potential locations, epochs, and initial conditions for planet formation.

**Key words.** hydrodynamics – turbulence – methods: numerical – stars: formation – dust, extinction

## 1. Introduction

Dust grains provide the solid material to form planets. However, the questions of how and when the first planetesimals form are still a matter of discussion (Drażkowska et al. 2023). Star formation and planetary formation are intimately connected by dust evolution. Even though dust represents a small fraction of the mass (about 1 percent), its dynamics strongly affects the coupling of the bulk mass (i.e., the gas) to the magnetic field, which is a prime mechanism in regulating angular momentum during a protostellar collapse (Maury et al. 2022; Tsukamoto et al. 2023). Indeed, dust grains carry most of the charges and, thus, they control the magnetic resistivities (Marchand et al. 2016; Zhao et al. 2016). Moreover, because the dust opacity mainly controls the optical thickness, it determines how much energy is radiated away during the collapse (Gaustad 1963; Larson 1969). Finally, dust grains allow for recombination and chemistry on their surfaces, making them a preferential channel for formation for abundant molecules such as H₂ (Gould & Salpeter 1963).

The multiple roles of the dust depend on the local size distribution and dust-to-gas ratio, which define, for instance, the available surface for chemistry and the mass for dynamical instabilities, such as the polydisperse streaming instability (Krapp et al. 2019). Mathis, Rumpl, and Nordsieck (hereafter MRN)

found, by fitting the interstellar extinction, a power-law size distribution ranging from 5 to 250 nm (Mathis et al. 1977), whereby large grains carry most of the dust mass and the small grains vastly outnumber large grains and provide most of the surface. Such conclusions rely on the optical and thermal properties of dust grains, resulting from the grain size, structure (compact or aggregate), and composition (grains are thought to be mainly carbonaceous and silicated). Efforts have to been made to better identify and trace the dust distribution, for example, by including the spectral energy distribution from Herschel and Planck (Compiègne et al. 2011) as well as laboratory measurements (Ysard et al. 2024). However, it remains poorly constrained during stellar formation.

In most astrophysical fluid codes, dust is accounted for in the bulk mass. We can assume, for instance, that the gas and dust mixture is perfectly coupled with constant dust-to-gas ratio and dust distribution to compute the opacity and the evolution of the radiation field (Muley et al. 2023) or the evolution of the magnetic field via magnetic resistivities (Wurster et al. 2016). The interest in the dust dynamics in the context of protoplanetary disk evolution and planet formation has greatly motivated dedicated numerical developments. The most frequently used codes include PHANTOM (Price et al. 2017), FARGO3D (Benítez-Llambay & Masset 2016), Athena++ (Stone et al. 2020), and Idefix (Lesur et al. 2023). They can rely on different approaches such as Lagrangian particles in SPH codes

<footnote n="1">¹ Université Paris Cité, Université Paris-Saclay, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author: gabriel.verrier@cea.fr</footnote>
as implemented in PHANTOM, or as a multifluid in grid-based codes as implemented in FARGO3D (Benítez-Llambay et al. 2019) and Athena++ (Huang & Bai 2022). More specifically, Benítez-Llambay et al. (2019) implemented a first-order implicit drag scheme whose direct inversion for noncolliding dust species is presented in Krapp & Benítez-Llambay (2020). Subsequently, Huang & Bai (2022) implemented two second-order schemes that remain stable in stiff regimes. Developing higher-order schemes coupling the hydrodynamics and the drag source term is an active area of research (Keppens et al. 2023; Krapp et al. 2024). Moreover, it can benefit high-performance from adaptive mesh refinement (AMR) in Stone et al. (2020); Keppens et al. (2023) and GPU parallelization in Benítez-Llambay & Masset (2016); Lesur et al. (2023).

In this work, we address the implementation of a multifluid solver in the RAMSES code (Teyssier 2002). RAMSES is a 3D finite-volume AMR code. The adaptive mesh refinement (AMR), together with subcycling, are essential to performing stellar formation simulations. Moreover, many physics modules have been successively implemented. Among the most useful ones in star formation simulations are radiative transfer (Commerçon et al. 2011; Rosdahl & Teyssier 2015; González et al. 2015), ideal magnetohydrodynamics (Fromang et al. 2006), and standard non-ideal magnetohydrodynamics (MHD), namely, Ohm diffusion and ambipolar drift (Masson et al. 2012), as well as the Hall effect (Marchand et al. 2018). For these reasons, RAMSES is a powerful tool to deal with the whole range of densities and scales in which dust plays a role (see, e.g., Ahmad et al. (2023) for the dynamical range to reach the birth of a protostar from a protostellar collapse). In particular, this could capture multi-scale and multi-physics processes such as infall, turbulence, and other transport mechanisms.

The first implementation of the dust dynamics in RAMSES was done by Lebreuilly et al. (2019). The dust dynamics is computed in the terminal velocity approximation, whose validity fails for strongly decoupled grains, which is likely to occur for (dynamically) large grains or charged grains. Moreover, this approximation fails in shocks (Lovascio & Paardekooper 2019).

Recently, charged dust has been introduced as massive super-particles (particle-in-cell method) coupled to the ideal MHD equations (Moseley et al. 2023).

The paper is organized as follows. We present the methods in Sect. 2 and the validation tests in Sect. 3. In particular, we include tests with multiple dust species and various coupling regimes. We compare to the terminal velocity approximation in protostellar collapses of dense cores and assess the limitations of the approaches in Sects. 4.1 and 4.2. In Sect. 4.3, we evaluate dust enrichment in the first hydrostatic core and the envelope as a function of the grain size and the initial turbulence. Section 5 presents our conclusions. More details on the numerical results are given in the appendix.

## 2. Methods

### 2.1. Equations

We consider a multifluid of $N$ dust species, each one referred by the index $d$:

$$ \partial_t \rho_d + \nabla \cdot (\rho_d \mathbf{V}_d) = 0, \qquad (1) $$

$$ \partial_t (\rho_d \mathbf{V}_d) + \nabla \cdot (\rho_d \mathbf{V}_d \mathbf{V}_d) = -\rho_d \nabla \phi + \mathbf{f}_{g \to d}. \qquad (2) $$

Each dust species is coupled to the gas via an individual drag force $\mathbf{f}_{g \to d}$. Then, $\phi$ is the gravitational potential. Contrary to

the dust fluids, the gas is pressure-supported and, thus, the gas energy equation is required. The equations of the gas dynamics are

$$ \partial_t \rho_g + \nabla \cdot (\rho_g \mathbf{V}_g) = 0, \qquad (3) $$

$$ \partial_t (\rho_g \mathbf{V}_g) + \nabla \cdot (\rho_g \mathbf{V}_g \mathbf{V}_g + P_g \mathbf{1}) = -\rho_g \nabla \phi + \sum_d \mathbf{f}_{d \to g}, \qquad (4) $$

$$ \partial_t E_g + \nabla \cdot \left( (E_g + P_g) \mathbf{V}_g \right) = -\rho_g \nabla \phi \cdot \mathbf{V}_g + \sum_d \mathbf{f}_{d \to g} \cdot \mathbf{V}_g + Q. \qquad (5) $$

Each dust fluid back-reacts on the gas thus $\mathbf{f}_{g \to d} = -\mathbf{f}_{d \to g}$. We can express the drag force as

$$ \mathbf{f}_{g \to d} = \frac{\rho_d}{t_{s,d}} (\mathbf{V}_g - \mathbf{V}_d), \qquad (6) $$

defining $t_{s,d}$ as the stopping time of the dust grain in the gas fluid. In the Epstein regime (Epstein 1924), it is

$$ t_{s,d} = \sqrt{\frac{\pi \gamma \rho_{\text{grain,d}}}{8} \frac{s_{\text{grain,d}}}{\rho_g} \frac{c_s}{c_s}}, \qquad (7) $$

where $c_s$ is the sound speed of the gas, $\gamma$ its adiabatic index and $s_{\text{grain,d}}$ is the size (radius) of the grain, while $\rho_{\text{grain,d}}$ its intrinsic density.

Here, $E_g$ and $P_g$ are the total energy and the pressure of the gas, respectively, while $Q$ is the energy deposit in the gas due to frictional heating. When it is not set to zero, to conserve the energy of the gas-dust multifluid system, we can express it as

$$ Q = - \sum_d (\mathbf{f}_{g \to d} \cdot \mathbf{V}_d + \mathbf{f}_{d \to g} \cdot \mathbf{V}_g) = \sum_d \frac{\rho_d}{t_{s,d}} (\mathbf{V}_g - \mathbf{V}_d)^2. \qquad (8) $$

Finally, the gravitational potential, which both the gas and the dust multifluid contribute to, can be obtained using the Poisson equation

$$ \Delta \phi = 4\pi G \left( \rho_g + \sum_d \rho_d \right). \qquad (9) $$

### 2.2. Operator splitting

Equations (1)–(5) can be written as follows:

$$ \partial_t \mathbf{U}_g + \nabla \cdot \mathbf{F}_g (\mathbf{U}_g) = \mathbf{S}_{\text{grav,g}} + \mathbf{S}_{\text{drag,g}}, \qquad (10) $$

$$ \partial_t \mathbf{U}_d + \nabla \cdot \mathbf{F}_d (\mathbf{U}_d) = \mathbf{S}_{\text{grav,d}} + \mathbf{S}_{\text{drag,d}}, \quad d \in \{1, \dots, N\}, \qquad (11) $$

where $\mathbf{U}_g$ and $(\mathbf{U}_d)_{d \in \{1, \dots, N\}}$ are the conserved variables of the gas and of the dust multifluid. Then, $\mathbf{F}_g$ and $(\mathbf{F}_d)_{d \in \{1, \dots, N\}}$ are the fluxes corresponding to the hydrodynamics of the conserved variables, $\mathbf{S}_{\text{grav,g}}$ and $(\mathbf{S}_{\text{grav,d}})_{d \in \{1, \dots, N\}}$ the gravity source terms, and $\mathbf{S}_{\text{drag,g}}$ and $(\mathbf{S}_{\text{drag,d}})_{d \in \{1, \dots, N\}}$ the drag source terms.

The hydrodynamical scheme and the gravity source term are already implemented in RAMSES (Teyssier 2002). We include the drag source term using a Lie-Trotter splitting, as explained hereafter, carrying out the computation of the whole step over $\Delta t$ as two successive steps over $\Delta t$.

The hydrodynamical step is performed by an unsplit second-order Godunov integrator using the MUSCL-HANCOCK
scheme with various slope limiters, including minmod (Roe 1986), superbee (Roe 1986), and Van-Leer (van Leer 1974). When using multiple levels of mesh refinement, multifluid variables can be restricted from fine levels to coarse levels by averaging down and they can be prolonged from coarse levels to fine levels using different linear interpolation strategies (minmod slope, monotonized central slope, unlimited central slope). This prolongation can be applied either to the conservative variables or to the primitive variables, which are ($ρ_d, ρ_dV_d$) and ($ρ_d, V_d$) respectively. Various Riemann solvers are implemented, as described in the next section. The gravity source term is added to the hydrodynamical term, following a second-order midpoint scheme.

After computing the hydrodynamical and gravity steps, the drag step is performed using the first-order implicit scheme addressed in Krapp & Benítez-Llambay (2020) on the velocity of the gas and the dust species. In the absence of frictional heating, the kinetic energy of the gas is simply updated. To account for the frictional heating, we remove the increment of the kinetic energy in the dust multifluid to the total energy of the gas. We test the order of the scheme in Sect. 3.2.

## 2.3. Riemann solver for a gas and dust multifluid mixture

The Riemann solver of the multifluid requires special care, depending on the coupling of the dust with the gas. Usually, multifluid Riemann solvers split into independent problems for each fluid. This strategy is entirely valid in the absence of interactions between the fluids. However, when dealing with infinitely low Stokes grains (dynamically small grains), we could not recover dust-to-gas ratio variations vanishing to zero (see Sect. 4.1), while the drag solver was correctly adapting the dust velocity to the one of the gas. Therefore, this numerical difference only originates from the advection of the dust density, which comes from the computation of Riemann fluxes. As far as we know, this problem and alternative strategies, such as using individual Riemann solvers or one common Riemann solver, have not been discussed in previous dust multifluid implementations. Our new solving strategy for interacting gas and dust is based on the propagation of waves in the mixture.

In the absence of source terms, for a pressureless dust fluid, the Riemann problem leads to solutions of delta shocks and vacuum states (LeVeque 2004). We implemented the Riemann solver of Huang & Bai (2022), based on upwind fluxes, as a reference solver used in other codes. We implemented a local Lax-Friedrichs (LLF) solver, denoted as LLFd and presented in Appendix A.

When the dust is well-coupled to the gas, the mixture behaves as a single fluid and, thus, the same waves propagate in the dust and in the gas. In the Riemann solver presented hereafter, we consider the same wave fan for the gas and for the dust fluid, based on the HLL approach, to model this strong coupling situation. It can switch to the individual local Lax-Friedrichs solver for a specific dust fluid if a specific decoupling criterion is met. Such a criterion is based on the dust fluid kinematics and its purpose is to inject the right amount of numerical viscosity depending on the situation. Considering more appropriate wave speeds in the Riemann solver should stabilize the hydrodynamical scheme.

We go on to define the solvers denoted by HLLgd and LLFgd (when the HLL solver or the LLF solver, respectively, is used for the gas). For the Riemann problem at the interface $i + 1/2$, we consider the left and right states of the dust multifluid, $U_{d,L}$ and $U_{d,R}$. The normal components of the dust velocity are $V_{d,L}$

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Modeling of the HLLgd solver of each dust fluid depending on the Riemann problem at the interface: Riemann problem on the gas in lines and kinematic coupling situation in columns. In blue, the gas wave fan defined by $S_{g,L}$ and $S_{g,R}$. With the dotted green lines, we show the dust normal velocities, $V_{d,L}$ and $V_{d,R}$. The HLLgd flux for the dust species d is described in red by its wave fan.</figcaption></figure>

and $V_{d,R}$ and the fluxes are denoted by $F_{d,L} = F_d(U_{d,L})$ and $F_{d,R} = F_d(U_{d,R})$. The Riemann problem for the gas, neglecting the interaction with the multifluid, is modeled by the two wave speeds, $S_{g,L}$ and $S_{g,R}$, for the HLL solver (Harten et al. 1983), such that $S_{g,L} < S_{g,R}$. We define the flux at the interface for each dust fluid, $F_{d,i+1/2}$, as follows.

If $S_{g,L} > 0$ (upwind wave speeds of the gas from the left), then

$$F_{d,i+1/2} = \begin{cases} F_{d,L} & \text{if } V_{d,L}, V_{d,R} > 0 \text{ (dust upwind flow),} \\ F_{LLF,d} & \text{otherwise (switch).} \end{cases} \quad (12)$$

Symmetrically, if $S_{g,R} < 0$ (upwind wave speeds of the gas from the right), then

$$F_{d,i+1/2} = \begin{cases} F_{d,R} & \text{if } V_{d,L}, V_{d,R} < 0 \text{ (dust upwind flow),} \\ F_{LLF,d} & \text{otherwise (switch).} \end{cases} \quad (13)$$

Otherwise $S_{g,L} \le 0 \le S_{g,R}$, the solution at the interface can be modeled as

$$F_{d,i+1/2} = \begin{cases} F_{HLL,d}(S_{g,L}, S_{g,R}) & \text{if } S_{g,L} < V_{d,L}, V_{d,R} < S_{g,R}, \\ F_{LLF,d} & \text{otherwise (switch),} \end{cases} \quad (14)$$

depending on whether the left and right dust velocities belong to the gas wave fan (first subcase) or not (switch).

We define the HLL flux as

$$F_{HLL,d}(s_L, s_R) = \frac{s_R F_{d,L} - s_L F_{d,R} + s_R s_L (U_{d,R} - U_{d,L})}{s_R - s_L}. \quad (15)$$

LLFgd is the particular case of HLLgd where $S_{g,R} = -S_{L,R} > 0$. The different Riemann problems with the corresponding solutions of the HLLgd solver are illustrated in Fig. 1.

We have decided to switch to a local Lax-Friedrichs solver (LLFd, see Appendix A) because of its simplicity, following the HLL approach, and robustness, and also because some situations could be difficult to model properly, for example if $V_{d,L}, V_{d,R} > 0$ and $S_{g,R} < 0$, which is when the gas and the dust flows cross each other.
We propose an interpretation of the switch criterion in the linear regime. For small perturbations in the gas and dust velocities, $\delta v_g$ and $\delta v_d$, and for a compressible wave propagating at the wave speed $c_\phi$, the equations of mass conservation lead to

$$ \frac{\delta \theta_d}{\theta_d} = \frac{\delta v_d - \delta v_g}{c_\phi}, \qquad (16) $$

where $\theta_d = \rho_d/\rho_g$ is the dust-to-gas ratio. More generally, we could expect a significant variation of the dust-to-gas ratio if $|\delta v_d - \delta v_g| > c_\phi$. This switch criterion also means that the dust species is not in the influence area of the gas, modeled by its wave fan. We note that the wave speed values used for the gas wave fan overestimate the wave speeds expected for a perfect monofluid. Indeed, the usual MHD wave speeds scales with $1/\sqrt{\rho}$, where $\rho = \rho_g + \sum_d \rho_d$ is the bulk mass carrying the wave; for example, the effective sound speed is $c_{s,eff} = c_s/\sqrt{1 + \sum_d \rho_d/\rho_g < c_s}$ (Laibe & Price 2011). This overestimation is weak for a small dust-to-gas ratio and it provides a small and safe boost of numerical viscosity.

As far as we know, adapting the Riemann solver of the multifluid to model the dynamical coupling with the gas has never been considered and thus we have evaluated the impact of the choice of the Riemann solver for studying the dust decoupling in the conditions of protostellar collapses (Sect. 4). The source terms are not explicitly included in the approximation of the solution of the Riemann problem. Jump conditions on fluxes are no longer satisfied and thus the development of an extension of the HLLgd method to a HLLD wave fan, following Miyoshi & Kusano (2005), would require a dedicated work. Because fractional step schemes can fail in capturing states that are close to equilibrium, LeVeque & Bale (1999) incorporate the source terms in the Riemann problem. However, this approach applies to the specific situation of quasi-steady states, which is not likely to occur for a vast range of grain sizes. Developing a more general and sophisticated Riemann solver is beyond the scope of this work.

## 2.4. Time-stepping

One advantage of the implemented implicit drag scheme is to be asymptotically and unconditionally stable (Krapp & Benítez-Llambay 2020). An explicit scheme would require to limit the time step to the stopping times of the dust species, which could be problematic during a protostellar collapse simulation for which the stopping times vary a lot because of the size distribution and the local density.

The Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1928) should guarantee that any information cannot leave the cell during one time step. We therefore add the classical condition for each dust species d and for each cell sharing a common subcycled domain,

$$ \Delta t < C_{CFL} \frac{\Delta x}{\max_j(|(V_{d,x})_j| + |(V_{d,y})_j| + |(V_{d,z})_j|)}, \qquad (17) $$

where $C_{CFL} < 1$ is a safety factor. In the case of a highly coupled mixture, the sound speed of the monofluid should be considered. This should not be problematic because the CFL condition on the hydrodynamical step of the gas should be sufficient, as discussed in Sect. 2.3. We have not found it necessary to include the acceleration due the drag force because this force tends to reduce the velocity drifts between the fluids.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Exact solution (solid lines) and two numerical solutions (circles and squares, respectively) of the dustybox at $x = 0.25$, relaxing towards the barycenter velocity (dotted line) over six damping times ($t = 0.06$). Equivalently, $t \approx 3t_{s,d}$, where $t_{s,d} = \rho_d/K_d \approx 0.02$ is the stopping time. Numerical solutions of the first-order implicit drag solver are computed in the stiff regime ($\Delta t = t_{s,d} = 0.02$, squares) and in the non-stiff regime ($\Delta t = 0.05t_{s,d} = 10^{-3}$, circles).</figcaption></figure>

## 3. Validation tests

We test the coupling between the hydrodynamical scheme, the drag solver, and the self-gravity solver. We test the drag solver alone with the dustybox test (Sect. 3.1) and its coupling with the hydrodynamical solver with the dustywave test (Sect. 3.2). Both tests were developed in Laibe & Price (2011) and tested in Benítez-Llambay et al. (2019). We go further in the dustywave test, by checking the time convergence (Sect. 3.2.2), the static mesh refinement, and the time subcycling (Appendix C.2). We test the self-gravity with the dustyjeanswave test (Appendix D), similarly to the test in Krapp et al. (2024). We recover the damping mode in a distribution of five dust fluids. We perform the disk settling test with ten dust fluids (Sect. 3.3), as in Hutchison et al. (2018), and the shock test (Appendix E), following Benítez-Llambay et al. (2019).

In the dustybox, dustywave, and dustyjeanswave tests, the initial conditions are

$$ \rho_g(x, 0) = (\rho_g)_0 + \delta(\rho_g)_0 \sin(kx + \phi_g), \qquad (18) $$

$$ V_g(x, 0) = (v_g)_0 \sin(kx + \psi_g). \qquad (19) $$

A similar expression stands for the dust species. The tests are isothermal, setting $P_g = C_s^2\rho_g$. The numerical domain is a periodic 1D box of a length $L = 1$. We denote $M_x = L/\Delta x = 2^l$ as the number of cells of a uniform grid of refinement level $l$ and $n = t/\Delta t$ as the number of time steps at the coarsest level of mesh refinement.

### 3.1. Dustybox

The dustybox test consists in the relaxation towards the velocity barycenter of colliding fluids. It tests the drag solver only. In
other words, the flux terms and the gravity terms in Eqs. (10) and (11) are ignored. The evolution of the conservative variables or, more precisely, the momentum of the fluids, are governed by the drag force, expressed as

$$f_{g \to d} = K_d(V_g - V_d), \qquad (20)$$

where $K_d$ is a constant drag coefficient. For two fluids, the evolution of the state $U = (\rho_g v_g, \rho_d v_d)^T$ is given by the system

$$\partial_t \begin{pmatrix} \rho_g(x)v_g(x, t) \\ \rho_d(x)v_d(x, t) \end{pmatrix} = M(x) \begin{pmatrix} \rho_g(x)v_g(x, t) \\ \rho_d(x)v_d(x, t) \end{pmatrix}, \qquad (21)$$

$$M(x) = K_d \begin{pmatrix} -1/\rho_g(x) & 1/\rho_d(x) \\ 1/\rho_g(x) & -1/\rho_d(x) \end{pmatrix}. \qquad (22)$$

The drag source term is colocated with the conservative variables, so the system Eq. (21) to solve is an ordinary differential equation (ODE) in time. The exact solution is

$$U(x, t) = \exp(M(x)t)U(x, 0). \qquad (23)$$

For only one dust species, it is easy to compute explicitly (Laibe & Price 2011). Therefore, Eq. (23) provides, over one time step $\Delta t$, an exact drag solver and we will refer to it later.

We tested the solver with the initial parameters $K_d = 50$, $k = 2\pi$, $(\rho_g)_0 = (\rho_d)_0 = 1$, $\delta(\rho_g)_0 = \delta(\rho_d)_0 = 10^{-4}$, $(v_g)_0 = (v_d)_0 = 10^{-4}$, $\phi_g = \phi_d = 0$, $\psi_g = 0$, and $\psi_d = -\pi$. We recovered the relaxation of gas and dust velocities towards the barycenter, even in stiff regimes (Fig. B.1), in agreement with the original work (Benítez-Llambay et al. 2019). In Appendix B, we go further by proving that the drag solver computes the time sequence of the first-order implicit Euler scheme and we present the solution within the entire box domain.

## 3.2. Dustywave

### 3.2.1. Linearized equations

We test the scheme resulting from the coupling between the drag solver to the hydrodynamical solver. The linearization of the equations Eqs. (1)–(4) provides the evolution of the perturbations:

$$d_t \begin{pmatrix} \delta \rho_g \\ \delta \rho_d \\ v_g \\ v_d \end{pmatrix} = \begin{pmatrix} 0 & 0 & -ik\rho_g & 0 \\ 0 & 0 & 0 & -ik\rho_d \\ -ikC_s^2/\rho_g & 0 & -K_d/\rho_g & K_d/\rho_g \\ 0 & 0 & K_d/\rho_d & -K_d/\rho_d \end{pmatrix} \begin{pmatrix} \delta \rho_g \\ \delta \rho_d \\ v_g \\ v_d \end{pmatrix}. \qquad (24)$$

### 3.2.2. Convergence

We found it necessary to test the convergence of the resulting scheme both in space and in time. Indeed, we identify potential sources of spatial error when coupling the hydrodynamical solver with the drag solver in Appendix C.1. Moreover, we check the error due to the operator splitting (in time). To do so, we define the convergence error at time $t$ on the field $U$ as

$$\text{err}(t, U, \Delta t, \Delta x) = \frac{1}{M_x} \sum_{k \in [1, M_x]} |U_{\text{num},k}^n - U_{\text{ref}}(x_{k+1/2}, t)|. \qquad (25)$$

Here, $\text{err}(t, U, \Delta t, \Delta x)$, as a function of $\Delta t$, is the $l^1$ convergence error in time and $\text{err}(t, U, \Delta t, \Delta x)$, as a function of $\Delta x$, is the Riemann sum associated to the $L^1$ spatial error. The numerical solution $U_{\text{num},k}^n$ is compared to a reference solution $U_{\text{ref}}$,

<figure><img src="image_4.png" /><figcaption>Fig. 3. Error (in time: upper panel and in space: lower panel) of the global scheme on the dust density variable with individual local Lax-Friedrichs solvers on the gas and on the dust (LLFg-LLFd). Test in the strong coupling regime (the setup parameters are the same as for the dustybox test (Sect. 3.1), expecting $\psi_g = \psi_d = 0$).</figcaption></figure>

which should be the exact solution. We use the solution of the linear system Eq. (24) for $U_{\text{ref}}$.

As illustrated in Fig. 3, as long as the time error dominates the convergence error, it scales with $\Delta t$ (upper panel), and as long as the spatial error dominates, it scales with $\Delta x^2$ (lower panel). This demonstrates the second-order in space and first-order in time accuracy of the combined scheme.

## 3.3. Disk settling

We assume a disk of gas at isothermal hydrostatic equilibrium with an analytical gravity acceleration. Dust grains can decouple from the gas and thus settle towards the mid-plane depending on their Stokes number, with the stopping time given by Eq. (7). We assume an infinitely thin radial disk slice and thus the test is in 1D. We sample a size distribution with ten dust species, ranging from 100 nm ($d = 1$) to 1 mm ($d = 10$). This is the same setup as Hutchison et al. (2018) and Lebreuilly et al. (2019), where the terminal velocity approach was used, as well as Lebreuilly et al. (2023), where a multifluid approach was used. We also used the reference solution from Hutchison et al. (2018). The HLLgdt solver provides a satisfying solution for the whole range of dynamical coupling regimes (Fig. 4).
<figure><img src="image_5.png" /><figcaption>Fig. 4. Settling test after ten orbital periods, performed with the HLLgd solver at level 11 of mesh refinement (2048 cells). The time step provided by the CFL condition is divided by 5 to reach convergence ($C_{CFL} = 0.2$). Colored dots are the numerical solutions for each dust fluid, whereas the small dots are the analytical solution.</figcaption></figure>

## 4. Protostellar collapses

The conversion of dust grains to small planetesimals could happen before the class II disk is formed (Manara et al. 2018). Consequently, we should not exclude early planet formation scenarios, and we should explore the most embedded phases, from the protostellar collapse to class 0/I disks with their envelopes. By performing numerical simulations of a protostellar collapse, we are able to provide self-consistent initial conditions for young disks. However, some ingredients remain unknown, such as the initial dust distribution and the initial level of turbulence. The role of these parameters are highly nonlinear and, thus, we vary these two parameters in simulations. In particular, we consider which grains remain coupled to the gas and where. Indeed, for a given grain size, the degree of coupling with the gas varies depending on the local density. The stopping time is shorter in high gas density regions. The multifluid implementation allows us to treat these changes of dynamical regimes self-consistently for any dust population. Ultimately, we want to explore how grains respond to turbulence, since it is a key parameter in star formation and disk formation. Thus, we investigate the effect of turbulence on the dust enrichment conditions when forming the first hydrostatic core.

We model the initial dense core as a sphere of uniform density in solid body rotation. The gas is non-ideally coupled to the magnetic field. Here, dust represents $\theta_{d,0} = 1\%$ of the initial mass of the gas. All dust mass is represented by only one fluid characterized by one grain size. The details of the numerical setup can be found in Appendix F.

The degree of coupling of a grain with the gas can be parametrized by the Stokes number, which we can define in the collapse conditions as the ratio of the stopping time of the grain to the free-fall time. Following our numerical setup, it can be evaluated as

$$St_{d,\text{ff}} \approx \frac{s_{\text{grain,d}}}{1 \text{ mm}} \left(\frac{\rho_g}{\rho_0}\right)^{-1/2}, \qquad (26)$$

where $\rho_0 \approx 2.5 \times 10^{-18} \text{ g.cm}^{-3}$ is the density of the initial core.

We only vary the grain size of the dust fluid and the initial Mach number of the turbulence. We also test for different numerical solvers.

We split our study as follows. In Sects. 4.1 and 4.2, we compare the multifluid implementation with the terminal velocity approach, respectively for small and large grains. In Sect. 4.3, we investigate the effect of the initial turbulence of the dense core and the grain size on the dynamics of the dust, and thus the dust enrichment within the first hydrostatic core and the envelope.

### 4.1. Submicron grain dynamics

When the stopping time is small compared to the dynamical time (which can be estimated in collapsing regions by the free-fall time), dust grains adapt their velocity to the gas dynamics. Consequently, when the gas and the dust behave as a single tightly coupled mixture, the terminal velocity approximation provides a first-order approximation of the velocity drift as a function of the stopping time and the difference in the accelerations between the gas and the dust induced by the force balance (Youdin & Goodman 2005; Laibe & Price 2014). From Eqs. (1)–(7) with one dust species d and by adding the Lorentz force $\mathbf{J} \times \mathbf{B}$ to the gas dynamics (MHD setup), the terminal velocity can be derived as

$$\mathbf{V}_d - \mathbf{V}_g = t_{s,d} \frac{\nabla P_g - \mathbf{J} \times \mathbf{B}}{\rho_g + \rho_d}. \qquad (27)$$

This approach allows us to reduce the number of equations we have to solve. Indeed, the continuity equation becomes sufficient to fully account for the dust dynamics. The implementation of this simplified treatment of the dust dynamics in RAMSES was presented by Lebreuilly et al. (2019). It was later used to model the dynamics of dust grains in nonturbulent protostellar collapses for multiple dust species in Lebreuilly et al. (2020). They found that grains decouple from the gas for sizes larger than a few 100 microns.

In this section, we compare the terminal velocity approximation with the multifluid implementation. In a smooth flow and for 100 nm grains (typical size to probe the dynamics of the MRN distribution, corresponding to a Stokes number of $10^{-4}$ according to Eq. (26)), the terminal velocity approximation provides a reference solution presented in Fig. 5. The dust remains very well-coupled to the gas and dust-to-gas variations are very weak (less than one percent), as underlined by the choice of the colorbar scale. In Fig. 6, we present the results from the multifluid with different Riemann solvers.

When using the HLLD-H&B22 solver (Miyoshi & Kusano 2005; Huang & Bai 2022), the dust-to-gas ratio profile is quite different from the one obtained using the terminal velocity approximation. It is spread much more around the initial 1%. Even when considering unphysical small grains (i.e., subnm grains, which should trace the gas perfectly), the dust-to-gas ratio profile does not change. Moreover, this does not change when increasing the time resolution, unless the spatial resolution is increased as well, as presented by the first column of Fig. 6. We obtain similar results when using individual local Lax-Friedrichs solvers for the gas and for the dust fluid (LLFg-LLFd, second column). This example demonstrates that this decoupling is artificial and it occurs because density fluxes are computed differently (LLFg and LLFd use independent wave speeds), even when the fluid velocities are the same. This artifact depends on the spatial resolution and we expect the total error to be dominated by the spatial error in the runs with individual Riemann solvers. Indeed, we expect the diffusion part of
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Collapse simulations with 100 nm grains using the terminal velocity at the formation of the first hydrostatic core (~60 kyr). The HLLD solver (Miyoshi & Kusano 2005) is used for the (gas-dust) monofluid. Gas density is on the left (slice), which is a reference for other simulation runs as long as the feedback of the dust remains weak. Mesh-refinement levels are displayed with the contours. White arrows indicate the gas velocity. Corresponding dust-to-gas ratio map on the right, as predicted by the terminal velocity implementation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Collapse simulations with the multifluid implementation for different Riemann solvers (columns) and for different spatial and time resolution (lines), expressed by the safety factor of the CFL condition and the number of Jeans length per cells for the mesh refinement. HLLD is from Miyoshi & Kusano (2005), H&B22 stands for the Riemann solver from Huang & Bai (2022), HLLDg-HLLgd means that we use the HLLD solver for the gas and the HLLgd solver only for the dust multifluid. We choose the same colorbar scale and box size as Fig. 5 in order to compare to the terminal velocity approximation. However, for the three first solvers, some regions are saturated: for the low-resolution runs, the dust-to-gas ratio (divided by $10^{-2}$) in these regions vary from 0.9 to 1.2 for HLLD-H&B22, from 0.75 to 1.3 for the LLFg-LLFd, and from 0.9 to 1.5 for HLLDg-HLLgd.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Direct comparison of the dust-to-gas ratio maps between the terminal velocity simulation (Fig. 5, right) and the multifluid simulation using the HLLgd solver with the CFL condition $\Delta t < 0.08t_{CFL}$ (Fig. 6, fourth column, second line). The colorbar scale is adapted from previous figures to emphasize the very small deviations to the initial dust-to-gas ratio of $\theta_d = 1\%$.</figcaption></figure>

the local Lax-Friedrichs flux, which depends on the considered wave speed (Appendix A), to vanish at high spatial resolution. The fluxes of LLFg-LLFd and LLFgd are asymptotically the same and therefore these two solvers are in agreement at high resolution.

The problem is solved when using a common Riemann solver (HLLgd, last column). We can converge in time with the same spatial resolution and we can achieve a similar result to the prediction with the terminal velocity (Fig. 5, with a more contrasted map in Fig. 7). Indeed, for the HLLgd solver, the dust-to-gas ratio relative variations are of the order of 2% for $\Delta t < 0.8\Delta t_{CFL}$ and of the order of 0.2% for $\Delta t < 0.08\Delta t_{CFL}$, which becomes negligible compared to the variations of physical interest. We recover the advection of the low Stokes grains with the gas. In the context of protostellar collapse simulations, once the artificial decoupling is fixed by the novel Riemann solver, the multifluid code would benefit from higher order-in-time schemes as currently developed in the literature. We note that using the HLLD solver for the gas and the HLLgd flux as defined in Sect. 2.3 only for the dust fluid does not provide a satisfying solution at all (third column). This is probably due to the details of the HLLD wave fan.

The resolution level required to obtain the agreement between the different Riemann solvers is difficult to estimate and to test. It depends on the performance of the Riemann solvers. One possible criterion consists in refining the mesh enough to properly follow the coupling between the dust and the gas. Typically, this could be $\Delta x < L_d \sim c_s t_{s,d}$, similarly to what was found in the context of SPH simulations of the dustywave (Laibe & Price 2012) and where $L_d$ corresponds to the recoupling length after a shock (Lovascio & Paardekooper 2019). In the context of a protostellar collapse, this length scales as $L_d \propto s_{\text{grain,d}}/\rho$, which can be more stringent than the refinement criterion on the Jeans length $L_J \sim c_s/\sqrt{G\rho}$. For 100 nm grains, this criterion is not satisfied within the whole collapse region, whereas for 1 mm grains (next section), this mainly concerns the first hydrostatic core and the forming disk (Fig. G.1). This observation depends on the mesh refinement (numerical setup in Appendix F).

## 4.2. Millimeter grain dynamics

Protostellar collapses could already host large dust grains, even though grain models remain too poorly constrained to make any firm conclusions on the size distribution. The micrometric emission and the spectroscopy of the dense regions of the molecular clouds (see Pagani et al. (2010) for the coreshine effect and

Dartois et al. (2024) for the JWST spectroscopy) suggest that dust grains can already grow and reach micron sizes before the protostellar phase, which are above the typical sizes of the MRN distribution (5 to 250 nm) measured in the diffuse interstellar medium. In the envelope of young protostars, polarized dust emission indicates the presence of grain sizes above 10 µm, found by modeling the grain alignment mechanism (Valdivia et al. 2019). The dust emissivity index is compatible with very large grains (up to (sub)millimeter sizes), but these conclusions are still under discussion, especially with respect to the optical properties of the grains, the dust growth mechanisms, and the dynamical origins and scales to reach such sizes (Galametz et al. 2019; Cacciapuoti et al. 2023). Multifluid simulations accounting for the full environment (envelope, infall, core, disk, outflow) can greatly help in clarifying the latest aspect.

By selecting only 1 mm grains for protostellar collapse simulations, we explore the regime where accounting for the full inertia of dust grains matters. Indeed, the Stokes number is initially close to unity according to Eq. (26), meaning that the stopping time of a 1 mm grain is comparable to the free-fall time (~42 kyr).

As presented in Fig. 8, the velocity drift between the dust and the gas is comparable to the velocity of the gas in the lowest density regions, which breaks the terminal velocity approximation. The dust recouples to the gas in the regions of high density only (the disk and the first hydrostatic core). This is properly captured by the HLLgd solver (Appendix G, with comparison with the H&B22 solver). The results of the terminal velocity approximation lack of robustness and depend on the control parameters of the run (Stokes number and velocity drift limiters to avoid unrealistic dust ratios; here $St_{max} = 0.3$ and $||V_d - V_g|| < 5 \times 10^5$ cm.s$^{-1}$, respectively). Therefore, the dust-to-gas ratio features are different between the two runs in the lowest density regions. More precisely, the terminal velocity approximation produces dust-to-gas ratio waves escaping the disk (Fig. 9, lower panels) and a different rotating and settling envelope (Fig. 9, upper panels). Still, as captured by the terminal velocity approximation, the pressure gradients (and the Lorentz force) seem to be the main mechanism for dust enrichment.

In these strong decoupling regions, the velocity drift is of the same order of magnitude as the sound speed. We note that correcting the stopping time, following the works of Kwok (1975) and Draine & Salpeter (1979), could partially limit such large drifts. Moreover, when using only one pressureless fluid, we cannot track a large velocity dispersion between grains, which is likely to occur for high Stokes grains in the turbulent motion of the gas. For example, the model of Ormel & Cuzzi (2007) predicts that the velocity dispersion is of the order of the gas velocity for grains whose Stokes number is close to unity.

## 4.3. Dust in turbulent protostellar collapses

Measurements of molecular lines suggest that there are internal turbulent motions within protostellar dense cores which are typically subsonic (Barranco & Goodman 1998; André et al. 2007). Here, we investigate the impact of the initial turbulence on local dust enrichment during the protostellar collapse. In addition to the initial solid body rotation, we apply a velocity field generated in the Fourier space by random phases according to the $-5/3$ power law. We varied the initial turbulent Mach $M_i \in \{0, 0.5, 1\}$, defined by the velocity dispersion divided by the sound speed. We performed simulations for one dust fluid representing grain sizes of $s_{\text{grain}} \in \{10 \, \mu\text{m}, 100 \, \mu\text{m}, 1 \, \text{mm}\}$. We set the intrinsic density of the grains to $\rho_{\text{grain}} = 1 \, \text{g.cm}^{-3}$ in this paper, but more
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Multifluid simulation (HLLgd) of 1 mm grains. Gas density profile with arrows representing the gas velocities (left) and relative velocity drift $||V_d - V_g||/||V_g||$ (right).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Comparison at a similar epoch (59.2 kyr, i.e. 0.8 kyr after the formation of the first hydrostatic core) between the terminal velocity approximation (left) and the multifluid (right). The snapshot of the multifluid is the same as in Fig. 8. Dust-to gas ratio for 1 mm grains.</figcaption></figure>

compact grains reach typically $\rho_{\text{grain}} = 3 \text{ g.cm}^{-3}$. We should keep in mind that more compact grains are equivalent to larger grains because of the expression of the stopping time in Eq. (7).

We found out that grains smaller than 10 µm remain highly coupled to the gas at all the tested levels of turbulence (Fig. 10). The 100 µm grains are decoupled from the gas, but the terminal velocity approximation still stands. For larger grains, the velocity drifts are such that a multifluid approach is required. This can be illustrated by modeling the evolution of the dust ratio $\epsilon_d = \rho_d/(\rho_g + \rho_d)$ as a function of the grain size (Fig. 10). We can see the similarity between the three panels, for each grain size. Increasing the grain size by one order of magnitude increases

<table><caption>**Table 1.** Dust-to-gas ratio within the first hydrostatic core as a function of grain size and initial turbulent Mach number.</caption><thead><tr><th></th><th>s<sub>grain</sub> = 10 µm</th><th>s<sub>grain</sub> = 100 µm</th><th>s<sub>grain</sub> = 1 mm</th></tr></thead><tbody><tr><td>M<sub>i</sub> = 0</td><td>0.01041</td><td>0.0128</td><td>0.051</td></tr><tr><td>M<sub>i</sub> = 0.5</td><td>0.01040</td><td>0.0128</td><td>0.07</td></tr><tr><td>M<sub>i</sub> = 1</td><td>0.01056</td><td>0.0146</td><td>0.16</td></tr></tbody></table>

**Notes.** Grain sizes $s_{\text{grain}}$ are indicated in columns and initial turbulent Mach numbers $M_i$ are indicated in lines. Dust-to-gas ratios $\theta_d$ are estimated at the formation of the FHSC, that is the first time the gas density reaches $\rho_g = 10^{-11} \text{ g.cm}^{-3}$. For each run, the initial dust-to-gas ratio is $\theta_{d,0} = 0.01$.

$\log(\epsilon_d/\epsilon_{d,0})$ by about one order of magnitude. This is compatible with the model from Lebreuilly et al. (2020), which predicts that $\epsilon_d/\epsilon_{d,0} = \exp(s_{\text{grain}}/s_{\text{ref}})$, where $s_{\text{ref}}$ is a function of the density profile. Thus, $s_{\text{ref}}$ depends on the evolution stage and the initial Mach number. The dust enrichment trend is a bit more shallow than this exponential trend. However, rough estimates from 10 µm grains and 100 µm grains can be computed at the formation of the first hydrostatic core (FHSC) from values in Table 1. We found $s_{\text{ref}} = 250-410 \text{ µm}$ at $M_i = 0$ and 0.5 and $s_{\text{ref}} = 180-270 \text{ µm}$ at $M_i = 1$. These sizes mean that initial turbulence promotes dust enrichment in the FHSC as if the effective size of the grains is lowered. The dust enrichment model is no longer valid for millimeter grains, most likely because most of the assumptions required to obtain this relationship, such as strong coupling and weak back-reaction, cannot be satisfied. However, the dust enrichment remains an increasing function of the size and initial turbulence. The initial turbulence delays the formation of the FHSC, thereby giving more time to the dust to fall and to settle. This is underlined in Fig. 10, where the dust-to-gas ratio reaches a maximum before the gas starts to form a dense hydrostatic core. Then the dust-to-gas ratio decreases as less dust-enriched material starts to fall and to mix in the inner region.

We extended our analysis to lower density regions. We computed the probability density function (PDF) of the dust-to-gas
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Evolution of the dust enrichment during the collapse. Logarithm of the normalized dust ratio at the maximum gas density, depending on the size of the grains (individual panels) and the initial turbulent Mach (green for $M_i = 0$, blue for $M_i = 0.5$, red for $M_i = 1$). Diamonds on the dotted line indicate the interpolated value to get the dust enrichment within the first hydrostatic core (Table 1).</figcaption></figure>

ratio in different regions of the collapse, within a radius of 2500 AU. We select cells belonging to the first hydrostatic core (FHSC), the disk, the pseudo-disk, the outflows, and the envelope following the criteria used in Lebreuilly et al. (2020), with the references therein. We caution that the selection criterion for the cells belonging to the FHSC is $\rho_g > 10^{-12.5}$ g.cm$^{-3}$. The PDF of the dust-to-gas ratio for 10 µm grains and 100 µm grains is narrow and peaks around the enrichment within the FHSC (details in Appendix H), as expected; whereas the PDF for 1 mm grains is particularly extended (Fig. 11). Therefore, we dedicate the rest of the section to the description of the dust enrichment and depletion of 1 mm grains. We present in Fig. 11 the gas density profile (the feedback from dust grains may not be negligible) and the corresponding dust-to-gas ratio map.

Without any initial turbulence, the main mechanism for dust enrichment in low density regions seems to be the settling due to the pressure gradient created by the stable dense regions of the pseudo-disk. On the contrary, initial turbulence form dust-rich structures in the envelope, where grains are highly decoupled. These structures seem to be decorrelated from the gas density profile in the envelope at the end of the first protostellar collapse phase. Indeed, even though the resulting pseudo-disk is distorted, the gas profile is smooth (upper panels); this is most likely because thermal (and magnetic) pressure provides support against compression, which is not the case for the dust fluid. The probability density function provides the typical maximum dust-to-gas ratios in the low-density regions (envelope and pseudo-disk). We obtained $\theta_d = 10\%, 25\%, 40\%$ for $M_i = 0, 0.5, 1$, respectively. This very promising enrichment of large grains in some locations of the envelope and the pseudo-disk also means that most of these regions becomes highly depleted in dust.

With such enrichment in the envelope and the hydrostatic core, the feedback of the dust on the gas, via gravity and drag, cannot be neglected and it strongly affects the formation of the disk. Indeed, it modifies the balance between the rotational and the gravitational energy, from the gas and the dust, and the thermal energy, only from the gas. These simulations are possibly extreme, because of the initial dust distribution and resulting decoupling and feedback. However, it demonstrates the versatility and robustness of the current multifluid implementation regarding the vast range of dynamical scales.

In dense regions (the first hydrostatic core and the disk) or in very dust-enriched regions (potentially the envelope), the grain coagulation time becomes shorter than the dynamical time scales and, thus, it becomes a dominant mechanism for the evolution of the dust distribution. In protoplanetary disks, the dust continuum emission indicates that grains reach submillimeter sizes (Kataoka et al. 2015). We should therefore account for dust growth by using, for instance, the method of Lombart & Laibe (2021) based on the Eulerian multifluid approach. This requires accurate hydrodynamical velocity drifts between the dust species. In that respect, the common Riemann solver should also help in preventing unphysical velocity drifts when advecting the velocity components. Coupling the hydrodynamics and the growth of the dust distribution should offer a better hint at the conditions of early planet formation.

## 5. Conclusion

We implemented a multifluid method in the RAMSES code for the multiscale physics of dusty flows. We tested our implementation, including up to ten dust fluids in various coupling regimes (Sect. 3.3 and Appendix D). We emphasize the difficulty in capturing the coupling regimes between the dust and the gas, particularly during a protostellar collapse. We present a novel Riemann solver based on the HLL approach to deal with the coupling regimes consistently.

- Current multifluid methods use individual Riemann solvers for each fluid with corresponding truncation errors. In the strong coupling regime, this leads to unbalanced density advection steps and, thus, unphysical dust-to-gas ratio variations.

- When using individual solvers, we found out that the dust enrichment within the first hydrostatic core cannot be studied properly. Indeed, with this strategy, the enrichment of MRN grains, which are well-coupled to the gas, is dominated by spatial errors (Sect. 4.1), where as for millimeter grains, they cannot properly recouple to the gas within the first hydrostatic core (Appendix G).

- Instead, we use a common Riemann solver for the gas and the dust multifluid, which is based on the HLL wave fan of the gas. In the strong coupling regime, this solution bypasses known problems with Eulerian methods for
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Collapse simulations with 1 mm grains and three initial levels of turbulence. They are presented in columns, indicated by the initial turbulent Mach number $M_i$ and the snapshot time, corresponding to the formation of the first hydrostatic core ($\rho_{g,max} = 10^{-11}$ g.cm$^{-3}$). In lines, the gas density map (slice), the dust-to-gas ratio map (slice), and the dust-to-gas ratio probability density function (PDF) in different regions of the collapse (FHSC, disk, outflow, pseudo-disk, envelope, if detected). The histogram is weighted by the gas mass within the cell and we use 100 log-spaced bins. This means, for example for $M_i = 1$, that 70% of the gas mass contains a dust enrichment lower than the initial one ($\theta_d < 10^{-2}$), and about 7% of the gas mass is very dust-enriched with $\theta_d > 10^{-1}$.</figcaption></figure>

pressureless fluids, such as situations of converging flows (Appendix C.2).

- This novel solver eliminates the spatial truncation error observed for a tightly coupled gas and dust mixture. It can reproduce the results of the terminal velocity approximation with reasonable time resolution. This solver allows us to go beyond the terminal velocity approximation and to study situations where accounting for the inertia of dust grains self-consistently is necessary; for instance, in shocks (Appendix E).

- Here, we use a solver switch to deal with weak coupling regimes. The switch criterion is purely based on the fluid kinematics or, more precisely, whether the dust velocities belong to the influence area of the gas, modeled by the wave fan (Sect. 2.3). This strategy limits numerical diffusion. We have not identified any issue with this switch at this time in the context of protostellar collapses.

- The common Riemann solver could provide unnecessarily more numerical diffusion than individual Riemann solvers, for instance, in the case of weakly coupled grains whose
velocity drift remains small. Individual Riemann solvers could be a preferential solution in situations of weak coupling (without changes of regimes), at high spatial resolution compared to the coupling scales, or when the dust-to-gas ratio is not a quantity of interest (compared to the precision on the gas dynamics with HLLD for instance). Even though these conditions are not met in protostellar collapses, this needs to be investigated in future work. In particular, the dust enrichment we found in some regions of the collapse could favor the development of dynamical instabilities, such as the streaming instability (Youdin & Goodman 2005) and the resonant drag instabilities (Squire & Hopkins 2018). We highlight the fact that small-scale dust enrichment, turbulence, and dynamical instabilities are not resolved in the simulations we present in this paper. Investigating small-scale dynamics using the conditions set by the multi-scale collapse simulations is a direct perspective of the code. Moreover, it would allow us to test the limits of the solvers in another context.

In this paper, we emphasize the importance of accounting for the physics of a dust and gas mixture to model the hydrodynamical solver, even when using a fractional step method as we do. This is of key importance to limit errors on the dust-to-gas ratio, a central parameter for triggering dynamical instabilities such as the streaming instability and the resonant drag instabilities, which are potential paths to planetesimal formation. While current works are mainly focused on developing modular and sophisticated high-order-in-time schemes, these multifluid implementations could suffer from truncation errors on the dust-to-gas ratio dominated by the spatial resolution.

Our multifluid method is well-suited to simulations of turbulent protostellar collapses, which are multi-scale and multi-physics. We found that the dust enrichment within the first hydrostatic core is an increasing function of the grain size and the initial turbulence. Grains under 100 µm remain well-coupled to the gas while 100 µm grains are enriched within the first hydrostatic core (between 20 and 50% of dust-to-gas ratio variations). This is in agreement with the terminal velocity approximation. This underlines the role of the force balance between the dust and the gas. Millimeter grains significantly drift relative to the gas within the envelope. Therefore, modeling their inertia thanks to the multifluid approach is necessary. In the presence of turbulence in the initial core, dust can spread in the envelope and form enriched filaments. Millimeter grains also fall faster than the gas and enrich the inner region very early on, that is, prior to the formation of the first hydrostatic core and later fed by less enriched material. This dust mass back-reacts on the gas and affects disk formation as a result. This is a promising avenue for testing early planet formation scenarios.

**Acknowledgements.** We thank the referee for their comments, which contribute in strengthening the quality of the paper. We thank Leodasce Sewanou, Benoît Commerçon, Guillaume Laibe, Geoffroy Lesur and Maxime Lombart for discussions on numerical methods for dust evolution. We also thank Adnan Ali Ahmad for his advice and help in post-processing the simulations. This research has received funding from the European Research Council synergy grant ECOGAL (Grant: 855130).

## References

Ahmad, A., González, M., Hennebelle, P., & Commerçon, B. 2023, A&A, 680, A23

André, P., Belloche, A., Motte, F., & Peretto, N. 2007, A&A, 472, 519

Barranco, J. A., & Goodman, A. A. 1998, ApJ, 504, 207

Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11

Benítez-Llambay, P., Krapp, L., & Pessah, M. E. 2019, ApJS, 241, 25

Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289

Cacciapuoti, L., Macias, E., Maury, A. J., et al. 2023, A&A, 676, A4

Commerçon, B., Teyssier, R., Audit, E., Hennebelle, P., & Chabrier, G. 2011, A&A, 529, A35

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

Courant, R., Friedrichs, K., & Lewy, H. 1928, Mathematische Annalen, 100, 32

Dartois, E., Noble, J. A., Caselli, P., et al. 2024, Nat. Astron., 8, 359

Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 77

Drążkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, ASP Conf. Ser., 534, 717

Epstein, P. S. 1924, Phys. Rev., 23, 710

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632, A5

Gaustad, J. E. 1963, ApJ, 138, 1050

González, M., Vaytet, N., Commerçon, B., & Masson, J. 2015, A&A, 578, A12

Gould, R. J., & Salpeter, E. E. 1963, ApJ, 138, 393

Harten, A., D. Lax, P., & van Leer, B. 1983, SIAM Rev., 25, 35

Huang, P., & Bai, X.-N. 2022, ApJS, 262, 11

Hutchison, M., Price, D. J., & Laibe, G. 2018, MNRAS, 476, 2186

Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78

Keppens, R., Popescu Braileanu, B., Zhou, Y., et al. 2023, A&A, 673, A66

Krapp, L., & Benítez-Llambay, P. 2020, Res. Notes Am. Astron. Soc., 4, 198

Krapp, L., Benítez-Llambay, P., Gressel, O., & Pessah, M. E. 2019, ApJ, 878, L30

Krapp, L., Garrido-Deutelmoser, J., Benítez-Llambay, P., & Kratter, K. M. 2024, ApJS, 271, 7

Kwok, S. 1975, ApJ, 198, 583

Laibe, G., & Price, D. J. 2011, MNRAS, 418, 1491

Laibe, G., & Price, D. J. 2012, MNRAS, 420, 2345

Laibe, G., & Price, D. J. 2014, MNRAS, 440, 2136

Larson, R. B. 1969, MNRAS, 145, 271

Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96

Lebreuilly, U., Commerçon, B., & Laibe, G. 2020, A&A, 641, A112

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lesur, G. R. J., Baghdadadi, S., Wafflard-Fernandez, G., et al. 2023, A&A, 677, A9

LeVeque, R. J. 2004, J. Hyperbolic Differential Eq., 01, 315

LeVeque, R. J., & Bale, D. S. 1999, in Hyperbolic Problems: Theory, Numerics, Applications, eds. R. Jeltsch, & M. Fey (Basel: Birkhäuser Basel), 609

Lombart, M., & Laibe, G. 2021, MNRAS, 501, 4298

Lovascio, F., & Paardekooper, S.-J. 2019, MNRAS, 488, 5290

Manara, C. F., Morbidelli, A., & Guillot, T. 2018, A&A, 618, L3

Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18

Marchand, P., on, B., & Chabrier, G. 2018, A&A

Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24

Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425

Maury, A., Hennebelle, P., & Girart, J. M. 2022, Front. Astron. Space Sci., 9, 949223

Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315

Moseley, E. R., Teyssier, R., & Draine, B. T. 2023, MNRAS, 518, 2825

Mouschovias, T. C., & Spitzer, L. J. 1976, ApJ, 210, 326

Muley, D., Melon Fuksman, J. D., & Klahr, H. 2023, A&A, 678, A162

Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413

Pagani, L., Steinacker, J., Bacmann, A., Stutz, A., & Henning, T. 2010, Science, 329, 1622

Price, D. J., Wurster, J., Nixon, C., et al. 2017, Astrophysics Source Code Library [record ascl:1709.002]

Roe, P. L. 1986, Ann. Rev. Fluid Mech., 18, 337

Rosdahl, J., & Teyssier, R. 2015, MNRAS, 449, 4380

Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011

Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4

Teyssier, R. 2002, A&A, 385, 337

Tsukamoto, Y., Maury, A., Commerçon, B., et al. 2023, ASP Conf. Ser., 534, 317

Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897

van Leer, B. 1974, J. Comput. Phys., 14, 361

Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90

Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

Ysard, N., Jones, A. P., Guillet, V., et al. 2024, A&A, 684, A34

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050
<figure><img src="image_13.png" /><figcaption>Fig. B.1: Numerical solution of the dust velocity within the box domain, where points belonging the level 8 of mesh refinement are in blue, and points belonging to the level 9 are in green. The final time $t = 0.103$, corresponds to $n = 33$ steps at the level 8 of mesh refinement.</figcaption></figure>

## Appendix A: Local Lax-Friedrichs solver for dust fluids

We define hereafter the local Lax-Friedrichs solver for individual dust fluids, denoted as LLFd. The local Lax-Friedrichs solver (Rusanov flux) was already presented and implemented in the context of pressureless fluids for example in Krapp et al. (2024). We consider for the Riemann problem, the left and right states of the dust fluid $U_{d,L}$ and $U_{d,R}$. The normal components of the dust velocity are $V_{d,L}$ and $V_{d,R}$. and the fluxes are denoted by $F_{d,L} = F_d(U_{d,L})$ and $F_{d,R} = F_d(U_{d,R})$. The fluxes for each dust species $d \in \{1, N\}$ are formally defined as

$$F_{LLF,d} := \frac{1}{2}(F_{d,L} + F_{d,R}) - \frac{S_{LLF,d}}{2}(U_{d,R} - U_{d,L}). \quad (A.1)$$

The term involving $S_{LLF,d}$ corresponds to the diffusion part of the solver. $S_{LLF,d}$ should be an estimate of the maximum wave speed in the fluid (see discussion in Sect. 2.3). The answer is not straight-forward since the homogeneous system is not hyperbolic. Considering an isothermal pressure in the dust fluid with corresponding sound speed $c_d$ leads to a hyperbolic system characterized by its eigenvalues $\{V_d - c_d, V_d, +c_d\}$. Therefore, considering the limit $c_d \to 0$, we should use the normal velocity $V_d$ as the typical wave speed. We define the speed associated to the local Lax-Friedrichs solver for each dust fluid as

$$S_{LLF,d} = \max(|V_{d,L}|, |V_{d,R}|). \quad (A.2)$$

## Appendix B: Dustybox: Time sequence of the drag scheme

In this section, we derive the sequence produced by two drag schemes at each time step. The exact drag scheme, following Eq. (23), provides, over one time step $\Delta t = t_{n+1} - t_n$,

$$U(x, t_{n+1}) = \exp(M(x)\Delta t)U(x, t_n). \quad (B.1)$$

Obviously, the sequence leads to the exact solution whatever the choice of the time step

$$U(x, n\Delta t) = (\exp(M(x)\Delta t))^n U(x, t_n) = \exp(M(x)n\Delta t)U(x, 0). \quad (B.2)$$

On the other hand, the first-order implicit Euler scheme is defined over $\Delta t$ as

$$U_1(x, t_{n+1}) = U_1(x, t_n) + M(x)U_1(x, t_{n+1})\Delta t, \quad (B.3)$$

and thus the scheme is formally

$$U_1(x, t_{n+1}) = (\mathbf{1} - M(x)\Delta t)^{-1}U_1(x, t_n). \quad (B.4)$$

Consequently, it produces the sequence

$$U_1(x, n\Delta t) = ((\mathbf{1} - M(x)\Delta t)^{-1})^n U(x, 0). \quad (B.5)$$

In mesh-refined levels, we can decide to subcycle in time. When doing so, a time step $\Delta t$ is split into 2 steps of $\Delta t/2$. Thus, for these points, the first-order implicit Euler scheme over one time step $\Delta t$ produces

$$U_2(x, t_{n+1}) = (\mathbf{1} - M(x)\Delta t/2)^{-1} \times (\mathbf{1} - M(x)\Delta t/2)^{-1} U_2(x, t_n), \quad (B.6)$$

and thus the sequence

$$U_2(x, n\Delta t) = ((\mathbf{1} - M(x)\Delta t/2)^{-1})^{2n} U(x, 0). \quad (B.7)$$

We can compute these matrices and sequences, and compare to the result of the implemented drag scheme. The setup is the one presented in Sect. 3.1 and here we explicitly present the solution within the box domain with two levels of mesh refinement (level 8 in green and level 9 in blue in Fig. B.1). The final time $t = 0.103$ corresponds to $n = 33$ steps at the level 8 of mesh refinement and $2n = 66$ steps at the level 9, thanks to time subcycling. This is the CFL condition at the level 8 of mesh refinement if the hydrodynamics was activated, with $C_{CFL} = 0.8$ and $C_s = 1$ (see the dustywave test, Sect. 3.2).

In Fig. B.1, we can see the impact of the small local density perturbation on the velocity profile, because the final time is longer than the local damping time, that is when $K_d(1/\rho_g(x) + 1/\rho_d(x))t > 1$. We recover the result of the first-order implicit Euler scheme with $((\mathbf{1} - M(x)\Delta t)^{-1})^n U(x, 0)$ for the coarse level and $((\mathbf{1} - M(x)\Delta t/2)^{-1})^{2n} U(x, 0)$ for the subcycled level. Consequently, a discontinuity appears between subcycled levels, unless the exact solver (exponential operator) is used (Eq. (B.2)).

This test is stronger than a time convergence test since the expected truncation error is also tested.

## Appendix C: Dustywave appendices

### C.1. Spatial error when coupling a finite volume method with an ODE solver for the drag step

A finite volume method lies on the value of conserved quantities, averaged within a cell. Thus, we illustrate by volume-averaging the drag step in Eq. (2), and by using the form of the drag force in Eq. (20).

$$\langle \rho_d v_d \rangle_I(t_2) = \langle \rho_d v_d \rangle_I(t_1) + \int_{t_1}^{t_2} \langle K_d(v_g - v_d) \rangle_I(t) dt \quad (C.1)$$
If $K_d$ is constant, the source term that needs to be time-integrated is

$$ \langle F_{g \to d} \rangle_I = K_d \langle (v_g - v_d) \rangle_I = K_d \left( \frac{\langle \rho_g v_g \rangle_I}{\langle \rho_g \rangle_I} - \frac{\langle \rho_d v_d \rangle_I}{\langle \rho_d \rangle_I} \right) + K_d \left( \frac{1}{\rho_d} \frac{\partial \rho_d}{\partial x} \frac{\partial v_d}{\partial x} (0) - \frac{1}{\rho_g} \frac{\partial \rho_g}{\partial x} \frac{\partial v_g}{\partial x} (0) \right) \frac{\Delta x^2}{12} + o(\Delta x^2). \quad (C.2) $$

The latter expression shows that estimating the velocity fields by diving the momentum by the density contributes to a second-order spatial error on the drag term. Another common source of spatial error is made when the volume-averaged value is considered to be the value at the center of the cell (initialization of the fields, comparison with a reference solution as done in Eq. (25)), which contributes to a second-order spatial error according to the midpoint rule.

### C.2. Mesh refinement and time subcycling

In Fig. C.1, we present three snapshots of the dustywave test with two subcycled levels of refinement. A glitch on the dust density appears in the cells in the neighborhood of the level change. As shown by the different snapshots, it oscillates on a short time scale as the dustywave propagates. By adding an ad hoc pressure to the dust fluid and by vanishing the pressure in the gas, we found that this behavior is not a bug, but a limit of the scheme to deal with pressureless fluids. We emphasize that, in this linear test, the dust density has no significant feedback on the system. Thus, the dust density is advected according to the velocity field. Velocity changes can generate density peaks but pressure support could help in limiting the formation of these peaks. In the case of the pressureless dust, these small scales density structures can persist and propagate at its group velocity in the domain, which could make them difficult to track. Even when increasing the spatial resolution, this glitch is present, similarly to the case of Gibbs phenomena. When such errors pass the level boundary, an interpolation error is made. This glitch is easily seen when activating the time subcycling, but it cannot be solved by the exact drag solver (Fig. C.2, in orange). This underlines the fact that velocity discontinuities at the interface, as presented in the dustybox test, are not the main mechanism that generates this glitch. Indeed, the same discontinuity is present in the gas velocity.

We observe a linear decrease of the glitch with the time step. Moreover, using the common Riemann solver HLLgd, whose wave fan is here mainly given by the gas sound speed, softens the glitch but at the cost of numerical diffusion (Fig. C.2, in red). However, thanks to this numerical viscosity, improving the precision of the drag solver helps in reaching a satisfying solution (Fig. C.2, in blue).

We face a similar situation in the disk settling test (Sect. 3.3). A shock feature appears at the mid-plane of the disk when using an individual Riemann solver for the dust, which is particularly visible at low resolution. This is because the mid-plane separate two regions with the velocity change of sign $v_{d,x}(x < 0) > 0$ and $v_{d,x}(x > 0) < 0$. Dust should not form a shock in this high density region, because it is very well coupled to the static gas (the

<table><caption>Table D.1. Initial conditions of the dustyjeanswave.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>(δρ<sub>g</sub>)<sub>0</sub></td><td>2.7314798795353203 × 10<sup>-5</sup></td></tr><tr><td>φ<sub>g</sub></td><td>0.05058203158724116</td></tr><tr><td>(δρ<sub>1</sub>)<sub>0</sub></td><td>1.7343406399117649 × 10<sup>-6</sup></td></tr><tr><td>(δρ<sub>2</sub>)<sub>0</sub></td><td>1.689443433957737 × 10<sup>-6</sup></td></tr><tr><td>(δρ<sub>3</sub>)<sub>0</sub></td><td>2.2962419848802065 × 10<sup>-6</sup></td></tr><tr><td>(δρ<sub>4</sub>)<sub>0</sub></td><td>1.189737412609774 × 10<sup>-5</sup></td></tr><tr><td>(δρ<sub>5</sub>)<sub>0</sub></td><td>1.3702893874930756 × 10<sup>-5</sup></td></tr><tr><td>φ<sub>1</sub></td><td>-2.6633323209942565</td></tr><tr><td>φ<sub>2</sub></td><td>-2.784967044742389</td></tr><tr><td>φ<sub>3</sub></td><td>2.4001985938943213</td></tr><tr><td>φ<sub>4</sub></td><td>0.7389309844855019</td></tr><tr><td>φ<sub>5</sub></td><td>0.12295937958860421</td></tr><tr><td>(v<sub>g</sub>)<sub>0</sub></td><td>5.616883358274461 × 10<sup>-5</sup></td></tr><tr><td>ψ<sub>g</sub></td><td>-0.0723773480013628</td></tr><tr><td>(v<sub>1</sub>)<sub>0</sub></td><td>7.132828728400922 × 10<sup>-6</sup></td></tr><tr><td>(v<sub>2</sub>)<sub>0</sub></td><td>6.948179834704889 × 10<sup>-6</sup></td></tr><tr><td>(v<sub>3</sub>)<sub>0</sub></td><td>9.443762326844544 × 10<sup>-6</sup></td></tr><tr><td>(v<sub>4</sub>)<sub>0</sub></td><td>4.893037158114604 × 10<sup>-5</sup></td></tr><tr><td>(v<sub>5</sub>)<sub>0</sub></td><td>5.6355938876176875 × 10<sup>-5</sup></td></tr><tr><td>ψ<sub>1</sub></td><td>-2.7862917005828582</td></tr><tr><td>ψ<sub>2</sub></td><td>-2.907926424331034</td></tr><tr><td>ψ<sub>3</sub></td><td>2.2772392143057014</td></tr><tr><td>ψ<sub>4</sub></td><td>0.6159716048968983</td></tr><tr><td>ψ<sub>5</sub></td><td>0.0</td></tr></tbody></table>

**Notes.** These parameter values describe the eigenmode corresponding to the eigenvalue $-1.584688946477482 - 12.82288981648163i$ of matrix (D.1).

terminal velocity approximation stands), but an individual Riemann solver treats it as a problem of converging dust flows at the interface.

## Appendix D: Dustyjeanswave

In this test, several dust fluids and the gas interact via gravity and drag, using the completeness of the equations Eqs. (1)-(5). For one dust species, the evolution of the state $W = (\delta\rho_g, \delta\rho_d, v_g, v_d)^T$ is given by the system $d_t W = M_j W$ with

$$ M_j = \begin{pmatrix} 0 & 0 & -ik\rho_g & 0 \\ 0 & 0 & 0 & -ik\rho_d \\ -ikC_s^2/\rho_g - 4\pi G/(ik) & -4\pi G/(ik) & -K_d/\rho_g & K_d/\rho_g \\ -4\pi G/(ik) & -4\pi G/(ik) & K_d/\rho_d & -K_d/\rho_d \end{pmatrix}. \quad (D.1) $$

Contrary to the dustywave test (Sect. 3.2) for which the density of the dust is mostly advected according to its velocity field, here the dust density can strongly feed back the dynamics of the system. We replicate the linear analysis of the Jeans instability extended to a mixture of gas and several dust fluids. In particular, we simulate one of the damping mode of the system with 5 dust fluids. We chose $(\rho_g)_0 = 1$ for the gas and $(\rho_1)_0 = (\rho_2)_0 = (\rho_3)_0 = (\rho_4)_0 = (\rho_5)_0 = 0.5$ for the dust fluids such that the contribution of dust to the gravity is high. We chose $K_1 = 0.01$, $K_2 = 0.1$, $K_3 = 1$, $K_4 = 10$, and $K_5 = 100$ for the drag coefficients of each dust species such that the test covers a vast
<figure><img src="image_14.png" /><figcaption>Fig. C.1: Dustywave test with two subcycled levels of mesh refinement. Dust density at t = 0.7, 1.5, 4.5. The Riemann solvers are HLL for the gas and LLF for the dust (HLLg-LLFd).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. C.2: Zoom-in on the dust density at t = 0.7 for different solvers, represented by different colors (with a darker shade for the level 9 of mesh refinement). In the legend, BK19 refers to the drag solver from Benítez-Llambay et al. (2019) and pred means that we include the drag in the predictor step of the hydro solver.</figcaption></figure>

range of Stokes numbers simultaneously. We chose $C_s = 3$ to get the propagation of a compressible wave. We set $k = 2\pi$ and $G = 1$.

We selected the propagating and damped eigenmode of $M_J$ (extended to 5 dust species) described in Table D.1. We successfully recovered the dynamics of the multifluid mixture, as presented in Fig. D.1.

## Appendix E: Shock in a dust and gas mixture

We followed the shock setup and the solution from Benítez-Llambay et al. (2019) for one dust fluid. We set $K_d = 1$ for the drag coefficient, $C_s = 1$ for the sound speed, and $t = 500$ for the final state. We demonstrate in Fig. E.1 the ability of HLLg to capture the shock and the recoupling region.

## Appendix F: Numerical setup of protostellar collapses

The collapse setup follows the test of Boss & Bodenheimer (1979). The mass of the core is $M_0 = 1M_\odot$. We recall, in the context of the multifluid model, the ratio between the thermal

<figure><img src="image_16.png" /><figcaption>Fig. D.1: Densities at t = 0.2 (which corresponds to ~ 0.5 crossing time of the wave) of the 5 dust species and of the gas (divided by 2) along with the eigenmode described in Table D.1. Dots are the numerical solutions. Continuous lines are the linear solutions. The time step is divided by a factor of 10 ($C_{CFL} = 0.1$) to achieve the convergence in time for this resolution (level 6, 64 points).</figcaption></figure>

energy in the gas and the gravitational energy of the total mass

$$ \alpha = \frac{5}{2} \frac{1}{1 + \sum_d \theta_d} \frac{R_0 k_B T_g}{G M_0 \mu_g m_p}, \qquad (F.1) $$

where $\theta_d$ is the dust-to-gas ratio, $R_0$ is the radius of the initial core, $T_g$ is the temperature of the gas and $\mu_g$ is the mean molecular weight. We define the ratio between the rotational energy and the gravitational energy as

$$ \beta = \frac{1}{3} \frac{R_0^3 \Omega_0^2}{G M_0}, \qquad (F.2) $$

where $\Omega_0$ is the angular velocity of the core. We define the mass-to-flux-to-critical-mass-to-flux ratio and, thus, the strength of the initial magnetic field, as

$$ \mu = \frac{M_0 / \Phi_B}{(M_0 / \Phi_B)_c}, \qquad (F.3) $$

where $\Phi_B$ is the magnetic flux, and $(M_0 / \Phi_B)_c = \frac{0.53}{3\pi} \sqrt{5/G}$ (Mouschovias & Spitzer 1976). We set the initial conditions to
<figure><img src="image_17.png" /><figcaption>Fig. E.1: Dustyshock test with the HLLgd Riemann solver, using 1024 points (level 10 of mesh refinement). Numerical solution in dots and analytical solution in black dashed lines.</figcaption></figure>

$\alpha = 0.4, \beta = 0.04, \mu = 0.3, \delta\rho/\rho = 0.1$ for the initial azimuthal density perturbation, and $\phi_{\text{mag}} = 30^{\circ}$ for the angle between the magnetic field (z-axis in figures) and the rotation axis. The gas is coupled to the magnetic field (non-ideal MHD with ambipolar diffusion only). The reference table of ambipolar resistivities is from Marchand et al. (2016). The equation of state $P_g(\rho_g)$

follows the barotropic law of Vaytet et al. (2013). We consider the Epstein regime for the drag force (Epstein 1924), for which the stopping time is recalled by Eq. 7. We use $\gamma = 5/3$ for the adiabatic index of the gas. The size of the grains $s_{\text{grain,d}}$ are typically between few nanometers to millimeters. We set the density of the grain to $\rho_{\text{grain,d}} = 1 \text{ g.cm}^{-3}$.

We carry out an adaptive refinement with 15 cells per local Jeans length for initially nonturbulent collapses if not specified in Sects. 4.1 and 4.2, and with 40 cells per local Jeans length for turbulent collapses in Sect. 4.3. The minimum (and initial) mesh refinement level is $l_{\text{min}} = 7$ and the maximum level is $l_{\text{max}} = 14$. They correspond respectively to a cell size of 124.38 au and 0.972 au.

## Appendix G: Performance of the Riemann solvers for millimeter grains

In this section, we compare the performances of the HLLgd Riemann solver and the Riemann solver from Huang & Bai (2022) denoted by H&B22. Contrary to the rest of the paper, in this section we deactivate the dust feedback on the gas. Indeed, with such dust enrichment and velocity drift, the feedback via the gravity and the drag force is strong. We also use the same Riemann solver for the gas, that is the HLL solver, so that the gas dynamics is identical in the two simulations. It makes the comparison easier.

In Fig. G.1, dust-to-gas ratio maps are different between the two solvers. Dust settles more in the pseudo-disk for the H&B22 solver, and more generally, it produces more contrasted dust-to-gas ratio profiles, certainly because it is less diffusive than HLLgd. The main differences take place in regions where the recoupling length of the dust is not resolved (in white and in green in the last column of Fig. G.1, that is in the first hydrostatic core, the disk, and the pseudo-disk). This issue worsens as the gas density increases. Indeed, we recall that $\Delta x \propto 1/\sqrt{\rho_g}$ (refinement based on the Jeans length) and $c_s t_{s,d} \propto s_{\text{grain}}/\rho_g$, thus $\Delta x/(c_s t_{s,d}) \propto \sqrt{\rho_g}/s_{\text{grain}}$. Moreover, in high density regions, the velocity drift vanishes and thus the dust velocity enters in the gas wave fan, and thus Riemann solvers differ.

More importantly, the dust recouples within the first hydrostatic core and thus the dust-to-gas ratio is homogenized. This is correctly reproduces by the HLLgd solver but not by the H&B22 solver (Fig. G.2). Here, the H&B22 solver lacks of robustness and can even produce negative dust densities with the chosen numerical parameters. Resolving the recoupling length within the first hydrostatic core, even for millimeter grains, requires to increase the space resolution by a factor of 100. This corresponds to a fraction of 0.01 au, or even a fraction of the solar radius for smaller grains. Time integration of the disk formation could not be achieved with such stringent refinement.

## Appendix H: Spatial enrichment of super-micron dust grains as a function of the initial level of turbulence

We complete the analysis of the dust enrichment in low density regions in Sect. 4.3 for 10 µm grains (Fig. H.1) and 100 µm grains (Fig. H.2). Contrary to 1 mm grains, in the envelope and the pseudo-disk, the PDF of the dust-to-gas ratio of smaller grains is peaked around the enrichment within the FHSC. The dust enrichment better corresponds to the gas density profile. The main rotation plane is mostly enriched (Figs. H.1 and H.2, face-on view). Some regions outside the rotation plane are depleted,
<figure><img src="image_18.png" /><figcaption>Fig. G.1: Comparison between two Riemann solvers for the dust fluid: HLLgd and H&B22, the solver from Huang & Bai (2022). Gas density map in the first column (HLL solver) at $t = 60.6$ kyr, dust-to-gas ratio maps from the two Riemann solvers (HLLgd in the second column and H&B22 in the third column), and resolution of the recoupling length expressed as $\Delta x/(c_s t_{s,1mm})$, with mesh-refinement levels indicated by contours, in the last column. Zoom-in of the collapse region (upper panels) to the disk scale (lower panels). Dust feedback has been deactivated to ease the comparison between the two solvers.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. G.2: Distribution of the dust-to-gas ratio as a function of the gas density (above $10^{-13}$ g.cm$^{-3}$), at $t = 60.8$ kyr. The colorbar indicates the number of cells in the bin of the histogram.</figcaption></figure>

leading to the tail of the dust-to-gas ratio PDF. The dust depletion in the envelope is more important as the initial turbulent Mach $M_i$ increases. One possible explanation is that turbulence delays the FHSC formation time, providing more time for dust grains in the envelope to settle in the pseudo-disk.
<figure><img src="image_20.png" /><figcaption>Fig. H.1: Same figure as Fig. 11, but for 10 µm grains.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

Fig. H.2: Same figure as Fig. 11, but for 100 µm grains. For $M_i = 1$, the PDF continues to decrease for lower dust-to-gas ratio values until $\theta_d \approx 4 \times 10^{-8}$ where the PDF is about $2 \times 10^{-4}$.