Protostellar collapse: the conditions to form dust rich protoplanetary disks

Dust plays a key role during star, disk and planet formation. Yet, its dynamics during the protostellar collapse remains a poorly investigated field. Recent studies seem to indicate that dust may decouple efficiently from the gas during these early stages. We aim to understand how much and in which regions dust grains concentrate during the early phases of the protostellar collapse, and see how it depends on the properties of the initial cloud and of the solid particles. We use the multiple species dust dynamics solver of the grid-based code RAMSES to perform various simulations of dusty collapses. We perform hydrodynamical and MHD simulations where we vary the maximum grain size, the thermal-to-gravitational energy ratio and the magnetic properties of the cloud. We simulate the simultaneous evolution of ten neutral dust grains species with grain sizes varying from a few nm to a few hundredth of microns. We obtain a significant decoupling between the gas and the dust for grains of typical sizes a few 10 microns. This decoupling strongly depends on the thermal-to-gravitational energy ratio, the grain sizes or 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. 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 and is depleted in low density regions. The size at which grains decouple from the gas depends on the initial properties of the clouds. Since dust can not necessarily be used as a proxy for gas during the collapse, we emphasize on the necessity of including the treatment of its dynamics in collapse simulations.


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 2 at present days (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 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 Hennebelle et al. 2020) and the fragmentation process .
Until recently, one paradigm was that dust of the interstellar medium (ISM) is usually composed of grains with sizes up to ∼ 0.1 µm with a typical size distribution well modelled by the Mathis-Rumpl-Nordsieck distribution (MRN, 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. 2015Pohl et al. 2016;Sadavoy et al. 2018aSadavoy et al. ,b, 2019Valdivia et al. 2019). Galametz et al. (2019) has also shown 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 longstanding problem in star formation. Indeed, the specific angular momentum of prestellar cores differs to those of young stars by more than three orders of magnitude (Bodenheimer 1995;Belloche 2013). In numerous studies, the 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;Masson et al. 2016). State-of-the-art simulations account now for the effect of magnetic fields both in a 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 investigated the dynamics of dust during the 3D protostellar collapse (dustycollapse). They report that ∼ 100 µm grains can significantly decouple from the gas lead- ing to large increase of dust-to-gas ratio in the disk and the first Larson core. In 2D,  have studied the gas and dust decoupling in gravitoviscous protoplanetry disks including dust growth and also report strong variations of dust-to-gas ratio. So far, no 3D dustycollapse simulation has been performed in a MHD/non ideal MHD context or with multiple dust species.
Multifluid-multi-dust species simulations are prohibitive as they require to solve two additional equations per dust species (mass and momentum conservation). Recent theoretical (Laibe & Price 2014a,b) and numerical  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 that 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 organised 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 have 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.

Dusty hydrodynamics for the protostellar collapse
A gas and dust mixture with N small grains 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 ρ, flows at its barycenter velocity v. The k-th dust phase, of density ρ 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 ∂ρ ∂t + ∇ · ρv = 0, ∂ρ k ∂t + ∇ · ρ k (v + w k ) = 0, ∀k ∈ [1, N] , ∂ρv ∂t where P g is the thermal pressure of the gas and I the identity matrix. The gravitational potential φ is set by the Poisson equation where 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 ρ 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 The gas density ρ g is ρ g = ρ − k ρ k . Regions of low densities are isothermal and have for sound speed 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 ρ 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 where ρ is the total density of the gas and dust mixture, c s is the sound speed of the gas and γ its adiabatic index. If the differential velocity ∆v k ≡ 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) where M d = |∆v k | c s is the differential velocity Mach number. In the remaining of this paper, unless specified, we consider this correction.
In the terminal velocity approximation, the differential velocity of the phase k is and the gas and dust velocities, v g and v k are given by For later purposes, we define the dust ratio k ≡ ρ k ρ , the total dust ratio ≡ N k k , and the dust-to-gas ratio θ d ≡ N k ρ k ρ g . For any quantity A, we defineĀ ≡ A/A 0 where A 0 is its initial value.¯ andθ are called 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). 1 distinct from the dust density ρ k 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, i.e. 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) write where B is the magnetic field, η A is the ambipolar resistivity and c is the speed of light. We note that η A c 2 4π|B| 2 [(∇ × B) × B] is the differential velocity between the ions and the neutrals in the gas phase (Masson et al. 2012. We then must correct the neutrals velocity by accounting for the differential velocity between the gas and the barycenter. The term − k ρ k ρ−ρ k w k appears in the induction equation to take that into account. In addition from the pressure force, a magnetic force now applies on the plasma. In a previous work, Laibe & Price (2014c) have given the expression of the differential velocity for a general force budget. Using this formula we find that we note that this expression is very similar to what Fromang & Papaloizou (2006) found for single dust species mixtures. To further simplify Eqs (8), we assume in this paper that the plasma velocity v is the barycenter velocity which is valid when k ||w k || ||w k || |v|. In this case the induction equation writes

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 proficient for problems that require a treatment of magnetohydrodynamics Masson et al. 2012;Marchand et al. 2018;Marchand et al. 2019), radiation hydrodynamics Rosdahl et al. 2013;Commerçon et al. 2014;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).

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 0 , the temperature of the gas T g and a mean molecular weight µ g . The ratio between the thermal and the gravitational energy α is and sets the initial radius of the cloud R 0 and its density ρ 0 . In addition, we impose an initial solid body rotation around the zaxis at the angular velocity Ω 0 by setting the ratio between the rotational and the gravitational energy β given by Eventually, we apply an initial azimuthal density perturbation according to 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 an uniform magnetic field using the mass-toflux-to-critical-mass-to-flux-ratio the critical mass-to-flux ratio being given by (M/Φ) c = 0.53 3π √ 5/G (Mouschovias & Spitzer 1976). We set an angle φ mag between the magnetic fields and the rotation axis to reduce the efficiency of the magnetic braking.

