© The Authors 2023

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 signiflcantly 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 ρ > 10⁻¹³ g cm⁻³ 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 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 χ 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: $$\text{d}\chi {=n}_{\text{H}}^{\frac{3}{4}}{T}^{-\frac{1}{4}}\text{d}t,$$(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 χ 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: $$\frac{\partial \chi }{\partial t}+u\cdot \nabla \chi =n_{\text{H}}^{\frac{3}{4}}{T}^{-\frac{1}{4}},$$(2)

where u is the velocity of the gas. We can combine Eq. 2 with the mass conservation equation: $$\frac{\partial \chi }{\partial t}+\nabla \cdot \left[\rho u\right]=0,$$(3)

where ρ is the gas mass density. This yields: $$\frac{\partial \rho \chi }{\partial t}+\nabla \cdot \left[\rho \chi u\right]=\rho n_{\text{H}}^{\frac{3}{4}}{T}^{-\frac{1}{4}},$$(4)

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 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 χ values in a large relevant interval (between χ = 10¹³ and χ = 10¹⁹ 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: $$\frac{\text{d}n}{\text{da}}\propto a^{-3.5}.$$(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 χ. Below χ = 10¹⁷ 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¹² cm⁻³ (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.

Fig. 1 State of the grain size distribution for values of χ, between 1015 cgs and 1019 cgs. The points represent the fractional abundance of the size-bin as function of the effective radius of the bin.

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, µ_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⁻¹⁷ s⁻¹.

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.

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 \left[\rho u\right]=0,$$(6) $$\frac{\partial \rho u}{\partial t}+\nabla \cdot \left[\rho uu+\left(P+\frac{B^2}{2}\right)\mathbb{I}-BB\right]=-\rho \nabla \text{Φ,}$$(7) $$\frac{\partial B}{\partial t}-\nabla \times \left[u\times B-{\eta }_{\text{Ω}}J-{\eta }_{\text{H}}\frac{J\times B}{B}+{\eta }_{\text{AD}}\frac{\left(J\times B\right)\times B}{B^2}\right]=0,$$(8) $$\nabla \cdot B=0,$$(9)

where u is the velocity of the gas, P its pressure, B the magnetic field, J = ∇ × B the current, $\mathbb{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: $$T=T_0\left(1+{\left[\frac{\rho }{{10}^{-13}\text{g\hspace{0.17em}}{\text{cm}}^{-3}}\right]}^{\gamma -1}\right),$$(10)

with T₀ = 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 $$\rho \left(\phi \right)={\rho }_0\left(1+{\delta }_{\rho }\sin \phi \right),$$(11)

where ρ₀ 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 ρ₀/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): $${\mu }_{\text{B}}=\frac{M/{\text{Φ}}_{\text{B}}}{{\left(M/{\text{Φ}}_{\text{B}}\right)}_{\text{crit}}},$$(12)

with $${\left(\frac{M}{{\text{Φ}}_{\text{B}}}\right)}_{\text{crit}}=\frac{0.53}{3\pi }\sqrt{\frac{5}{G}.}$$(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 µ_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.

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

Fig. 2 Evolution of the coagulation variable, χ, with increasing density, nH during the isothermal collapse. The solid line represents an initial density of ρ0 = 3.8 × 10−20 g cm−3 is and the dashed line represents ρ0 = 3.8 × 10−18 g cm−3.

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₀ can be described as (Flower et al. 2005): $$\frac{\text{d}x}{\text{d}t}=-\frac{\pi }{2{\tau }_{\text{ff}}}\sqrt{\frac{1}{x}-1,}$$(14)

where τ_ff is the free-fall time, given by: $${\tau }_{\text{ff}}=\sqrt{\frac{3\pi }{32G{\rho }_0},}$$(15)

where ρ₀ 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) = ρ₀(R₀/R[t])³. The gas density then evolves as: $$\frac{1}{n_{\text{H}}}\frac{\text{d}{n}_{\text{H}}}{\text{d}t}=\frac{1}{\rho }\frac{\text{d}\rho }{\text{d}t}=-\frac{3}{x}\frac{\text{d}x}{\text{d}t}.$$(16)

We used a second-order Runge-Kutta scheme to numerically integrate this equation from the beginning of the collapse at ρ = ρ₀ until the formation of the first Larson core at ρ = 10⁻¹³ g cm⁻³. The evolution of χ with ρ is plotted in Fig. 2, assuming T = 10 K. The solid and dashed lines represent different initial densities, ρ₀ = 3.8 × 10⁻²⁰ g cm⁻³ and ρ₀ = 3.8 × 10⁻¹⁸ g cm⁻³, respectively. At densities n_H > 10⁻¹⁶ g cm⁻³, the value of χ is independent from ρ₀ and increases as χ ∝ ρ^1/4. This evolution is expected as χ ~ ρ^3/4t and t ~ τ_ff ~ ρ^−1/2. At ρ = 10⁻¹³ g cm⁻³, the coagulation variable reaches χ ≈ 5.6 × 10¹⁷ 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⁻¹³ g cm⁻³, 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.

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.

3.2.1 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 isother-mally collapsing envelope for ρ < 10¹⁵ g cm⁻³. 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.

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⁻¹⁶ g cm⁻³). 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³ 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 µ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.

Fig. 4 Evolution of the coagulation variable χ with increasing density nH 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.

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 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²/(ρη_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⁴ 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 ρ > 10⁻¹² g cm⁻³, Lebreuilly et al. (2023) showed that the replenishment of small grains by fragmentation would lead to an increase in both Ohmic and ambipolar resistivities.

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. 6 Slices of the C-3 simulation with coagulation showing the mode of the grain size distribution amax, 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 amax = 1, 5 and 10 µm. Bottom panel: amax = 0.5 and 1 µm.

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 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 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 cŕrcumstellar disks could predict an early onset for this process.

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.

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.

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

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.

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¹³ cm⁻³ (see their Fig. 7), while in our simulations, the peak exceeds 100 µm at a lower maximum density of n_H ≈ 3 × 10¹² cm⁻³ (or ρ ≈ 10⁻¹¹ g cm⁻³; 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⁸, while we used (Ormel et al. 2009; Guillet et al. 2020) $$\mathrm{Re}=6.7\times {10}^7{\left(\frac{n_{\text{H}}}{{10}^5{\text{cm}}^{-3}}\right)}^{\frac{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 pseudo-disk is an over-density, typically denser than ρ ≈ 10⁻¹⁵ g cm⁻³, 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 over-estimation 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 (ρ > 10⁻¹⁵ g cm⁻³) 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³ 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.

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

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-{ϵ}^2{\text{Θ}}^2}{1+ϵ\text{Θ}{\alpha }_k+{ϵ}^2{\text{Θ}}^2},$$(A.1) $$\langle \tilde{J}\left({\tau }_k\right)\rangle =\left(1-\psi \right)+\frac{\frac{2}{{\tau }_k}\left[{ϵ}^2{\text{Θ}}^2+ϵ\text{Θ}\right]}{1+ϵ\text{Θ}{\alpha }_k+{ϵ}^2{\text{Θ}}^2},$$(A.2) $$ϵ=\frac{1-\psi }{{\text{Θe}}^{\psi }},$$(A.3) $$n_i=-\frac{1}{1-ϵ}{\displaystyle \sum _k{n}_k{Z}_k,}$$(A.4) $$f\left(\psi \right)=\frac{{\langle \sigma \upsilon \rangle }_{\text{ie}}ϵn_{\text{i}}^2}{\zeta n_{\text{H}}}+\frac{n_{\text{i}}{\upsilon }_{\text{i}}}{\zeta n_{\text{H}}}{\displaystyle \sum _k{n}_k\pi a_k^2\langle \tilde{J}\left({\tau }_k\right)\rangle -1=0.}$$(A.5)

Here, ϵ = n_e/n_i < 1 is the ratio between the number of electrons and ions, ψ 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, ζ is the CR ionization rate, and v_i = [8k_BT/πµ_im_H]^1/2 is the thermal velocity of ions, with T the temperature, µ_i the average atomic mass of ions, and m_H the mass of the proton. We also have Θ = s_e[µ_im_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}\left({\tau }_k\right)\rangle$ represents the enhancement factor for ion recombination on grains, τ_k = a_kk_BT/e² is the reduced temperature of grains (Draine & Sutin 1987), and α_k = [8/(πτ_k)]^1/2. Finally, 〈σv〉_ie = 2 × 10⁻⁷[T/300]^−1/2 is the collision rate between ions and electrons.

We solve the system of relevant equations (A.1-A.5) for ψ 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 }_s{}_{{\text{H}}_2}=\frac{1}{{\text{a}}_s\text{He}}\frac{m_s+m_{{\text{H}}_2}}{m_{{\text{H}}_2}}\frac{1}{n_{{\text{H}}_2}{\langle \sigma \omega \rangle }_s}.$$(A.6)

Here, a_sHe 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_{{\text{H}}_2}$ are the mass of the species, s, and the H₂ molecule; $n_{{\text{H}}_2}$ represents the number density of H₂ (roughly equal to the density of the gas); 〈σω〉_s is the collision rate, taken from Pinto & Galli (2009) for electrons and ions: $${\langle \sigma \omega \rangle }_{\text{e}}=3.16\times {10}^{-11}{\upsilon }_{\text{rms,e}}^{1.3},$$(A.7) $${\langle \sigma \omega \rangle }_{\text{i}}=2.4\times {10}^{-9}{\upsilon }_{\text{rms,i}}^{0.6},$$(A.8)

where $${\upsilon }_{\text{rms,s}}={\left(\frac{8k_{\text{B}}T}{\pi {\mu }_{s,{\text{H}}_2}}\right)}^{\frac{1}{2}},$$(A.9)

is in km s⁻¹, with ${\mu }_{s,{\text{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={\left(\frac{8k_{\text{B}}T}{\pi m_{{\text{H}}_2}}\right)}^{\frac{1}{2}}\left(1+{\left(\frac{\pi }{2{\tau }_k}\right)}^{\frac{1}{2}}\right).$$(A.10)

We also define the conductivity, σ_s, and cyclotron frequency, ω_s, of species s: $${\sigma }_s=\frac{n_s{q}_s^2{\tau }_{s{\text{H}}_2}}{m_s},$$(A.11) $${\omega }_s=\frac{q_{s}B}{c{m}_s},$$(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 }_{\Vert }={\displaystyle \sum _s{\sigma }_s,}$$(A.13) $${\sigma }_{\perp }={\displaystyle \sum _s\frac{{\sigma }_s}{1+\left({\omega }_s{\tau }_{s{\text{H}}_2}\right)},}$$(A.14) $${\sigma }_{\text{H}}={\displaystyle \sum _s\frac{{\sigma }_s{\omega }_s{\tau }_{s{\text{H}}_2}}{1+{\left({\omega }_s{\tau }_{s{\text{H}}_2}\right)}^2}.}$$(A.15), (A.16)

Finally, the Ohmic, Hall, and ambipolar resisitivities are defined as: $${\eta }_{\text{O}}=\frac{1}{{\sigma }_{\parallel }},$$(A.17) $${\eta }_{\text{H}}=\frac{{\sigma }_{\text{H}}}{{\sigma }_{\perp }^2+{\sigma }_{\text{H}}^2},$$(A.18) $${\eta }_{\text{A}}=\frac{{\sigma }_{\text{⊥}}}{{\sigma }_{\perp }^2+{\sigma }_{\text{H}}^2}-\frac{1}{{\sigma }_{\parallel }}.$$(A.19), (A.20)

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 χ, is (Paper I, Equation 5): $$\frac{\text{d}X\left(a,t\right)}{\text{d}t}=C_{g_{\text{local}}{n}_{\text{H}}}I\left(a,X,t\right),$$(B.1)

where C is a constant, g_local is a function of the local properties of the gas (density, temperature…), and $$\begin{array}{l}I\left(a,X,t\right)=-{\displaystyle {\int }_0^{\infty }h\left(m,m^{\prime }\right)X\left(m,t\right)X\left(m^{\prime },t\right)\text{d}{m}^{\prime }}\hfill \\ +\frac{1}{2}{\displaystyle {\int }_0^{m}h\left(m-m^{\prime },m^{\prime }\right)X\left(m-m^{\prime },t\right)X\left(m^{\prime },t\right)\text{d}{m}^{\prime }}.\hfill \end{array}$$(B.2)

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_gx, where d_g is the dust-to-gas ratio, and x is the normalized relative number of grains. We can then rewrite Equation B.1: $$\frac{\text{d}x\left(a,t\right)}{\text{d}t}=C_{g_{\text{local}}{n}_{\text{H}}}{d}_g{I}^{\prime }\left(a,x,t\right),$$(B.3)

and include d_g in the definition of our coagulation variable, giving it the following form: $$\text{d}{\chi }^{\prime }=g_{\text{local}}{n}_{\text{H}}{d}_g\text{d}t.$$(B.4)

In our case, for the Ormel & Cuzzi (2007) kernel, $$\text{d}{\chi }^{\prime }=n_{\text{H}}^{\frac{3}{4}}{T}^{-\frac{1}{4}}{d}_g\text{d}t.$$(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 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: $$E_{\text{kin}}>5N_{\text{tot}}{E}_{\text{br}},$$(C.1)

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_{\text{br}}=A_{\text{br}}{\gamma }_{\text{grain}}^{5/3}\frac{{\left(a_0/2\right)}^{4/3}}{{\epsilon }^{2/3}},$$(C.2)

with A_br = 2.8 × 10³, γ_grain as the surface energy density of the material, a₀ 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.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⁻² and ε = 2.8 × 10¹¹ dyn cm⁻². The kinetic energy of two grains of mass m and m′ is expressed as: $$E_{\text{kin}}=\frac{1}{2}\frac{m{m}^{\prime }}{m+m^{\prime }}\Delta {\upsilon }^2,$$(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 \upsilon ={\left(\frac{3}{\sqrt{8}}{z}_0{\left[k_{\text{B}}G\right]}^{\frac{1}{2}}\frac{\gamma {\rho }_{\text{s}}}{\mu m_{\text{H}}}\right)}^{\frac{1}{2}}{n}_{\text{H}}^{-\frac{1}{4}}{T}^{-\frac{1}{4}}{a}^{\frac{1}{2}},$$(C.4)

where z₀ = 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⁻³ 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 $$m=\frac{4}{3}\pi {\rho }_{\text{s}}{a}_0^3{N}_{\text{tot}},$$(C.5)

Equation C.1 becomes $$a>\frac{15E_{\text{br}}}{\pi {\rho }_{\text{s}}{a}_0^3}{\left(\frac{3}{\sqrt{8}}{z}_0{\left[k_{\text{B}}G\right]}^{\frac{1}{2}}\frac{\gamma {\rho }_{\text{s}}}{\mu m_{\text{H}}}\right)}^{-1}{n}_{\text{H}}^{\frac{1}{2}}{T}^{\frac{1}{2}}.$$(C.6)

Replacing the variables with numeric values, we find $$a>826\mu \text{m}{\left(\frac{n_{\text{H}}}{{10}^{10}{\text{cm}}^{-3}}\right)}^{\frac{1}{2}}{\left(\frac{T}{10\text{K}}\right)}^{\frac{1}{2}}.$$(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 $E_{\text{br}}~{\gamma }_{\text{grain}}^{5/3}$). The value derived in Equation C.7 should then be considered very conservative.

References

All Tables

All Figures

	Fig. 1 State of the grain size distribution for values of χ, between 1015 cgs and 1019 cgs. The points represent the fractional abundance of the size-bin as function of the effective radius of the bin.
In the text

	Fig. 2 Evolution of the coagulation variable, χ, with increasing density, nH during the isothermal collapse. The solid line represents an initial density of ρ0 = 3.8 × 10−20 g cm−3 is and the dashed line represents ρ0 = 3.8 × 10−18 g cm−3.
In the text

	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.
In the text

	Fig. 4 Evolution of the coagulation variable χ with increasing density nH 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.
In the text

	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.
In the text

	Fig. 6 Slices of the C-3 simulation with coagulation showing the mode of the grain size distribution amax, 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 amax = 1, 5 and 10 µm. Bottom panel: amax = 0.5 and 1 µm.
In the text

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.

In the text

	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.
In the text

	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.
In the text