Issue 
A&A
Volume 670, February 2023



Article Number  A61  
Number of page(s)  11  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202244291  
Published online  08 February 2023 
Fast methods for tracking grain coagulation and ionization
III. Protostellar collapse with nonideal MHD
^{1}
Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, UT3PS, CNRS, CNES,
9 av. du Colonel Roche,
31028
Toulouse Cedex 4, France
^{2}
Department of Astrophysics, American Museum of Natural History,
200 Central Park West,
New York, NY
10024, USA
email: pierre.marchand.astr@gmail.com
^{3}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot,
Sorbonne Paris Cité,
91191
GifsurYvette, France
^{4}
Université ParisSaclay, CNRS, Institut d’Astrophysique spatiale,
rue Jean Teillac,
91405
Orsay, France
^{5}
Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3,
CC 72,
Place Eugène Bataillon,
34095
Montpellier Cedex 5, France
Received:
17
June
2022
Accepted:
23
December
2022
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 nonideal magnetohydrodynamics effects selfconsistently and at low numerical cost. We find that the coagulation efficiency is mostly affected by the time spent in highdensity regions. Starting from submicron 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)
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Grains play a major role during star formation. First, they are the seeds of planet formation. While their characteristic size is submicron 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 0I 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 nonideal 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). Nonideal 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 selfconsistently 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 preprocessed or postprocessed with no selfconsistent 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 selfconsistent way that is particularly suited for modeling star formation. Here, we apply this method to nonideal 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 nonideal MHD resistivities. These 3D simulations are the first to include a selfconsistent grain growth with direct feedback to the dynamics through the selfconsistent 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 nonideal MHD simulations with the RAMSES code (Teyssier 2002). RAMSES is an Eulerian gas dynamics code with adaptive mesh refinement (AMR) and selfgravity. It includes a monofluid treatment of nonideal 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 onthefly in a selfconsistent 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 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: (1)
where n_{H} is the number density of the gas, T its temperature, and t the time. In threedimensional (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: (2)
where u is the velocity of the gas. We can combine Eq. 2 with the mass conservation equation: (3)
where ρ is the gas mass density. This yields: (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 selfconsistently in all cells at each timestep, as a massweighted average, ignoring any diffusion of dust through the gas. We used the Ishinisan code (Marchand et al. 2021) to precalculate 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: (5)
The total quantity of grains is determined by the dusttogas 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.
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 sizebin 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 gasgrain 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 cosmicray (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 NewtonRaphson 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.
Parameters of the protostellar collapse simulations.
2.2 Star formation simulations
2.2.1 Model
We performed four numerical simulations using the RAMSES code. We solved the following MHD equations: (6) (7) (8) (9)
where u is the velocity of the gas, P its pressure, B the magnetic field, J = ∇ × B the current, 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: (10)
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 (11)
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 masstoflux ratio over the critical value (Mouschovias & Spitzer 1976): (12)
Observations show that dense cores are slightly supercritical (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 C3. We changed the value of α to 0.4, and the initial mass, M, to 5 M_{⊙} for two other cases referred to as C4 and C3M5, respectively, to investigate the influence of the collapse time on the grain coagulation. The final run, named NC3, is similar to C3 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 C3 and NC3 (resolution of 1.4 au), 14 for C4 (resolution of 0.96 au), and 16 for C3M5 (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 nonmagnetic 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 onezone model in Sect. 3.1, then in 3D MHD simulations in Sect. 3.2.
Fig. 2 Evolution of the coagulation variable, χ, with increasing density, n_{H} 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_{0} can be described as (Flower et al. 2005): (14)
where τ_{ff} is the freefall time, given by: (15)
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: (16)
We used a secondorder RungeKutta 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.
3.2 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 C3, a small circumstellar disk with a radius of ≈20 au forms and a disk wind is launched by magnetocentrifugal acceleration (Blandford & Payne 1982). Figure 3 shows faceon and edgeon slices of density at the final timestep of the simulation, with arrows indicating the gas velocity.
Fig. 3 Density slices of the C3 simulation, with grain coagulation. Top panel is a faceon 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 timestep, 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 timestep, for runs C3, C4, and C3M5. 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 pseudodisk, 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 C4 needing 50% more time to collapse to the first Larson core stage, and C3M5 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 pseudodisk and the early protoplanetary disk, where grains grow by a factor 100 from submicron 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 C3, C4, and C3M5, 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 C3 is displayed in Fig. 6. Size distributions shifting significantly from the initial MRN distribution are indeed confined to the midplane in the disk and pseudodisk. The bottom panel also shows moderately larger grains in the outflow, as they only traveled through the upper layers of the pseudodisk before being ejected (Marchand et al. 2020). If they had been in the midplane 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 n_{H} in the numerical collapse models C3 (purple), C3M5 (light blue), and C4 (green). Each point corresponds to a simulation cell at the final timestep. 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 nonideal 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 volumeweighted average resistivities as a function of density, at the final timestep, for simulations C3 and NC3. A nonevolving 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 C3 and NC3, azimuthally averaged in the midplane. The higher resistivity in run NC3 results in Am < 1 in the inner ~12 au, indicating active ambipolar diffusion and magnetic field dissipation, while Am ≳ 10^{4} for C3 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 NC3. A second effect of the weaker magnetic forces from the stronger ambipolar diffusion in run NC3 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 largemassend 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 truncatedMRN 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 nonideal 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.
Fig. 5 Mode of the coagulated grain size distribution as a function of density in simulations (purple), C3M5 (light blue), and C4 (green). The discrete values of the sizes are due to the binning of the size distribution. 
Fig. 6 Slices of the C3 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 µm. Bottom panel: a_{max} = 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 Nonideal MHD effects in the C3 and NC3 simulations. Top panel: abundances of ions (purple) and electrons (green) as a function of density in models with coagulation (C3; solid lines) or without it (NC3; dashed lines). Middle panel: volumeaveraged Ohmic (purple), ambipolar (green), and Hall (blue) resistivities for the same models. Bottom panel: average ambipolar Elsasser number Am in the midplane as a function of radius. The thin grey line represents Am = 1. 
Fig. 8 Density slices of the NC3 simulation, without grain coagulation. Top panel is a faceon 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 timestep, 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 mmsize grains in low numbers. The timescale 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 dusttogas 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 millimetersize) 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 largeend of the size distribution. In our case, the disk wind is fueled by the upper layers of the pseudodisk, in which dust grows only moderately. However, dusttogas 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 C3 with coagulation, and dashed lines are simulation NC3 without coagulation. 
4.3 Coagulation in the pseudodisk
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 used (Ormel et al. 2009; Guillet et al. 2020) (17)
Hence, in the central regions, our Reynolds number can be larger by three orders of magnitude. Their grains are then stuck in the tightlycoupled regime where relative grain velocities are lower than in the intermediate couplingregime, 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 pseudodisk that forms in our simulations. The pseudodisk is an overdensity, 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 pseudodisk and the overall efficiency of coagulation is mainly affected by the time spent in highdensity regions such as the pseudodisk. This is what appears in the bottom panel of Fig. 6, where there is a ~30 authick layer of large grains in the midplane. 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 lesscoagulated 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 wellknown 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 steadystate Kolmogorov turbulence spectrum and that the injectionscale 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 selfconsistent calculation of nonideal 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 freefall times. What matters is the time spent in highdensity regions (ρ > 10^{−15} g cm^{−3}) such as the pseudodisk. Fragmentation can also be ignored in the cloudcollapse 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 nonideal MHD resistivities. Coagulated grains result in resistivities up to four orders of magnitude lower than noncoagulated 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 AST1815461. 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): (A.1) (A.2) (A.3) (A.4) (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_{B}T/πµ_{i}m_{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}[µ_{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, represents the enhancement factor for ion recombination on grains, τ_{k} = a_{k}k_{B}T/e^{2} is the reduced temperature of grains (Draine & Sutin 1987), and α_{k} = [8/(πτ_{k})]^{1/2}. Finally, 〈σv〉_{ie} = 2 × 10^{−7}[T/300]^{−1/2} is the collision rate between ions and electrons.
We solve the system of relevant equations (A.1A.5) for ψ and find n_{i}, n_{e}, and Z_{k} for all bins. A NewtonRaphson algorithm is described in details in the Appendix A of Paper I.
A.2 The resistivities
The collision timescale of species s with neutral H_{2}, the most abundant species in the gas, is given by: (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 are the mass of the species, s, and the H_{2} molecule; represents the number density of H_{2} (roughly equal to the density of the gas); 〈σω〉_{s} is the collision rate, taken from Pinto & Galli (2009) for electrons and ions: (A.7) (A.8)
is in km s^{−1}, with the reduced mass of species s and H_{2}. These velocities have been calculated in a threefluid 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): (A.10)
We also define the conductivity, σ_{s}, and cyclotron frequency, ω_{s}, of species s: (A.11) (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: (A.13) (A.14) (A.15), (A.16)
Finally, the Ohmic, Hall, and ambipolar resisitivities are defined as: (A.17) (A.18) (A.19), (A.20)
Appendix B Coagulation and dusttogas ratio
Our coagulation model assumes that grains are perfectly coupled with the gas and wellmixed, so that the dusttogas 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 dusttogas mass ratio during the protostellar collapse. The disk and first core tend to be dustenriched while the envelope becomes dustdepleted. 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): (B.1)
where C is a constant, g_{local} is a function of the local properties of the gas (density, temperature…), and (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_{g}x, where d_{g} is the dusttogas ratio, and x is the normalized relative number of grains. We can then rewrite Equation B.1: (B.3)
and include d_{g} in the definition of our coagulation variable, giving it the following form: (B.4)
In our case, for the Ormel & Cuzzi (2007) kernel, (B.5)
This expression means that the dusttogas ratio does not change the evolution path of the size distribution, only the coagulation speed. A dusttogas 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 dusttogas 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 selfconsistent treatment of dust grains in hydrodynamics simulations. This method will lead to an overestimate of the smalltolarger grain ratio, but this is probably a decent approximation to get the total dusttogas 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: (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): (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: (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: (C.4)
where 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 (C.5)
Equation C.1 becomes (C.6)
Replacing the variables with numeric values, we find (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 ). The value derived in Equation C.7 should then be considered very conservative.
References
 Balsara, D. S. 2012, J. Comput. Phys., 231, 7476 [CrossRef] [MathSciNet] [Google Scholar]
 Bate, M. R. 2022, MNRAS, 514, 2145 [CrossRef] [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Crutcher, R. M. 1999, ApJ, 520, 706 [NASA ADS] [CrossRef] [Google Scholar]
 Desch, S. J., & Mouschovias, T. C. 2001, ApJ, 550, 314 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Tielens, A. G. G. M. 1997, ApJ, 480, 647 [Google Scholar]
 Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Duffin, D. F., & Pudritz, R. E. 2008, MNRAS, 391, 1659 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [CrossRef] [EDP Sciences] [Google Scholar]
 Epstein, P. S. 1924, Phys. Rev., 23, 710 [Google Scholar]
 Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2005, A&A, 436, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Galli, D., & Shu, F. H. 1993, ApJ, 417, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Guillet, T., & Teyssier, R. 2011, J. Comput. Phys., 230, 4756 [Google Scholar]
 Guillet, V., Hennebelle, P., Pineau des Forêts, G., et al. 2020, A&A, 643, A17 [EDP Sciences] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [Google Scholar]
 Kawasaki, Y., Koga, S., & Machida, M. N. 2022, MNRAS, 515, 2072 [NASA ADS] [CrossRef] [Google Scholar]
 Kunz, M. W., & Mouschovias, T. C. 2009, ApJ, 693, 1895 [Google Scholar]
 Kuznetsova, A., Hartmann, L., & Heitsch, F. 2020, ApJ, 893, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Kwon, W., Looney, L. W., Mundy, L. G., Chiang, H.F., & Kemball, A. J. 2009, ApJ, 696, 841 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [Google Scholar]
 Lebreuilly, U., Commerçon, B., & Laibe, G. 2019, A&A, 626, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lebreuilly, U., Commerçon, B., & Laibe, G. 2020, A&A, 641, A112 [EDP Sciences] [Google Scholar]
 Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Lebreuilly, U., VallucciGoy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326 [Google Scholar]
 Le Gouellec, V. J. M., Hull, C. L. H., Maury, A. J., et al. 2019, ApJ, 885, 106 [NASA ADS] [CrossRef] [Google Scholar]
 Li, Z.Y., Krasnopolsky, R., & Shang, H. 2011, ApJ, 738, 180 [Google Scholar]
 Lombart, M., & Laibe, G. 2021, MNRAS, 501, 4298 [Google Scholar]
 Machida, M. N., Inutsuka, S., & Matsumoto, T. 2006, ApJ, 647, L151 [Google Scholar]
 Marchand, P., Masson, J., Chabrier, G., et al. 2016, A&A, 592, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchand, P., Commerçon, B., & Chabrier, G. 2018, A&A, 619, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchand, P., Tomida, K., Tanaka, K. E. I., Commerçon, B., & Chabrier, G. 2020, ApJ, 900, 180 [CrossRef] [Google Scholar]
 Marchand, P., Guillet, V., Lebreuilly, U., & Mac Low, M. M. 2021, A&A, 649, A50 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masson, J., Teyssier, R., MuletMarquis, C., Hennebelle, P., & Chabrier, G. 2012, Astrophys. J. Suppl. Ser., 201, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, J., Chabrier, G., Hennebelle, P., Vaytet, N., & Commercon, B. 2016, A&A, 587, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [Google Scholar]
 Mellon, R. R., & Li, Z. 2009, ApJ, 698, 922 [Google Scholar]
 Miotello, A., Testi, L., Lodato, G., et al. 2014, A&A, 567, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204 [NASA ADS] [CrossRef] [Google Scholar]
 Mouschovias, T. C., & Spitzer, L., Jr. 1976, ApJ, 210, 326 [NASA ADS] [CrossRef] [Google Scholar]
 Ormel, C. W., & Cuzzi, J. N. 2007, A&A, 466, 413 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. T. 2009, A&A, 502, 845 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pinto, C., & Galli, D. 2009, ThreeFluid Magnetohydrodynamics in Star Formation, eds. K. Tsinganos, T. Ray, & M. Stute, 597 [Google Scholar]
 Rossi, S. C. F., BenevidesSoares, P., & Barbuy, B. B. 1991, A&A, 251, 587 [Google Scholar]
 Silsbee, K., Ivlev, A. V., Sipilä, O., Caselli, P., & Zhao, B. 2020, A&A, 641, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Silsbee, K., Akimkin, V., Ivlev, A. V., et al. 2022, ApJ, 940, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Smoluchowski, M. V. 1916, Zeitsch. Physik, 17, 557 [NASA ADS] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
 Tomida, K., Okuzumi, S., & Machida, M. N. 2015, ApJ, 801, 117 [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [CrossRef] [Google Scholar]
 Tsukamoto, Y., Machida, M. N., & Inutsuka, S.I. 2021, ApJ, 920, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Valdivia, V., Maury, A., Brauer, R., et al. 2019, MNRAS, 488, 4897 [NASA ADS] [CrossRef] [Google Scholar]
 Vaytet, N., & Haugbølle, T. 2017, A&A, 598, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vaytet, N., Commerçon, B., Masson, J., González, M., & Chabrier, G. 2018, A&A, 615, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vericel, A., Gonzalez, J.F., Price, D. J., Laibe, G., & Pinte, C. 2021, MNRAS, 507, 2318 [NASA ADS] [CrossRef] [Google Scholar]
 Vorobyov, E. I., Skliarevskii, A. M., Molyarova, T., et al. 2022, A&A, 658, A191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wurster, J., Price, D. J., & Bate, M. R. 2016, MNRAS, 457, 1037 [Google Scholar]
 Yang, C.C., Johansen, A., & Carrera, D. 2017, A&A, 606, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
 Zhao, B., Caselli, P., Li, Z.Y., et al. 2016, MNRAS, 460, 2050 [Google Scholar]
 Zhao, B., Caselli, P., Li, Z.Y., & Krasnopolsky, R. 2018, MNRAS, 473, 4868 [Google Scholar]
All Tables
All Figures
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 sizebin as function of the effective radius of the bin. 

In the text 
Fig. 2 Evolution of the coagulation variable, χ, with increasing density, n_{H} 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 C3 simulation, with grain coagulation. Top panel is a faceon 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 timestep, 1000 yr after the formation of the first Larson core. 

In the text 
Fig. 4 Evolution of the coagulation variable χ with increasing density n_{H} in the numerical collapse models C3 (purple), C3M5 (light blue), and C4 (green). Each point corresponds to a simulation cell at the final timestep. 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), C3M5 (light blue), and C4 (green). The discrete values of the sizes are due to the binning of the size distribution. 

In the text 
Fig. 6 Slices of the C3 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 µm. Bottom panel: a_{max} = 0.5 and 1 µm. 

In the text 
Fig. 7 Nonideal MHD effects in the C3 and NC3 simulations. Top panel: abundances of ions (purple) and electrons (green) as a function of density in models with coagulation (C3; solid lines) or without it (NC3; dashed lines). Middle panel: volumeaveraged Ohmic (purple), ambipolar (green), and Hall (blue) resistivities for the same models. Bottom panel: average ambipolar Elsasser number Am in the midplane as a function of radius. The thin grey line represents Am = 1. 

In the text 
Fig. 8 Density slices of the NC3 simulation, without grain coagulation. Top panel is a faceon 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 timestep, 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 C3 with coagulation, and dashed lines are simulation NC3 without coagulation. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.