Dust grain size distributions
In our multigrain simulations, N > 1 dust bins are considered. The dust ratio of each bins is set from power-law distributions with 0 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 The typical grain size of a bin k required to compute the stopping time is 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

Cloud setup
We impose initial conditions that are typical of the first protostellar collapse (Larson 1969) with T g = 10 K, µ g = 2.31 and a solar mass cloud. We also set γ = 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 Eqs. 1 with ρ 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 10 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 ρ grain = 1 g cm −3 . Finally, we impose an uniform initial total dust-to-gas ratio of θ d,0 = 0.01 in all the models.
The two magnetic models have been computed with µ = 5 and φ mag = 40 • . For the non-ideal MHD model the ambipolar resistivity is computed similarly to the case of reference of Marchand et al. (2016).

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 4 point per Jeans length, Truelove et al. 1997) must be satisfied to avoid artificial clump formation. We therefore enforce a refinement criteria that imposes at least 15 points per local Jeans length. The grid is initialized to the level min = 5 and allows refinement up to a maximum level max = 16 (which gives a resolution between 32 3 and 65536 3 cells).

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/the fragments are any object of density larger than 10 −12.5 g cm −3 , as in Lebreuilly et al. (2019). The FHSC or F 0 corresponds to the central fragment. F 1 and F 2 are the secondary fragments/FHSC. -The disk D is the region that satisfies Joos et al. (2012) criterion. For the analysis, we place ourselves in cylindrical coordinates (r, φ, z). A region is identified as a disk if it is Keplerian (v φ > f thre v r ), in hydrostatic equilibrium (v φ > f thre v z ), rotationally supported ( 1 2 ρv 2 φ > f thre P g ) and dense ρ > 3.9 × 10 −15 g cm −3 . As in Joos et al. (2012) , we choose f thre = 2.
-The pseudo-disk P (Galli & Shu 1993), for magnetic runs, is defined as the regions with r < 5000 AU and densities above 3.9 × 10 −17 gcm −3 that are not in the disk and fragments. The criterion is similar to Hincelin et al. (2016). 5000 AU is an arbitrary distance that is sufficiently larger than the central objects while being smaller than the initial cloud. -Jets/outflows O correspond to any region with r < 5000 AU with v · r |r| > 0.2 km sec −1 . The criterion is also similar as Hincelin et al. (2016).
-The envelope E encompasses the regions with r < 5000 AU that exclude the fragments, the disk/pseudo-disk and the jets/outflows.
We consider two different weights for averaging a quantity A over a volume V in the computational box. Volume averaging is computed according to Mass averaging is computed according to where ρ i and ∆x i are the total density and length of individual cells i in the averaged volume. Volume averages emphasize on regions of large spatial extension, i.e. the envelope, while mass averages emphasize on regions of high density, i.e. the core+disk and the denser regions of the envelope.

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 timesteps 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 timestep to estimate the differential velocity mach number.
Finally, we set the drift velocity to zero at densities lower than the ones of the initial cloud, i.e. the background. This is a way to ensure that the regions where the terminal velocity is not valid do not affect significantly the calculation.

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 freefall 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 s grain 0.05 cm The dynamics of grains smaller than 0.05 cm can be simulated using the diffusion approximation. Note that the Stokes number varies as ∝ 1 √ ρ (since t ff ∝ 1 √ ρ and t s ∝ 1 ρ ) and hence can increases with a decreasing density. We show in Fig 1 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 external regions of the collapse for all our models and it 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 α decreases, implying that the diffusion approximation remains valid for small values of α.
In the induction equation, we consider for simplicity that the plasma is moving at the barycenter velocity. This approximation is valid when 1, i.e. when the back-reaction from the dust onto the gas is negligible. Note that we investigate the impact of the back-reaction in Appendix B.

