<figure><img src="image_1.png" /></figure>

HAL
open science

# Protostellar collapse: the conditions to form dust-rich protoplanetary disks

U. Lebreuilly, B. Commerçon, G. Laibe

► To cite this version:

U. Lebreuilly, B. Commerçon, G. Laibe. Protostellar collapse: the conditions to form dust-rich protoplanetary disks. Astronomy and Astrophysics - A&A, 2020, 641, pp.A112. 10.1051/0004-6361/202038174. hal-02942468

HAL Id: hal-02942468

https://hal.science/hal-02942468

Submitted on 17 Sep 2020

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Protostellar collapse: the conditions to form dust-rich protoplanetary disks

U. Lebreuilly, B. Commerçon, and G. Laibe

École normale supérieure de Lyon, CRAL, UMR CNRS 5574, Université de Lyon, 46 allée d Italie, 69364 Lyon Cedex 07, France
e-mail: ugo.lebreuilly@ens-lyon.fr

Received 15 April 2020 / Accepted 12 July 2020

## ABSTRACT

**Context.** Dust plays a key role during star, disk, and planet formation. Yet, its dynamics during the protostellar collapse remain a poorly investigated field. Recent studies seem to indicate that dust may decouple efficiently from the gas during these early stages.

**Aims.** We aim to understand how much and in which regions dust grains concentrate during the early phases of the protostellar collapse, and to see how this depends on the properties of the initial cloud and of the solid particles.

**Methods.** We used the multiple species dust dynamics MULTIGRAIN solver of the grid-based code RAMSES to perform various simulations of dusty collapses. We performed hydrodynamical and magnetohydrodynamical simulations where we varied the maximum size of the dust distribution, the thermal-to-gravitational energy ratio, and the magnetic properties of the cloud. We simulated the simultaneous evolution of ten neutral dust grain species with grain sizes varying from a few nanometers to a few hundreds of microns.

**Results.** We obtain a significant decoupling between the gas and the dust for grains of typical sizes of a few tens of microns. This decoupling strongly depends on the thermal-to-gravitational energy ratio, the grain sizes, and the inclusion of a magnetic field. With a semi-analytic model calibrated on our results, we show that the dust ratio mostly varies exponentially with the initial Stokes number at a rate that depends on the local cloud properties.

**Conclusions.** We find that larger grains tend to settle and drift efficiently in the first-core and in the newly formed disk. This can produce dust-to-gas ratios of several times the initial value. Dust concentrates in high-density regions (cores, disk, and pseudo-disk) and is depleted in low-density regions (envelope and outflows). The size at which grains decouple from the gas depends on the initial properties of the clouds. Since dust cannot necessarily be used as a proxy for gas during the collapse, we emphasize the necessity of including the treatment of its dynamics in protostellar collapse simulations.

**Key words.** ISM: kinematics and dynamics – hydrodynamics – stars: formation – methods: numerical

## 1. Introduction

Small dust grains are essential ingredients of star, disk, and planet formation. They regulate the thermal budget of star forming regions through their opacity and thermal emission (McKee & Ostriker 2007; Draine 2004). In addition, they are thought to be the main formation site of H₂ at present (Gould & Salpeter 1963). It is widely accepted that planet formation is induced by dust growth within protoplanetary disks (see the recent review by Birnstiel et al. 2016). Finally, the dust grains are significant charge carriers (Marchand et al. 2016; Wurster et al. 2016; Zhao et al. 2016) and therefore regulate the evolution of magnetic fields during the protostellar collapse which can affect, among others, the disk formation (Masson et al. 2016; Hennebelle et al. 2020) and the fragmentation processes (Commerçon et al. 2011a).

Until recently, the accepted paradigm was that the dust of the interstellar medium (ISM) is usually composed of grains with sizes up to ~0.1 µm with a typical size distribution well modeled by the Mathis-Rumpl-Nordsieck (MRN) distribution (Mathis et al. 1977). Recent observations seem to indicate that larger grains could exist in the denser regions of the ISM. Pagani et al. (2010) proposed that over-bright envelopes of prestellar cores (coreshine) could be explained by the presence of micrometer grains. In addition, it was suggested that recent observations with

ALMA of the polarised light at (sub)millimeter wavelengths in Class 0 and I objects could be interpreted as the presence of grains up to ≈100 µm (Kataoka et al. 2015, 2016; Pohl et al. 2016; Sadavoy et al. 2018a,b, 2019; Valdivia et al. 2019). Galametz et al. (2019) also showed that the low values of the dust emissivity in Class 0 objects could indicate the presence of these large grains in their envelope. Finally, Tychoniec et al. (2020) estimated that the mass of solids in Class 0 disks is sufficient to grow planets only if large grains are included in the opacity models, which might indicate dust growth in the early phases of protostar formation.

Over the past few years, significant improvements have been made in numerical models to better understand the early phases of the protostellar collapse that leads to the first Larson core formation (Larson 1969). The angular momentum budget is a long-standing problem in star formation. Indeed, the specific angular momentum of prestellar cores differs from those of young stars by more than three orders of magnitude (Bodenheimer 1995; Belloche 2013). In numerous studies, magnetic braking has been investigated as one of the possible solutions to address this issue (Allen et al. 2003; Price & Bate 2007; Hennebelle & Fromang 2008; Commerçon et al. 2011a; Masson et al. 2016). State-of-the-art simulations account now for the effect of magnetic fields both in an ideal (Commerçon et al. 2010) and nonideal (Tomida et al. 2015; Vaytet et al. 2018; Wurster et al. 2019)
magnetohydrodynamics (MHD) framework, radiative feedback (Commerçon et al. 2010; González et al. 2015; Tomida et al. 2015), and other physical mechanisms. Only Bate & Lorén-Aguilar (2017) have so far investigated the dynamics of dust during the 3D protostellar collapse (DUSTYCOLLAPSE). These latter authors report that ~100 µm grains can significantly decouple from the gas leading to a large increase of the dust-to-gas ratio in the disk and the first Larson core. In 2D, Vorobyov & Elbakyan (2019) studied the gas and dust decoupling in gravitoviscous protoplanetry disks including dust growth, and also report strong variations in dust-to-gas ratio. So far, no 3D DUSTYCOLLAPSE simulation has been performed in a MHD/nonideal MHD context or with multiple dust species.

Multifluid-multi-dust species simulations are prohibitive as they require that two additional equations be solved per dust species (mass and momentum conservation). Recent theoretical (Laibe & Price 2014a,b) and numerical (Hutchison et al. 2018) developments have shown that, for strongly coupled dust and gas mixtures, only the mass conservation equation needs to be solved as the differential velocity between gas and dust can be directly determined by the force budget on both phases. This is the so-called one-fluid or monofluid approach in the diffusion limit (Laibe & Price 2014a). In a recent study (Lebreuilly et al. 2019), we presented a fast, accurate, and robust implementation of dust dynamics for strongly coupled gas and dust mixtures, which allows an efficient treatment of multiple grain species in the adaptive-mesh-refinement (AMR, Berger & Oliger 1984) finite-volume code RAMSES (Teyssier 2002). In this study, we extend this one-fluid formalism to neutral dust grains in a partially ionized plasma. We present the first simulations of protostellar collapse of gas and dust mixtures with multiple grain species (MULTIGRAIN hereafter). In particular, we investigate how the maximum grain size of the dust distribution, the ratio between the thermal and gravitational energy, i.e., the thermal support, and the magnetic fields may affect the decoupling between the gas and neutral dust grains.

This paper is organized as follows. We recall the general framework of gas and dust mixtures dynamics and extend it to neutral grains in the presence of a magnetic field in Sect. 2 and the methods in Sect. 3. Section 4 gives an overview of the different models considered. We summarize the important features of the dusty-collapses obtained in our simulations in Sect. 5. We show that dust grains present very distinct behaviors in the core and the fragments; the disk and the pseudo-disk; the outflow and the envelope. In Sect. 6, we propose a semi-analytical and a simplified analytical model to estimate the central dust enrichment during the collapse. We discuss the implications and caveats of this work in Sect. 7. Finally we present our conclusions and perspectives in Sect. 8.

## 2. Framework

### 2.1. Dusty hydrodynamics for the protostellar collapse

A gas and dust mixture with $N$ small grain species can be modeled as a monofluid in the diffusion approximation (Laibe & Price 2014c; Price & Laibe 2015; Hutchison et al. 2018; Lebreuilly et al. 2019). This fluid, of density $\rho$, flows at its barycenter velocity $v$. The $k$th dust phase, of density $\rho_k$, has the specific velocity $v + w_k$, where $w_k$ is the differential velocity between the dust and the barycenter. In the context of the protostellar collapse and in the absence of magnetic fields, this

mixture is well described by the following set of equations

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot [\rho v] &= 0, \\
\frac{\partial \rho_k}{\partial t} + \nabla \cdot [\rho_k (v + w_k)] &= 0, \quad \forall k \in [1, N], \\
\frac{\partial \rho v}{\partial t} + \nabla \cdot [P_g \mathbb{I} + \rho (v \otimes v)] &= -\rho \nabla \phi,
\end{aligned}
\quad (1) $$

where $P_g$ is the thermal pressure of the gas and $\mathbb{I}$ the identity matrix. The gravitational potential $\phi$ is set by the Poisson equation

$$ \Delta\phi = 4\pi\mathcal{G}\rho, \quad (2) $$

where $\mathcal{G}$ denotes the gravitational constant.

The previous equations are closed using a barotropic law that reproduces both the isothermal regime at low density and the adiabatic regime when the density reaches the critical value $\rho_{ad}$ which corresponds to the density at which dust becomes opaque to its own radiation (Larson 1969). Similarly to Commerçon et al. (2008), we express the gas pressure as

$$ P_g = \rho_g c_{s,iso}^2 \left[ 1 + \left( \frac{\rho_g}{\rho_{ad}} \right)^{\gamma-1} \right]. \quad (3) $$

The gas density $\rho_g$ is $\rho_g = \rho - \sum_k \rho_k$. Regions of low density are isothermal and have a sound speed of $c_{s,iso}$.

As in Lebreuilly et al. (2019), we model a single dust grain $k$ as a small compact and homogeneous sphere of radius $s_{grain,k}$ and intrinsic density $\rho_{grain,k}^1$. When the grain is smaller than the mean free path of the gas (the so-called Epstein drag regime, Epstein 1924), the drag stopping time $t_{s,k}$ is given by

$$ t_{s,k} \equiv \sqrt{\frac{\pi\gamma}{8} \frac{\rho_{grain,k}}{\rho} \frac{s_{grain,k}}{c_s}}, \quad (4) $$

where $\rho$ is the total density of the gas and dust mixture, $c_s$ is the sound speed of the gas, and $\gamma$ its adiabatic index.

If the differential velocity $\Delta v_k \equiv v_k - v_g$ between the gas and the dust is supersonic, a correction in the drag regime must be applied. In this case the stopping time is given by (Kwok 1975)

$$ t_{s,k} \equiv \sqrt{\frac{\pi\gamma}{8} \frac{\rho_{grain,k}}{\rho} \frac{s_{grain,k}}{c_s} \left(1 + \frac{9}{128\pi} M_d^2\right)^{-1/2}}, \quad (5) $$

where $M_d = |\Delta v_k|/c_s$ is the differential velocity Mach number. Hereafter, unless specified, we consider this correction.

In the terminal velocity approximation, the differential velocity of the phase $k$ is

$$ w_k = \left[ \frac{\rho}{\rho - \rho_k} t_{s,k} - \sum_{l=1}^{N} \frac{\rho_l}{\rho - \rho_l} t_{s,l} \right] \frac{\nabla P_g}{\rho}, \quad (6) $$

and the gas and dust velocities, $v_g$ and $v_k$ are given by

$$ 
\begin{aligned}
v_g &= v - \sum_k \frac{\rho_k}{\rho - \rho_k} w_k, \\
v_k &= v + w_k.
\end{aligned}
\quad (7) $$

<footnote n="1">¹ Distinct from the dust density $\rho_k$.</footnote>
For later purposes, we define the dust ratio $\epsilon_k \equiv \frac{\rho_k}{\rho}$, the total dust ratio $\epsilon \equiv \sum_k^N \epsilon_k$, and the dust-to-gas ratio $\theta_d \equiv \frac{\sum_k^N \rho_k}{\rho_g}$. For any quantity A, we define $\bar{A} \equiv A/A_0$ where $A_0$ is its initial value. Here, $\bar{\epsilon}$ and $\bar{\theta}$ are the dust-ratio and dust-to-gas-ratio enrichment, respectively. Further details about the monofluid formalism and its implementation in RAMSES can be found in Lebreuilly et al. (2019).

## 2.2. Dusty MHD with neutral grains

We extend the above formalism to neutral grains embedded in a weakly ionized plasma. Here, we only consider the resistive effect of ambipolar diffusion, that is, the drift between ions and neutrals other than the dust. This formalism can be straightforwardly extended to more general Ohm's laws. In this context, the equations of dusty-magnetohydrodynamics with neutral grains (Ndusty-MHD) can be written as follows.

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot [\rho \mathbf{v}] &= 0, \\
\frac{\partial \rho_k}{\partial t} + \nabla \cdot [\rho_k (\mathbf{v} + \mathbf{w}_k)] &= 0, \quad \forall k \in [1, N], \\
\frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot \left[ \left( P_g + \frac{\mathbf{B}^2}{2} \right) \mathbf{I} + \rho (\mathbf{v} \otimes \mathbf{v}) - \mathbf{B} \otimes \mathbf{B} \right] &= -\rho \nabla \phi, \\
\frac{\partial \mathbf{B}}{\partial t} - \nabla \times \left[ \left( \mathbf{v} - \sum_k \frac{\rho_k}{\rho - \rho_k} \mathbf{w}_k \right) \times \mathbf{B} \right] \\
+ \nabla \times \left[ \frac{\eta_A c^2}{4\pi|\mathbf{B}|^2} [(\nabla \times \mathbf{B}) \times \mathbf{B}] \times \mathbf{B} \right] &= 0, \\
\nabla \cdot \mathbf{B} &= 0,
\end{aligned}
\quad (8) $$

