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

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 1000yr 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.


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(Marchand et al. , 2020;;Zhao et al. 2016Zhao et al. , 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 ρ > 10 −13 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 postprocessed 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.

Methods
We performed non-ideal MHD simulations with the RAMSES code (Teyssier 2002).RAMSES is an Eulerian gas dynamics 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.

Coagulation
In Paper I, we demonstrate that the coagulation process as described by the Smoluchowski (1916) equation is a onedimensional (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 χ that encompasses the whole history of the physical conditions seen by the grains.At a given χ 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 χ can be derived from integrating: 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 χ is needed to track the coagulation of grains.Equation ( 1) is a Lagrangian derivative of χ with respect to time along the path of the grain.We can transform it into a partial (Eulerian) derivative: where u is the velocity of the gas.We can combine Eq. 2 with the mass conservation equation: where ρ is the gas mass density.This yields: which means that the quantity ρχ 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 χ is therefore calculated self-consistently in all cells at each timestep, 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 χ values in a large relevant interval (between χ = 10 13 and χ = 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 χ.
Our initial distribution in this paper is a Mathis et al. (1977) distribution (i.e., MRN).The minimum and maximum radii are 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: 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 χ.Below χ = 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 ρ > 10 12 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.

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, µ i , and the sticking probability of electrons on grains, s e .Here, we assume µ i = 28, which corresponds to the ion HCO + , and s e = 0.6 as in Marchand et al. (2016).We also assume ζ = 5 × 10 −17 s −1 .The density and the temperature, are taken from the hydrodynamic 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.

Model
We performed four numerical simulations using the RAMSES code.We solved the following MHD equations: where u is the velocity of the gas, P its pressure, B the magnetic field, J = ∇ × B the current, I the identity matrix, Φ the gravitational potential, and η Ω , η H , and η AD are the ohmic, Hall, and ambipolar resistivities, respectively.The temperature evolution is prescribed by the barotropic equation: with T 0 = 10 K and γ = 5/3 as the adiabatic index.The initial condition is a sphere of gas of 1 M ⊙ .The radius of the sphere is controlled by the thermal over gravitational energy ratio α.We chose α = 0.3, which sets a radius of R = 2946 au.The density distribution has an m = 2 azimuthal perturbation of where ρ 0 is the density of a uniform sphere of same mass and radius, and φ the azimuthal angle.We choose δ ρ = 0.05.The computational domain outside of the sphere is filled with gas of density ρ 0 /100.The sphere undergoes a solid rotation, with a ratio of rotational to gravitational energy of β = 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): with 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 µ B = 5 as a fiducial value.Those parameters are used in our reference case C-3.We changed the value of α to 0.4, and the initial mass, M, to 5 M ⊙ 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.

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.

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.

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): where τ ff is the free-fall time, given by: where ρ 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 ρ(t) = ρ 0 (R 0 /R[t]) 3 .The gas density then evolves as: 1 We used a second-order Runge-Kutta scheme to numerically integrate this equation from the beginning of the collapse at ρ = ρ 0 until the formation of the first Larson core at ρ = 10 −13 g cm −3 .The evolution of χ with ρ is plotted in Fig. 2, assuming T = 10 K.The solid and dashed lines represent different initial densities, ρ 0 = 3.8 × 10 −20 g cm −3 and ρ 0 = 3.8 × 10 −18 g cm −3 , respectively.At densities n H > 10 −16 g cm −3 , the value of χ is independent from ρ 0 and increases as χ ∝ ρ 1/4 .This evolution is expected as χ ∼ ρ 3/4 t and t ∼ τ ff ∼ ρ −1/2 .At ρ = 10 −13 g cm −3 , the coagulation variable reaches χ ≈ 5.6 × 10 17 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.

Numerical collapse
The numerical simulations were run as described in Sect. 2 until 1000 yr after the density reaches 10 −13 g cm −3 , the formation of 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.

Grain growth
We show in Fig. 4 the value of χ 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 ρ < 10 15 g cm −3 .Beyond this density, there is a large spread of χ 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.
A61, page 4 of 11 There is no significant difference in the χ 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 χ 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 (ρ < 10 −16 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 µm.The third region is the first Larson core, which has even larger grains that reach 400 µ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 midplane 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 µ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 µm and a max = 10 µm.

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 removing the very small grains (a < 0.1 µm) and redistributing their mass to the larger end of the distribution (Zhao et al. 2016(Zhao et al. , 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 timestep, 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 /(ρη AD Ω).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 ≳ 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 largemass-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 ρ > 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.

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 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.

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 (>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.

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 cm −3 (see their Fig.7), while in our simulations, the peak exceeds 100 µm at a lower maximum density of n H ≈ 3 × 10 12 cm −3 (or ρ ≈ 10 −11 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 A61, page 7 of 11 A&A 670, A61 (2023) used (Ormel et al. 2009;Guillet et al. 2020) Re = 6.7 × 10 7 n H 10 5 cm −3 1 2 . (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 pseudodisk is an over-density, typically denser than ρ ≈ 10 −15 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 pseudodisk 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.

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.

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 (ρ > 10 −15 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 pseudodisk 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 ≥100 µm in radius in disks.The effects of grain growth on chemistry and radiative transfer will be explored in future papers.
where C is a constant, g local is a function of the local properties of the gas (density, temperature...), and The relative number of grains of size, a, at a given time, t, is X(a, t) = n grain /n H where n grain is the number density of grains.We can then write X = d g x, where d g is the dust-to-gas ratio, and x is the normalized relative number of grains.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 overestimate 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 differential 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 roll , that determines the energy needed to restructure the grain.In Section 4.2 of paper I, we misinterpreted their results by assuming the fragmentation would occur for a kinetic energy E kin = 5E roll .The actual criterion is: where E kin is the kinetic energy, E br is the breaking energy and N tot is the total numbers of monomers composing the two colliding grains; E br is defined by (Dominik & Tielens 1997): E br = A br γ 5/3 grain (a 0 /2) 4/3 ε 2/3 , (C.2) with A br = 2.8 × 10 3 , γ grain as the surface energy density of the material, a 0 as the size of the monomers composing the grains, and ε as the reduced elastic modulus.As in Ormel et al. (2009), we adopt a 0 = 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, γ grain = 25 erg cm −2 and ε = 2.8 × 10 11 dyn cm −2 .The kinetic energy of two grains of mass m and m ′ is expressed as: 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: (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 γ grain we use is taken from Ormel et al. (2009), but the surface energies of different materials can reach 10 to 100 times this value, resulting in a much higher estimate of the fragmentation threshold (as

Fig. 1 .
Fig. 1.State of the grain size distribution for values of χ, 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.

Fig. 3 .
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 edgeon 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.

PFig. 4 .
Fig. 4. Evolution of the coagulation variable χ 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 ρ 0 = 3.8 × 10 −20 g cm −3 .

Fig. 5 .Fig. 6 .
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.

Fig. 7 .Fig. 8 .
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

Fig. 9 .
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.
z 0 = 2.97, k B and G are the Boltzmann and gravitational constants, γ = 5/3 the adiabatic index of the gas, ρ s = 2.3 g cm −3 the bulk density of the grains, µ = 2.3 the average atomic mass of the gas, m H the proton mass, and a the radius of the larger of the two grains.Assuming two identical grains and

Table 1 .
Parameters of the protostellar collapse simulations.Columns are: name of the simulation, thermal over gravitational energy ratio α, radius R, mass M, initial density, ρ 0 , and initial magnetic field B of the initial cloud, formation time of the first hydrostatic core, t fc , and use of coagulation.