Models
The models presented in this section are referenced in Table 1. All of them have been evolved up to at least 2 kyr after the formation of the first core. Our fiducial case mmMRN has been run over a longer time.

Fiducial simulation
In this section, we present our fiducial case mmMRN where α = 0.5 and β = 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 ∼ 160 µm. The total initial dust-to-gas ratio is θ 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 kyr 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 these strong enrichment in the mid-plane and depletion in the envelope have already occurred at t core + 1 kyr, and continues 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 enrich significantly the fragments and the disk in large grains and the dust-to-gas ratio in the disk even decreases. We note that, the latter is still very enriched by the end of the calculation. Most of the enrichment of dust-to-gas ratio in the central objects is indeed actually taking 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 (PDF) of the dust ratio enrichment, denoted log(¯ ). It compares the distributions of the two most coupled and the three least coupled dust species (colored dashed lines) at t core + 2 kyr. It also indicates the PDF integrated over the grain size distribution (black line). These PDF 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. Small 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 Table 1. Syllabus of the different simulations, with the thermal-to-gravitational energy ratios α, maximum grain sizes S max . The initial mass-to-flux ratio µ as well as the tilt between the magnetic field and the rotation axis φ mag are given for simulations with magnetic field. Additionally, we also provide the formation time of the FHSC t core and and the initial Stokes number of the largest grains St 0,10 , the mass of the initial core and the number of dust bins. dust-ratio enrichment¯ 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 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 ≈ v g . Figure 4 shows the dust ratio for the total dust distribution (left) and for the 10 th species (right) in the midplane 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 essential of the 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 with large grains. Finally, we note that 4 kyr after the first core formation, the average value for the total dust-to-gas ratio is roughly unchanged but generally increasing (of about 1%) in the core and fragments since their formation. For the disk, we note a decrease of dust ratio of ∼ 30% for the largest grains (∼ 22% in total) in the disk since 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 threshold 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 ρ > 10 −13 g cm −3 delimited by the brown line indeed happens exactly when fragments form. This explains why these fragment tend to be more dust-rich than the first hydrostatic core. Interestingly, the dustto-gas ratio does not vary significantly in the volume where ρ > 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 the time of its formation by almost a factor of two.