where $\mathbf{B}$ is the magnetic field, $\eta_A$ is the ambipolar resistivity, and $c$ is the speed of light. We note that $\frac{\eta_A c^2}{4\pi|\mathbf{B}|^2}[(\nabla \times \mathbf{B}) \times \mathbf{B}]$ is the differential velocity between the ions and the neutrals in the gas phase (Masson et al. 2012, 2016). We must then correct the velocity of the neutrals by accounting for the differential velocity between the gas and the barycenter. The term $-\sum_k \frac{\rho_k}{\rho-\rho_k} \mathbf{w}_k$ appears in the induction equation to take this into account.

In addition to the pressure force, a magnetic force now applies on the plasma. In a previous study, Laibe & Price (2014c) provided the expression of the differential velocity for a general force budget. Using this formula we find that

$$ \mathbf{w}_k = \left[ \frac{\rho}{\rho - \rho_k} t_{s,k} - \sum_{l=1}^{N} \frac{\rho_l}{\rho - \rho_l} t_{s,l} \right] \frac{\nabla P_g - (\nabla \times \mathbf{B}) \times \mathbf{B}}{\rho}, \quad (9) $$

and we note that this expression is very similar to what Fromang & Papaloizou (2006) found for single dust species mixtures.

To further simplify Eq. (8), we assume in this paper that the plasma velocity $\mathbf{v}$ is the barycenter velocity which is valid when $\epsilon_k ||\mathbf{w}_k|| \ll ||\mathbf{w}_k|| \ll |\mathbf{v}|$. In this case, the induction equation writes

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times [\mathbf{v} \times \mathbf{B}] + \nabla \times \left[ \frac{\eta_A c^2}{4\pi|\mathbf{B}|^2} [(\nabla \times \mathbf{B}) \times \mathbf{B}] \times \mathbf{B} \right] = 0. \quad (10) $$

# 3. Method

## 3.1. RAMSES

For this work, we take advantage of the RAMSES code (Teyssier 2002). This finite-volume Eulerian code solves the hydrodynamics equations using a second-order Godunov method (Godunov

1959) on an adaptive-mesh-refinement grid (Berger & Oliger 1984). With a proper refinement criteria, the AMR grid is a powerful tool to study multi-scale problems such as the protostellar collapse of a dense core. The RAMSES code is very effective for problems that require the treatment of magnetohydrodynamics (Teyssier et al. 2006; Fromang et al. 2006; Masson et al. 2012; Marchand et al. 2018, 2019), radiation hydrodynamics (Commerçon et al. 2011b, 2014; Rosdahl et al. 2013; Rosdahl & Teyssier 2015; González et al. 2015; Mignon-Risse et al. 2020) or cosmic rays (Dubois & Commerçon 2016; Dubois et al. 2019).

We extended the code to the treatment of dust dynamics with multiple species in the diffusion approximation and terminal velocity regime (Lebreuilly et al. 2019). The method, based on an operator splitting technique has been extensively presented and tested. It uses a predictor-corrector MUSCL scheme (van Leer 1974) and is second-order accurate in space. The solver can be used to simultaneously and efficiently model several dust species (MULTIGRAIN).

## 3.2. Boss and Bodenheimer test

We perform Boss and Bodenheimer tests (Boss & Bodenheimer 1979) to follow the dynamics of the dust during the first collapse and first core formation. The parameters of the setup are the initial mass of the dense core (or prestellar core) $M_0$, the total dust ratio $\epsilon_0$, the temperature of the gas $T_g$, and a mean molecular weight $\mu_g$. The ratio between the thermal and the gravitational energy $\alpha$ is

$$ \alpha = \frac{5(1 - \epsilon_0)R_0 k_B T_g}{2 \mathcal{G} M_0 \mu_g m_H}, \quad (11) $$

and sets the initial radius of the cloud $R_0$ and its density $\rho_0$. In addition, we impose an initial solid body rotation around the z-axis at the angular velocity $\Omega_0$ by setting the ratio between the rotational and the gravitational energy $\beta$ given by

$$ \beta = \frac{1}{3} \frac{R_0^3 \Omega_0^2}{\mathcal{G} M_0}. \quad (12) $$

Eventually, we apply an initial azimuthal density perturbation according to

$$ \rho = \rho_0 [1 + A \cos(m\theta)]. \quad (13) $$

In this paper, we aim to investigate the impact of magnetic fields on the dynamics of neutral dust grains in two simulations: one with ideal MHD and one with ambipolar diffusion. For these runs, we impose a uniform magnetic field using the mass-to-flux-to-critical-mass-to-flux-ratio

$$ \mu = \left( \frac{M_0}{\Phi} \right) / \left( \frac{M}{\Phi} \right)_c, \quad (14) $$

the critical mass-to-flux ratio being given by $(M/\Phi)_c = \frac{0.53}{3\pi} \sqrt{5/\mathcal{G}}$ (Mouschovias & Spitzer 1976). We set an angle $\phi_{mag}$ between the magnetic fields and the rotation axis to reduce the efficiency of the magnetic braking.

## 3.3. Dust grain size distributions

In our MULTIGRAIN simulations, $N > 1$ dust bins are considered. The dust ratio of each bin is set from power-law distributions

$$ \frac{d\epsilon}{ds} = \frac{\epsilon_0}{\int_{S_{\min}}^{S_{\max}} s^{3-m} ds} s^{3-m}, \quad (15) $$
with $\epsilon_0$ being the total initial dust ratio, and $S_{min}$ and $S_{max}$ being the minimum and maximum sizes of the grains present in the medium, respectively. For the standard MRN distribution, $S_{min} = 5$ nm, $S_{max} = 250$ nm, and $m = 3.5$.

The method described in Hutchison et al. (2018) is used to compute the initial dust ratio and typical grain size of each bin. A logarithmic grid is used to determine the edges $S_k$ of the bins

$$ \log(S_k) = \log\left(\frac{S_{\max}}{S_{\min}}\right) \frac{k}{N} + \log(S_{\min}). \quad (16) $$

The typical grain size of a bin $k$ required to compute the stopping time is

$$ s_k = \sqrt{S_k S_{k+1}}. \quad (17) $$

We note that $S_{min}$ and $S_{max}$ are the edges of the distribution and must not be confused with the minimum and maximum bin size that are averaged quantities. The initial dust ratios in each bin are computed according to

$$ \epsilon_{0k} = \epsilon_0 \left[ \frac{S_{k+1}^{4-m} - S_k^{4-m}}{S_{\max}^{4-m} - S_{\min}^{4-m}} \right]. \quad (18) $$

## 3.4. Setup

### 3.4.1. Cloud setup

We impose initial conditions that are typical of the first protostellar collapse (Larson 1969) with $T_g = 10$ K, $\mu_g = 2.31$, and a solar mass cloud. We also set $\gamma = 5/3$ since molecular hydrogen behaves as a monoatomic gas at low temperatures (Whitworth & Clarke 1997). As explained in Sect. 2.1, a barotropic law is used to close Eq. (1) with $\rho_{ad} = 10^{-13}$ g cm$^{-3}$ (Larson 1969). Finally, we always set $m = 2$ and $A = 0.1$ to favor fragmentation and the formation of two spiral arms.

For the dust, we always consider ten bins with grain sizes distributed according to Sect. 3.3. In all the models, $S_{min} = 5$ nm and $m = 3.5$, and $S_{max}$ is specified individually. In all our models we set $\rho_{grain} = 1$ g cm$^{-3}$. Finally, we impose a uniform initial total dust-to-gas ratio of $\theta_{d,0} = 0.01$ in all the models.

The two magnetic models have been computed with $\mu = 5$ and $\phi_{mag} = 40^\circ$. For the nonideal MHD model the ambipolar resistivity is computed similarly to the case of reference of Marchand et al. (2016).

### 3.4.2. Numerical setup

We use the HLLD Riemmann solver (Miyoshi & Kusano 2005) for the barycenter part of the conservation equations with a MINMOD slope limiter (Roe 1986) for both the gas and the dust. The Truelove criterion (at least four point per Jeans length, Truelove et al. 1997) must be satisfied to avoid artificial clump formation. We therefore enforce a refinement criterion that imposes at least 15 points per local Jeans length. The grid is initialized to the level $\ell_{min} = 5$ and allows refinement up to a maximum level $\ell_{max} = 16$ (which gives a resolution between 32$^3$ and 65536$^3$ cells).

### 3.4.3. Analysis of the models

We consider that the first hydrostatic core (FHSC) is fully formed when the peak density reaches $10^{-11}$ g cm$^{-3}$ for the first time. We denote the corresponding time $t_{core}$, and use this definition to compare our models at similar evolutionary stages. We present

here the different objects that are observed in our model and their definition in this work.

- The first hydrostatic core and the fragments are any object of density larger than $10^{-12.5}$ g cm$^{-3}$, as in Lebreuilly et al. (2019). The FHSC or $\mathcal{F}_0$ corresponds to the central fragment. $\mathcal{F}_1$ and $\mathcal{F}_2$ are the secondary fragments or the secondary FHSC.

- The disk $\mathcal{D}$ is the region that satisfies Joos et al. (2012) criterion. For the analysis, we place ourselves in cylindrical coordinates $(r, \phi, z)$. A region is identified as a disk if it is Keplerian ($v_\phi > f_{thre}v_r$), in hydrostatic equilibrium ($v_\phi > f_{thre}v_z$), rotationally supported ($\frac{1}{2}\rho v_\phi^2 > f_{thre}P_g$), and dense $\rho > 3.9 \times 10^{-15}$ g cm$^{-3}$. As in Joos et al. (2012), we choose $f_{thre} = 2$.

- The pseudo-disk $\mathcal{P}$ (Galli & Shu 1993), for magnetic runs, is defined as the regions with $r < 5000$ AU and densities above $3.9 \times 10^{-17}$ g cm$^{-3}$ that are not in the disk or fragments. The criterion is similar to that used in Hincelin et al. (2016). Here, 5000 AU is an arbitrary distance that is sufficiently larger than the central objects while being smaller than the initial cloud.

- Jets/outflows $\mathcal{O}$ correspond to any region with $r < 5000$ AU with $v \cdot \frac{r}{|r|} > 0.2$ km s$^{-1}$. This criterion is also similar to the one used in Hincelin et al. (2016).

- The envelope $\mathcal{E}$ encompasses the regions with $r < 5000$ AU that exclude the fragments, the disk, the pseudo-disk, the jets and the outflows.

We consider two different weights for averaging a quantity $A$ over a volume $\mathcal{V}$ in the computational box. Volume averaging is computed according to

$$ \langle A \rangle_v = \frac{\sum_{i \in \mathcal{V}} \Delta x_i^3 A_i}{\sum_{i \in \mathcal{V}} \Delta x_i^3}. \quad (19) $$

Mass averaging is computed according to

$$ \langle A \rangle_m = \frac{\sum_{i \in \mathcal{V}} \rho_i \Delta x_i^3 A_i}{\sum_{i \in \mathcal{V}} \rho_i \Delta x_i^3}, \quad (20) $$

where $\rho_i$ and $\Delta x_i$ are the total density and length of individual cells $i$ in the averaged volume. Volume averages emphasize regions of large spatial extension, such as the envelope, while mass averages emphasize regions of high density, such as the core+disk and the denser regions of the envelope.

## 3.5. Regularization of the differential velocity and dust density

The terminal velocity approximation is unrealistic in low-density regions or in shocked regions where the pressure is discontinuous (Lovascio & Paardekooper 2019). We cap the differential velocities to $w_{cap}$ in our models to avoid prohibitively small time-steps and unrealistically large variations in the dust ratio in strong shock fronts. We impose $w_{cap} = 1$ km s$^{-1}$. To verify that this does not impact the results, we ran extra models with $w_{cap} = 0.1$ km s$^{-1}$, $w_{cap} = 0.5$ km s$^{-1}$, and $w_{cap} = 2$ km s$^{-1}$. A comparison between these models and our fiducial is given in Appendix A.

In our models, the drift velocity can in some regions be supersonic. To account for the correction presented in Eq. (5), we use the drift velocity estimated at the previous time-step to estimate the differential velocity mach number.

Finally, we set the drift velocity to zero at densities lower than those of the initial cloud, i.e., the background. This is a way
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Logarithm of the maximum (10th bin) Stokes number as a function of the spherical radius for the four models with the largest Stokes numbers at $t_{core} + 2$ kyr. Top-left: MMMRN, top-right: MMMRNA0.25, bottom-left: MMMRNmhd, bottom-right: MMMRNnimhd.</figcaption></figure>

to ensure that the regions where the terminal velocity is not valid do not significantly affect the calculation.

### 3.6. Validity of the model

The diffusion approximation is valid as long as the ratio between the stopping time and the dynamical timescale of the gas is small compared to unity (Laibe & Price 2014c). This ratio is called the Stokes number St.

During the first collapse, the dynamical timescale is the free-fall time $t_{ff}$. We have shown in Lebreuilly et al. (2019) that for an initial dust-to-gas ratio of 0.01 and a temperature of 10 K, the initial Stokes number $St_0$ of a spherical collapse is given by

$$St_0 \sim 0.038 \left(\frac{M_0}{1 M_\odot}\right) \left(\frac{\rho_{\text{grain}}}{1 \text{ g cm}^{-3}}\right) \left(\frac{s_{\text{grain}}}{0.05 \text{ cm}}\right) \left(\frac{\alpha}{0.5}\right)^{3/2} < 1. \quad (21)$$

The dynamics of grains smaller than 0.05 cm can be simulated using the diffusion approximation. We note that the Stokes number varies as $\propto \frac{1}{\sqrt{\rho}}$ (since $t_{ff} \propto \frac{1}{\sqrt{\rho}}$ and $t_s \propto \frac{1}{\rho}$) and hence can increase with a decreasing density. Figure 1 shows the values of the maximum Stokes number as a function of the radius for the four models with the least coupled dust (see Sect. 4 for a description of the models). The maximum value for St is smaller than ~0.15 in the external regions of the collapse for all our models and this value is typically smaller than 0.05 inside the collapsing regions of the models. Rotation provides an additional support to the collapse, which causes an increase of the free-fall timescale. Hence, as the initial angular velocity increases, the initial Stokes number decreases and the diffusion approximation is even more accurate. Similarly, magnetic fields increase the duration of the collapse which broadens the validity domain of the diffusion approximation. Finally, the Stokes number also decreases when $\alpha$ decreases, implying that the diffusion approximation remains valid for small values of $\alpha$.

In the induction equation, we consider for simplicity that the plasma is moving at the barycenter velocity. This approximation is valid when $\epsilon \ll 1$, that is, when the back-reaction from the dust onto the gas is negligible. We note that we investigate the impact of the back-reaction in Appendix B.

