# Fast methods to track grain coagulation and ionization. III. Protostellar collapse with non-ideal MHD

Pierre Marchand, Ugo Lebreuilly, Mordecai-Mark Mac Low, V. Guillet

► To cite this version:

Pierre Marchand, Ugo Lebreuilly, Mordecai-Mark Mac Low, V. Guillet. Fast methods to track grain coagulation and ionization. III. Protostellar collapse with non-ideal MHD. Astronomy and Astrophysics - A&A, 2023, 670, pp.A61. 10.1051/0004-6361/202244291. hal-03926347

HAL Id: hal-03926347

https://hal.science/hal-03926347

Submitted on 13 Oct 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Fast methods for tracking grain coagulation and ionization

## III. Protostellar collapse with non-ideal MHD

P. Marchand<fnref n="1" />,<fnref n="2" /> ID, U. Lebreuilly<fnref n="3" /> ID, M.-M. Mac Low<fnref n="2" /> ID, and V. Guillet<fnref n="4" />,<fnref n="5" /> ID

Received 17 June 2022 / Accepted 23 December 2022

## ABSTRACT

Dust grains influence many aspects of star formation, including planet formation and the opacities for radiative transfer, chemistry, and the magnetic field via Ohmic, Hall, as well as ambipolar diffusion. The size distribution of the dust grains is the primary characteristic influencing all these aspects. Grain size increases by coagulation throughout the star formation process. In this work, we describe numerical simulations of protostellar collapse using methods described in earlier papers of this series. We compute the evolution of the grain size distribution from coagulation and the non-ideal magnetohydrodynamics effects self-consistently and at low numerical cost. We find that the coagulation efficiency is mostly affected by the time spent in high-density regions. Starting from sub-micron radii, grain sizes reach more than 100 µm in an inner protoplanetary disk that is only 1000 yr old. We also show that the growth of grains significantly affects the resistivities, while also having an indirect effect on the dynamics and angular momentum of the disk.

**Key words.** stars: formation – dust, extinction – magnetohydrodynamics (MHD)

## 1. Introduction

Grains play a major role during star formation. First, they are the seeds of planet formation. While their characteristic size is sub-micron in the interstellar medium (ISM; Mathis et al. 1977), they grow by coagulation during the collapse and can reach sizes larger than 10 µm in the early stages of protostellar collapse (Guillet et al. 2020; Silsbee et al. 2020; Tsukamoto et al. 2021; Vorobyov et al. 2022; Bate 2022). Observations also suggest they reach these sizes (and possibly larger ones) in the envelopes of Class 0-I objects (Kwon et al. 2009; Miotello et al. 2014; Le Gouellec et al. 2019; Galametz et al. 2019). Their growth then continues in protoplanetary disks until they eventually become planetesimals. Variations in their size also significantly impact non-ideal magnetohydrodynamical (MHD) effects through their ionization and their chemical interactions with the gas, with direct feedback to the dynamics of the gas (Marchand et al. 2016, 2020; Zhao et al. 2016, 2018; Guillet et al. 2020). Non-ideal MHD effects have been shown to be critical for the regulation of magnetic field and angular momentum during the protostellar collapse and the protoplanetary disk evolution (Mouschovias & Paleologou 1979; Machida et al. 2006; Duffin & Pudritz 2008; Mellon & Li 2009; Li et al. 2011; Tomida et al. 2015; Wurster et al. 2016; Masson et al. 2016; Vaytet et al. 2018; Marchand et al. 2020; Lebreuilly et al. 2021). Grains are also the main source of opacity in protostellar environments, affecting the cooling of the gas and the observations made of those systems. Their high optical depth at densities of $\rho > 10^{-13} \text{ g cm}^{-3}$ leads to the formation of the first hydrostatic core (Larson 1969).

However, numerical simulations usually do not account for grain coagulation self-consistently due to the great cost of computing a coagulation algorithm on the fly (although new methods are being developed; see the recent work by Lombart & Laibe 2021). The dust evolution used to be pre-processed or post-processed with no self-consistent feedback on the dynamics (Rossi et al. 1991; Dullemond & Dominik 2005; Zhao et al. 2016; Marchand et al. 2020). Recently, more and more studies include the growth of grains in their hydrodynamics simulations (Tsukamoto et al. 2021; Vericel et al. 2021; Vorobyov et al. 2022; Bate 2022). In Marchand et al. (2021, hereafter Paper I), we presented a simple and fast method to track coagulation in a self-consistent way that is particularly suited for modeling star formation. Here, we apply this method to non-ideal MHD protostellar collapse simulations. It is coupled with the second method presented in Paper I: a fast calculation of the ionization and grain charge to obtain non-ideal MHD resistivities. These 3D simulations are the first to include a self-consistent grain growth with direct feedback to the dynamics through the self-consistent calculation of MHD resistivities.

The paper is organized as follows. In Sect. 2, we describe the methods used in our study. Section 3 presents the results of our calculations, both analytical in Sect. 3.1 and numerical in Sect. 3.2. We compare our results to other works and discuss the caveats in Sect. 4. We present our conclusions in Sect. 5.

## 2. Methods

We performed non-ideal MHD simulations with the RAMSES code (Teyssier 2002). RAMSES is an Eulerian gas dynamics

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UT3-PS, CNRS, CNES, 9 av. du Colonel Roche, 31028 Toulouse Cedex 4, France</footnote>

<footnote n="2">² Department of Astrophysics, American Museum of Natural History, 200 Central Park West, New York, NY 10024, USA
e-mail: pierre.marchand.astr@gmail.com</footnote>

<footnote n="3">³ AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique spatiale, rue Jean Teillac, 91405 Orsay, France</footnote>

<footnote n="5">⁵ Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France</footnote>
code with adaptive mesh refinement (AMR) and self-gravity. It includes a monofluid treatment of non-ideal MHD effects (Masson et al. 2012; Marchand et al. 2018). We have implemented the methods presented in Paper I to calculate the coagulation and ionization of grains on-the-fly in a self-consistent manner.

## 2.1. Grain coagulation and ionization

### 2.1.1. Coagulation

In Paper I, we demonstrate that the coagulation process as described by the Smoluchowski (1916) equation is a one-dimensional (1D) process with certain types of coagulation kernels. Consequently, the size distribution of the coagulated grains depends only on the initial size distribution and a variable $\chi$ that encompasses the whole history of the physical conditions seen by the grains. At a given $\chi$ that is integrated along the path of the grains, the coagulated distribution is always the same, independently of the actual path taken. This method works for every coagulation kernel for which the dependence on the gas variables, such as density and temperature, can be separated from the grain properties such as size and mass. In Paper I and in the present paper, we use the turbulent kernel derived by Ormel & Cuzzi (2007) in the intermediate coupling regime, which is suited for star formation conditions. We show that in this case $\chi$ can be derived from integrating:

$$d\chi = n_H^{3/4} T^{-1/4} dt, \quad (1)$$

where $n_H$ is the number density of the gas, $T$ its temperature, and $t$ the time. In three-dimensional (3D) hydrodynamical simulations, only the knowledge of $\chi$ is needed to track the coagulation of grains.

Equation (1) is a Lagrangian derivative of $\chi$ with respect to time along the path of the grain. We can transform it into a partial (Eulerian) derivative:

$$\frac{\partial \chi}{\partial t} + \mathbf{u} \cdot \nabla \chi = n_H^{3/4} T^{-1/4}, \quad (2)$$

where $\mathbf{u}$ is the velocity of the gas. We can combine Eq. 2 with the mass conservation equation:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot [\rho \mathbf{u}] = 0, \quad (3)$$

where $\rho$ is the gas mass density. This yields:

$$\frac{\partial \rho \chi}{\partial t} + \nabla \cdot (\rho \chi \mathbf{u}) = \rho n_H^{3/4} T^{-1/4}, \quad (4)$$

which means that the quantity $\rho\chi$ can be treated in an Eulerian framework as a passive scalar with a source term. We exploit this property and implement it as such in RAMSES. The value of $\chi$ is therefore calculated self-consistently in all cells at each time-step, as a mass-weighted average, ignoring any diffusion of dust through the gas. We used the Ishinisan code (Marchand et al. 2021) to pre-calculate a table containing the grain size distribution for a large number of $\chi$ values in a large relevant interval (between $\chi = 10^{13}$ and $\chi = 10^{19}$ in cgs units). When the size distribution is needed during the hydrodynamical simulation, it is interpolated from the table based on the value of $\chi$.

Our initial distribution in this paper is a Mathis et al. (1977) distribution (i.e., MRN). The minimum and maximum radii are

<figure><img src="image_2.png" /><figcaption>Fig. 1. State of the grain size distribution for values of $\chi$, between $10^{15}$ cgs and $10^{19}$ cgs. The points represent the fractional abundance of the size-bin as function of the effective radius of the bin.</figcaption></figure>

$a_{\min} = 5$ nm and $a_{\max} = 250$ nm, and the slope of the distribution is -3.5, so that the number density of grains, $n$, follows the following variation:

$$\frac{dn}{da} \propto a^{-3.5}. \quad (5)$$

The total quantity of grains is determined by the dust-to-gas mass ratio that we assume to be 0.01. In this work, we sample the distribution with 60 bins of size logarithmically spaced between 5 nm and 5000 µm. In Fig. 1, we present the coagulated MRN size distribution at various $\chi$. Below $\chi = 10^{17}$ cgs, the shift in the size distribution is negligible. For higher values, the maximum size and the peak of the distribution are located at larger and larger radii, while the slope of the distribution remains similar. In all cases, the small grains are more abundant while the large grains hold more mass. The mode of the size distribution $a_{\max}$ is located near the largest relevant grain size.

As grains grow, they may experience fragmentation when the kinetic energy of the collision is high. Contrarily to what we derived in Paper I, we find that fragmentation does not occur in early stages of the disk, but it is present, rather, only at very high densities $\rho > 10^{12} \text{ cm}^{-3}$ (Lebreuilly et al. 2023) for very large grains with $a > 0.08$ cm. We demonstrate this result in Appendix C. We therefore neglect fragmentation in this work. We also do not account for the grain drift with respect to the gas in this work. We discuss the possible consequences in Sect. 4.2. Drift is, however, compatible with our coagulation model, and we detail the method in Appendix B. Other limitations are discussed in Sect. 4.4.

### 2.1.2. Ionization and resistivities

In Paper I, we also presented a fast method to calculate the ionization of the gas-grain mixture. For an arbitrary size distribution, we can calculate the average electric charge of each grain size, the number of ions, and the number of electrons, provided the cosmic-ray (CR) ionization rate, the density and temperature of the gas, the average atomic mass of ions, $\mu_i$, and the sticking probability of electrons on grains, $s_e$. Here, we assume $\mu_i = 28$, which corresponds to the ion HCO$^+$, and $s_e = 0.6$ as in Marchand et al. (2016). We also assume $\zeta = 5 \times 10^{-17} \text{ s}^{-1}$.
<table><caption>**Table 1.** Parameters of the protostellar collapse simulations.</caption><thead><tr><th>Name</th><th>α</th><th>R (au)</th><th>M (M<sub>☉</sub>)</th><th>ρ<sub>0</sub> (g cm<sup>-3</sup>)</th><th>B (μG)</th><th>t<sub>fc</sub> (kyr)</th><th>Coagulation</th></tr></thead><tbody><tr><td>C-3</td><td>0.3</td><td>2946</td><td>1</td><td>5.49 × 10<sup>-18</sup></td><td>133</td><td>~30</td><td>Yes</td></tr><tr><td>C-3-M5</td><td>0.3</td><td>14738</td><td>5</td><td>2.19 × 10<sup>-19</sup></td><td>27</td><td>~150</td><td>Yes</td></tr><tr><td>C-4</td><td>0.4</td><td>3930</td><td>1</td><td>2.31 × 10<sup>-18</sup></td><td>75</td><td>~47</td><td>Yes</td></tr><tr><td>NC-3</td><td>0.3</td><td>2946</td><td>1</td><td>5.49 × 10<sup>-18</sup></td><td>133</td><td>~30</td><td>No</td></tr></tbody></table>

**Notes.** Columns are: name of the simulation, thermal over gravitational energy ratio $\alpha$, radius $R$, mass $M$, initial density, $\rho_0$, and initial magnetic field $B$ of the initial cloud, formation time of the first hydrostatic core, $t_{fc}$, and use of coagulation.

The density and the temperature, are taken from the hydrodynamical simulation. The calculation is performed by the Newton-Raphson scheme described in Appendix A of Paper I. The resistivities are computed using a similar method to Marchand et al. (2016), with one difference. For each grain size, they sum over the contributions of the whole charge distribution (between -1 to +1 in their case). Instead, we average the contributions using the mean electric charge. We explain the method in greater detail in Appendix A.

## 2.2. Star formation simulations

### 2.2.1. Model