Maximum grain size
As seen in Sect. 4.1 for the mmMRN model, the differential dynamics between gas and dust during the protostellar collapse depends critically on the grain sizes. Therefore, we perform two simulations 100micMRN and MRN with the same set of parameters as in mmMRN, but where we vary the maximum grain size. We choose S max = 100 µm for the 100micMRN model (which yields s 10 ∼ 22.6 µm) and S max = 250 nm for the MRN model (which yields s 10 ∼ 139 nm).
Let us first consider MRN, which is the model that has the smallest S max . Figure 7 shows the density of the gas (left) of the least coupled dust species (right) at t 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 ≈ 10 −5 cm which corresponds to an initial Stokes number St 0,10 ∼ 1.06 × 10 −5 1(see Sect 6 for a theoretical justification). As a result, dust is a excellent tracer of the gas in this model. This is illustrated by Fig. 8 that shows the probability density function of the dust ratio enrichment log(¯ ) in the core and fragments, the disk and the envelope at t 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 ≈ 1% in the core and the fragments, the disk and the envelope. The standard deviation for the dust ratio enrichment ranges between 2 × 10 −4 % (in the core) and 1.4 × 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 that a standard MRN distribution appears 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(¯ ) in the core and fragments, the disk and the envelope at t 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 the envelope. The typical standard deviation for the dust-ratio enrichment are 7 × 10 −3 in the core, 2.2 × 10 −2 in the disk and 7.8 × 10 −2 in Fig. 2. mmMRN test at t core + 1 kyr and t core + 2 kyr (t 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 right respectively). Values of the gas density are indicated by the colorbar on 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 0 ρ d < min(ρ g ) have voluntarily not been displayed to highlight the enriched regions. These depleted regions are delimited by the dashed grey lines. This choice of colors applies for all density maps in this study. Gas and dust are clearly not perfectly coupled. the envelope. This confirms that the larger the grains are, the more significant the decoupling with the gas is. We note that, for 100micMRN, it is reasonable to infer the gas density from the dust.

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 α = 0.25. This parameter is expected to affect strongly the dust dynamics. A lower value of α produces faster protostellar collapses prior to the first core formation due to the smaller thermal support. It thus develops faster high densities regions where dust strongly couples. In addition, a cloud with a smaller α 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 than in mmMRN in virtue of a smaller initial disk radius. Smaller disks with a smaller initial value of α 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 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 probability density function of the dust ratio enrichment log(¯ ) in the core and the fragments, the disk and the envelope at t 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 × 10 −2 in the FHSC and fragments, 3 × 10 −2 in the disk and 0.1 in the envelope. This was actually expected as the initial Stokes number scales as α 2/3 and is thus ≈ 0.6 times smaller in mmMRNa0.25 than in mmMRN.

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 µ = 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 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 (∝ J × B see Eq. 9) those are negligible compared to the hydrodynamical terms (∝ ∇P g ) in the decoupling of gas and dust in our magnetised collapse models.
We show fig. 14 and 15 the PDF of the dust ratio enrichment for the different objects at t 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 to the same conclusion that is a peaked distribution with Fig. 3. mmMRN test at t core + 2 kyr. Probability density function (PDF) of the dust ratio enrichment log(¯ ) 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. Fig. 4. mmMRN test ∼ 1 kyr (top), ∼ 2 kyr (middle) and ∼ 4 kyr (bottom) after the first core formation (t core = 72.84 kyr). Midplane 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 grains (right) dust-to-gas ratio is at its initial value, which can also be regarded as a dust enrichment line.  5. mmMRN test at t core + 1 kyr (t 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 . 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 and the core are 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 regions. We note that these pseudo-disks have a very broad PDF and show enriched and depleted regions in both the ideal and non-ideal cases. A similar behavior is observed in the magnetically driven outflow for the ideal case, that are dust-rich in dense regions and depleted at low densities. In the mmMRNnimhd, 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 seems 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), 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 at the benefit of an enrichment of high density regions. Cores and fragments at ρ > 10 −11 g cm −3 are formed in a dust rich environment. The dust-to-gas ratio is almost constant for ρ > 10 −11 g cm −3 because the increase of temperature due to the adiabatic contraction strengthen the coupling between the gas and the dust. 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.
the disk, the first hydrostatic core. and the high density regions of the outflow. Fig. 8. MRN test at t core + 2 kyr. Probability density function (PDF) of the total dust ratio enrichment log(¯ ) in the core and the fragments, the disk and the envelope. Dust is a very good tracer for the mass here and the dust distribution is almost uniform in the objects considered.

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 above it, it is enriched in dust during the collapse. If it lies under, it is dust depleted. This information is collected in Table D.2. Fig. 9. 100micMRN test at t core + 2 kyr. Probability density function (PDF) of the total dust ratio enrichment log(¯ ) in the core and the fragments, the disk and the envelope. Dust is a relatively good tracer for the gas although significant variations of the dustto-gas are observed. Fig. 10. mmMRNa0.25 test at t core + 2 kyr (t core = 23 kyr). Edgeon (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.

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 of Θ d,k≤6 m that remains unchanged for all the cores and fragments. Indeed, small grains have short stopping times Fig. 11. mmMRNa0.25 test at t core + 2 kyr. Probability density function (PDF) of the total dust ratio enrichment log(¯ ) in the core+fragments, the disk and the envelope. Dust is a relatively good tracer for the gas in the dense object although a notable depletion is observed in the envelope. 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 concentrate in the high density regions such as the core, the disk, the pseudo-disk and the inner regions of the outflow. 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 re- Fig. 13. mmMRNnimhd 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). As in the ideal case, dust is significantly decoupled from the gas and concentrate in the high density regions such as the core, the disk, the pseudo-disk and the inner regions of the outflow . Fig. 14. mmMRNmhd test at t core +2 kyr. Probability density function (PDF) of the dust ratio enrichment log(¯ ) in the core (blue), the disk (green) the pseudo-disk (purple), the outflow (orange) and the envelope (red). Dust is not a good tracer for the gas here and the dust distribution is not uniform in the considered objects. mains 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 frag- Fig. 15. mmMRNnimhd test at t core + 2 kyr. Probability density function (PDF) of the dust ratio enrichment log(¯ ) in the core (blue), the disk (green) the pseudo-disk (purple), the outflow (orange) and the envelope (red). Dust is not a good tracer for the gas here and the dust distribution is not uniform in the considered objects. Fig. 16. All the models at t core + 2 kyr. The dust-to-gas ratio enrichment averaged in mass is shown as a function of the grain size for all the objects. Grain with sizes smaller than 10 −4 cm are almost always perfectly coupled with the gas. For larger 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 F or the disk D and pseudo-disk P are enriched in dust. Low density objects such as the envelope E or the outflows O are depleted in dust. Magnetized models exhibit the stronger decoupling between the gas and the dust. ments have larger enrichment in large grains (k > 6, see Tab. 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 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 before, 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 that can be explained by two mechanisms. First, fragments form after the central object and thus stay a longer time in the isothermal phase where dust is less coupled since the temperature is smaller. Second, the dust-rich spiral arms developing through the envelope (see Figs. 4) are mainly accreted by the fragments (see Fig. 6). This provides an additional channel to enrich the fragments in solids.

Disks
We review the disk properties of all the models at t core + 2 kyr . Essentially, the values are similar to what is measured in the cores. This is essentially 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 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 non-uniform 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 core and t 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 J × B ∼ ∇P g for magnetic runs -dust cannot always be used as a direct proxy to trace the gas density. Although, as shown in Fig 5, the sub-structures seen in the dust originate from the ones in the gas, 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, longer integration time is required for the disk to grow significantly (Hennebelle et al. 2020).

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 ∼ 2 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 are more easily drifting 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. 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 of more than two orders of magnitude are not displayed (black background).

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 (9th and 10th 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 Θ d,10 m ∼ 0.22 for mmMRNmhd and Θ d,10 m ∼ 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 then inner regions and rather collects the grains from the envelope. This explains why the enrichment measured in the outflow are 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.

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 ≤ 6) do not experience significant dustto-gas ratio variations. Larger grains can however show significant dust-to-gas ratio variations. We measure a value of Θ d,10 m as low as 0.38 for mmMRNnimhd and mmMRNmhd, and 0.73 for mmMRN. This values are averaged over all the 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 is experiencing the larger dust-to-gas ratio variations (see for example Fig. 3 and 15). This is expected since the density in the envelope is lower than in the other objects. Typically, the depletion in large grain of the envelope increases with a decreasing density and have a larger dust content in their inner regions (see Fig.17 for both behavior).
We have shown that dust is not necessarily a good tracer for the gas density, i.e. we find important local variations in the dust distribution in mmMRN, mmMRNmhd, mmMRNnimhd. This has strong consequences for observations, since dust continuum radiation fluxes depend on densities integrated along the line of sight. It is therefore interesting to estimate error that arises when the total column density is estimated from the mass of a single dust bin k. We compute this error Err k from the total column density and the dust column density Σ d,k according to The total error Err is defined 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 by 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. Fig. 18. mmMRNmhd at t 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).

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.