## 4. Models

The models presented in this section are referenced in Table 1. All of them were evolved up to at least 2 kyr after the formation

of the first core. Our fiducial case MMMRN was run over a longer time.

### 4.1. Fiducial simulation

In this section, we present our fiducial case MMMRN where $\alpha = 0.5$ and $\beta = 0.03$. The grain size distribution is extended up to $S_{max} = 1$ mm. The value of $S_{max}$ leads to an average size of the last and largest bin $s_{10} \sim 160$ µm. The total initial dust-to-gas ratio is $\theta_{d,0} = 1\%$.

Figure 2 shows edge-on (four top figures) and mid-plane (four bottom figures) density cuts of the gas (left) and of the dust fluid with the largest grain size (160 µm, right) at 1 and 2 kyr after the formation of the first Larson core ($t_{core} = 72.84$ kyr for this model), respectively. The density distributions obtained for the gas and the dust are clearly different. This discrepancy originates from an imperfect coupling between the two phases which causes a drift of the dust toward the inner regions of the collapse. This general trend can be explained by a simple force budget on the gas and the dust. Although the gas is partially supported by pressure, dust grains are only subjected to gravity and gas drag. As such, the dust fluid collapses essentially faster than the gas. It therefore enriches the first core and the disk at the cost of a depletion of solids in the envelope. Figure 2 shows that this strong enrichment in the mid-plane and depletion in the envelope have already occurred at $t_{core} + 1$ kyr, and continue for more than 1 kyr. In the mid-plane, the envelope is enriched in large dust grains close to the central object and depleted further away. In the vertical direction, it is mostly depleted in large grains. In short, after the first core formation these grains are concentrated in a very thin layer of 10–100 AU above/under the mid-plane. At this stage, the envelope is mostly a reservoir of low dust densities for the large grains. Hence accretion of dust arising from the envelope does not significantly enrich the fragments or the disk with large grains and the dust-to-gas ratio in the disk even decreases. We note that the disk is still very enriched by the end of the calculation. Most of the enrichment of the dust-to-gas ratio in the central objects indeed takes place during the initial phases of the collapse, when the densities are low everywhere and the coupling between the gas and the dust is the weakest.

Figure 3 shows probability density functions (PDFs) of the dust ratio enrichment, denoted $\log(\bar{\epsilon})$. Here we compare the distributions of the two most coupled and the three least coupled dust species (colored dashed lines) at $t_{core} + 2$ kyr. We also indicate the PDFs integrated over the grain size distribution (black line). These PDFs are displayed in three different regions, namely the core and the fragments (left), the disk (middle), and the envelope (right). Figure 3 shows that the dust enrichment in the inner regions is size dependent. Smaller grains experience larger drag that reduces their differential velocity with respect to the gas.

Here, grains with sizes smaller than a few microns remain very well coupled with the gas in all the considered regions, whereas for larger grains the PDF of the dust ratio is broad. For 160 µm grains, the dust ratio increases by one order of magnitude in some regions of the disk and up to two orders of magnitude in the envelope. On average, the dust ratio is 0.018 in the core, 0.0175 in the disk, and ~0.0086 in the envelope. In addition, the dust has experienced a strong and local dynamical sorting. We indeed measure a typical standard deviation for the dust-ratio enrichment $\bar{\epsilon}$ of 0.072 in the core, 0.14 in the disk, and 0.23 in the envelope. The standard deviation is the largest in the envelope, a region which is depleted in dust in the outer regions and enriched close to the disk and fragments. The disk experiences larger
<table><caption>**Table 1.** Syllabus of the different simulations, with the thermal-to-gravitational energy ratios $\alpha$, and the maximum grain sizes $S_{\text{max}}$.</caption><thead><tr><th>Model</th><th>α</th><th>S<sub>max</sub> (cm)</th><th>μ</th><th>φ<sub>mag</sub> (°)</th><th>Ambipolar</th><th>t<sub>core</sub> (kyr)</th><th>St<sub>0,10</sub></th></tr></thead><tbody><tr><td>MMMRN</td><td>0.5</td><td>0.1</td><td>—</td><td>—</td><td>—</td><td>72.84</td><td>1.22 × 10<sup>-2</sup></td></tr><tr><td>MRN</td><td>0.5</td><td>2.5 × 10<sup>-5</sup></td><td>—</td><td>—</td><td>—</td><td>73.6</td><td>1.06 × 10<sup>-5</sup></td></tr><tr><td>100MICMRN</td><td>0.5</td><td>0.01</td><td>—</td><td>—</td><td>—</td><td>72.9</td><td>1.72 × 10<sup>-3</sup></td></tr><tr><td>MMMRNA0.25</td><td>0.25</td><td>0.1</td><td>—</td><td>—</td><td>—</td><td>23</td><td>4.31 × 10<sup>-3</sup></td></tr><tr><td>MMMRNMHD</td><td>0.5</td><td>0.1</td><td>5</td><td>40</td><td>NO</td><td>81.1</td><td>1.22 × 10<sup>-2</sup></td></tr><tr><td>MMMRNNIMHD</td><td>0.5</td><td>0.1</td><td>5</td><td>40</td><td>YES</td><td>81.1</td><td>1.22 × 10<sup>-2</sup></td></tr><tr><td>Initial cloud mass<br>1 M<sub>☉</sub></td><td>Number of dust bins<br>10</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

**Notes.** The initial mass-to-flux ratio $\mu$ as well as the tilt between the magnetic field and the rotation axis $\phi_{\text{mag}}$ are given for simulations with magnetic field. Additionally, we provide the formation time of the FHSC $t_{\text{core}}$, the initial Stokes number of the largest grains $St_{0,10}$, the mass of the initial core, and the number of dust bins.

variations of dust-to-gas ratio compared to the core. Indeed, the latter is in adiabatic contraction. This induces high temperatures, which in return causes a strong decrease of the Stokes number. Hence, dust is essentially frozen with the gas in the core and $v_d \approx v_g$. Figure 4 shows the dust ratio for the total dust distribution (left) and for the tenth species (right) in the mid-plane of the collapse, at 1 kyr (top), 2 kyr (middle), and 4 kyr (bottom) after the formation of the first Larson core, respectively. The maps on the left and on the right are very similar as most of the evolution of the dust ratio is due to the dynamics of the least coupled species, which represents a large fraction of the dust mass (see also the PDF). The structures in the dust ratio observed in Fig. 4 can be interpreted by looking at the thermal pressure distribution shown in Fig. 5. Dust grains tend to drift toward local pressure maxima (see Fig. 4, top panel) where the differential velocity is zeroed (see Eq. (6)). The greatest variations of the total dust ratio are due to the largest grains, since they represent most of the dust mass and have the largest drift velocities. Hence, a significant fraction of the dust mass in the inner regions is composed of large grains. Finally, we note that 4 kyr after the first core formation, the average value for the total dust-to-gas ratio remains roughly unchanged but generally increases (of about 1%) in the core and fragments after their formation. For the disk, we note a decrease in the dust ratio of ~30% for the largest grains (~22% in total) in the disk following the first core formation. In addition, the total dust-to-gas ratio continues to diminish in the envelope at this time as settling goes on, with a final average value of ~0.0083.

Figure 6 shows the evolution of the dust-to-gas ratio enrichment averaged in volume for different density thresholds in the regions where $R < 5000$ AU. The dust depletion at large scales in the envelope at low densities is a relatively slow process, which occurs during the entire collapse. Once regions with larger densities are formed, they usually experience a relatively quick enrichment from the dust content of the low-density regions, and then a quick depletion in favor of even denser regions. Figure 6 shows that fragmentation occurs in a dust-rich environment. A strong enrichment of the volume where $\rho > 10^{-13}$ g cm$^{-3}$, delimited by the brown line, indeed happens exactly when fragments form. This explains why these fragments tend to be more dust-rich than the first hydrostatic core. Interestingly, the dust-to-gas ratio does not vary significantly in the volume where $\rho > 10^{-11}$ g cm$^{-3}$ which is in adiabatic contraction. Again, the temperatures in this region are high and therefore the Stokes numbers are very low, which significantly slows down the differential

dynamics of the gas and the dust. We note that this volume is already dust enriched by almost a factor of two by the time of its formation.

## 4.2. Parameter exploration

### 4.2.1. Maximum grain size

As seen in Sect. 4.1 for the MMMRN model, the differential dynamics between gas and dust during the protostellar collapse depend critically on grain size. Therefore, we perform two simulations 100MICMRN and MRN with the same set of parameters as in MMMRN, but with a different maximum grain size. We choose $S_{\text{max}} = 100$ µm for the 100MICMRN model (which yields $s_{10} \sim 22.6$ µm) and $S_{\text{max}} = 250$ nm for the MRN model (which yields $s_{10} \sim 139$ nm).

Let us first consider MRN, which is the model that has the smallest $S_{\text{max}}$. Figure 7 shows the density of the gas (left) of the least coupled dust species (right) at $t_{\text{core}} + 2$ kyr. Because the coupling between gas and dust is almost perfect, the two maps are indistinguishable by eye. This is expected because the maximum grain size is $\approx 10^{-5}$ cm which corresponds to an initial Stokes number $St_{0,10} \sim 1.06 \times 10^{-5} \ll 1$ (see Sect. 6 for a theoretical justification). As a result, dust is an excellent tracer of gas in this model. This is illustrated by Fig. 8, which shows the probability density function of the dust ratio enrichment $\log(\bar{\epsilon})$ in the core and fragments, the disk, and the envelope at $t_{\text{core}} + 2$ kyr. Contrary to the MMMRN model, these PDFs are strongly peaked. The average dust-to-gas ratio integrated over the total dust distribution is $\approx 1\%$ in the core and the fragments, the disk, and the envelope. The standard deviation for the dust ratio enrichment ranges between $2 \times 10^{-4}\%$ (in the core) and $1.4 \times 10^{-2}\%$ (in the envelope). For this model, the variations of the dust ratio are very small. Therefore, in absence of coagulation, one may expect a standard MRN distribution to remain extremely well preserved during the protostellar collapse at all scales.

To investigate an intermediate scenario, we now focus on the 100MICMRN model. We do not show the density maps in this case due to their strong resemblance with the MRN case. The PDFs of the dust ratio enrichment $\log(\bar{\epsilon})$ in the core and fragments, the disk, and the envelope at $t_{\text{core}} + 2$ kyr are shown in Fig. 9. In this case, the variations of dust ratio are more significant than in MRN. However, compared to the MMMRN case, these variations still remain quite small. The average dust-to-gas ratio is 0.0106 in the core and the disk and 0.0099 in
<figure><img src="image_4.png" /></figure>

**Fig. 2.** MMMRN test at $t_{\text{core}} + 1$ kyr and $t_{\text{core}} + 2$ kyr ($t_{\text{core}} = 72.84$ kyr). Edge-on and mid-plane cuts of the gas and the dust densities for the least coupled species are provided (left and rights respectively). Values of the gas density are indicated by the colorbar at the bottom. Dust densities have been multiplied by a factor 100 to be represented on the same scale. Hence, colors of the gas and the dust maps match when the dust-to-gas ratio equals its initial value 0.01. Dust density variations in regions where $\epsilon_0\rho_d < \min(\rho_g)$ have voluntarily not been displayed to highlight the enriched regions. These depleted regions are delimited by the dashed gray lines. This choice of colors applies for all density maps in this study. Gas and dust are clearly not perfectly coupled.

the envelope. The typical standard deviation for the dust-ratio enrichment are $7 \times 10^{-3}$ in the core, $2.2 \times 10^{-2}$ in the disk, and $7.8 \times 10^{-2}$ in the envelope. This confirms that the larger the grains are, the more significant the decoupling with the gas. We note that, for 100MICMRN, it is reasonable to infer the gas density from the dust.

### 4.2.2. Thermal-to-gravitational energy ratio

The free-fall timescale depends on the ratio between the thermal and the gravitational energy. We therefore present in this section MMMRNA0.25, a model similar to the reference case but with $\alpha = 0.25$. This parameter is expected to strongly affect the dust dynamics. A lower value of $\alpha$ produces faster protostellar collapses prior to the first core formation due to the smaller thermal support. It therefore develops high-density regions more quickly, where dust strongly couples. In addition, a cloud with a smaller $\alpha$ has a smaller initial Stokes number, which means that the dust is also initially better coupled with the gas. The post-core evolution of MMMRNA0.25 is different from that in MMMRN in virtue of a smaller initial disk radius. Smaller disks with a smaller initial value of $\alpha$ are expected as shown by Hennebelle et al. (2016; their Eq. (14)).

Figure 10 shows the densities of the gas (left) and the least coupled dust species (right) at $t_{\text{core}} + 2$ kyr for MMMRNA0.25. Apart from the outer regions that are quite depleted in dust content, we do not see any significant difference between gas and dust. Indeed, the core forms quickly, leaving no time for the differential motion between gas and dust to develop. For a comparison with the fiducial case, we show in Fig. 11 the PDF of the dust ratio enrichment $\log(\bar{\epsilon})$ in the core and the fragments, the disk, and the envelope at $t_{\text{core}} + 2$ kyr. The PDFs are much more peaked in MMMRNA0.25 than in MMMRN. The values of the standard deviation of the dust-ratio enrichment are $2 \times 10^{-2}$ in the FHSC and fragments, $3 \times 10^{-2}$ in the disk, and 0.1 in the envelope. This was actually expected as the initial Stokes number scales as $\alpha^{2/3}$ and is thus $\approx 0.6$ times smaller in MMMRNA0.25 than in MMMRN.

### 4.2.3. Magnetic fields

We now consider the dynamics of neutral grains in collapsing magnetized clouds. The two models are performed with the same parameters as MMMRN but with an initial magnetic field given by $\mu = 5$ and a tilt of 40°. For MMMRNMHD we use an ideal MHD solver and for MMMRNNIMHD we consider ambipolar diffusion.