We performed four numerical simulations using the RAMSES code. We solved the following MHD equations:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot [\rho \mathbf{u}] = 0, \qquad (6) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot \left[ \rho \mathbf{u} \mathbf{u} + \left( P + \frac{B^2}{2} \right) \mathbf{I} - \mathbf{B} \mathbf{B} \right] = -\rho \nabla \Phi, \qquad (7) $$

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times \left[ \mathbf{u} \times \mathbf{B} - \eta_{\Omega} \mathbf{J} - \eta_{H} \frac{\mathbf{J} \times \mathbf{B}}{B} + \eta_{AD} \frac{(\mathbf{J} \times \mathbf{B}) \times \mathbf{B}}{B^2} \right] = 0, \qquad (8) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (9) $$

where $\mathbf{u}$ is the velocity of the gas, $P$ its pressure, $\mathbf{B}$ the magnetic field, $\mathbf{J} = \nabla \times \mathbf{B}$ the current, $\mathbf{I}$ the identity matrix, $\Phi$ the gravitational potential, and $\eta_{\Omega}$, $\eta_{H}$, and $\eta_{AD}$ are the ohmic, Hall, and ambipolar resistivities, respectively. The temperature evolution is prescribed by the barotropic equation:

$$ T = T_0 \left( 1 + \left[ \frac{\rho}{10^{-13} \text{ g cm}^{-3}} \right]^{\gamma-1} \right), \qquad (10) $$

with $T_0 = 10$ K and $\gamma = 5/3$ as the adiabatic index. The initial condition is a sphere of gas of 1 $M_\odot$. The radius of the sphere is controlled by the thermal over gravitational energy ratio $\alpha$. We chose $\alpha = 0.3$, which sets a radius of $R = 2946$ au. The density distribution has an $m = 2$ azimuthal perturbation of

$$ \rho(\varphi) = \rho_0(1 + \delta_{\rho} \sin \varphi), \qquad (11) $$

where $\rho_0$ is the density of a uniform sphere of same mass and radius, and $\varphi$ the azimuthal angle. We choose $\delta_{\rho} = 0.05$. The computational domain outside of the sphere is filled with gas of density $\rho_0/100$. The sphere undergoes a solid rotation, with a ratio of rotational to gravitational energy of $\beta = 0.02$. The magnetic field is initially uniform and parallel to the rotation axis.

It is defined using the mass-to-flux ratio over the critical value (Mouschovias & Spitzer 1976):

$$ \mu_B = \frac{M/\Phi_B}{(M/\Phi_B)_{\text{crit}}}, \qquad (12) $$

with

$$ \left( \frac{M}{\Phi_B} \right)_{\text{crit}} = \frac{0.53}{3\pi} \sqrt{\frac{5}{G}}. \qquad (13) $$

Observations show that dense cores are slightly super-critical (Crutcher 1999), although recent numerical simulations indicate that observations may overestimate the actual strength of the magnetic field due to projection effects (Kuznetsova et al. 2020). We chose $\mu_B = 5$ as a fiducial value. Those parameters are used in our reference case C-3. We changed the value of $\alpha$ to 0.4, and the initial mass, $M$, to 5 $M_\odot$ for two other cases referred to as C-4 and C-3-M5, respectively, to investigate the influence of the collapse time on the grain coagulation. The final run, named NC-3, is similar to C-3 without coagulation. The grain distribution and ionization evolve as described in Sects. 2.1 and 2.1.2, and in Paper I. The four cases are summarized in Table 1.

### 2.2.2. Grid and algorithm

The simulation box is a cube that is four times as large as the radius of the sphere, with periodic boundary conditions. The initial grid is composed of $32^3$ cells (level 5 of AMR) and is refined to ensure at least ten points per Jeans length, strongly satisfying the Truelove et al. (1997) criterion. The maximum AMR refinement level is 13 for C-3 and NC-3 (resolution of 1.4 au), 14 for C-4 (resolution of 0.96 au), and 16 for C-3-M5 (resolution of 0.90 au).

Simulations were performed with a 3D unsplit slope limiter to avoid overshooting of the magnetic field at shock boundaries, while keeping the second order convergence for the Hall effect. We used the HLLD Riemann solver (Miyoshi & Kusano 2005) for non-magnetic variables and the 2D HLL Riemann solver for the magnetic field and the Hall effect (Balsara 2012; Marchand et al. 2018). The Poisson equation is solved using the multigrid method of Guillet & Teyssier (2011) in our periodic domain.

## 3. Application to star formation

In this section, we present the results of our calculations. We first applied our coagulation method to an analytical one-zone model in Sect. 3.1, then in 3D MHD simulations in Sect. 3.2.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of the coagulation variable, $\chi$, with increasing density, $n_H$, during the isothermal collapse. The solid line represents an initial density of $\rho_0 = 3.8 \times 10^{-20} \text{ g cm}^{-3}$ is and the dashed line represents $\rho_0 = 3.8 \times 10^{-18} \text{ g cm}^{-3}$.</figcaption></figure>

### 3.1. Analytical collapse

During spherical protostellar collapse, the time evolution of the contraction ratio of a gas cloud compared to its original radius $x = R(t)/R_0$ can be described as (Flower et al. 2005):

$$ \frac{dx}{dt} = -\frac{\pi}{2\tau_{ff}} \sqrt{\frac{1}{x} - 1}, \qquad (14) $$

where $\tau_{ff}$ is the free-fall time, given by:

$$ \tau_{ff} = \sqrt{\frac{3\pi}{32G\rho_0}}, \qquad (15) $$

where $\rho_0$ is the initial density. Guillet et al. (2020) showed that assuming a uniform compression of the gas nicely reproduces the isothermal phase of the collapse, particularly when comparing the dust size distribution at the same gas density. In this case, the density of mass scales as $\rho(t) = \rho_0(R_0/R[t])^3$. The gas density then evolves as:

$$ \frac{1}{n_H} \frac{dn_H}{dt} = \frac{1}{\rho} \frac{d\rho}{dt} = -\frac{3}{x} \frac{dx}{dt}. \qquad (16) $$

We used a second-order Runge-Kutta scheme to numerically integrate this equation from the beginning of the collapse at $\rho = \rho_0$ until the formation of the first Larson core at $\rho = 10^{-13} \text{ g cm}^{-3}$. The evolution of $\chi$ with $\rho$ is plotted in Fig. 2, assuming $T = 10 \text{ K}$. The solid and dashed lines represent different initial densities, $\rho_0 = 3.8 \times 10^{-20} \text{ g cm}^{-3}$ and $\rho_0 = 3.8 \times 10^{-18} \text{ g cm}^{-3}$, respectively. At densities $n_H > 10^{-16} \text{ g cm}^{-3}$, the value of $\chi$ is independent from $\rho_0$ and increases as $\chi \propto \rho^{1/4}$. This evolution is expected as $\chi \sim \rho^{3/4}t$ and $t \sim \tau_{ff} \sim \rho^{-1/2}$. At $\rho = 10^{-13} \text{ g cm}^{-3}$, the coagulation variable reaches $\chi \approx 5.6 \times 10^{17} \text{ cgs}$, which corresponds to a peak of the size distribution of ~10 µm, indicating a significant grain growth. That is consistent with the results of Guillet et al. (2020), who use the same collapse model, but solve coagulation on the fly.