Enrichment equation
In the terminal velocity approximation, and making the assumption that the collapse is is isothermal and purely hydrodynamical, the evolution of the dust ratio for a species k is given by 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¯ = 0 and are constrained by the same equation, i.e.¯ does not depend on 0 . The dust ratio enrichment is described by the equation where d dt ≡ ∂ ∂t + v · ∇. We now use the dimensionless variables whereρ = ρ ρ 0 . The termρ −1 appears in the divergence owing to t s,k ∝ 1 ρ .

Semi-analytical model
Neglecting local variations of¯ in comparison with local density variations yields We then obtain where χ ≡ e − τ 0 ∇ x (ρ −1 ∇ xρ )/ρdτ 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.

Estimate in the core
We can roughly approximate χ in the core and after a free-fall time t ff,0 = 3π 32 τ ff,0 . An order of magnitude estimate provides 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 χ 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.
where r var is the typical length at which variations of density become significant. Above ρ ad , the temperature are high and the gas and dust differential dynamics is negligible. A reasonable choice forρ is thereforeρ = ρ ad ρ 0 . In the typical condition of a protostellar collapse, one obtains |ln (χ) | ≈ 126 r var 1 AU −2 T gas 10 K .
We note that taking r var = 1 AU seems reasonable as it is about a tenth of the first-core radius in our models. For r var ≈ 5 AU, we find that |ln (χ) | ≈ 5. The value of χ strongly depends on the steepness of the pressure gradients, hence on r var . This model only provides an rough estimate of¯ and should not replace either a numerical treatment of the dust or a proper estimate of χ during the collapse. 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 χ. Finally, the blue area represents the range of dust enrichment obtained with the two extreme values of χ estimated in the previous section. We do not show the dust enrichment for MRN as it is clearly negligible (see Fig. 16).