Figures 12 and 13 show the densities for the gas (left) and the least coupled dust species (right) at $t_{\text{core}} + 2$ kyr for the two models MMMRNMHD and MMMRNNIMHD, respectively. For both models, the dust is significantly decoupled from the gas and the settling in the core/disk/pseudo-disk is very efficient. As in MMMRN, dense regions (disk, core, pseudo-disk, and high-density regions of the outflow) are prone to be enriched in solid particles while low-density regions are depleted (envelope and low-density regions of the outflow). We note that the decoupling is slightly more efficient in these models than in MMMRN. This is mostly due to the 10 kyr difference in free-fall timescale. Although additional decoupling terms due to the magnetic field appear in the dust differential velocity ($\propto J \times B$ see Eq. (9)) those are negligible compared to the hydrodynamical terms ($\propto \nabla P_g$) in the decoupling of gas and dust in our magnetised collapse models.

Figures 14 and 15 show the PDFs of the dust ratio enrichment for the different objects at $t_{\text{core}} + 2$ kyr for the MMMRNMHD and MMMRNNIMHD runs, respectively. Although the shapes of the distributions are different from our fiducial case, we essentially reach the same conclusion, that is, a peaked distribution with a significantly large average, indicating a strong initial enrichment in the [core+disk] system and a broad distribution in the envelope. We note that the average dust-to-gas ratio in the disk
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** MMMRN test at $t_{\text{core}} + 2$ kyr: probability density function (PDF) of the dust ratio enrichment $\log(\bar{\epsilon})$ for the two most coupled and three least coupled dust species (colored dashed lines) and for all the dust (black line) in the core+fragments (**left**), the disk (**middle**), and the envelope (**right**). Dust is not a good tracer for the gas here and the dust distribution is not uniform in the considered objects.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** MMMRN test ~1 kyr (**top**), ~2 kyr (**middle**), and ~4 kyr (**bottom**) after the first core formation ($t_{\text{core}} = 72.84$ kyr). Mid-plane view of the total dust ratio (**left**) and the dust ratio of the least coupled species (**right**). The colorbar is the same for both figures. The dotted white lines represent the regions where the total (**left**) or 160 µm grain (**right**) dust-to-gas ratio is at its initial value, which can also be regarded as a dust enrichment line.</figcaption></figure>

<caption>**Fig. 5.** MMMRN test at $t_{\text{core}} + 1$ kyr ($t_{\text{core}} = 72.84$ kyr). Mid-plane view of the gas pressure (up-close). The white arrows represent the direction of the differential velocity $w_{10}$.</caption>

and the core is higher in these two models than in the fiducial case. We indeed measure an average dust-to-gas ratio of ~0.022–0.023 in the disk and the first hydrostatic core for these models. In these magnetic runs, the pinching of the magnetic

field lines during the collapse produces a pseudo-disk, which is a dense but not rotationally supported region. We note that these pseudo-disks have a very broad PDF and show enriched and depleted regions in both the ideal and nonideal cases. Similar behavior is observed in the magnetically driven outflow for the ideal case, which is dust-rich in dense regions and depleted at low densities. In the MMMRN NIMHD, the outflow is less evolved and is mostly dust-depleted, similarly to the envelope.

In summary, neutral dust dynamics in the presence of magnetic fields seem to follow the same general trend as in the hydrodynamical case. Dust collapses faster than the gas and is enriched in the inner regions of the collapse a few thousand years after the first core formation. This enrichment is mainly located in the pseudo-disk (only observed in the magnetized models), the disk, the first hydrostatic core, and the high-density regions of the outflow.

## 5. Features of dusty collapses

We summarize here the properties of the dusty collapse in its different regions. Figure 16 shows the dust-to-gas ratio enrichment averaged in mass as a function of the grain size for all the models and all the objects defined in Sect. 3.4.3. We refer to the dashed horizontal line as the enrichment line; if an object lies
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** MMMRN test. Volume averaged total dust-to-gas ratio enrichment as a function of time for different density thresholds. The FHSC formation can be identified by the dotted vertical blue line while fragmentation occurs during a time delimited by the green area. We observe a slow decrease of the dust-to-gas ratio at low density which gives rise to an enrichment of high-density regions. Cores and fragments at $\rho > 10^{-11}$ g cm$^{-3}$ are formed in a dust-rich environment. The dust-to-gas ratio is almost constant for $\rho > 10^{-11}$ g cm$^{-3}$ because the increase of temperature due to the adiabatic contraction strengthens the coupling between the gas and the dust.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** MRN test at $t_{core} + 2$ kyr ($t_{core} = 73.6$ kyr). Edge-on (top) and mid-plane (bottom) cuts of the gas density (left) and dust density of the least coupled species (right). The two maps are almost indistinguishable due to the very strong coupling between gas and all dust species.</figcaption></figure>

above this line, it is enriched in dust during the collapse; if it lies underneath, it is dust depleted. This information is collected in Table D.2.

## 5.1. Core and fragments

Here we detail the dust and gas properties in the first hydrostatic core and the fragments. In every model, we observe a value

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** MRN test at $t_{core} + 2$ kyr. Probability density function (PDF) of the total dust ratio enrichment $\log(\bar{\epsilon})$ in the core and the fragments, the disk, and the envelope. Dust is a very good tracer of mass here and the dust distribution is almost uniform in the objects considered.</figcaption></figure>

of $\langle \Theta_{d,k \le 6} \rangle_m$ that remains unchanged for all the cores and fragments. Indeed, small grains have short stopping times and remain very well coupled to the gas all along the collapse. Simulations with the largest maximum grain size ($S_{max} = 0.1$ cm) have the largest dust enrichment. On the contrary, for the MRN model, the dust-to-gas ratio preserves its initial value in all the fragments. Moreover, the dust distribution itself remains extremely well preserved. In the MMMRNA0.25 case, the enrichment of the largest dust grains is much less efficient than in MMMRN. This is explained by a shorter free-fall timescale and higher initial densities, which implies smaller initial Stokes numbers. In MMMRN, MMMRNA0.25, and 100MICMRN, the fragments have greater enrichment in larger grains ($k > 6$, see Table D.1 for the corresponding grain sizes). For example, in the MMMRN case, the dust-to-gas ratio of the 160 µm grains is enriched by a
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** 100MICMRN test at $t_{core} + 2$ kyr. Probability density function of the total dust ratio enrichment $log(\bar{\epsilon})$ in the core and the fragments, the disk, and the envelope. Dust is a relatively good tracer of gas although significant variations of the dust-to-gas ratio are observed.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** MMMRNA0.25 test at $t_{core} + 2$ kyr. Probability density function of the total dust ratio enrichment $log(\bar{\epsilon})$ in the core+fragments, the disk, and the envelope. Dust is a relatively good tracer for gas in the dense object although a notable depletion is observed in the envelope.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** MMMRNA0.25 test at $t_{core} + 2$ kyr ($t_{core} = 23$ kyr). Edge-on (top) and mid-plane (bottom) cuts of the gas density (left) and dust density of the least coupled species (right). Dust has less time to significantly decouple from the gas than in the MMMRN case. A strong dust depletion is observed in the envelope.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** MMMRNMHD at $t_{core} + 2$ kyr ($t_{core} = 81.1$ kyr). Edge-on (top) and mid-plane (bottom) cuts of the gas density (left) and dust density of the least coupled species (right). Dust is significantly decoupled from the gas and concentrates in the high-density regions such as the core, the disk, the pseudo-disk, and the inner regions of the outflow.</figcaption></figure>

factor of ≈2.6 in the central object, and ≈3 in the fragments. We note that the dust-to-gas ratio enrichment of the first hydrostatic core is stronger in simulations that include magnetic fields than in MMMRN. As explained above, magnetic fields bring an additional decoupling in the case of neutral grains which explains in part why the dust enrichment is even stronger in MMMRNMHD and MMMRNNIMHD than in MMMRN. More importantly, the collapse is longer for these models due to the magnetic support. This leaves more time for dust grains to enrich the central

regions. For the fragmenting cases, we note a preferred concentration of dust in the fragments, which can be explained by two mechanisms. First, fragments form after the central object and thus remain for longer in the isothermal phase where dust is less coupled because the temperature is lower. Second, the dust-rich spiral arms developing through the envelope (see Fig. 4) are mainly accreted by the fragments (see Fig. 6). This provides an additional channel to enrich the fragments in solids.
<figure><img src="image_14.png" /></figure>

**Fig. 13.** MMMRNNIMHD at $t_{\text{core}} + 2$ kyr ($t_{\text{core}} = 81.1$ kyr). Edge-on (top) and mid-plane (bottom) cuts of the gas density (left) and dust density of the least coupled species (right). As in the ideal case, dust is significantly decoupled from the gas and concentrates in the high-density regions such as the core, the disk, the pseudo-disk, and the inner regions of the outflow.

<figure><img src="image_15.png" /></figure>

**Fig. 14.** MMMRNMHD test at $t_{\text{core}} + 2$ kyr. Probability density function of the dust ratio enrichment $\log(\bar{\epsilon})$ in the core (blue), the disk (green), the pseudo-disk (purple), the outflow (orange) and the envelope (red). Dust is not a good tracer of gas here and the dust distribution is not uniform in the considered objects.

## 5.2. Disks

We review the disk properties of all the models at $t_{\text{core}} + 2$ kyr. In general, the values are similar to what is measured in the cores. This is mostly caused by the fact that the dust enrichment happens prior to core formation at low densities (see Fig. 6). We measure a total dust-to-gas ratio enrichment

<figure><img src="image_16.png" /></figure>

**Fig. 15.** MMMRNNIMHD test at $t_{\text{core}} + 2$ kyr. Probability density function of the dust ratio enrichment $\log(\bar{\epsilon})$ in the core (blue), the disk (green), the pseudo-disk (purple), the outflow (orange) and the envelope (red). Dust is not a good tracer of gas here and the dust distribution is not uniform in the considered objects.

<figure><img src="image_17.png" /></figure>

**Fig. 16.** All the models at $t_{\text{core}} + 2$ kyr. The dust-to-gas ratio enrichment averaged in mass is shown as a function of grain size for all the objects. Grains with sizes smaller than $10^{-4}$ cm are almost always perfectly coupled with the gas. For larger grain sizes, the enrichment is model dependent. Grains with typical sizes larger than $10^{-3}$ cm decouple from the gas. Dense objects such as the fragments $\mathcal{F}$ or the disk $\mathcal{D}$ and pseudo-disk $\mathcal{P}$ are enriched in dust. Low-density objects such as the envelope $\mathcal{E}$ or the outflows $\mathcal{O}$ are depleted in dust. Magnetized models exhibit the strongest decoupling between the gas and the dust.

of ~1.75 in MMMRN, ~1 in MRN, ~1.1 in 100MICMRN, ~1.1 in MMMRNA0.25, ~2.3 in MMMRNMHD, and ~2.1 in MMMRNNIMHD. We emphasize once again that the decoupling between the gas and the dust depends strongly on the initial properties of the cloud. We note that in MMMRN the dust ratio is highly nonuniform in the disk (see Figs. 3 and 4 for the MMMRN case) or even constant (see Fig. 6). As explained in Sect. 4.1, there is a decrease of ~22% of dust-to-gas ratio in MMMRN between $t_{\text{core}}$ and $t_{\text{core}} + 4$ kyr. This is most likely due to the fact that the disk is accreting dust-depleted material from the envelope. In addition, since dust drifts toward pressure bumps – or regions where $\mathbf{J} \times \mathbf{B} \sim \nabla P_g$ for magnetic runs – dust cannot always be used as a direct proxy to trace the gas density. Although the substructures seen in the dust originate from
those in the gas, as shown in Fig 5, they may have very distinct morphologies. This is similar to what is found for T Tauri disks, where gaps could be opened in the dust only (Dipierro & Laibe 2017). Once again, simulations with magnetic fields show a more significant dust enrichment because of the longer timescale of the collapse. We note that the disk masses are however much smaller in the two models where magnetic braking occurs. Indeed, a longer integration time is required for the disk to grow significantly (Hennebelle et al. 2020).

### 5.3. Pseudo-disks

In this section, we describe the principal features of the dusty pseudo-disks that are observed in the two magnetic runs MMMRNMHD and MMMRNNIMHD. These pseudo-disks are strongly enriched with a total dust-to-gas ratio enrichment of approximately two for both cases (see values in Table D.2). Interestingly, the pseudo-disk is generally more enriched in smaller grains (up to 47 µm grains) than the other objects. This is due to two effects. First, the pseudo-disk is less dense than the central regions of the collapse, namely the disk and the core; this explains why smaller grains drift more easily towards it. Second, the strong pressure gradient orthogonal to the pseudo-disk generates a drift from the envelope towards it, even for small grains. Once these grains have reached the pseudo-disk, they couple strongly to the gas while larger grains are able to drift to even deeper regions such as the disk and the core.

### 5.4. Outflows

We now describe the major features of the dusty outflows that can be observed in the two magnetic runs MMMRNMHD and MMMRNNIMHD. For MMMRNMHD, Fig. 17 shows the relative variations of the dust ratio at three different times, for the 47 µm (left) and 160 µm (right) grains (ninth and tenth bins), respectively. The magenta arrows represent the differential velocity with the barycenter. These two dust species have a completely different evolution. Indeed, the outflow does not carry a significant quantity of 160 µm grains at $t_{core} + 2$ kyr because they are already strongly depleted in low-density regions. Subsequently, the outflow is strongly depleted in these species with $\langle \bar{\Theta}_{d,10} \rangle_m \sim 0.22$ for MMMRNMHD and $\langle \bar{\Theta}_{d,10} \rangle_m \sim 0.44$ for MMMRNNIMHD. On the contrary, 47 µm grains are significantly enriched by a factor 1.04–1.28 at that time. Initially, the outflow is not powerful enough to eject matter from the inner regions and rather collects the grains from the envelope. This explains why the enrichment measured in the outflow is similar to those measured in the envelope. Interestingly, at $t_{core} + 2$ kyr, we see that the outflow is well established and starts to carry the inner regions that are denser and more enriched in 160 µm grains. This indicates that outflows provide a channel to re-enrich the envelope in large grains.

### 5.5. Envelope

In general, the envelope is dust depleted owing to mass conservation since the inner parts are enriched. For all the models, the small dust grains (with $k \le 6$) do not experience significant dust-to-gas ratio variations. Larger grains can however show significant dust-to-gas ratio variations. We measure a value of $\langle \Theta_{d,10} \rangle_m$ as low as 0.38 for MMMRNNIMHD and MMMRNMHD, and 0.73 for MMMRN. These values are averaged over the whole envelope, but contrary to the core and fragments, the envelope is