### 3.2. Numerical collapse

The numerical simulations were run as described in Sect. 2 until 1000 yr after the density reaches $10^{-13} \text{ g cm}^{-3}$, the formation of

<figure><img src="image_4.png" /><figcaption>Fig. 3. Density slices of the C-3 simulation, with grain coagulation. Top panel is a face-on slice of the plane $z = 0$, and the bottom panel an edge-on slice of the plane $y = 0$. White arrows represent the direction of the gas velocity. The snapshot is taken at the final time-step, 1000 yr after the formation of the first Larson core.</figcaption></figure>

the first hydrostatic core at ~30 kyr (Table 1). In reference simulation C-3, a small circumstellar disk with a radius of ≈20 au forms and a disk wind is launched by magneto-centrifugal acceleration (Blandford & Payne 1982). Figure 3 shows face-on and edge-on slices of density at the final time-step of the simulation, with arrows indicating the gas velocity.

#### 3.2.1. Grain growth

We show in Fig. 4 the value of $\chi$ as a function of the density in simulation cells at the final time-step, for runs C-3, C-4, and C-3-M5. The increase is quasi unidimensional in the isothermally collapsing envelope for $\rho < 10^{15} \text{ g cm}^{-3}$. Beyond this density, there is a large spread of $\chi$ values of over an order of magnitude. This spread most likely occurs in gas falling in the pseudo-disk, then the disk and the first Larson core, or the outflow, over different timescales, spending unequal times in a given density range. The overall trend agrees well with the analytical calculation.
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Evolution of the coagulation variable $\chi$ with increasing density $n_H$ in the numerical collapse models C-3 (purple), C-3-M5 (light blue), and C-4 (green). Each point corresponds to a simulation cell at the final time-step. The red line is the analytical collapse solution for $\rho_0 = 3.8 \times 10^{-20} \text{ g cm}^{-3}$.

There is no significant difference in the $\chi$ values, and thus dust size distributions, between the three runs, despite run C-4 needing 50% more time to collapse to the first Larson core stage, and C-3-M5 needing 400% more time. Coagulation happening in the isothermally collapsing envelope is therefore hardly impacted by the initial conditions, as growth accelerates with increasing density.

Figure 5 shows the mode of the size distribution as a function of density, corresponding to the distribution of $\chi$ shown in Fig. 4. Three regions are clearly demarcated, the first being the envelope, in which grain coagulation is not efficient enough to form large grains ($\rho < 10^{-16} \text{ g cm}^{-3}$). The second comprises the pseudo-disk and the early protoplanetary disk, where grains grow by a factor 100 from sub-micron sizes to several tens of $\mu$m. The third region is the first Larson core, which has even larger grains that reach 400 $\mu$m within a mere $10^3$ yr after its formation. That value is in line with similar recent studies (Kawasaki et al. 2022; Lebreuilly et al. 2023). There is little difference between runs C-3, C-4, and C-3-M5, confirming that coagulation in the envelope does not impact large grains, as found by previous studies (for example Silsbee et al. 2022). We discuss observations of large grains in the envelope in Sect. 4.2.

The spatial distribution of those grain sizes for run C-3 is displayed in Fig. 6. Size distributions shifting significantly from the initial MRN distribution are indeed confined to the mid-plane in the disk and pseudo-disk. The bottom panel also shows moderately larger grains in the outflow, as they only traveled through the upper layers of the pseudo-disk before being ejected (Marchand et al. 2020). If they had been in the mid-plane of the disk, they would have grown much more, as coagulation is irreversible in our model. The upper panel shows that grains reach a radius larger than 1 $\mu$m in the outskirts of the disk, within 100 au of the center. Growth then occurs rapidly as density increases in the inner 15 au, which is shown by the almost overlapping contours delimiting $a_{max} = 5 \mu$m and $a_{max} = 10 \mu$m.

### 3.2.2. Resistivities and gas dynamics

We describe here the impact of grain coagulation on non-ideal MHD resistivities and their macroscopic effects on gas dynamics. Previous studies emulated the coagulation of grains by

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

**Fig. 5.** Mode of the coagulated grain size distribution as a function of density in simulations (purple), C-3-M5 (light blue), and C-4 (green). The discrete values of the sizes are due to the binning of the size distribution.

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

**Fig. 6.** Slices of the C-3 simulation with coagulation showing the mode of the grain size distribution $a_{max}$, 1000 yr after the first core formation. The color scale is different for each panel to yield a better contrast. Top panel: black contours indicate $a_{max} = 1, 5$ and $10 \mu$m. Bottom panel: $a_{max} = 0.5$ and $1 \mu$m.
removing the very small grains ($a < 0.1$ µm) and redistributing their mass to the larger end of the distribution (Zhao et al. 2016, 2018; Marchand et al. 2020). This method leads to an increase in resistivities, in particular the ambipolar resistivity, resulting in weaker coupling between the magnetic field and the gas, hence weaker magnetic braking and larger, more unstable disks. However, what we observe here is the exact opposite behavior.

The middle panel of Fig. 7 presents the volume-weighted average resistivities as a function of density, at the final time-step, for simulations C-3 and NC-3. A non-evolving size distribution produces resistivities relatively similar in the envelope, but two to four orders of magnitude larger at disk densities, particularly the Ohmic and ambipolar resistivities. Consequently, the magnetic braking is weakened, and the gas retains more angular momentum without the grain coagulation, forming a larger disk, as shown in Fig. 8.

This difference in regime can be quantified by the ambipolar Elsasser number $Am = B^2/(\rho\eta_{AD}\Omega)$. In regions where $Am < 1$, the ambipolar diffusion has a significant impact on the dynamics of the gas. The bottom panel of Fig. 7 shows the radial profile of the ambipolar Elsasser number in runs C-3 and NC-3, azimuthally averaged in the mid-plane. The higher resistivity in run NC-3 results in $Am < 1$ in the inner ~12 au, indicating active ambipolar diffusion and magnetic field dissipation, while $Am \gtrsim 10^4$ for C-3 over the same radial range, indicating weak ambipolar diffusion. Figure 9 compares the disk size and angular momentum in both simulations, further confirming the lower magnetic braking in run NC-3. A second effect of the weaker magnetic forces from the stronger ambipolar diffusion in run NC-3 is the absence of outflow at this stage of evolution, as shown in the lower panel of Fig. 8.