Comparison with the models
In addition, we do not display the enrichment in the secondary fragments for the sake of readability. A 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 non-linear behaviour of Eq. 23, either due to local variations of 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 variation of dust-to-gas ratio, e.g 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 value estimated with our toy-model. This model being quite crude, we acknowledge that it cannot compete with an eventual model based on an accurate estimate of χ.

Critical Stokes number
Assuming that the value of χ is known, it can be used to determine the critical Stokes number St crit,¯ above which some regions can reach a given enrichment¯ . Equation 27 can indeed by inverted as Using the value of χ obtained by fitting our models (Fig.  19), we can estimate the typical Stokes numbers needed to get a dust enrichment by a factor of 2 in the core and the disk is approximately St crit,2 ∼ 0.01 − 0.027. With our toy model, we find St crit,2 ∼ 0.006 − 0.13. It 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 crit,1/2 ∼ 0.027 − 0.11. In short, it is easier to enrich the disk and the core than it is to deplete the envelope.

Summary of the models
We investigated the effect of several parameters on the dust dynamics during the protostellar collapse such as the thermal-togravitational energy ratio α, 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 a 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 2 − 3, even in the very first thousands of years of the protostellar collapse (see mmMRN, mmMRNmhd or mmMRNnimhd). Hence, in order to understand the initial dust and gas content of protoplanetary disks, it is crucial to measure accurately the dust size distribution in early prestellar cores. This conclusion is reinforced by recent observations (Galametz et al. 2019) or synthetic observations ) that seem to probe the existence of ∼ 100 µm in Class 0 objects.
For small initial values of α, e.g. in mmMRNa0.25, the collapse is fast and large grains do not have the time to significantly 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 densities regions (see Fig.6) and on the range of densities experienced during the collapse (see

Dust depleted
Dust enriched Schematic view of a dusty protostellar collapse 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 is 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.
Sect. 6). Hence, although 100micMRN has a smaller maximum grain size as 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 initially properties of the cloud appear to be extremely important to quantify the evolution of the dust distribution during the protostellar collapse. It would be therefore interesting to study the dust collapse of a Bonor-Ebert sphere, since its free-fall timescale is usually longer than the one of the Boss and Bodenheimer test (Machida et al. 2014). We leave this proper comparison to further works. With mmMRNmhd and mmMRNnimhd, we investigated the effect of a magnetic field on the dynamics of dust during the protostellar collapse. We qualitatively find similar results as in our fiducial case. Quantitatively, the decoupling between the gas and the dust does however produce more significant variations 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 time evolution.
For the sake of summarizing, we show in Fig. 20 a schematic view of a dusty protostellar collapse a few kyr 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 of 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 were we do not quantitatively show the local variations of the dust-to-gas ratio.

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 as 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 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 provides, less significant dust-to-gas ratio variations as in Bate & Lorén-Aguilar (2017). This effect was predicted in Bate & Lorén-Aguilar (2017). We note that¯ 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 enrichment is independent from the initial value of the dust-to-gas ratio. This allows us a more direct comparison with Bate & Lorén-Aguilar (2017). In their study, they observe an increase of dust-to-gas ratio of about one order of magnitude for the 100 µm grains, which is about 3 times larger than what we observe in mmMRN for example. Using the calibrated value of χ and Eq.27, we can estimate that the dust ratio enrichment of 100 µm for grains with ρ grain = 3 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 densities (∼ 10 −20 g cm −3 )and therefore larger initial Stokes numbers (∼ 1 for 100 µm grains). In 2D simulation of collapsing gravitoviscous protoplanetary disks, Vorobyov et al.  Elbakyan et al. (2020) have focused on the evolution of dust including grain growth and fragmentation. Similarly to our models, they observe local variations of the dust-to-gas ratio in the disks. They also found larger dust-to-gas ratios in the inner regions of the collapse and smaller dust-to-gas ratios a few hundreds of AU away from the core. 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,  observe particularly large increase of 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 simulations of non-rotating gas and dust mixture, considering only single dust species (1 µm, 10 µm and 100 µm). In this work, we already observed a significant decoupling occurring for 100 µm grains. However the increase of dust-to-gas ratio in the core was strong only in the outer regions of the collapse. As previously said, we can estimate that in mmMRN the dust-to-gas ratio enrichment of 100 µm with ρ grain = 3 g cm −3 would be about 2.27. It was only ∼ 1.2 in our first non-rotating spherical collapse calculation, although both models have the same initial α. 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 details since 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.