very contrasted in terms of density. Among all objects (pseudo-disk excluded), the envelope experiences the largest dust-to-gas ratio variations (see e.g., Figs. 3 and 15). This is expected since the density in the envelope is lower than in the other objects. Typically, the large-grain depletion in the envelope increases with decreasing density and there is a larger dust content in its inner regions (see Fig. 17 for both behavior).

<figure><img src="image_18.png" /><figcaption>Fig. 17. mmMRNmhd. Edge-on view of the relative variations of the dust ratio at four different times for the 47 µm (left) 160 µm grains (right). The magenta arrows represent the differential velocity with the barycenter. Regions that are dust depleted by more than two orders of magnitude are not displayed (black background).</figcaption></figure>
We have shown that dust is not necessarily a good tracer of gas density, that is, we find important local variations in the dust distribution in MMMRN, MMMRNMHD, and MMMRNNIMHD. This has significant consequences for observations, since dust continuum radiation fluxes depend on densities integrated along the line of sight. It is therefore interesting to estimate the error that arises when the total column density is estimated from the mass of a single dust bin k. We compute this error Errk from the total column density and the dust column density Σd,k according to

$$Err_k = \frac{\Sigma - \Sigma_{d,k} / \epsilon_{k,0}}{\Sigma}. \qquad (22)$$

The total error Err is defined in the same way but using the total dust column density and initial dust ratio.

Figure 18 shows an edge-on view of the total column density Σ (top), the total error Err (middle), and the error estimated using the largest grains only – 160 µm in this case (bottom) for the MMMRNNIMHD model 2 kyr after the formation of the first core. The error is large when considering either all the grains (middle) or only the largest ones (bottom). We note that the total column density inferred from the total dust mass would be underestimated in the upper layers and overestimated in the inner regions because dust drifts toward the center of the collapse. The effect is maximal for the largest grains, where the error can reach values as high as ~250% in the inner envelope.

# 6. Estimate of the dust enrichment

Here, we provide a semi-analytic estimate of the dust enrichment occurring during the protostellar collapse. We use it to infer the typical minimal Stokes number above which a given dust enrichment can be reached.

## 6.1. Enrichment equation

In the terminal velocity approximation, and making the assumption that the collapse is isothermal and purely hydrodynamical, the evolution of the dust ratio for a species k is given by

$$\frac{d\epsilon_k}{dt} = -\frac{1}{\rho}\nabla \cdot \left[ \epsilon_k \left( \frac{t_{s,k}}{1-\epsilon_k} - \sum_{i=1}^{N} \frac{\epsilon_i}{1-\epsilon_i} \frac{t_{s,i}}{1-\epsilon_i} \right) c_s^2 \nabla \rho \left( 1 - \sum_{j=1}^{N} \epsilon_j \right) \right]. \quad (23)$$

To provide analytical estimates of dust-ratio enrichment during the collapse, we now neglect the cumulative back-reaction of dust onto the gas. When the cumulative back-reaction is negligible, the evolution of $\bar{\epsilon} = \epsilon/\epsilon_0$ and $\epsilon$ are constrained by the same equation, that is, $\bar{\epsilon}$ does not depend on $\epsilon_0$. The dust ratio enrichment is described by the equation

$$\frac{d\bar{\epsilon}_k}{dt} = -\frac{1}{\rho}\nabla \cdot \left[ \bar{\epsilon}_k t_s c_s^2 \nabla \rho \right], \qquad (24)$$

where $\frac{d}{dt} \equiv \frac{\partial}{\partial t} + \nu \cdot \nabla$. We now use the dimensionless variables

$$\tau = \frac{t}{\tau_{\text{ff,0}}} \quad \text{and} \quad x = \frac{r}{\lambda_{J,0}},$$

where $\tau_{\text{ff,0}} = \sqrt{\frac{1}{g\bar{\rho}_0}}$ and $\lambda_{J,0} = c_s \tau_{\text{ff,0}}$. Equation (24) becomes

$$\frac{d\bar{\epsilon}_k}{d\tau} = -\text{St}_{k,0} \frac{1}{\bar{\rho}} \nabla_x \cdot \left[ \bar{\epsilon}_k \bar{\rho}^{-1} \nabla_x \bar{\rho} \right], \qquad (25)$$

where $\bar{\rho} = \rho/\rho_0$. The term $\bar{\rho}^{-1}$ appears in the divergence owing to $t_{s,k} \propto \frac{1}{\bar{\rho}}$.

<figure><img src="image_19.png" /><figcaption>Fig. 18. MMMRNMHD at $t_{\text{core}} + 2$ kyr. Edge-on view of the total column density log(Σ) (top), the total error Err (middle) and the error when using the largest grains only – 160 µm in this case (bottom).</figcaption></figure>

## 6.2. Semi-analytical model

Neglecting local variations of $\bar{\epsilon}$ in comparison with local density variations yields

$$\frac{d\bar{\epsilon}_k}{d\tau} \approx -\text{St}_{k,0} \bar{\epsilon}_k \frac{1}{\bar{\rho}} \nabla_x \cdot \left[ \bar{\rho}^{-1} \nabla_x \bar{\rho} \right]. \qquad (26)$$
<figure><img src="image_20.png" /><figcaption>Fig. 19. Semi-analytic and measured enrichment in our hydrodynamical models (MRN excluded) against the initial Stokes number. The lines represent the semi-analytical development using a best fit of $\chi$ for the FHSC (dotted), disk (solid), and envelope (dashed), respectively. The extreme values of the enrichment given by our toy model are delimited by the blue areas. The color coding and choice of markers is the same as in Fig. 16.</figcaption></figure>

We then obtain

$$ \bar{\epsilon}(x, \tau) = \chi^{\text{St}_{k,0}}, \qquad (27) $$

where $\chi \equiv e^{-\int_0^\tau \nabla_s(\bar{\rho}^{-1}\nabla_s\bar{\rho})/\bar{\rho}d\tau}$ is independent of the dust properties. Hence, at a given time and position, the dust enrichment varies essentially exponentially with the initial Stokes number. A proper mathematical estimate of the integral quantity is beyond the scope of this study.

### 6.3. Estimate in the core

We can roughly approximate $\chi$ in the core and after a free-fall time $t_{\text{ff},0} = \sqrt{\frac{3\pi}{32}\tau_{\text{ff},0}}$. An order of magnitude estimate provides

$$ |\ln(\chi)| \approx \sqrt{\frac{3\pi}{32} \frac{c_s^2}{G \rho_{\text{ad}} r_{\text{var}}^2}}, \qquad (28) $$

where $r_{\text{var}}$ is the typical length at which variations of density become significant. Above $\rho_{\text{ad}}$, the temperature is high and the gas and dust differential dynamics is negligible. A reasonable choice for $\bar{\rho}$ is therefore $\bar{\rho} = \frac{\rho_{\text{ad}}}{\rho_0}$. In the typical condition of a protostellar collapse, one obtains

$$ |\ln(\chi)| \approx 126 \left( \frac{r_{\text{var}}}{1 \text{ AU}} \right)^{-2} \left( \frac{T_{\text{gas}}}{10 \text{ K}} \right). \qquad (29) $$

We note that taking $r_{\text{var}} = 1$ AU seems reasonable as it is about a tenth of the first-core radius in our models. For $r_{\text{var}} \approx 5$ AU, we find that $|\ln(\chi)| \approx 5$. The value of $\chi$ strongly depends on the steepness of the pressure gradients, and therefore on $r_{\text{var}}$. This model only provides a rough estimate of $\bar{\epsilon}$ and should not replace either a numerical treatment of the dust or a proper estimate of $\chi$ during the collapse.

### 6.4. Comparison with the models

Figure 19 shows the dust ratio enrichment as a function of the initial Stokes number for MMMRN, MMMRNA0.25, and

100MICMRN for the first core, the disk, and the envelope. The dashed, dotted, and solid lines represent the values obtained with Eq. (27) by fitting the values of $\chi$. Finally, the blue area represents the range of dust enrichment obtained with the two extreme values of $\chi$ estimated in the previous section. We do not show the dust enrichment for MRN as it is clearly negligible (see Fig. 16). In addition, we do not display the enrichment in the secondary fragments for the sake of readability.

Fairly good agreement between the fits and the measured dust ratio enrichment is observed in all the regions, especially in the cores. This suggests that the enrichment indeed mostly varies exponentially with the initial Stokes number. We do observe small deviations in the disk and the envelope for MMMRN. The nonlinear behavior of Eq. (23), either due to local variations of $\epsilon$ or the cumulative back-reaction of dust on the gas, is therefore not completely negligible. In Appendix B, we show that the dust-to-gas variations induced by the dust back-reaction are almost negligible and that the discrepancy with the exponential increase of the dust ratio observed is most likely caused by local variations in the dust-to-gas ratio; for example the mixing between depleted material of the envelope and dust-rich material from the disk.

Finally, we emphasize that the dust enrichment in all the cores is comprised between the lowest and highest values estimated with our toy model. As this model is relatively crude, we acknowledge that it cannot compete with an eventual model based on an accurate estimate of $\chi$.

### 6.5. Critical Stokes number

Assuming that the value of $\chi$ is known, it can be used to determine the critical Stokes number $St_{\text{crit},\bar{\epsilon}}$ above which some regions can reach a given enrichment $\bar{\epsilon}$. Equation (27) can indeed by inverted as

$$ St_{\text{crit},\bar{\epsilon}} = \frac{\ln \bar{\epsilon}}{\ln \chi}. \qquad (30) $$

Using the value of $\chi$ obtained by fitting our models (Fig. 19), we can estimate the typical Stokes numbers needed to get a dust enrichment by a factor of two in the core, and the disk is approximately $St_{\text{crit},2} \sim 0.01-0.027$. With our toy model, we find $St_{\text{crit},2} \sim 0.006-0.13$. This latter is a significantly wider range but it contains what is measured with our models. Similarly, we can estimate that to get a dust ratio depletion of 50%, the grains must typically have $St_{\text{crit},1/2} \sim 0.027-0.11$. In short, it is easier to enrich the disk and the core than it is to deplete the envelope.

## 7. Discussion

### 7.1. Summary of the models

We investigated the effect of several parameters on the dust dynamics during the protostellar collapse such as the thermal-to-gravitational energy ratio $\alpha$, the maximum grain size of the dust distribution, and the presence of magnetic fields. It appears that the first two parameters are critical for the development of significant differential dynamics between the gas and the dust during the early phases of the collapse.

From MMMRN, MRN, and 100MICMRN, we established that the maximum grain size is a critical parameter for the differential gas and dust dynamics. Typically, if the maximum grain size in the core is smaller than a few microns (see Fig. 16), we do not observe significant variations of the dust ratio. On the contrary, when larger grains are considered, the total dust ratio can
increase by a factor of two to three, even in the very first few
thousands of years of the protostellar collapse (see MMMRN,
MMMRN**MHD or MMMRN**NIMHD). Hence, in order to under-
stand the initial dust and gas content of protoplanetary disks, it
is crucial to accurately measure the dust size distribution in early
prestellar cores. This conclusion is reinforced by recent observa-
tions (Galametz et al. 2019) or synthetic observations (Valdivia
et al. 2019) that seem to probe the existence of ~100 µm in Class
0 objects.

For small initial values of α, for example in MMMRNA0.25,
the collapse is fast and large grains do not have the time to sig-
nificantly enrich the core and the disk in one free-fall timescale.
Besides, as MMMRNA0.25 is set with a higher initial density,
dust is initially more coupled with the gas in this particular
model. We point out that the efficiency of the dust enrichment
relies strongly on the lifetime of low-density regions (see Fig. 6)
and on the range of densities experienced during the collapse
(see Sect. 6). Therefore, although 100MICMRN has a smaller
maximum grain size than MMMRNA0.25, both models have a
similar total dust content by the end of the simulation as the free-
fall timescale in 100MICMRN is longer. The initial properties of
the cloud appear to be extremely important to quantify the evo-
lution of the dust distribution during the protostellar collapse.
It would therefore be interesting to study the dust collapse of a
Bonor-Ebert sphere, since its free-fall timescale is usually longer
than that of the Boss and Bodenheimer test (Machida et al. 2014).
We leave this proper comparison to future studies.

With MMMRN**MHD and MMMRN**NIMHD, we investigated
the effect of a magnetic field on the dynamics of dust during
the protostellar collapse. We find qualitatively similar results to
those in our fiducial case. Quantitatively, the decoupling between
the gas and the dust does however produce more significant vari-
ations of the dust-to-gas ratio in the magnetic case. The presence
of a dense and stratified pseudo-disk strengthens the envelope
and the outflow depletion. This pseudo-disk is consequently
strongly enriched in solids. It is in fact almost as enriched as
the disk and the first hydrostatic core. However, it is much more
massive than the disk by the end of the calculation. Therefore,
understanding how the pseudo-disk is accreted by the core and
the disk is of particular interest and future studies should focus
on its long-term evolution.

To summarize, Fig. 20 shows a schematic view of a dusty
protostellar collapse a few thousand years after the formation of
the first hydrostatic core. The blue areas depict the dust-depleted
regions (low-density regions of the envelope and outflow) and
the red areas represent the regions enriched in dust (cores, disk,
high-density regions of the envelope and outflow, and pseudo-
disk). The intensity of the gas and dust decoupling depends
naturally on the parameters that we presented earlier in this
article. We emphasize that this cartoon illustration is only a
simplified picture of a dusty protostellar collapse that does not
account for the variability between the models and we do not
quantitatively show the local variations of the dust-to-gas ratio.

## 7.2. Comparison with previous works

Our results are in qualitative agreement with the previous study
of Bate & Lorén-Aguilar (2017), where a decoupling between
gas and dust for grains larger than ≈100 µm was also identified.
A main difference is that we do not obtain equally large dust-to-
gas ratio enhancements. Indeed, in Bate & Lorén-Aguilar (2017),
the dust mass is distributed in a single bin of dust with a mass
of 1% of the mass of the gas. In our MULTIGRAIN simulations,
only a fraction of the dust mass lies in the largest grains, which

<figure><img src="image_21.png" /><figcaption>Schematic view of a dusty protostellar collapse</figcaption></figure>