The discrepancy of resistivity values between actual coagulation and methods simply redistributing the mass to the large-mass-end of the distribution originates from the lack of very large grains (>10 µm) in the latter. In both cases, removing the small grains decreases the electron absorption by dust, and therefore decreases the Ohmic resistivity. However, the ambipolar resistivity is controlled by the relative abundance of ions and charged grains in the gas. Although the small grain removal method barely changes the abundance of ions, the dominant positive charge carriers, it reduces significantly the charged grain population, leading to an increase in resistivity at low density. At high density, both with a standard MRN and a truncated-MRN distribution, the grains become the dominant charge carriers and the abundance of ions decreases (see the dashed lines in the top panel of Fig. 7). That does not happen with coagulation because the number of grains decreases significantly, hence, leading to a higher number of ions and a lower resistivity than without proper coagulation. This kind of method without larger grain creation is therefore inappropriate for emulating coagulation alone at a low cost for non-ideal MHD calculations. At later times, at densities $\rho > 10^{-12}$ g cm$^{-3}$, Lebreuilly et al. (2023) showed that the replenishment of small grains by fragmentation would lead to an increase in both Ohmic and ambipolar resistivities.

# 4. Discussion and caveats

## 4.1. Grain growth

As explained in previous sections and displayed in Figs. 4 and 5, initial gas conditions have little to no influence on grain coagulation for later stages and coagulation is ineffective in the envelope for growing large grains. This reinforces the idea that the system forgets the initial conditions at the formation of the first

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

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

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

**Fig. 7.** Non-ideal MHD effects in the C-3 and NC-3 simulations. Top panel: abundances of ions (purple) and electrons (green) as a function of density in models with coagulation (C-3; solid lines) or without it (NC-3; dashed lines). Middle panel: volume-averaged Ohmic (purple), ambipolar (green), and Hall (blue) resistivities for the same models. Bottom panel: average ambipolar Elsasser number $Am$ in the mid-plane as a function of radius. The thin grey line represents $Am = 1$.

hydrostatic core and the disk (Vaytet & Haugbølle 2017). Consequently, calculations such as the ones presented in this work may provide standard initial dust grain size distributions for studies of protoplanetary disk evolution. Although other coagulation kernels affect the size distribution in different ways (see Sect. 4.4), it is certain that even young protoplanetary disks contain grains
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Density slices of the NC-3 simulation, without grain coagulation. *Top panel* is a face-on slice of the plane $z = 0$, and the *bottom panel* an edge-on slice of the plane $y = 0$. White arrows represent the direction of the gas velocity. The snapshot is taken at the final time-step, 1000 yr after the formation of the first Larson core.</figcaption></figure>

that are significantly larger than those in the classical MRN distribution. This has important implications for the dynamics of grains in the disk. Larger grains couple differently to the gas and may trigger the streaming instability (Youdin & Goodman 2005; Johansen & Youdin 2007; Yang et al. 2017), which is an early step toward planet formation. Although this regime is not reached in our simulations, the fast growth of the grains in circumstellar disks could predict an early onset for this process.

## 4.2. Large grains in envelopes and dust diffusion

Although grain coagulation is negligible in the envelope of our simulations, large grain signatures in envelopes have been observed. Galametz et al. (2019), for example, report “low and varying dust emissivity indices” at the envelope scale for some Class 0 and Class I protostars. This could be due to the presence of mm-size grains in low numbers. The time-scale to form such large grains in the envelope and cold ISM cores is large

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Disk size (purple, left axis) and angular momentum (green, right axis) as a function of time after the formation of the first Larson core. Solid lines represent simulation C-3 with coagulation, and dashed lines are simulation NC-3 without coagulation.</figcaption></figure>

(>100 Myr), so we exclude the possibility of early coagulation. Similarly, Valdivia et al. (2019) reported that their synthetic polarization observations of young protostellar envelopes in RAMSES calculations require grains larger than 10 µm in order to be consistent with observations, which is also inconsistent with our findings and those of recent studies (Silsbee et al. 2022; Lebreuilly et al. 2023, for the latest ones).

The aerodynamic properties of the grains may cause differential velocities between grain populations and the gas, leading to varying dust-to-gas ratio throughout the cloud (Lebreuilly et al. 2020). We note, however, that (as these authors have found) dust diffusion will start to play a significant role only for very large grains of a few hundred microns. Generally, large grains tend to accumulate in higher density regions. Tsukamoto et al. (2021) found that this dust diffusion can lead to what they call an ash fall phenomenon, in which large coagulated grains (up to millimeter-size) from the disk are ejected by an outflow, then decouple from the gas and fall back in the envelope. This process may explain the observations of low spectral index in Class 0 envelopes, as the outflow fuels the envelope with large grains formed in the central region. Eventually, those ejected grains circle back to the outer edge of the disk, enriching the large-end of the size distribution. In our case, the disk wind is fueled by the upper layers of the pseudo-disk, in which dust grows only moderately. However, dust-to-gas ratio enhancement in this region due to the grain differential velocities may lead to an accelerated growth. Further studies are needed to investigate this discrepancy about the size of grains in envelopes between observations and theory.

## 4.3. Coagulation in the pseudo-disk

Grains in our simulations seem to grow faster than in Bate (2022) despite their use of the same turbulent kernel as in our work. The peaks of the distributions in that work reach a few microns at a maximum density of $n_H = 10^{13} \text{ cm}^{-3}$ (see their Fig. 7), while in our simulations, the peak exceeds 100 µm at a lower maximum density of $n_H \approx 3 \times 10^{12} \text{ cm}^{-3}$ (or $\rho \approx 10^{-11} \text{ g cm}^{-3}$; see our Fig. 5). That discrepancy is mainly due to a different modeling of the Reynolds number to calculate the turbulent coagulation kernel. They assumed a constant value of Re = $10^8$, while we
used (Ormel et al. 2009; Guillet et al. 2020)

$$Re = 6.7 \times 10^7 \left( \frac{n_H}{10^5 \text{ cm}^{-3}} \right)^{\frac{1}{3}}. \quad (17)$$

Hence, in the central regions, our Reynolds number can be larger by three orders of magnitude. Their grains are then stuck in the tightly-coupled regime where relative grain velocities are lower than in the intermediate coupling-regime, which is more adapted to this situation (as demonstrated in Sect. 4.1 of Paper I). The lower relative velocities result in lower coagulation rates and, therefore, a slower growth rate. The poor constraints on the value of the Reynolds number in protostellar environments therefore represents a source of uncertainty for grain growth by coagulation.