Possible implications for planet formation
The simulations presented in this study consolidate the idea that protostellar collapses may form protoplanetary disks containing 0 2 − 3% of their mass under the form of solids. When the dust ratio 0 is larger than the square of the aspect ratio of the disk H r 2 -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 2 0.01 (Andrews et al. 2010). The former condition is strengthened by the fact that 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 one, dust may exert strong drag onto the gas and powers 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 0 2 − 3%. In any case, large dust contents favour the formation of planetesimals through the development of the streaming instability (e.g. Johansen et al. 2007Johansen et al. , 2009Draż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 , dust-rich density clumps, similar to the secondary cores in our fragmenting models, could be a favored locus of giant planet formation. They indeed noted both the piling-up of large grains and important growth in these clumps.
In short, a larger initial dust content always favour planet formation in disks, and this may be in a dramatic manner. A quantitative knowledge of the differential dynamics of gas and dust during the protostellar collapse and the initial dust size distribution in prestellar cores therefore appears essential to understand the early stages of planet formation. In that perspective, an extensive study of dust dynamics and coagulation/fragmentation should be done from the scales of molecular clouds and, through the protostellar collapse, up to protoplanetary disks.

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 where Z k is the number of charges on the dust grain k, e = 1.6 × 10 −20 Abc the electron charge and E is the electric field. The 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 magnetic and electric fields evolution. To compare the magnetic and neutral drags strength, it is interesting to compare the stopping-to-gyration time ratio also known as Hall factor Γ k . In the Epstein regime, this ratio writes Z k e|B| ρc s s 2 grain,k . (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 γ = 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 where B 0 is the initial magnetic field, given by µ and ρ 0 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); Vaytet et al. (2018). One can then show that Fig. 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 Γ k is larger than unity for very small grain (s grain,k 10 −6 cm). These grains are very well coupled the the gas in the neutral case but, when charged, could experience a strong decoupling with the neutrals. This would occur if the latter are decoupled from the magnetic field e.g., in the non-ideal regime. The Lorentz drag might play a crucial role 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 collapses in future studies. We emphasize that the electromotive term in the Lorentz force applies only on 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 Γ k < 1.

Caveat: Coagulation/fragmentation during the collapse
Dust coagulation has been neglected during this study. It may however affect strongly dust evolution during the collapse since dust decoupling depend on grain sizes. Following Draine (1985), one can estimate the coagulation timescale t coag,i, j between two dust phases i and j due to their relative motion within the collapsing clouds as Assuming the same density for all the grains, neglecting cumulative back-reaction effects and magnetic fields, and considering an isothermal collapse, differential velocities can be estimated from the diffusion approximation as As said in Sect. 6, the density profile of the free-falling material can be approximated as a power law with an exponent ζ = −2 (Larson 1969). In this case, one obtains at a distance r from the central region where q i, j ≡ s grain,j s grain,i is the ratio between the grain sizes. In the limit q i, j 1, we find Growth induced by the relative dynamics between dust grains is therefore expected to be not very efficient during protostellar collapse away from the core but could be non-negligible in the inner regions. At r = 10 AU, t coag,i, j ∼ 17 kyr, which is quite smaller than the typical free-fall timescale of a prestellar core. We note that growth could also be enhanced by the 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) hence, assuming again that q i, j 1 the coagulation timescale writes as We now consider a region of density 10 −12 g cm −3 and temperature of 10 K. Assuming s grain, j = 0.1 µm and s grain,i = 100 µm, we get t coag,i, j ∼ 240 kyr. We note that, in the case of the Brownian motions, the coagulation timescale depends on the grains size of both species. If we now consider s grain, j = 0.01 µm, t coag,i, j ∼ 7 kyr. In short, we are tempted to say that in the presence of large grains, very small grain could be efficiently removed during the collapse in high density regions. Coagulation should therefore be included in future studies. We admit however 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, e.g up to a few ≈ 0.1km s −1 in the envelope of mmMRN in the case of the two least coupled species. Typically, the velocity above which fragmentation can occur is thought to be about a few 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 coagulation timescale (e.g. Gonzalez et al. 2017). Fragmentation could simply not have the time to occur during the collapse, especially considering that large grains quickly drift toward regions of high density where their drift velocity is typically around a few meter per seconds. Let us now estimate the typical fragmentation radius r frag , which corresponds to the distance from the center at which the fragmentation velocity v frag is equal to the dust drift velocity. It can be obtained by solving ρ grain s grain c s |∇ρ| Let us now assume first core formation with a core profile of the type ρ = ρ F 1+ r r F 2 .We can show that Assuming ρ F = 10 −11 g cm −3 , r F = 10 AU, ρ grain = 1g cm −3 , s grain = 160 µm and T g = 10 K r frag = 2.5 × 10 4 AU r F 10 AU 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 frag = 2.5 × 10 2 AU. Hence, fragmentation may be relevant for small cores and should be investigated in further studies.