**Fig. 20.** Schematic view of a dusty protostellar collapse. Blue regions are dust depleted and red regions are dust enriched. Typically, the outer regions of the envelope are depleted. The outflow, only observed in magnetic runs, is enriched on its surface and depleted elsewhere. Dense regions such as the core and fragments *F* (the fragments are observed only in the hydrodynamical case), the pseudo-disk *P* (only in magnetic runs), and the disk *D* are generally enriched. The strength of dust decoupling depends on the initial choice of parameters such as the maximum grain size, the thermal-to-gravitational energy ratio, or the presence of a magnetic field. This view of a dusty protostellar collapse is simplified and provides only a global sketch of the evolution.

provides less significant dust-to-gas ratio variations than in Bate
& Lorén-Aguilar (2017). This effect was predicted in Bate &
Lorén-Aguilar (2017). We note that $\bar{\epsilon}$ only depends on the initial
dust content via the cumulative back-reaction of the dust on the
gas. If this back-reaction is neglected as in Eq. (24), the enrich-
ment is independent of the initial value of the dust-to-gas ratio.
This allows us to make a more direct comparison with Bate &
Lorén-Aguilar (2017). In their study, these latter authors observe
an increase in the dust-to-gas ratio of about one order of mag-
nitude for the 100 µm grains, which is about three times larger
than what we observe in MMMRN for example. Using the cal-
ibrated value of $\chi$ and Eq. (27), we can estimate that the dust
ratio enrichment of 100 µm for grains with $\rho_{\text{grain}} = 3 \text{ g cm}^{-3}$
would be ~2.27 in the first hydrostatic core. The difference with
Bate & Lorén-Aguilar (2017) is likely due to their use of Bonor-
Ebert spheres as initial conditions that have a longer free-fall
timescale (~120 kyr) and because they have lower initial den-
sities (~$10^{-20}$ g cm$^{-3}$) and therefore larger initial Stokes numbers
(~1 for 100 µm grains). In a 2D simulation of collapsing gravito-
viscous protoplanetary disks, Vorobyov et al. (2019); Vorobyov
& Elbakyan (2019); Elbakyan et al. (2020) focus on the evolu-
tion of dust including grain growth and fragmentation. Similarly
to our models, these latter authors observe local variations of
the dust-to-gas ratio in the disks. They also find larger dust-
to-gas ratios in the inner regions of the collapse and smaller
dust-to-gas ratios a few hundred astronomical units away from
the core. Furthermore, in their high-density clumps, they find
dust-to-gas ratios between 1.7 and 2.3% which is quite similar
to our findings. Locally, Vorobyov & Elbakyan (2019) observe
a particularly large increase in the dust-to-gas ratio in density
clumps, which is typically what we observe in our secondary
fragments.

In Lebreuilly et al. (2019), we performed three collapse sim-
ulations of nonrotating gas and dust mixtures, considering only
single dust species (1, 10 and 100 µm). We already observed a significant decoupling occurring for grains as little as 100 µm. However, the increase of dust-to-gas ratio in the core was strong only in the outer regions of the collapse. As mentioned above, we can estimate that in MMMRN the dust-to-gas ratio enrichment of 100 µm with $\rho_{\text{grain}} = 3 \text{ g cm}^{-3}$ would be about 2.27. This value was only ~1.2 in our first nonrotating spherical collapse calculation, although both models have the same initial $\alpha$. We interpret this result as an effect of rotation. Firstly, because it slows down the collapse (by a factor ~0.87 here), which leaves more time for the central regions to be enriched in dust. Secondly, because it generates steeper vertical pressure gradients which allows a more efficient settling of the dust grains. We do not aim to investigate the effect of the initial angular velocity in detail because this was done by Bate & Lorén-Aguilar (2017). The initial angular velocity was found to simply enhance the differential gas and dust dynamics similarly to what the thermal-to-gravitational energy ratio would do. We choose not to explore the impact of grain density because the dependence of the Stokes number in this quantity is the same as for the grain size.

### 7.3. Possible implications for planet formation

The simulations presented in this study consolidate the idea that protostellar collapses may form protoplanetary disks containing $\epsilon_0 \ge 2-3\%$ of their mass under the form of solids. When the dust ratio $\epsilon_0$ is larger than the square of the aspect ratio of the disk ($H/r$)² – even by a tiny amount – grain growth is expected to occur so efficiently that pebbles can decouple from the gas before drifting and falling onto the central star (Laibe 2014). This condition is likely to be fulfilled as protoplanetary disks typically have ($H/r$)² ≈ 0.01 (Andrews et al. 2010). The former condition is strengthened by the fact that the back-reaction may also inhibit radial drift and vertical settling (e.g., Kanagawa et al. 2017; Dipierro et al. 2018; Lin 2019), and holds until grains fragment. Two scenarios have been debated when fragmentation occurs. In the first scenario, grains may fall onto the central star if the disk does not contain a pressure trap (e.g., Brauer et al. 2008; Birnstiel et al. 2009). In the second, dust may exert strong drag on the gas and power the development of self-induced dust traps (Gonzalez et al. 2017). The formation of these traps occurs when back-reaction dominates locally over gas viscosity, and may be extremely effective for $\epsilon_0 \ge 2-3\%$. In any case, large dust contents favor the formation of planetesimals through the development of the streaming instability (e.g., Johansen et al. 2007, 2009; Drążkowska et al. 2016). This instability may be even more effective when it develops through unstable epicyclic modes (Jaupart & Laibe 2020), although peculiar dust distributions may quench it (Krapp et al. 2019). As argued by Vorobyov & Elbakyan (2019), dust-rich density clumps, similar to the secondary cores in our fragmenting models, could be a favored locus of giant planet formation. These latter authors indeed noted both the piling-up of large grains and important growth in these clumps.

In short, a larger initial dust content always favors planet formation in disks, and this may be in a dramatic manner. Quantitative information on the differential dynamics of gas and dust during the protostellar collapse and the initial dust size distribution in prestellar cores therefore appears to be essential to understanding the early stages of planet formation. In that perspective, an extensive study of dust dynamics and coagulation and fragmentation should be done from the scales of molecular

clouds, through the protostellar collapse, up to protoplanetary disks.

### 7.4. Neutral grains approximation

Grains are likely to be significant or even the main charge carriers during the protostellar collapse (Umebayashi & Nakano 1990; Marchand et al. 2016). Charged dust fluids feel the Lorentz drag $f_{L,k}$ in addition to the Epstein drag. The expression of this force is

$$f_{L,k} = \frac{3Z_k e}{4\pi\rho_{\text{grain},k} s_{\text{grain},k}^3} (\mathbf{E} + \mathbf{v}_k \times \mathbf{B}), \quad (31)$$

where $Z_k$ is the number of charges on the dust grain $k$, $e = 1.6 \times 10^{-20}$ Abc the electron charge and $E$ is the electric field. The gyration time of a grain is expressed as

$$t_{\text{gyr},k} = \frac{4\pi\rho_{\text{grain},k} s_{\text{grain},k}^3}{3Z_k e |\mathbf{B}|}, \quad (32)$$

A grain subject to a significant Lorentz drag could preferentially couple to the magnetic field rather than the gas. In addition, by controlling the Ohmic, ambipolar, and Hall resistivities (Kunz & Mouschovias 2009), charged grains most likely affect the evolution of the magnetic and electric fields. To compare the magnetic and neutral drag strength, it is interesting to compare the stopping-to-gyration time ratio also known as the Hall factor $\Gamma_k$. In the Epstein regime, this ratio writes

$$\Gamma_k = \sqrt{\frac{9\gamma}{128\pi} \frac{Z_k e |\mathbf{B}|}{\rho c_s s_{\text{grain},k}^2}}. \quad (33)$$

To model the Hall factor in a collapsing core, we consider a gas with a barotropic equation of state with $T_g = 10$ K, an adiabatic index $\gamma = 5/3$, and negatively charged grains. The grain charge computation follows Wurster et al. (2016) and is detailed in Appendix C. We use simple assumptions for the magnetic field, stating that

$$|\mathbf{B}| = \min \left( B_0 \left( \frac{\rho}{\rho_0} \right)^{2/3}, 0.1 \text{G} \right), \quad (34)$$

where $B_0$ is the initial magnetic field given by $\mu$, and $\rho_0$ is the initial core density. The magnetic threshold at 0.1 G is imposed to reproduce the plateau systematically observed when considering ambipolar diffusion Masson et al. (2016), Hennebelle et al. (2016) and Vaytet et al. (2018). One can then show that

$$\Gamma_k = \sqrt{\frac{9\gamma}{128\pi} \frac{Z_k e}{\rho c_s s_{\text{grain},k}^2} \min \left( \frac{M_0^{1/3}}{\mu \left( \frac{M_0}{\Phi} \right)_c} \left( \frac{4\pi\rho}{3} \right)^{2/3}, 0.1 \text{G} \right)} \quad (35)$$

Figure 21 shows the absolute value of the Hall factor as a function of the density and the grain size. For a wide range of grain sizes, the Hall factor is always much smaller than unity. These grains can therefore be considered neutral at least dynamically. We note that $\Gamma_k$ is larger than unity for very small grains ($s_{\text{grain},k} \lesssim 10^{-6}$ cm). These grains are very well coupled to the gas in the neutral case but, when charged, could experience a strong decoupling with the neutrals. This would occur if the latter were decoupled from the magnetic field, for example, in the nonideal regime. The Lorentz drag might play a crucial role
<figure><img src="image_22.png" /><figcaption>Fig. 21. Hall factor as a function of the density and the grain size for $\mu = 5$ for negatively charged grains. The dashed line denotes the equality between the gyration and the stopping time.</figcaption></figure>

for the dynamics of very small grains in star formation but is probably not very important for the large grains that we observe to decouple in our models. Efforts to study the dynamics of such grains have been made in the past (Guillet et al. 2007; Hopkins & Squire 2018) and should be extended to dusty col- lapses in future studies. We emphasize that the electromotive term in the Lorentz force applies only to the charged species and not the barycenter and might play a very important role in the decoupling between the charged grains and the neutrals (barycenter) even when $\Gamma_k < 1$.

### 7.5. Caveat: coagulation and fragmentation during the collapse

Dust coagulation has been neglected during this study. However, it may significantly affect dust evolution during the collapse as dust decoupling depends on grain size. Following Draine (1985), one can estimate the coagulation timescale $t_{\text{coag},i,j}$ between two dust phases i and j due to their relative motion within the collapsing clouds as

$$t_{\text{coag},i,j} = \left[ \frac{3}{4} \rho_i \frac{(s_{\text{grain},i} + s_{\text{grain},j})^2}{\rho_{\text{grain},i} s_{\text{grain},i}^3} |\mathbf{v}_i - \mathbf{v}_j| \right]^{-1}. \quad (36)$$

Assuming the same density for all the grains, neglecting cumu- lative back-reaction effects and magnetic fields, and considering an isothermal collapse, differential velocities can be estimated from the diffusion approximation as

$$|\mathbf{v}_i - \mathbf{v}_j| \sim \rho_{\text{grain}} |s_{\text{grain},i} - s_{\text{grain},j}| c_s \frac{|\nabla \rho|}{\rho^2}. \quad (37)$$

As said in Sect. 6, the density profile of the free-falling mate- rial can be approximated as a power law with an exponent $\zeta = -2$ (Larson 1969). In this case, one obtains at a distance $r$ from the central region

$$t_{\text{coag},i,j} = \left[ \frac{3|\zeta|}{4r} \epsilon_i (1 + q_{i,j})(1 - q_{i,j}^2) c_s \right]^{-1}, \quad (38)$$

where $q_{i,j} \equiv \frac{s_{\text{grain},i}}{s_{\text{grain},i}}$ is the ratio between the grain sizes. In the limit $q_{i,j} \ll 1$, we find

$$t_{\text{coag},i,j} \sim 170 \text{ kyr} \left( \frac{\epsilon_j}{10^{-2}} \right)^{-1} \left( \frac{r}{100 \text{ AU}} \right) \left( \frac{c_s}{0.19 \text{ km s}^{-1}} \right)^{-1}. \quad (39)$$

Growth induced by the relative dynamics between dust grains is therefore expected to be inefficient during protostellar col- lapse away from the core but could be non-negligible in the inner regions. At $r = 10$ AU, $t_{\text{coag},i,j} \sim 17$ kyr, which is slightly smaller than the typical free-fall timescale of a prestellar core.

We note that growth could also be enhanced by turbulence, the Brownian motions of dust grains, or focalization due to grain charges (Blum & Wurm 2008). In the case of Brownian motions, for example, assuming that the grain temperature is equal to the gas temperature, the differential velocity for two grains of different mass can be expressed as (Birnstiel et al. 2016)