An additional explanation may involve an excess of growth in the pseudo-disk that forms in our simulations. The pseudo-disk is an over-density, typically denser than $\rho \approx 10^{-15} \text{ g cm}^{-3}$, created by the convergence of gas flowing along magnetic field lines (Galli & Shu 1993) that takes the shape of a disk perpendicular to the magnetic field, but is not supported against gravity. It is shown in the bottom panels of Figs. 3 and 8.

Gas in the rotationally supported disk generally comes directly from the pseudo-disk and the overall efficiency of coagulation is mainly affected by the time spent in high-density regions such as the pseudo-disk. This is what appears in the bottom panel of Fig. 6, where there is a ~30 au-thick layer of large grains in the mid-plane. A passage of grains through the pseudo-disk would therefore provide an acceleration of coagulation in the early stage of star formation, even before arrival in the disk. That does not happen in Bate (2022) since they do not consider magnetic fields, resulting in less-coagulated size distributions as grains enter the disk.

### 4.4. Coagulation kernel

Our coagulation methods works for every coagulation kernel where the environment variables and grain variables can be separated. One example of this is the well-known turbulent kernel derived by Ormel & Cuzzi (2007) that we use in this work. This kernel is appropriate to calculate the growth of large grains, but it has several limitations. We assume a steady-state Kolmogorov turbulence spectrum and that the injection-scale of the turbulence is equal to the Jeans length corresponding to the local density (see Sect. 2.2 of Paper I), which may lead to an overestimation of the grain collision rate. Other kernels may have different effects on the size distribution of grains. Guillet et al. (2020) showed that this turbulent kernel would increase the maximum size of the distribution, while a kernel derived from ambipolar diffusion, which creates a drift between charged and neutral grains, is efficient at removing small grains. This also happens in Bate (2022) and Lebreuilly et al. (2023), where combine several processes are combined, including Brownian motion and pressure gradients that generate drift velocities between grains and rapidly remove the smaller grains. These processes are not, however, efficient at producing larger grains without the help of the turbulent kernel. Contrary to fragmentation, the addition of those kernels would steepen the distribution at its lower end.

### 5. Conclusions

We present here the results of numerical simulations of protostellar collapse with dust coagulation and self-consistent calculation

of non-ideal MHD resistivities, using the methods detailed in Marchand et al. (2021). We performed four simulations, including three with coagulation and different collapse times and one without coagulation for reference. Here, are our main results:

- Coagulated size distributions retain some characteristics of the initial MRN distribution, in particular, the dominance of small grains in number and large grains in mass.

- Dust coagulation is inefficient at growing larger grains in the envelope, even with long free-fall times. What matters is the time spent in high-density regions ($\rho > 10^{-15} \text{ g cm}^{-3}$) such as the pseudo-disk. Fragmentation can also be ignored in the cloud-collapse phase.

- Grain growth is extremely rapid in the disk. The peak of the size distribution in mass, which is also the largest relevant grain size of the distribution, reaches 1 µm in the pseudo-disk and more than 100 µm in the inner disk only $10^3$ yr after its formation.

- Grain sizes have a significant impact on non-ideal MHD resistivities. Coagulated grains result in resistivities up to four orders of magnitude lower than non-coagulated grains, with a significant impact on the dynamics of the disk. Simple redistribution approximations fail to capture this effect, as it occurs because of the growth of the largest grains.

Accounting for grain coagulation is therefore necessary in star formation and protoplanetary disk simulations, as grains grow rapidly to $\ge 100 \text{ µm}$ in radius in disks. The effects of grain growth on chemistry and radiative transfer will be explored in future papers.

*Acknowledgements.* P.M. acknowledges the financial support of the Kathryn W. Davis Postdoctoral Fellowship of the American Museum of Natural History and of the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (ERC Starting Grant "Chemtrip", grant agreement No 949278). M.-M.M.L. was partly supported by US NSF grant AST18-15461. U.L. acknowledges the financial support of the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130).

### References

Balsara, D. S. 2012, J. Comput. Phys., 231, 7476

Bate, M. R. 2022, MNRAS, 514, 2145

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883

Crutcher, R. M. 1999, ApJ, 520, 706

Desch, S. J., & Mouschovias, T. C. 2001, ApJ, 550, 314

Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

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

Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2005, A&A, 436, 933

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

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

Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756

Guillet, V., Hennebelle, P., Pineau des Forêts, G., et al. 2020, A&A, 643, A17

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627

Kawasaki, Y., Koga, S., & Machida, M. N. 2022, MNRAS, 515, 2072

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

Kuznetsova, A., Hartmann, L., & Heitsch, F. 2020, ApJ, 893, 73

Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.-F., & Kemball, A. J. 2009, ApJ, 696, 841

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

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

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

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

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

Le Gouellec, V. J. M., Hull, C. L. H., Maury, A. J., et al. 2019, ApJ, 885, 106

Li, Z.-Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180

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

Machida, M. N., Inutsuka, S., & Matsumoto, T. 2006, ApJ, 647, L151
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., Tanaka, K. E. I., Commerçon, B., & Chabrier, G. 2020, ApJ, 900, 180

Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50 (Paper I)

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

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

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

Mellon, R. R., & Li, Z. 2009, ApJ, 698, 922

Miotello, A., Testi, L., Lodato, G., et al. 2014, A&A, 567, A32

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

Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204

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

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

Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. T. 2009, A&A, 502, 845

Pinto, C., & Galli, D. 2009, Three-Fluid Magnetohydrodynamics in Star Formation, eds. K., Tsinganos, T., Ray, & M., Stute, 597

Rossi, S. C. F., Benevides-Soares, P., & Barbuy, B. B. 1991, A&A, 251, 587

Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39

Silsbee, K., Akimkin, V., Ivlev, A. V., et al. 2022, ApJ, 940, 188

Smoluchowski, M. V. 1916, Zeitsch. Physik, 17, 557

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

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

Tsukamoto, Y., Machida, M. N., & Inutsuka, S.-i. 2021, ApJ, 920, L35

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

Vaytet, N., & Haugbølle, T. 2017, A&A, 598, A116

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

Vericel, A., Gonzalez, J.-F., Price, D. J., Laibe, G., & Pinte, C. 2021, MNRAS, 507, 2318

Vorobyov, E. I., Skliarevskii, A. M., Molyarova, T., et al. 2022, A&A, 658, A191

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

Yang, C.-C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80

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

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050

Zhao, B., Caselli, P., Li, Z.-Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868
# Appendix A: Calculating the resistivities

## A.1: Solving for the ionization

Let us assume a size distribution of grains divided into N bins. We then have the following system of equations (see Paper I):

$$Z_k = \psi \tau_k + \frac{1 - \epsilon^2 \Theta^2}{1 + \epsilon \Theta \alpha_k + \epsilon^2 \Theta^2}, \qquad (A.1)$$