Conclusion and perspective
In this study, we presented the first multigrain and non-ideal MHD simulations of dusty protostellar collapses using the new dust dynamics solver of RAMSES. We presented six dustycollapse simulations with a simultaneous treatment of 10 dust species with different sizes. In these simulations, we investigate the impact of the maximum grain size, the thermal-togravitational 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 less than a few 10 µm 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 250nm is however extremely well preserved during the protostellar collapse in absence of coagulation. 5. Dust dynamics is strongly affected by the initial cloud properties. Variations of the dust-to-gas ratio reach the largest values 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 have shown that it can be expressed as 0 χ St 0 where χ is a dimensionless function of the time independent on the dust properties. Fitting the values of χ gives a very good agreement between the semi-analytical model and the measured dust-ratios in the range of Stokes number considered in our simulations. 7. Using the calibration of χ 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 2. Similarly we show that grains with St 0 0.0027 can potentially be depleted by a factor of 2 in the envelope after the first core formation.
Dust evolution during the protostellar collapse could have serious consequence 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, since the Lorentz drag cannot necessarily be neglected for small grains. Coagulation and fragmentation of dust grains should also be considered to investigate more realistically dust evolution during the star formation process.

Appendix A: Impact of velocity regularization
As explained in Sect. 3.5, the simulations of our fiducial model have been 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 4 models. With w cap = 0.5 − 2 km s −1 we essentially find the same results as in the fiducial case that has w cap = 1 km s −1 . Having w cap = 0.1 km s −1 however appears to be a too 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: Non-linear 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 backreaction or not, we perform an additional simulation similar to the fiducial one, but with θ d,0 = 10 −7 (mmMRNeps1e-7). For mmMRNeps1e-7, we expect the back-reaction to be definitely negligible as θ d,0 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 for mmMRN (grey). 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 χ for the FHSC (dotted), disk (solid) and envelope (dashed), respectively. We display the same information for mmMRN at t core in the disk (grey).
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, i.e. a non-negligible ∇ , are the source of discrepancy between a pure exponential enrichment and the values measured in our models in the disk. We indeed observe a 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-in-size grains, with a dust-to-gas ratio of 1%. Local electroneutrality ensures that n ions − n e + Z d n d = 0, (C.1) n ions , n e and n d being the number density of the ions, electrons and dust, the latter being given by where m grain is the average mass of a dust grain. Two additional constrains are provided by the evolution equations for the ions and electrons charge density (Umebayashi & Nakano 1980;Fujii et al. 2011). Assuming steady state and only considering charge capture by the grains (Keith & Wardle Table D.1. Initial dust distributions, rounded quantities for the three maximum grain size used in our models (the exact calculation can be made using the method presented in 3.3) . For each S max , the first line represent 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.