$$|\mathbf{v}_i - \mathbf{v}_j| = \sqrt{\frac{8k_B T_g}{\pi}} \sqrt{\frac{m_{\text{grain},i} + m_{\text{grain},j}}{m_{\text{grain},i} m_{\text{grain},j}}. \quad (40)$$

Therefore, assuming again that $q_{i,j} \ll 1$ the coagulation timescale can be written as

$$t_{\text{coag},i,j} = \left[ \frac{3}{4} \rho_i \frac{1}{\rho_{\text{grain},i} s_{\text{grain},i}} \sqrt{\frac{8k_B T_g}{\pi m_{\text{grain},j}}} \right]^{-1}. \quad (41)$$

We now consider a region of density $10^{-12}$ g cm$^{-3}$ and temper- ature of 10 K. Assuming $s_{\text{grain},j} = 0.1$ µm and $s_{\text{grain},i} = 100$ µm, we get $t_{\text{coag},i,j} \sim 240$ kyr. We note that, in the case of the Brow- nian motions, the coagulation timescale depends on the grains size of both species. If we now consider $s_{\text{grain},j} = 0.01$ µm, $t_{\text{coag},i,j} \sim 7$ kyr.

In short, we are tempted to say that in the presence of large grains, very small grains could be efficiently removed during the collapse in high-density regions. Coagulation should there- fore be included in future studies. However, we admit that the presence of such large grains during the early phases of the protostellar collapse is still under debate. It is indeed unclear how these large dust grains can overcome the fragmentation barrier, as the differential velocities between two dust species is typically quite large; for example, up to a few $\approx 0.1$ km s$^{-1}$ in the envelope of MMMRN in the case of the two least cou- pled species. Typically, the velocity above which fragmentation can occur is thought to be about a few times $10$ m s$^{-1}$ (Blum & Wurm 2008) although it could be higher (Yamamoto et al. 2014). Large grains could overcome the fragmentation barrier because the fragmentation timescale is typically equivalent to the coagulation timescale (e.g., Gonzalez et al. 2017). Frag- mentation could simply not have the time to occur during the collapse, especially considering the fact that large grains quickly drift toward regions of high density where their drift velocity is typically around a few meters per second. Let us now esti- mate the typical fragmentation radius $r_{\text{frag}}$, which corresponds to the distance from the center at which the fragmentation veloc- ity $v_{\text{frag}}$ is equal to the dust drift velocity; it can be obtained by solving

$$\rho_{\text{grain}} s_{\text{grain}} c_s \frac{|\nabla \rho|}{\rho^2} = v_{\text{frag}}. \quad (42)$$

Let us now assume first core formation with a core profile of the type $\rho = \frac{\rho_{\text{F}}}{1 + \left(\frac{\rho_{\text{F}}}{r_{\text{F}}}\right)^2}$. We can show that

$$r_{\text{frag}} = \frac{1}{2} \frac{r_{\text{F}}^2}{s_{\text{grain}}} \frac{\rho_{\text{F}}}{\rho_{\text{grain}}} \frac{v_{\text{frag}}}{c_s}. \quad (43)$$
Assuming $\rho_{\mathcal{F}} = 10^{-11} \text{ g cm}^{-3}$, $r_{\mathcal{F}} = 10$ AU, $\rho_{\text{grain}} = 1 \text{ g cm}^{-3}$, $s_{\text{grain}} = 160 \text{ µm}$ and $T_g = 10 \text{ K}$

$$r_{\text{frag}} = 2.5 \times 10^4 \text{ AU} \left(\frac{r_{\mathcal{F}}}{10 \text{ AU}}\right)^2 \left(\frac{\rho_{\mathcal{F}}}{10^{-11} \text{ g cm}^{-3}}\right) \left(\frac{c_s}{0.19 \text{ km s}^{-1}}\right) \left(\frac{v_{\text{frag}}}{10 \text{ m s}^{-1}}\right) \left(\frac{s_{\text{grain}}}{160 \text{ µm}}\right)^{-1} \left(\frac{\rho_{\text{grain}}}{1 \text{ g cm}^{-3}}\right)^{-1} \quad (44)$$

The fragmentation radius is significantly larger than the scales that are involved during the first protostellar collapse. We note a quadratic dependency with respect to the first core radius. For a core of size 1 AU, one finds $r_{\text{frag}} = 2.5 \times 10^2 \text{ AU}$. Hence, fragmentation may be relevant for small cores and should be investigated in further studies.

## 8. Conclusion and perspective

In this study, we present the first MULTIGRAIN and nonideal MHD simulations of dusty protostellar collapses using the new dust dynamics solver of RAMSES. We present six DUSTYCOLLAPSE simulations with a simultaneous treatment of ten dust species with different sizes. In these simulations, we investigate the impact of the maximum grain size, the thermal-to-gravitational energy ratio, and the presence of magnetic fields on the dynamics of the DUSTYCOLLAPSE. We summarize below our principal findings:

1. Small grains with sizes of less than a few tens of microns are strongly coupled to the gas during the protostellar collapse. On the contrary, grains larger than ~100 µm tend to decouple significantly.

2. When the first hydrostatic core forms, high-density regions – the core, the fragments, the disk, and the pseudo-disk – are enriched in dust by a typical factor of two, whereas low-density regions – the envelope and the outflow – are strongly dust depleted.

3. Dust is not necessarily a proxy for gas during the collapse. Inferring gas densities from dust is found to potentially lead to extremely large errors (up to ~250%).

4. A standard MRN grain-size distribution with a maximum grain size of 250 nm is however extremely well preserved during the protostellar collapse in the absence of coagulation.

5. Dust dynamics is strongly affected by the initial cloud properties. Variations of the dust-to-gas ratio are greatest when the free-fall is long and the initial density is low. An additional decoupling occurs for neutral grains in the presence of magnetic fields because collapse proceeds over longer timescales.

6. With a semi-analytical model, we show that the dust-ratio varies exponentially with the initial Stokes number during the collapse. More precisely, we show that it can be expressed as $\epsilon_0 \chi^{\text{St}_0}$ where $\chi$ is a dimensionless function of time and is independent of the dust properties. Fitting the values of $\chi$ leads to very good agreement between the semi-analytical model and the measured dust ratios in the range of Stokes numbers considered in our simulations.

7. Using the calibration of $\chi$ with the results of our model, we show that a Stokes number of at least 0.01 is required to enrich the core and the disk in a dust species by a factor of two. Similarly we show that grains with $\text{St}_0 \gtrsim 0.027$ can potentially be depleted by a factor of two in the envelope after the first core formation.

Dust evolution during the protostellar collapse could have serious consequences on the initial state of protoplanetary disks and the further formation of the planets. In the future, substantial efforts should be made to include the dynamics of charged dust grains during the protostellar collapse because the Lorentz drag cannot necessarily be neglected for small grains. Coagulation and fragmentation of dust grains should also be considered to investigate dust evolution during the star formation process more realistically.

**Acknowledgements.** First, we thank the referee for providing useful comments and advice that helped us to ameliorate our manuscript. We acknowledge financial support from “Programme National de Physique Stellaire” (PNPS) of CNRS/INSU, CEA and CNES, France. This work was granted access to the HPC resources of CINES (Occigen) under the allocation DARI A0020407247 made by GENCI. Computations were also performed at the Common Computing Facility (CCF) of the LABEX Lyon Institute of Origins (ANR-10-LABX-0066). This work took part under the programs ISM3D and Core2disk of the PSI2 project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02. This project was partly supported by the IDEXLyon project (contract n ANR-16-IDEX-0005) under University of Lyon auspices. The plots were generated using the very efficient Osiris library developed by Neil Vaytet, Tommaso Grassi and Matthias González whom we thank. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 823823. We also thank Etienne Jaupart for useful discussions on the theoretical model for the dust enrichment during the collapse.

## References

Allen, A., Shu, F. H., & Li, Z.-Y. 2003, ApJ, 599, 351

Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C. P. 2010, ApJ, 723, 1241

Bate, M. R., & Lorén-Aguilar, P. 2017, MNRAS, 465, 1089

Belloche, A. 2013, EAS Pub. Ser., 62, 25

Berger, M. J., & Oliger, J. 1984, J. Comput. Phys., 53, 484

Birnstiel, T., Dullemond, C. P., & Brauer, F. 2009, A&A, 503, L5

Birnstiel, T., Fang, M., & Johansen, A. 2016, Space Sci. Rev., 205, 41

Blum, J., & Wurm, G. 2008, ARA&A, 46, 21

Bodenheimer, P. 1995, ARA&A, 33, 199

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

Brauer, F., Dullemond, C. P., & Henning, T. 2008, A&A, 480, 859

Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371

Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3

Commerçon, B., Hennebelle, P., & Henning, T. 2011a, ApJ, 742, L9

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

Commerçon, B., Debout, V., & Teyssier, R. 2014, A&A, 563, A11

Dipierro, G., & Laibe, G. 2017, MNRAS, 469, 1932

Dipierro, G., Laibe, G., Alexander, R., & Hutchison, M. 2018, MNRAS, 479, 4187

Draine, B. T. 1985, in Protostars and Planets II, eds. D. C. Black, & M. S. Matthews (Tucson, AZ: University of Arizona Press), 621

Draine, B. T. 2004, in The Cold Universe (Berlin: Springer), 213

Drążkowska, J., Alibert, Y., & Moore, B. 2016, A&A, 594, A105

Dubois, Y., & Commerçon, B. 2016, A&A, 585, A138

Dubois, Y., Commerçon, B., Marcowith, A. r., & Brahimi, L. 2019, A&A, 631, A121

Elbakyan, V. G., Johansen, A., Lambrechts, M., Akimkin, V., & Vorobyov, E. I. 2020, A&A, 637, A5

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

Fromang, S., & Papaloizou, J. 2006, A&A, 452, 751

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

Fujii, Y. I., Okuzumi, S., & Inutsuka, S.-i. 2011, ApJ, 743, 53

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

Galli, D., & Shu, F. H. 1993, ApJ, 417, 220

Godunov, S. K. 1959, Mat. Sb., 47, 271

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

González, J.-F., Laibe, G., & Maddison, S. T. 2017, MNRAS, 467, 1984

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

Guillet, V., Pineau Des Forêts, G., & Jones, A. P. 2007, A&A, 476, 263

Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9
Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8

Hennebelle, P., Commerçon, B., Lee, Y.-N., & Charnoz, S. 2020, A&A, 635, A67

Hincelin, U., Commerçon, B., Wakelam, V., et al. 2016, ApJ, 822, 12

Hopkins, P. F., & Squire, J. 2018, MNRAS, 479, 4681

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

Jaupart, E., & Laibe, G. 2020, MNRAS, 492, 4591

Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, Nature, 448, 1022

Johansen, A., Youdin, A., & Mac Low, M.-M. 2009, ApJ, 704, L75

Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128

Kanagawa, K. D., Ueda, T., Muto, T., & Okuzumi, S. 2017, ApJ, 844, 142

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

Kataoka, A., Tsukagoshi, T., Momose, M., et al. 2016, ApJ, 831, L12

Keith, S. L., & Wardle, M. 2014, MNRAS, 440, 89

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

Kunz, M. W., & Mouschovias, T. C. 2009, ApJ, 693, 1895

Kwok, S. 1975, ApJ, 198, 583

Laibe, G. 2014, MNRAS, 437, 3037

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

Laibe, G., & Price, D. J. 2014b, MNRAS, 440, 2147

Laibe, G., & Price, D. J. 2014c, MNRAS, 444, 1940

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

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

Lin, M.-K. 2019, MNRAS, 485, 5221

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

Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2014, MNRAS, 438, 2278

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

Marchand, P., Commerçon, B., & Chabrier, G. 2018, A&A, 619, A37

Marchand, P., Tomida, K., Commerçon, B., & Chabrier, G. 2019, A&A, 631, A66

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

Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commerçon, B. 2016, A&A, 587, A32

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

McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565

Mignon-Risse, R., González, M., Commerçon, B., & Rosdahl, J. 2020, A&A, 635, A42

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

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

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

Pohl, A., Kataoka, A., Pinilla, P., et al. 2016, A&A, 593, A12

Price, D. J., & Bate, M. R. 2007, Ap&SS, 311, 75

Price, D. J., & Laibe, G. 2015, MNRAS, 454, 2320

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

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

Rosdahl, J., Blaizot, J., Aubert, D., Stranex, T., & Teyssier, R. 2013, MNRAS, 436, 2188

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018a, ApJ, 859, 165

Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018b, ApJ, 869, 115

Sadavoy, S. I., Stephens, I. W., Myers, P. C., et al. 2019, ApJS, 245, 2

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

Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44

Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179

Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19

Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405

Umebayashi, T., & Nakano, T. 1990, MNRAS, 243, 103

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

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

Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5

Vorobyov, E. I., & Elbakyan, V. G. 2019, A&A, 631, A1

Vorobyov, E. I., Skliarevskii, A. M., Elbakyan, V. G., et al. 2019, A&A, 627, A154

Whitworth, A. P., & Clarke, C. J. 1997, MNRAS, 291, 578

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

Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719

Yamamoto, T., Kadono, T., & Wada, K. 2014, ApJ, 783, L36

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050
# Appendix A: Impact of velocity regularization

Total dust-ratio $\epsilon$

<figure><img src="image_23.png" /><figcaption>**Fig. A.1.** Fiducial model for various maximum dust differential velocity $w_{cap}$. Mid-plane view of the total dust ratio at the FHSC formation. (Top-left) $w_{cap} = 0.1$ km s$^{-1}$, (top-right) $w_{cap} = 0.5$ km s$^{-1}$, (bottom-left) $w_{cap} = 1$ km s$^{-1}$, (bottom-right) $w_{cap} = 2$ km s$^{-1}$.</figcaption></figure>

As explained in Sect. 3.5, the simulations of our fiducial model were performed with various values for the maximum gas and dust differential velocity $w_{cap}$. We investigated the effect of varying $w_{cap}$ by performing complementary simulations with $w_{cap} = 0.1, 0.5, 1$ and $2$ km s$^{-1}$. We attempted to simulate an additional model with $w_{cap} = 10$ km s$^{-1}$ but this led to numerical instabilities due to unrealistically large dust velocities at the accretion shock, where the diffusion approximation is not valid.

In Fig. A.1, we show a face-one view of the dust ratio at the time of the FHSC formation for these four models. With $w_{cap} = 0.5-2$ km s$^{-1}$ we generally find the same results as in the fiducial case that has $w_{cap} = 1$ km s$^{-1}$. However, having $w_{cap} = 0.1$ km s$^{-1}$ appears to be an overly extreme choice; it indeed suppresses most of the initial dust enrichment that is due to the decoupling between the gas and the dust in the envelope. Having $w_{cap} = 1$ km s$^{-1}$ appeared to be a reasonable choice to cope with time stepping and physical constraints.

# Appendix B: Nonlinear dust enrichment: neglecting back-reaction

As stated in Sect. 6, Eq. (25) is valid as long as local dust-ratio variations and the cumulative back-reaction of the dust can be neglected. In this approximation, the dust ratio enrichment increases exponentially with the Stokes number. As pointed out in Sect. 6, we observe deviations from the exponential law in the disk and the envelope in MMMRN, while it is very well verified in the core.

To ascertain whether these deviations are due to back-reaction or not, we perform an additional simulation similar to the fiducial one, but with $\theta_{d,0} = 10^{-7}$ (MMMRNEPS1E-7). For MMMRNEPS1E-7, we expect the back-reaction to be definitely negligible as $\theta_{d,0} \ll 1$.

Figure B.1 shows the dust ratio enrichment as a function of the initial Stokes number for MMMRN (red) and MMMRNEPS1E-7 (pink) in the core, disk, and envelope at $t_{core} + 2$ kyr. On top of this, the dust ratio in the disk but at $t_{core}$ is plotted

<figure><img src="image_24.png" /><figcaption>**Fig. B.1.** Semi-analytic and measured dust-ratio enrichment in MMMRN (red) and MMMRNEPS1E-7 (pink) in the core (stars), disk (triangles) and envelope (circles) at $t_{core} + 2$ kyr against the initial Stokes number. The lines represent the semi-analytical development using a best fit of $\chi$ for the FHSC (dotted), disk (solid), and envelope (dashed), respectively. We display the same information for MMMRN at $t_{core}$ in the disk (gray).</figcaption></figure>

for MMMRN (gray). As can be seen, the differences between the models at $t_{core} + 2$ kyr have similar dust-ratio enrichment. This indicates that the cumulative back-reaction of dust is negligible when $St_{k,0} < 10^{-2}$ and only corrective above that. We are also confident that the mixing between dust-enriched and dust-depleted content, that is, a non-negligible $\nabla\epsilon$, are the source of discrepancy between a pure exponential enrichment and the values measured in our models in the disk. We indeed observe very good agreement with the exponential law for MMMRN at $t_{core}$. At this time the disk has just formed from the dense and dust-enriched material of the inner envelope. Later at $t_{core} + 2$ kyr, the disk has been accreting material from the envelope that is significantly depleted in large grains which causes a diminution of the average dust enrichment.

# Appendix C: Charges on dust grains

One can crudely estimate the charges on the grains during the collapse, similarly to Wurster et al. (2016), to compare the relative intensities of the Lorentz and the Epstein drags. For simplicity, we assume that all the grains have the same negative charge $Z_d$ (for all $k, Z_k = Z_d$) and the ions to have a charge $Z_{ions} = 1$. We consider 100 dust bins distributed as an MRN extended to millimeter-sized grains with a dust-to-gas ratio of 1%. Local electroneutrality ensures that

$$n_{\text{ions}} - n_e + Z_d n_d = 0, \qquad (C.1)$$

$n_{\text{ions}}$, $n_e$ and $n_d$ being the number density of the ions, electrons and dust, the latter being given by

$$n_d \equiv \frac{\mu_g m_H}{m_{\text{grain}}} \epsilon n, \qquad (C.2)$$

where $m_{\text{grain}}$ is the average mass of a dust grain.

Two additional constraints are provided by the evolution equations for the ion and electron charge density (Umebayashi & Nakano 1980; Fujii et al. 2011). Assuming steady state and only considering charge capture by the grains (Keith & Wardle 2014), these equations are

$$\begin{aligned} n_{\text{ions}} &= \frac{\zeta n}{k_{\text{ions,grains}} n_d}, \\ n_e &= \frac{\zeta n}{k_{\text{e,grains}} n_d}, \end{aligned} \qquad (C.3)$$
where $\zeta$ is the cosmic-ray ionization rate. Similarly to Wurster et al. (2016), we adopt the typical value $\zeta = 10^{-17} s^{-1}$. We also express $k_{\text{ions,grains}}$ and $k_{\text{e,grains}}$, the charge capture rates on neutral grains, as in (Fujii et al. 2011) as

$$k_{\text{ions,grains}} \equiv \pi s_{\text{grain}}^2 \sqrt{\frac{8k_b T_g}{\pi m_{\text{ions}}}} \left(1 - \frac{e^2 Z_d}{s_{\text{grain}} k_b T_g}\right), \\ k_{\text{e,grains}} \equiv \pi s_{\text{grain}}^2 \sqrt{\frac{8k_b T_g}{\pi m_e}} \exp\left[\frac{e^2 Z_d}{s_{\text{grain}} k_b T_g}\right], \quad (C.4)$$

where $s_{\text{grain}}$ is the average grain size, $m_e$ is the electron mass, and $m_{\text{ions}}$ the ions mass, assumed to be 24.3 proton masses. Using the previous equations, one obtains

$$Z_d = \frac{\zeta}{n} \left( \frac{m_{\text{grain}}}{\epsilon \mu_g m_{\text{H}}} \right)^2 \left[ \frac{1}{k_{\text{ions,grains}}(Z_d)} - \frac{1}{k_{\text{e,grains}}(Z_d)} \right], \quad (C.5)$$

which we invert using the Newton-Raphson method to get $Z_d$.

## Appendix D: Distributions

In Table D.1, we provide the initial dust distributions, rounded quantities for the three maximum grain sizes used in our models (the exact calculation can be made using the method presented in Sect. 3.3). Table D.2 shows $\langle \Theta_{d,k} \rangle_m$ (in %), the dust-to-gas ratio averaged in mass, along with its corresponding dust enrichment, and the gas mass (in units of $M_{\odot}$) for all runs $t_{\text{core}} + 2$ kyr and all the different objects.

<table><caption>**Table D.1.** Initial dust distributions, rounded quantities for the three maximum grain sizes used in our models (the exact calculation can be made using the method presented in Sect. 3.3).</caption><thead><tr><th>S<sub>max</sub></th><th>s (cm), &theta;:</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td rowspan="2">2.5(-5)</td><td>s (cm)</td><td>4.11(-7)</td><td>6.08(-7)</td><td>8.99(-7)</td><td>1.33(-6)</td><td>1.97(-6)</td><td>2.91(-6)</td><td>4.3(-6)</td><td>6.36(-6)</td><td>9.4(-6)</td><td>1.39(-5)</td></tr><tr><td>&theta;</td><td>3.56(-4)</td><td>4.33(-4)</td><td>5.26(-4)</td><td>6.40(-4)</td><td>7.78(-4)</td><td>9.46(-4)</td><td>1.15(-3)</td><td>1.4(-3)</td><td>1.7(-3)</td><td>2.1(-3)</td></tr><tr><td rowspan="2">0.01</td><td>s (cm)</td><td>3.05(-7)</td><td>8.2(-7)</td><td>2.21(-6)</td><td>5.95(-6)</td><td>1.6(-5)</td><td>4.31(-5)</td><td>1.16(-4)</td><td>3.12(-4)</td><td>8.41(-4)</td><td>2.26(-3)</td></tr><tr><td>&theta;</td><td>4.56(-5)</td><td>7.49(-5)</td><td>1.223(-4)</td><td>2.02(-4)</td><td>3.31(-4)</td><td>5.43(-4)</td><td>8.9(-4)</td><td>1.46(-3)</td><td>2.40(-3)</td><td>3.93(-3)</td></tr><tr><td rowspan="2">0.1</td><td>s (cm)</td><td>2.72(-7)</td><td>9.2(-7)</td><td>3.12(-6)</td><td>1.05(-5)</td><td>3.58(-5)</td><td>1.21(-4)</td><td>4.12(-4)</td><td>1.4(-3)</td><td>4.73(-3)</td><td>1.6(-2)</td></tr><tr><td>&theta;</td><td>1.88(-5)</td><td>3.47(-5)</td><td>6.39(-5)</td><td>1.18(-4)</td><td>2.16(-4)</td><td>3.99(-4)</td><td>7.34(-4)</td><td>1.35(-3)</td><td>2.49(-3)</td><td>4.58(-3)</td></tr></tbody></table>

**Notes.** For each $S_{\text{max}}$, the first line represents the grain size in cm and the second the initial dust-to-gas ratio. The exponents are given by the parenthesis and the integers from 1 to 10 correspond to bins of increasing sizes.

<table><caption>**Table D.2.** Mass-averaged dust-to-gas ratio $\langle \Theta_{d,k} \rangle_m$ (in %) and gas mass (in units of $M_{\odot}$) for all runs $t_{\text{core}} + 2$ kyr and all the objects;</caption><thead><tr><th>Model</th><th>j</th><th>⟨Θ<sub>d,k</sub>≤6⟩<sub>m</sub></th><th>⟨Θ<sub>d,7</sub>⟩<sub>m</sub></th><th>⟨Θ<sub>d,8</sub>⟩<sub>m</sub></th><th>⟨Θ<sub>d,9</sub>⟩<sub>m</sub></th><th>⟨Θ<sub>d,10</sub>⟩<sub>m</sub></th><th>m<sub>g</sub></th></tr></thead><tbody><tr><td rowspan="5">MMMRN</td><td>F<sub>0</sub></td><td>0.085 [1]</td><td>0.075 [1.03]</td><td>0.14 [1.03]</td><td>0.31 [1.24]</td><td>1.17 [2.6]</td><td>0.082</td></tr><tr><td>F<sub>1</sub></td><td>0.085 [1]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.37 [1.54]</td><td>1.38 [3]</td><td>0.024</td></tr><tr><td>F<sub>2</sub></td><td>0.085 [1]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.37 [1.54]</td><td>1.37 [2.98]</td><td>0.024</td></tr><tr><td>D</td><td>0.085 [1]</td><td>0.075 [1.03]</td><td>0.15 [1.1]</td><td>0.34 [1.41]</td><td>1.1 [2.4]</td><td>0.096</td></tr><tr><td>E</td><td>0.085 [1]</td><td>0.073 [1]</td><td>0.13 [0.96]</td><td>0.23 [0.95]</td><td>0.34 [0.73]</td><td>1.0</td></tr><tr><td rowspan="6">MRN</td><td>F<sub>0</sub></td><td>0.37 [1]</td><td>0.12 [1]</td><td>0.14 [1]</td><td>0.17 [1]</td><td>0.21 [1]</td><td>0.088</td></tr><tr><td>F<sub>1</sub></td><td>0.37 [1]</td><td>0.12 [1]</td><td>0.14 [1]</td><td>0.17 [1]</td><td>0.21 [1]</td><td>0.033</td></tr><tr><td>F<sub>2</sub></td><td>0.37 [1]</td><td>0.12 [1]</td><td>0.14 [1]</td><td>0.17 [1]</td><td>0.21 [1]</td><td>0.033</td></tr><tr><td>D</td><td>0.37 [1]</td><td>0.12 [1]</td><td>0.14 [1]</td><td>0.17 [1]</td><td>0.21 [1]</td><td>0.091</td></tr><tr><td>E</td><td>0.37 [1]</td><td>0.12 [1]</td><td>0.14 [1]</td><td>0.17 [1]</td><td>0.21 [1]</td><td>0.98</td></tr><tr><td>F<sub>0</sub></td><td>0.13 [1]</td><td>0.09 [1]</td><td>0.15 [1]</td><td>0.25 [1.04]</td><td>0.43 [1.1]</td><td>0.086</td></tr><tr><td rowspan="5">100MICMRN</td><td>F<sub>1</sub></td><td>0.13 [1]</td><td>0.09 [1]</td><td>0.15 [1]</td><td>0.25 [1.04]</td><td>0.46 [1.17]</td><td>0.023</td></tr><tr><td>F<sub>2</sub></td><td>0.13 [1]</td><td>0.09 [1]</td><td>0.15 [1]</td><td>0.25 [1.04]</td><td>0.46 [1.17]</td><td>0.023</td></tr><tr><td>D</td><td>0.13 [1]</td><td>0.09 [1]</td><td>0.15 [1]</td><td>0.25 [1.04]</td><td>0.44 [1.13]</td><td>0.097</td></tr><tr><td>E</td><td>0.13 [1]</td><td>0.09 [1]</td><td>0.15 [1]</td><td>0.24 [1]</td><td>0.38 [0.97]</td><td>0.99</td></tr><tr><td>F<sub>0</sub></td><td>0.085 [1]</td><td>0.074 [1.01]</td><td>0.14 [1.03]</td><td>0.26 [1.04]</td><td>0.51 [1.12]</td><td>0.20</td></tr><tr><td rowspan="5">MMMRNA0.25</td><td>F<sub>1</sub></td><td>0.085 [1]</td><td>0.074 [1.01]</td><td>0.14 [1.03]</td><td>0.27 [1.08]</td><td>0.58 [1.29]</td><td>0.10</td></tr><tr><td>F<sub>2</sub></td><td>0.085 [1]</td><td>0.074 [1.01]</td><td>0.14 [1.03]</td><td>0.27 [1.08]</td><td>0.58 [1.29]</td><td>0.10</td></tr><tr><td>D</td><td>0.085 [1]</td><td>0.074 [1.01]</td><td>0.14 [1.03]</td><td>0.26 [1.04]</td><td>0.56 [1.24]</td><td>0.16</td></tr><tr><td>E</td><td>0.085 [1]</td><td>0.073 [1]</td><td>0.13 [0.96]</td><td>0.24 [0.96]</td><td>0.4 [0.89]</td><td>0.7</td></tr><tr><td rowspan="5">MMMRNMHD</td><td>F<sub>0</sub></td><td>0.085 [1]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.37 [1.48]</td><td>1.5 [3.33]</td><td>0.075</td></tr><tr><td>D</td><td>0.085 [1]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.39 [1.56]</td><td>1.58 [3.58]</td><td>0.0012</td></tr><tr><td>P</td><td>0.086 [1.01]</td><td>0.077 [1.05]</td><td>0.16 [1.19]</td><td>0.43 [1.72]</td><td>1.23 [2.73]</td><td>0.1</td></tr><tr><td>O</td><td>0.085 [1]</td><td>0.074 [1.01]</td><td>0.14 [1.03]</td><td>0.26 [1.04]</td><td>0.1 [0.22]</td><td>0.0082</td></tr><tr><td>E</td><td>0.085 [1]</td><td>0.072 [0.98]</td><td>0.13 [0.96]</td><td>0.2 [0.8]</td><td>0.17 [0.38]</td><td>0.5</td></tr><tr><td rowspan="5">MMMRNNIMHD</td><td>F<sub>0</sub></td><td>0.085 [1]</td><td>0.077 [1.5]</td><td>0.15 [1.1]</td><td>0.37 [1.48]</td><td>1.5 [3.33]</td><td>0.081</td></tr><tr><td>D</td><td>0.085 [1]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.37 [1.48]</td><td>1.44 [3.2]</td><td>0.011</td></tr><tr><td>P</td><td>0.086 [1.01]</td><td>0.077 [1.05]</td><td>0.16 [1.19]</td><td>0.43 [1.72]</td><td>1.22 [2.71]</td><td>0.1</td></tr><tr><td>O</td><td>0.086 [1.01]</td><td>0.076 [1.04]</td><td>0.15 [1.1]</td><td>0.32 [1.28]</td><td>0.2 [0.44]</td><td>0.0018</td></tr><tr><td>E</td><td>0.085 [1]</td><td>0.072 [0.98]</td><td>0.13 [0.96]</td><td>0.2 [0.8]</td><td>0.17 [0.38]</td><td>0.5</td></tr></tbody></table>

**Notes.** The number inside the brackets is the corresponding mass averaged dust enrichment. When $\langle \Theta_d \rangle_m > 1$, the enrichment is referenced in red while it is referenced in blue when $\langle \Theta_d \rangle_m \le 1$. We denote the different objects as follows: $F_j$ represent the fragments (F0 being the FHSC), $D$ the disks, $O$ the outflow, $P$ the pseudo-disks, and $E$ the envelope.