$$\langle \tilde{J}(\tau_k) \rangle = (1 - \psi) + \frac{\frac{2}{\tau_k} [\epsilon^2 \Theta^2 + \epsilon \Theta]}{1 + \epsilon \Theta \alpha_k + \epsilon^2 \Theta^2}, \qquad (A.2)$$

$$\epsilon = \frac{1 - \psi}{\Theta e^{\psi}}, \qquad (A.3)$$

$$n_i = - \frac{1}{1 - \epsilon} \sum_k n_k Z_k, \qquad (A.4)$$

$$f(\psi) = \frac{\langle \sigma v \rangle_{ie} \epsilon n_i^2}{\zeta n_H} + \frac{n_i v_i}{\zeta n_H} \sum_k n_k \pi a_k^2 \langle \tilde{J}(\tau_k) \rangle - 1 = 0. \qquad (A.5)$$

Here, $\epsilon = n_e/n_i < 1$ is the ratio between the number of electrons and ions, $\psi$ is the ratio between the electric potential of the grains and the kinetic energy of electrons, $n_k$, $Z_k$, and $a_k$ represent the number density, the average charge and the radius of grains in bin $k$, $\zeta$ is the CR ionization rate, and $v_i = [8k_BT/\pi\mu_i m_H]^{1/2}$ is the thermal velocity of ions, with $T$ the temperature, $\mu_i$ the average atomic mass of ions, and $m_H$ the mass of the proton. We also have $\Theta = s_e[\mu_i m_H / m_e]^{1/2}$, with $s_e$ as the sticking probability of electrons on grains and $m_e$ the electron mass. Then, $\langle \tilde{J}(\tau_k) \rangle$ represents the enhancement factor for ion recombination on grains, $\tau_k = a_k k_B T / e^2$ is the reduced temperature of grains (Draine & Sutin 1987), and $\alpha_k = [8/(\pi \tau_k)]^{1/2}$. Finally, $\langle \sigma v \rangle_{ie} = 2 \times 10^{-7} [T/300]^{-1/2}$ is the collision rate between ions and electrons.

We solve the system of relevant equations (A.1-A.5) for $\psi$ and find $n_i$, $n_e$, and $Z_k$ for all bins. A Newton-Raphson algorithm is described in details in the Appendix A of Paper I.

## A.2: The resistivities

The collision time-scale of species s with neutral H₂, the most abundant species in the gas, is given by:

$$\tau_{sH_2} = \frac{1}{a_s He} \frac{m_s + m_{H_2}}{m_{H_2}} \frac{1}{n_{H_2} \langle \sigma \omega \rangle_s}. \qquad (A.6)$$

Here, $a_s$He accounts for the collisions with He and is equal to 1.14 for ions, 1.16 for electrons, and 1.28 for grains (Desch & Mouschovias 2001); $m_s$ and $m_{H_2}$ are the mass of the species, s, and the H₂ molecule; $n_{H_2}$ represents the number density of H₂ (roughly equal to the density of the gas); $\langle \sigma \omega \rangle_s$ is the collision rate, taken from Pinto & Galli (2009) for electrons and ions:

$$\langle \sigma \omega \rangle_e = 3.16 \times 10^{-11} v_{rms,e}^{1.3}, \qquad (A.7)$$

$$\langle \sigma \omega \rangle_i = 2.4 \times 10^{-9} v_{rms,i}^{0.6}, \qquad (A.8)$$

where

$$v_{rms,s} = \left( \frac{8k_B T}{\pi \mu_{s,H_2}} \right)^{\frac{1}{2}}, \qquad (A.9)$$

is in km s⁻¹, with $\mu_{s,H_2}$ the reduced mass of species s and H₂. These velocities have been calculated in a three-fluid formalism

but we use them in our monofluid framework. See some related concerns in the appendix of Marchand et al. (2020). For grains, the collision rate is given by (Draine & Sutin 1987; Kunz & Mouschovias 2009):

$$\langle \sigma \omega \rangle_k = \pi a_k^2 \left( \frac{8k_B T}{\pi m_{H_2}} \right)^{\frac{1}{2}} \left( 1 + \left( \frac{\pi}{2\tau_k} \right)^{\frac{1}{2}} \right). \qquad (A.10)$$

We also define the conductivity, $\sigma_s$, and cyclotron frequency, $\omega_s$, of species s:

$$\sigma_s = \frac{n_s q_s^2 \tau_{sH_2}}{m_s}, \qquad (A.11)$$

$$\omega_s = \frac{q_s B}{c m_s}, \qquad (A.12)$$

where $q_s$ is the electric charge of species s, B the magnetic field strength and c the speed of light. The parallel, perpendicular, and Hall conductivities are expressed as:

$$\sigma_{\parallel} = \sum_{s} \sigma_{s}, \qquad (A.13)$$

$$\sigma_{\perp} = \sum_{s} \frac{\sigma_{s}}{1 + (\omega_{s} \tau_{sH_{2}})^{2}}, \qquad (A.14)$$

$$\sigma_H = \sum_s \frac{\sigma_s \omega_s \tau_{sH_2}}{1 + (\omega_s \tau_{sH_2})^2}. \qquad (A.15)$$

$$\sigma_{\parallel} = \sum_{s} \sigma_{s}, \qquad (A.16)$$

Finally, the Ohmic, Hall, and ambipolar resistivities are defined as:

$$\eta_O = \frac{1}{\sigma_{\parallel}}, \qquad (A.17)$$

$$\eta_H = \frac{\sigma_H}{\sigma_{\perp}^2 + \sigma_H^2}, \qquad (A.18)$$

$$\eta_A = \frac{\sigma_{\perp}}{\sigma_{\perp}^2 + \sigma_H^2} - \frac{1}{\sigma_{\parallel}}. \qquad (A.19)$$

$$\eta_A = \frac{1}{\sigma_{\parallel}}, \qquad (A.19)$$

# Appendix B: Coagulation and dust-to-gas ratio

Our coagulation model assumes that grains are perfectly coupled with the gas and well-mixed, so that the dust-to-gas mass ratio is constant. However, grains of different sizes have different aerodynamic properties and may experience a significant drift through the gas. This characteristic is usually determined by their Stokes number, that is, the ratio between the stopping time of the grain (Epstein 1924) and the dynamical timescale of the system. In particular, Lebreuilly et al. (2020) showed strong variations in the dust-to-gas mass ratio during the protostellar collapse. The disk and first core tend to be dust-enriched while the envelope becomes dust-depleted. This has strong implication for the coagulation, as a higher grain density promotes collisions and speeds up their growth (conversely for a lower density). In this section, we briefly explore a refinement to the coagulation model presented in Paper I.

The general expression of the Smoluchowski (1916) equation, from which we derive the expression of $\chi$, is (Paper I, Equation 5):

$$\frac{dX(a, t)}{dt} = C g_{\text{local}} n_{\text{H}} I(a, X, t), \qquad (B.1)$$
where C is a constant, g<sub>local</sub> is a function of the local properties of the gas (density, temperature...), and

$$
\begin{equation}
\begin{aligned}
I(a, X, t) = & - \int_{0}^{\infty} h(m, m')X(m, t)X(m', t)dm' \\
& + \frac{1}{2} \int_{0}^{m} h(m - m', m')X(m - m', t)X(m', t)dm'.
\end{aligned}
\tag{B.2}
\end{equation}
$$

The relative number of grains of size, a, at a given time, t, is
X(a, t) = n<sub>grain</sub>/n<sub>H</sub> where n<sub>grain</sub> is the number density of grains.
We can then write X = d<sub>g</sub>x, where d<sub>g</sub> is the dust-to-gas ratio,
and x is the normalized relative number of grains. We can then
rewrite Equation B.1:

$$
\frac{dx(a,t)}{dt} = C g_{\text{local}} n_H d_g I'(a, x, t), \quad (\text{B.3})
$$

and include $d_g$ in the definition of our coagulation variable,
giving it the following form:

$$
d\chi' = g_{\text{local}} n_{\text{H}} d_{\text{g}} dt. \quad (\text{B.4})
$$

In our case, for the Ormel & Cuzzi (2007) kernel,

$$
d\chi' = n_H^{\frac{3}{4}} T^{-\frac{1}{4}} d_g dt. \quad (B.5)
$$

This expression means that the dust-to-gas ratio does not change the evolution path of the size distribution, only the coagulation speed. A dust-to-gas ratio, N, times higher requires a χ value that is N times lower to reach the same coagulated state. This new definition of χ' can be used in an environment with varying dust-to-gas ratio. In hydrodynamical simulations, it is possible to couple it with the drift of one grain size close to the peak of the mass density distribution, as allowed by the method proposed by Lebreuilly et al. (2019).

Although it is not yet possible to associate this coagulation
method with the differential drift of grains of different sizes, this
is a first step towards a self-consistent treatment of dust grains
in hydrodynamics simulations. This method will lead to an over-
estimate of the small-to-larger grain ratio, but this is probably a
decent approximation to get the total dust-to-gas ratio. As shown
in Lebreuilly et al. (2020), as long as the large grains dominate
the mass in the distribution, they also dominate the differen-
tial gas and dust dynamics. In other words, most of the dust
enrichment comes from the dynamics of large grains.

**Appendix C: The fragmentation barrier**

Ormel et al. (2009) derived criteria for the fragmentation of dust
aggregates. They define the rolling energy, E<sub>roll</sub>, that determines
the energy needed to restructure the grain. In Section 4.2 of
paper I, we misinterpreted their results by assuming the frag-
mentation would occur for a kinetic energy E<sub>kin</sub> = 5E<sub>roll</sub>. The
actual criterion is:

$$
E_{\text{kin}} > 5N_{\text{tot}}E_{\text{br}}, \quad (C.1)
$$

where E<sub>kin</sub> is the kinetic energy, E<sub>br</sub> is the breaking energy
and N<sub>tot</sub> is the total numbers of monomers composing the two
colliding grains; E<sub>br</sub> is defined by (Dominik & Tielens 1997):

$$
E_{\text{br}} = A_{\text{br}} \gamma_{\text{grain}}^{5/3} \frac{(a_0/2)^{4/3}}{\varepsilon^{2/3}}, \quad (C.2)
$$

with A<sub>br</sub> = 2.8 × 10<sup>3</sup>, γ<sub>grain</sub> as the surface energy density of the
material, a<sub>0</sub> as the size of the monomers composing the grains,
and ε as the reduced elastic modulus. As in Ormel et al. (2009),
we adopt a<sub>0</sub> = 0.1 µm. Ice mantles on grains make them more
resistant to fragmentation. Therefore, we assume bare silicates
to obtain a lower limit for a fragmentation criterion. In this case,
γ<sub>grain</sub> = 25 erg cm<sup>-2</sup> and ε = 2.8 × 10<sup>11</sup> dyn cm<sup>-2</sup>. The kinetic
energy of two grains of mass m and m' is expressed as:

$$
E_{\text{kin}} = \frac{1}{2} \frac{mm'}{m + m'} \Delta v^2, \quad (\text{C.3})
$$

where Δv is the relative velocity between the two grains. This
velocity is given by the kernel of Ormel & Cuzzi (2007) that we
use in Paper I:

$$
\Delta v = \left( \frac{3}{\sqrt{8}} z_0 [k_B G]^{\frac{1}{2}} \frac{\gamma \rho_s}{\mu m_H} \right)^{\frac{1}{2}} n_H^{-\frac{1}{4}} T^{-\frac{1}{4}} a^{\frac{1}{2}}, \quad (\text{C.4})
$$

where z₀ = 2.97, k<sub>B</sub> and G are the Boltzmann and gravitational
constants, γ = 5/3 the adiabatic index of the gas, ρ<sub>s</sub> = 2.3 g cm<sup>-3</sup>
the bulk density of the grains, μ = 2.3 the average atomic mass
of the gas, m<sub>H</sub> the proton mass, and a the radius of the larger of
the two grains. Assuming two identical grains and

$$
m = \frac{4}{3} \pi \rho_s a_0^3 N_{\text{tot}}, \quad (C.5)
$$

Equation C.1 becomes

$$
a > \frac{15E_{\text{br}}}{\pi \rho_s a_0^3} \left( \frac{3}{\sqrt{8}} z_0 [k_B G]^{\frac{1}{2}} \frac{\gamma \rho_s}{\mu m_H} \right)^{-1} n_H^{\frac{1}{2}} T^{\frac{1}{2}}. \quad (C.6)
$$

Replacing the variables with numeric values, we find

$$
a > 826 \, \mu m \left( \frac{n_H}{10^{10} \, \text{cm}^{-3}} \right)^{\frac{1}{2}} \left( \frac{T}{10 \, \text{K}} \right)^{\frac{1}{2}} . \quad (C.7)
$$

This limit is much higher than the one derived in Paper I and
we do not reach it in our simulations. The value of γ<sub>grain</sub> we use
is taken from Ormel et al. (2009), but the surface energies of
different materials can reach 10 to 100 times this value, result-
ing in a much higher estimate of the fragmentation threshold (as
E<sub>br</sub> ~ γ<sub>grain</sub><sup>5/3</sup>). The value derived in Equation C.7 should then be
considered very conservative.