© The Authors 2026

1 Introduction

The formation and early evolution of dust grains in young stellar systems represent a critical yet poorly understood component of star and planet formation processes. While dust grains constitute only ~1% of the mass of the interstellar medium (ISM), they fundamentally influence protostellar evolution through their roles in radiative cooling, chemistry catalysis, and magnetic field coupling (Draine 2003; Zhao et al. 2016; Lebreuilly et al. 2023b; Vallucci-Goy et al. 2024).

Observations of the dust thermal emission at different frequencies allow us to measure the dust emissivity, β, a quantity that is key to computing dust masses from dust emission but is also widely used to measure the mean dust grain size in astrophysical environments and discuss, for example, dust evolution from submicronic solid particles toward pebbles and the formation of planetary systems (Testi et al. 2014; Birnstiel 2024). Some pilot observational studies suggested that young (<1 Myr old) disks show a dust emissivity evolution that is consistent with an increase in the mean dust grain size (Kwon et al. 2009; Miotello et al. 2014). More recently, observations have revealed compelling evidence that dust evolution is more complex than expected and may begin earlier than previously assumed (maybe even starting during the pre-stellar phase), which challenges traditional models that assume significant grain growth only occurred in disks. For example, Dartois et al. (2024) modeled the CO₂ ice band observed by JWST in two lines of sight in a dense cloud (N(H) of ~1–2 × 10²³ cm⁻², A_v > 60). They interpret the extended blue wing of this band as an excess of scattered light and conclude that, to explain this spectral feature, a fair fraction of the dust mass must be in grains of ~1 micrometer in size. At longer wavelengths, observational studies of the Taurus B211/B213 star-forming filament have shown a systematic decrease in the dust emissivity index (β) from β ≈ 2 in pre-stellar cores at scales of ~10 000 au to β < 1 toward protostellar envelopes (Ysard et al. 2013; Bracco et al. 2017). Low dust emissivities, β, have also been observed in the inner envelopes of embedded protostars in high-angular-resolution studies of the millimeter dust emission (Galametz et al. 2019; Cacciapuoti et al. 2023). All together, these studies reveal a variety of dust optical properties. However, the lowest dust emissivity (β) values have only so far been reproduced by including large dust grains (>100 μm), which suggests that efficient grain aggregation could begin in dense cores and envelopes long before dust grains have entered the disk environment.

Improving our modeling of the dust size distribution during the star-formation sequence requires us to address the magnetohydrodynamics (MHD) of gas and dust mixtures simultaneously with a full treatment of dust evolution. State-of-the-art dust evolution models employ distinct numerical strategies to balance computational feasibility with physical fidelity. Three-dimensional simulations are required to capture asymmetries in the evolution of the collapse. However, they are challenging and therefore often adopt a monodisperse or bi-disperse¹ approximation to resolve dust growth (Vorobyov & Elbakyan 2019; Elbakyan et al. 2020; Tsukamoto et al. 2021; Vorobyov et al. 2024, 2025) or do not include dust-gas differential dynamics (Bate 2022; Marchand et al. 2023). We note the very recent effort of Bate et al. (2025) to couple dust dynamics with a full treatment of dust evolution with the Smoluchowski equation. They employed methods for collapse calculation that resemble those presented in this article but using Smooth Particle Hydrodynamic (SPH) methods. Conversely, several codes have been developed to solve the Smoluchowski equation across a wide span of grains size with the gas and dust dynamics in 1D, such as DustPy (Stammler & Birnstiel 2022) and shark (Lebreuilly et al. 2023b; Vallucci-Goy et al. 2024), or in 2D, such as LA-COMPASS (Drążkowska et al. 2019; Li et al. 2020; Ho et al. 2024) and CUDISC (Robinson et al. 2024), with a significant numerical cost. Recent efforts like the TriPoD framework (Pfeil et al. 2024) aim to reduce this numerical cost by embedding parameterized size distributions in 2D hydrodynamic grids, although challenges persist in reconciling growth timescales with dynamical environments. However, incorporating accurate numerical solutions to the Smoluchowski equation in 3D multi-physics hydrodynamic simulations is out of reach with the mentioned algorithms. This numerical challenge requires a new algorithm to efficiently solve the Smoluchowski equation. Thanks to the newly developed COALA code, it is now possible to fully solve the Smoluchowski equation for a wide range of grain sizes on the fly in 3D simulations that also follow dust dynamics. In the present study, by coupling COALA to the adaptive mesh-refinement code RAMSES, we ran simulations beyond the current state of the art that follow nonideal MHD protostellar collapses with a full solving of the Smoluchowski equation together with a treatment of dust dynamics in the terminal velocity approximation.

2 Methods

In this work, we performed numerical 3D MHD simulations of protostellar collapses using the RAMSES code (Teyssier 2002). In particular, we employed its extension to nonideal MHD (Masson et al. 2012). We followed the evolution of a gravitationally unstable pre-stellar core of 2.5 M_⊙ with an initial radius of ~8000 au² and an initial turbulent Mach 0 or 1 (denoted ℳ₀ and ℳ₁). Initial conditions were similar to those in the models presented in Maury et al. (2019) and Valdivia et al. (2022). The temperature of the cloud was computed assuming a barotropic equation of state as per Lebreuilly et al. (2020), which assumes isothermality below 10⁻¹³ g cm⁻³ and adiabaticity with an adiabatic index of 5/3 above. We used the adaptive-mesh refinement capability of RAMSES to follow the gravitational collapse. The coarse grid had a resolution of 1100 au and was refined by imposing at least ten points per Jeans length to avoid artificial fragmentation (Truelove et al. 1997); it reached a minimum cell size of 4.3 au in the densest regions.

2.1 Dust growth with COALA

We incorporated a full treatment of the dust coagulation by solving the Smoluchowski equation for single-composition spherical grains in the conservative form on the fly (Tanaka et al. 1996) using the COALA code (Lombart & Laibe 2021). Assuming a mass range [m_min, m_max] for the dust grains, the equation is written as $$\begin{array}{rl}& \frac{\mathrm{\partial }g(m,t)}{\mathrm{\partial }t}+\frac{\mathrm{\partial }{F}_{\text{growth }}[g](m,t)}{\mathrm{\partial }m}=0,\\ & F_{\text{growth}}[g](m,t)\\ & \equiv {\int }_{m_{min}}^m{\int }_{m-m_1+m_{min}}^{m_{max}-m_1+m_{min}}\frac{K(m_1,m_2)}{m_2}g(m_1,t)g(m_2,t)\mathrm{d}{m}_2\mathrm{d}{m}_1,\end{array}$$(1)

where m ∈ [m_min, m_max] is the mass of dust grains, g(m, t) is the dust mass distribution, and K(m₁, m₂) is the collision kernel between the two grains of masses m₁ and m₂. We used the ballistic kernel K(m₁, m₂) = σ(m₁, m₂) Δv(m₁, m₂), where σ is the collisional cross section and Δv is the differential velocity between grains. In this study, the sources of grain-grain differential velocities are Brownian motion, hydrodynamical drift, and gas turbulence (see details in Appendix A). We used the high accuracy of the COALA code to follow the evolution of the dust size distribution with a low number of dust size bins. The time solver implemented in COALA is a strong stability preserving Runge-Kutta third order (Gottlieb et al. 2009). In each spatial cell, we used 𝒩 > 1 dust species for which the set of dust mass density {ρ_k}_1≤k≤𝒩 defines the dust size distribution. The dust mass distribution (g) in Eq. (1) is linked to ρ_k such that $${\rho }_k={\int }_{I_k}g(m,t)\mathrm{d}m,$$(2)

with I_i the size of bin k. In this study, the code COALA gives the evolution of the mean value of the function g, noted $\overline{g}$, which means that the piecewise constant approximation is used. Therefore, the relation between ρ_k and ${\overline{g}}_k$ is directly obtained, ${\overline{g}}_k={\rho }_k/\mathrm{\Delta }{m}_k$. The accuracy of the approximation of g depends on the number of dust bins, 𝒩. In Appendix B, COALA is benchmarked with exact solutions in the setup given in Sect. 2.4.

2.2 Dust dynamics with RAMSES

The dust dynamics was treated by the terminal velocity approximation solver of Lebreuilly et al. (2019). In short, this approach considers dust species as separate phases in a monofluid composed of the full gas-dust mixture with density ρ and velocity v. Each dust phase (k) of density ρ_k required us to solve a mass conservation equation with a velocity v_k = v + w_k, where w_k is given by the formula $${\mathit{w}}_k=[\frac{\rho }{\rho -{\rho }_k}{t}_{\mathrm{s},k}-\sum _{l=1}^{\mathcal{N}}\frac{{\rho }_l}{\rho -{\rho }_l}{t}_{\mathrm{s},l}]\frac{\mathrm{\nabla }{P}_{\mathrm{t}\mathrm{h}}-(\mathrm{\nabla }\times \mathit{B})\times \mathit{B}}{\rho }.$$(3)

Here P_th is the gas thermal pressure, B is the magnetic field, and t_s,k is the grain stopping time. In the Epstein (1924) regimes, for a grain of size S_grain,k and intrinsic density ρ_grain,k, it is written as $$t_{\mathrm{s},k}\equiv \frac{{\rho }_{\text{grain},k}}{\rho }\frac{S_{\text{grain},k}}{w_{\text{th}}}.$$(4)

This regime is valid for small grains of the ISM provided that their drift velocity is low compared to the gas thermal speed, $$w_{\mathrm{t}\mathrm{h}}=\sqrt{\frac{8k_{\mathrm{B}}T}{\pi {\mu }_{\mathrm{g}}{m}_{\mathrm{H}}}}.$$

2.3 Coupling between RAMSES and COALA

The dust growth equation was added as a source term in the continuity equation for each dust fluid (see Lebreuilly et al. 2023b) as $$\{\begin{array}{l}\frac{\mathrm{\partial }{\rho }_k}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left[{\rho }_k{\mathit{v}}_k\right]={\mathrm{\Lambda }}_{k,\text{growth}},\\ {{\mathrm{\Lambda }}_{k,\text{growth}}\equiv \frac{\mathrm{\partial }{\rho }_k}{\mathrm{\partial }t}|}_{\text{growth}}\approx \frac{\mathrm{d}{\overline{g}}_k}{\mathrm{d}t}\mathrm{\Delta }{m}_k,\end{array}$$(5)

where g is the unknown variable in the Smoluchowski equation (see Sect. 2.1). Then, an ordinary operator splitting was applied to solve these continuity equations. First, the dynamic of the dust fluids was solved with the code RAMSES, and then the dust growth source term was solved by the code COALA. This is illustrated using an Euler forward time scheme, $$\begin{array}{rl}& {\rho }_k^{\ast }={\rho }_k\left(t_n\right)-\mathrm{\Delta }t(\mathrm{\nabla }\cdot \left[{\rho }_k{\mathit{v}}_k\right])\to \text{RAMSES},\\ & {\rho }_k\left(t_{n+1}\right)={\rho }_k^{\ast }+\mathrm{\Delta }t{\mathrm{\Lambda }}_{k,\text{growth}}\to \text{COALA},\end{array}$$(6)

where Δt is the hydrodynamic time-step. In this study, magnetic resistivities of grains do not evolve with changes in the dust size distribution.

2.4 Dust initial conditions

For this work, the binning was done logarithmically with 𝒩 = 40 between s_min = 5 nm and s_max = 1 nm. We set up the dust-to-gas ratio of each bin assuming a Mathis–Rumpl–Nordsieck (MRN)-like (Mathis et al. 1977) power-law distribution between s_min = 5 nm and s_cut~270 nm, which is the right interface of the bin N_cut = 11, including the grain size 250 nm. For bins with i > N_cut, we set an initial dust ratio with a very low value of 10⁻²⁰. For each bin (k) between 1 and N_cut, we have an initial dust ratio of $${ϵ}_{k,0}={ϵ}_0\left[\frac{S_{k,+}^{4+m}-S_{k,-}^{4+m}}{S_{+,{\mathrm{N}}_{\mathrm{c}\mathrm{u}\mathrm{t}}}^{4+m}-S_{-,1}^{4+m}}\right],$$(7)

assuming an MRN distribution with m = −3.5. The S_k,− and S_k,+ are the grain size at the left and right interface of bin k. As explained in Lebreuilly et al. (2019), the grain size that is used to dynamically evolve the bin (k) was taken to be $\sqrt{S_{k,-}{S}_{k,+}}$. Note that we assumed in this study that all grains have the same intrinsic density, ρ_grain = 2.3 g cm⁻³. In this work, all the terms “grain size” refer to the radius of grains.

3 Results

The RAMSES model of a collapsing protostellar envelope with dust growth computed by COALA presented here is the first 3D MHD model of dust evolution during the disk formation stage. This model self-consistently accounts for both the gas and dust dynamics (multiple dust sizes) and the dust size evolution modeled by the Smoluchowski equation.

To exclude the impact of the grid boundaries, we carried out the analysis of the dust evolution in the inner half of the initial core extent, i.e., in the inner envelope at radii <4000 au (see Fig. C.1). We ran two models with and without turbulence, referred to as ℳ₀ and ℳ₁, and followed their evolution until the protostar accreted one-half of the total core mass (0.85 M_⊙in the sink). From these two models, we selected a snapshot after the formation of the star when M_⋆ = 0.4 M_⊙, i.e., t_{0.4 M⊙} = 118.1 kyr for model ℳ₁ and t_0.4M⊙ = 102.4 kyr for model ℳ₀. This evolutionary stage is representative of a young embedded protostar. In the following, the disk region is defined as the region where the gas number density is greater than 10⁹ cm⁻³. For these snapshots, the disk radius is ≲80 au: thus, radii >100 au exclusively probe the envelope material.

Figure 1 shows models ℳ₀ and ℳ₁ taken when the central stellar object has accreted 0.4 M_⊙. The models are seen edge-on, with the disk rotation axis on the plane of the sky. The left column shows the map of the total {gas+dust} column density, Σ, computed as the density integrated along the line of sight. The middle column shows the map of the grain size, a_peak, at the peak of the size distribution (in each cell, the grain size, a_peak, is computed as the average size in the size bin that shows the largest mass fraction in the dust grain size distribution), averaged along the line of sight and weighted by the local gas density, ρ (see Appendix D). The right column shows the dust enrichment, defined as the ratio between the evolved dust ratio, ϵ = ρ_dust/(ρ_gas + ρ_dust), and the initial dust ratio, ϵ₀, averaged along the line of sight and weighted by the gas density along the line of sight (see Appendix D). The boxy structures around the star for model ℳ₀ only appear when integrating the density along the line of sight at scales greater than 1000 au, and are not seen in individual model slices. Moreover, they only appear in the model without turbulence. These structures are not numerical artifacts from the simulation but are due to a visualization effect. Therefore, these patterns are not a problem for the physical interpretation of the results.

Both models show significant dust evolution in the disks but also in the protostellar envelopes, even at early times. For example, for model ℳ₁, which accounts for both disk and envelope contributions in a sphere with radius 4000 au around the star, the grains with a size between 0.188 and 0.270 μm initially account for ~20% of the total dust mass, while at M_⋆ = 0.4 M_⊙, after the collapse starts, these grains have been depleted to produce larger grains and represent only ~0.1% of the dust mass in the disk and ~13% of the dust mass in the envelope. This depletion of the small grains is associated with an increase in the number of large, micron-sized grains. For example, for the dust bin containing grains of 5 μm, the dust mass enrichment is ~10% in the disk and ~0.6% in the envelope. We note that the denser the material is, the larger the grains. This is expected since as the material becomes denser, the collision rate between grains increases. As shown in Appendix E, if we assume a simple spherical collapse, we can show that the peak of the dust size distribution scales $\propto \sqrt{\rho }$.

Our work in 3D shows that the dust evolution is highly anisotropic in both models ℳ₀ and ℳ₁, with preferred locations for producing large dust grains associated with the development of higher density structures either due to a turbulent infall (ℳ₁) or due to the development of outflow cavity walls (ℳ₀). Figure F.1 shows the anisotropy of the dust evolution via a comparison of the simulation for model ℳ₁ with and without dust growth. Moreover, our models show that implementing the size evolution of the dust grains during the collapse amplifies local variations in the dust ratio. Figure 2 compares the evolution of the dust enrichment for model ℳ₁ with and without growth. The histograms are obtained with the dust enrichment calculated in each cell in a box with sides of 2000 au around the star. With dust growth, the dust enrichment in the envelope can become more than twice that in the disk region. However, the dust enrichment has a high value in only a few cells: the variance of the histogram has a low value of ~0.01 for model ℳ₁ when M_⋆ = 0.4 M_⊙. This low variance explains the small variation observed in the map in Figs. 1 and 2.

This does not come as a surprise since the Stokes number, the ratio between the stopping time and the free-fall time, is proportional to the grain size. As shown in Lebreuilly et al. (2020), larger Stokes numbers produce large dust-to-gas ratio variations. It is interesting to point out that the Stokes number is roughly constant because of the scaling of the grain size with the density. This explains why variations in the dust-to-gas ratio remain relatively small, even when incorporating dust growth. These variations essentially depend on the initial choice of grain size, and the Stokes number of the peak of the MRN distribution is ≪1. The strongest dust-to-gas ratio variations are in the regions with strong electric currents and are therefore not caused by pressure gradients. This effect was previously observed near strong current sheets in the models of Lebreuilly et al. (2023a, see the cartoon illustration in their Fig. 3). We note that, if the variations in the dust-to-gas ratio are small, in our case an instability could enhance them either by slowing the collapse (more rotation, more turbulence, stronger magnetic fields) or accelerating coagulation. The latter might be possible when considering other types of turbulence such as that of the Ormel kernel, as shown in Gong et al. (2020, 2021).

Figure 3 shows the global execution time as a function of the mass of the star with model ℳ₀ for three simulations: gas only, gas with 40 dust sizes, and gas with dust growth. By using COALA in RAMSES (orange circles), the computational cost is slightly higher, by a factor of 1.7, in the global execution time compared to the same simulation without dust growth (blue circles). The factor of 6 between the gas-only simulation and that for gas with 40 dust sizes is coherent with the computational time scaling as $\sqrt{\mathcal{N}}$ as shown in Lebreuilly et al. (2019). The next goal is to take advantage of the performances of COALA to reduce the number of dust sizes to ~20, which will reduce the computational cost.

Fig. 1 Snapshots taken from the RAMSES MHD models (ℳ0 and ℳ1) of a collapsing protostellar envelope (initial parameters can be found in the text) with dust grain size evolution computed by COALA. The maps show a 2000 au zoom in on the model when the central star has accreted 0.4 M⊙from the initial 2.5 M⊙core. The red cross represents the position of the star. The left column shows the total column gas density, Σ. The middle column shows the average grain size at the peak of the size distribution, apeak, weighted by the gas density along the line of sight (see the main text for further details). Similarly, the right column shows the average dust enrichment (ratio between the evolved dust-to-gas ratio, ϵ, and the initial dust-to-gas ratio, ϵ0) weighted by the gas density along the line of sight. In the top row, the square patterns near the sink are not numerical artifacts but result from the weighted average of the quantities along the line of sight.

Fig. 2 Snapshots and histograms of the dust enrichment at two given times of model ℳ1 with (top row) and without (bottom row) dust growth. The histograms are obtained with the dust enrichment calculated in each cell in a box with sides of 2000 au around the star. The two times are when M⋆ = 0.06 M⊙ (orange) and M⋆ = 0.4 M⊙ (blue). Dust growth and dynamics result in larger local variations in the dust enrichment, by a factor of greater than 2, in the envelope compared to the dust dynamics alone. The black line shows the histograms for the ℳ1 model with dust growth and 60 dust bins, highlighting the numerical convergence reached with 40-bin simulations. The red cross represents the position of the star. The variance value of ~0.01 of the histogram explains the low variation observed on the map for model ℳ1 when M⋆ = 0.4 M⊙.

Fig. 3 Comparison of the computational times of the three simulations for model ℳ0: gas only (black), gas with 40 dust sizes (blue), and gas with 40 dust sizes and dust growth (orange). Using COALA in RAMSES only increases the global execution time by a factor of 1.7.

4 Discussion

Figures 4a and 4c show the radial profiles of the dust grain size at the peak of the distribution, a_peak, in the envelope (100 to 4000 au) for the two models (ℳ₀ and ℳ₁) at t_{0.4 M⊙}. The radial profile was built by computing the mean a_peak in concentric shells around the central star (see Appendix G). The mean value of a_peak was calculated for the cells where the grain–grain differential velocity of grains with size a_peak, noted δv_apeak, is lower than the fragmentation velocity threshold. Both profiles show a clear increase in a_peak, up to ~10 μm, with decreasing envelope radius. We used the fit function of the form f(x) ≡ a + bx^−c. The power-law coefficient, c, goes from 0.9 for model ℳ₀ to 1.1 for model ℳ₁. This power law is perfectly consistent with the fact that the grain size is $\propto \sqrt{\rho }$ since the density scales roughly as $\propto \frac{1}{r^2}$. We note that the slope is steeper in ℳ₁ because turbulence (1) slows the collapse and (2) generates a dense filamentary structure far from the center. We observe that both models produce micron-sized grains, while model ℳ₁ produces slightly larger grains (up to ~8 μm) in the envelope. We interpret this as being due to the turbulence slowing the collapse and generating dense regions away from the collapse center; this increases collision rates, which in turn favors dust growth.

Figures 4b and 4d show 2D histograms of the grain size a_peak according to the distance to the central star and δv_apeak, the grain–grain differential velocity of grains with size a_peak, from turbulence models (Ormel & Cuzzi 2007, see our Appendix A.1). In this model, the level of turbulence is reduced in the disk region (when n_H > 10⁹) in order to be consistent with the model of the α-viscosity disk (Birnstiel et al. 2010). Therefore, δv_apeak is reduced by four orders of magnitude. The disk region is defined for low values of radii and δv_apeak (cross marker). The point markers, located at higher values of δv_apeak, define the envelope region. The accuracy of grain size predictions, when only considering dust growth, depends on the fragmentation limit used to evaluate if grains can survive. Our modeling does not yet include fragmentation processes, but as grains grow in size these processes are expected to become limiting factors to the dust growth. To mitigate the impact of this lacking physics, we used prescriptions from the literature to evaluate whether the grains formed in our models are expected to fragment. The hypotheses made to identify the parameter space where grains from our models are expected to fragment (blue and orange hashed regions in Figs. 4b and 4d) are as follows: To compute the fragmentation of the grains in our models, we first assumed grains are composed of 0.1 μm monomers. Second, we estimated the fragmentation velocity of the grains by comparing the kinetic energy of the grain and the energy required to break a bond between two monomers as v_frag ~ 22 m s⁻¹ for silicates and v_frag ~ 300 m s⁻¹ for ice-coated grains (Ormel et al. 2009; Lebreuilly et al. 2023b). Observations of Class 0 disks and hot-corino regions in low-mass protostars suggest that a significant fraction of the mass at r < 100 au is subject to accretion and viscous heating, with temperatures often higher than 50–100 K (Zamponi et al. 2021). Therefore, while dust grains are expected to be ice-coated in the envelope at radii r > 100 au, we assumed they are probably bare in the region r < 100 au. The fragmentation velocities for the two types of grains (ice-coated and bare grains) are respectively reported as orange and blue hashed areas in Figs. 4b and 4d. All grains with differential dust velocities (δv_apeak; enhanced by gas turbulence from Ormel’s model; see our Appendix A.1) higher than these fragmentation velocities have a high probability of fragmenting and are thus ignored in the following discussion. Note that these differential velocities were computed by taking into account the collision between grains representing the bulk of the dust mass (e.g., with sizes close to a_peak).

We stress that the fragmentation of dust grains in astrophysical environments is not well constrained, as the velocity threshold over which the grain-grain collision will result in grain fragmentation(s) depends on the microphysics properties of a grain, such as its size, shape, and composition. These grain properties are poorly constrained, and there is no consensus as to which grains survive which collisions. For instance, an analysis of grains collected by the Rosetta mission has shown that 50% of the aggregates with sizes of a few microns may have fragmented due to collecting velocities of ≳3–7 m s⁻¹ (Kim et al. 2023). This velocity threshold is much lower than, for example, the theoretical value of 22 m s⁻¹ widely used for bare silicate grains with the same porosity (Ormel et al. 2009). To better constrain the growth of grains up to sizes where their fragmentation is likely in protostellar envelopes, where typical gas velocities are a fraction of a kilometer per second, one needs to make a hypothesis regarding the structure of these grown grains, especially their porosity and fractal dimensions. Such a level of modeling is beyond the current state of the art: future COALA implementations will consistently treat both dust coagulation and fragmentation using a more realistic hypothesis.

To conclude, despite the fact that fragmentation processes are not yet included in the COALA dust evolution model implemented in RAMSES, most of the predicted grains with sizes of a few microns (1–10) are expected to survive in the envelope’s physical conditions. The grains ranging from 4.9 to 7.1 μm, which do not enter the fragmentation regime in the envelope, represent 0.5% of the total dust mass in model ℳ₁, and 0.4% in model ℳ₀. The precise evolution at disk scales requires a dedicated study that includes dust growth and fragmentation and discusses the hypothesis made regarding the disk model that results in the low values of δv_apeak seen in Fig. 4. Interestingly, if grains grow easily to ~10 μm in the envelopes during the first stages of disk formation, before entering the disks and evolving further, this could explain some properties from the dust in comet 67P (Bentley et al. 2016). Indeed, dust particles collected by the Rosetta probe show that grains larger than 8 μm are more compact and less elongated than the smaller particles – as would be expected from compaction processes in a disk environment. However, micron-sized grains consist of porous agglomerates that are suggested to be hierarchical on a scale probably covering at least hundreds of nanometers to several micrometers, such as would be expected from the growth processes at work in envelope environments.

Figure 5 shows the evolution of the fraction of the total dust mass represented by each grain size bin in the envelope in model ℳ₁. Grains can grow easily and very quickly to about 10 μm after 100 kyr in dense envelopes experiencing protostellar collapse. The mass fraction evolution of the grains with a size between ~3 and 5 μm reaches a steady state after 100 kyr, with a value of ~1% of the total dust mass.

The evolution of the dust size distribution affects the nonideal MHD resistivities and, thus, the magnetic braking and other magnetic effects (Zhao et al. 2016; Guillet et al. 2020; Silsbee et al. 2020; Kawasaki et al. 2022; Lebreuilly et al. 2023a). Resistivity evolution is governed by the evolution of the population of small grains. Because our model does not yet include fragmentation, the fraction of very small grains, which drives the ambilopar diffusion resistivity, is not accurately described in these current models. In a future study, we will include fragmentation to accurately account for the evolution of MHD resistivities.

Note that in this study, we did not consider the disruption of grains due to radiative torque from strong radiation fields (Hoang 2019; Reissl et al. 2024). This phenomenon may play an important role in the depletion of large grains in the envelope. This fragmentation process will be added in COALA to quantify which grains are rotationally disrupted to the extent that they affect the dust size distribution in the protostellar envelope, and the results will be compared with those of a previous study (Le Gouellec et al. 2023).

Our work shows that the dust size distribution rapidly evolves from the initial MRN distribution in protostellar envelopes due to the coagulation process, as shown in several previous studies (Bate 2022; Lebreuilly et al. 2023b; Vorobyov et al. 2025). Therefore, in disks whose material comes from the parent envelope, the initial dust size distribution most likely does not follow an MRN distribution and includes large, micron-sized dust grains. Such a different dust size distribution may significantly affect the disk evolution, as it will be key in several physical processes, such as the disk thermal equilibrium and the chemistry at the grain surfaces, that potentially modify the disk’s temperature, scale height, and structure. Moreover, a change in the pristine dust size distribution may also result in a different dust size evolution and spatial distribution in the young disks, as well as the coupling of the disk material with the local magnetic field.

Fig. 4 Radial profiles of the dust grain size at the peak of the size distribution, apeak (see the main text for further details) and the grain-grain differential velocity of grains with size apeak from the turbulence models (see Appendix A.1), noted δvapeak, for the two snapshots shown in Fig. 1: models ℳ0 and ℳ1 at t0.4M⊙. Left panels: Radial profiles built by computing the mean apeak in concentric shells from the apeak map (for further details, see Fig. G.1 and the text). The mean dust grain size clearly increases with decreasing envelope radius, up to ~10 μm. The power law in model ℳ1 is higher than for model ℳ0. This shows that the grains grow faster in model ℳ1 since the gas turbulence tends to slow the collapse, leaving more time for grains to grow. Right panels: Dust grain populations (color-coded by apeak) at all envelope radii, plotted against the δvapeak they experience at this envelope location (in individual cells in the model). A vertical gray line represents the 100 au radius, approximating the ice line expected for such solar-type protostars. On each side of this line, orange and blue shaded areas represent the locations where the dust grains are expected to fragment (see text for further details).

Fig. 5 Time evolution of the fraction of the total dust mass in the envelope represented by each grain size bin. Initially, the grains from 0.19 to 0.27 μm represent ~20% of the total dust mass. After 100 kyr, these grains account for ~13% of the total dust mass. A part of the mass is represented by larger grains, such as grains with sizes of ~3–5 μm, which represent a steady-state value of ~1% of the total dust mass (white lines).

5 Conclusions

In this work we performed, for the first time, 3D nonideal MHD simulations of protostellar collapse that self-consistently account for dust dynamics and growth; we modeled the dust growth using the Smoluchowski equation and coupling the codes COALA and RAMSES. The new code COALA allowed us to include dust growth at a low numerical cost. This new numerical tool allowed us to run simulations beyond the current state of the art. Our goal was to predict the sizes of the grains expected to survive up until the fragmentation regime in a collapsing protostellar envelope during the disk formation stage. We ran two models, with and without an initial turbulent Mach, up to an evolutionary stage representative of a young embedded protostar. Our main findings are as follows:

• In the envelope environment, dust grains grow rapidly to about 1–10 μm in both models (ℳ₀ and ℳ₁) and survive after 100 kyr from an initial MRN distribution;

• In the disk region, the dust size distribution does not follow an MRN distribution and includes large grains of up to several tens of microns in size;

• Dust growth results in a higher dust enrichment compared to dust dynamics alone. More models with various kernels are required to see how influential this can be on the dust content of the disk;

• Dust evolution is highly anisotropic due to the nature of the collapse. There are preferred locations for the production of large grains, such as high-density structures from a turbulent infall (model ℳ₁) or outflow cavity walls (model ℳ₀).

Dust fragmentation processes from grain-grain collisions and radiative torque counterbalance dust growth and thus play key roles in shaping the dust size distribution. Therefore, it is essential to also self-consistently account for the dust fragmentation in 3D MHD simulations in order to precisely predict the local dust size distribution in the envelope and in the disk. An accurate prediction of the dust size distribution is key when studying several physical and chemical processes in the protostellar collapse. New 3D simulations of protostellar collapse with dust growth and fragmentation will be performed with the new tool {COALA + RAMSES} in a future study.

Acknowledgements

We thank the referee for providing very useful comments that helped us to improve our manuscript. This work was made possible thanks to the support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (Grant agreement No. 101098309 – PEBBLES). The simulations were carried out on the Alfven supercomputing cluster of the Commissariat à l’Énergie Atomique et aux énergies alternatives (CEA). Post-processing and data visualization was done using the open source Osyris package.

References

Appendix A Grain-grain differential velocities

In this study, the sources of the grain-grain differential velocities are the turbulence, the Brownian motion and the hydrodynamical drift.

Appendix A.1 Differential velocity from gas turbulence

The gas turbulence is a source of relative velocity between grains through the drag force. The strength of this turbulent motion depends on the grain size and the Reynolds number Re. Ormel & Cuzzi (2007) formulated the expression of the relative velocity between two grains i (the larger grain) and j (the smaller grain) in three regimes which depends on the Stokes number St_i ≡ t_s,i/t_dyn where t_dyn is the dynamical timescale of the system. As shown in Marchand et al. (2021), we can consider only the intermediate regime in the context of the collapse. In that case, the differential velocity writes $$\mathrm{\Delta }{v}_{\text{turb},i,j}\equiv \sqrt{\alpha {\beta }_{i,j}}{c}_{\mathrm{s}}\sqrt{{\mathrm{S}\mathrm{t}}_i},$$(A.1)

The term β_i,j is written as $${\beta }_{i,j}=3.2-(1+{\xi }_{i,j})+\frac{2}{1+{\xi }_{i,j}}(\frac{1}{2.6}+\frac{{\xi }_{i,j}^3}{1.6+{\xi }_{i,j}}),$$(A.2)

where ξ_i,j ≡ St_i/St_j. In this study, the α coefficient, describing the level of turbulence, has a value of 1.5 in the envelope and 10⁻⁴ in the disk region, defined with a criteria on the number density of hydrogen n_H > 10⁹ cm⁻³. As in Tsukamoto et al. (2021), we compute the dynamical time as $$t_{\mathrm{d}\mathrm{y}\mathrm{n}}=\frac{c_{\mathrm{s}}}{\left|{\mathit{a}}_{\mathrm{g}}\right|},$$(A.3)

a_g being the gravitational acceleration.

Appendix A.2 Brownian motion

The differential velocity between two grains of masses m_i and m_j is formulated as $$\mathrm{\Delta }{v}_{\mathrm{B}\mathrm{r},i,j}\equiv \sqrt{\frac{8k_{\mathrm{B}}T}{\pi }}\sqrt{\frac{1}{m_i}+\frac{1}{m_j}},$$(A.4)

where k_b is the Boltzmann constant and T the gas temperature. The differential velocity increases as the masses of the grains decrease. Therefore, the dust growth by Brownian motion is efficient for small grains.

Appendix A.3 Hydrodynamical drift

The dynamics of the grains depends on their size due to the drag force from the gas. Therefore, two grains of different sizes can collide. The relative velocity between grains sizes i and j is written as $$\mathrm{\Delta }{v}_{\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t},i,j}\equiv |{\mathit{v}}_i-{\mathit{v}}_j|.$$(A.5)

This velocity is named hydrodynamical drift velocity.

Fig. B.1 Benchmarks of the code COALA to solve the Smoluchowski coagulation equation (Eq. (1)) with the constant and the additive kernels. In the two plots, the y-axis are the dust grain mass distribution and the x-axis the grain masses in dimensionless unit. The mass range is chosen in order to have 19 orders of magnitude in mass similar to the mass range obtained from smin = 5 nm to smax = 1 nm. The black lines are the initial mass distribution. The blue lines are the exact solutions at a final time tfinal. The green circle are the numerical solutions for which the mass range is sampled in 40 mass bins. The numerical solutions follow accurately the exact solutions.

Appendix B Benchmarking COALA

The code COALA is benchmarked with solutions to the Smoluchowski equation, for the constant and additive kernel, in order to check the good accuracy in the setup with the coupling to the code RAMSES. The mass range spans 19 orders of magnitude, which is similar to the one obtained from the size range s_min = 5 nm to s_max = 1 nm. The mass range is discretized in 40 bins. The results are shown in Fig. B.1 where the numerical solutions are in green and the exact ones are in blue. In both tests, the analytical solution is well reproduced up to a final time where grains grown by 3 orders of magnitude in size. We note a slight numerical diffusion at the high mass tail for the additive kernel. This numerical problem is well known from most of the solvers for the Smoluchowski equation. The new code COALA is well designed to drastically reduce this numerical diffusion by taking advantage of high-order approximation (Lombart & Laibe 2021). In this study, only the piecewise constant approximation is used to keep low numerical costs for the 3D simulations.

Appendix C Selected region in the simulation for the analysis

The RAMSES simulations is around 32000 au side. As indicated in 3, we select 8000 au side cube in order to exclude the impact of the grid boundaries. Figure C.1 shows the full map and the selected region in the edge-on view.

Fig. C.1 Selected cube of 8000 au sides (grey square) in the simulations of around 32000 au side to perform the analysis in Fig. 4. Upper row: Total column gas density maps. Lower row: Grain size at the peak of the dust size distribution.

Appendix D Line-of-sight averaged quantities

The two average quantities shown in the maps and discussed throughout the paper are calculated as $$a_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}\equiv a_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k},i,j}=\frac{\sum _k{a}_{\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k},i,j,k}{\rho }_{i,j,k}}{\sum _k{\rho }_{i,j,k}},\frac{ϵ}{{ϵ}_0}\equiv \frac{{ϵ}_{i,j}}{{ϵ}_0}=\frac{\sum _k\frac{{ϵ}_{i,j,k}}{{ϵ}_0}{\rho }_{i,j,k}}{\sum _k{\rho }_{i,j,k}},$$(D.1)

where i, j and k are the indices for the spatial coordinates of the cells, for which k is the index of the direction along the line of sight. The quantities a_peak,i,j,k and ϵ_i,j,k are the grain size at the peak of the size distribution and the dust-to-gas ratio in a given cell.

Appendix E Scaling of the peak of the dust size distribution

We demonstrate here a simple relation between the typical grain size s and the gas density ρ. We assume a simple spherical collapse such as $$\frac{d\rho }{dt}=\frac{\rho }{t_{\mathrm{f}\mathrm{f}}},$$(E.1)

and a monodisperse formulation of the grain size evolution (Stepinski & Valageas 1996; Tsukamoto et al. 2021) $$\frac{ds}{dt}=\frac{s}{3t_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}},$$(E.2)

where the coagulation timescale writes $$t_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}=\frac{m_d}{4\pi s^2{\rho }_d\mathrm{\Delta }v},$$(E.3)

and the free-fall timescale writes $$t_{\mathrm{f}\mathrm{f}}=\sqrt{\frac{3\pi }{32G\rho }}.$$(E.4)

Combining the two equations yields $$\frac{ds}{d\rho }=\frac{s}{3\rho }\frac{t_{\mathrm{f}\mathrm{f}}}{t_{\mathrm{c}\mathrm{o}\mathrm{a}\mathrm{g}}}.$$(E.5)

Assuming coagulation by turbulence, we can write $$\mathrm{\Delta }v=A\sqrt{\frac{t_{\mathrm{s}}}{t_{\mathrm{f}\mathrm{f}}}}.$$(E.6)

Finally, injecting Eq. 4 with S_grain ← s and Eqs. E.3, E.4, and E.6 into Eq. E.5, we obtain $$\frac{ds}{\sqrt{s}}={\left(\frac{3\pi }{32G}\right)}^{1/4}A\frac{{ϵ}_{\mathrm{d}\mathrm{g}}}{\sqrt{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}}}\frac{1}{\sqrt{w_{\mathrm{t}\mathrm{h}}}}\frac{d\rho }{{\rho }^{3/4}},$$(E.7)

where ϵ_dg ≡ ρ_dust/ρ is the dust-to-gas ratio. By integration we finally get $$\sqrt{s(t)}-\sqrt{s_0}=2{\left(\frac{3\pi }{32G}\right)}^{1/4}A\frac{{ϵ}_{\mathrm{d}\mathrm{g}}}{\sqrt{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{i}\mathrm{n}}}}\frac{1}{\sqrt{w_{\mathrm{t}\mathrm{h}}}}[\rho (t{)}^{1/4}-{\rho }_0^{1/4}],$$(E.8)

or $$s(t)={(2{\left(\frac{3\pi }{32G}\right)}^{1/4}A\frac{{ϵ}_{\mathrm{d}\mathrm{g}}}{\sqrt{{\rho }_{\text{grain }}}}\frac{1}{\sqrt{w_{\text{th }}}}[\rho (t{)}^{1/4}-{\rho }_0^{1/4}]+\sqrt{s_0})}^2.$$(E.9)

Assuming a small initial size and density yields $$s(t)=4{\left(\frac{3\pi }{32G}\right)}^{1/2}{A}^2\frac{{ϵ}_{\mathrm{d}\mathrm{g}}^2}{{\rho }_{\text{grain }}}\frac{1}{w_{\mathrm{t}\mathrm{h}}}\sqrt{\rho }.$$(E.10)

We note that $A\simeq \sqrt{2\alpha }{c}_s$ and $w_{\text{th}}=\sqrt{\frac{8}{\pi \gamma }}{c}_s$. We finally get $$s(t)\simeq \sqrt{8\pi \gamma }{\left(\frac{3\pi }{32G}\right)}^{1/2}\alpha c_s\frac{{ϵ}_{\mathrm{d}\mathrm{g}}^2}{{\rho }_{\text{grain }}}\sqrt{\rho }$$(E.11)

Or $$\begin{array}{rl}s(t)\simeq & 0.004\mathrm{c}\mathrm{m}\left(\frac{\alpha }{1}\right){\left(\frac{T}{10K}\right)}^{1/2}{\left(\frac{{ϵ}_{\mathrm{d}\mathrm{g}}}{0.01}\right)}^2{\left(\frac{{\rho }_{\text{grain }}}{2.3\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}}\right)}^{-1}\\ & \left(\frac{\rho }{{10}^{-13}\mathrm{g}{\mathrm{c}\mathrm{m}}^{-3}}\right)\end{array}$$(E.12)

We point out that, as shown in Lebreuilly et al. (2020), the Stokes number St, which is the ratio between the stopping time and the free-fall time scales as $\propto \frac{s}{\sqrt{\rho }}$. This means that this quantity is roughly constant in the isothermal regions of the collapse when dust growth is included.

Appendix F Model ℳ₁ with and without dust growth

In Fig. F.1 we compare simulations of the model ℳ₁ with (bottom row) and without (top row) dust growth. The total column density Σ and a_peakare defined in Sect. 3. The dust growth highlights the anisotropy of the dust evolution, where grains grow faster in some specific regions, i.e. structures due to turbulent infall in the model ℳ₁. We note that the initial MRN distribution evolved locally in the simulation due to the dust growth.

Fig. F.1 Comparison without dust growth (upper row) and with dust growth (lower row) for the model ℳ1 with M⋆ = 0.4M⊙. The red cross gives the position of the star. The left and middle columns show the total column gas density Σ and the average grain size at the peak of the dust size distribution apeak similarly to Fig. 1. For the model without growth, the apeak is the maximum size of the MRN distribution (~250 nm). The right column shows the evolution of the dust size distribution in a cell in the simulation.

Appendix G Radial profile of the dust size distribution

The radial profile of the dust size distribution in the 3D simulation is built by computing the mean of the dust size distribution in selected concentric shells around the star. The width of the shells is defined as 10% of the radius in order to capture enough cells to compute the mean. Figure G.1 shows the shells with radius [100,200,300,400,500] (top row) with the corresponding average dust size distribution (bottom row). From these average dust size distributions, we compute the mean grain size a_peak at the peak of the distribution.

Fig. G.1 Radial dependence of the dust grain size distribution, from the model ℳ1 with M⋆ = 0.4 M⊙. Concentric shells of varying radii from 100 to 4000 au (as shown in the upper panels) are used to compute the mean dust grain size distributions shown in the lower panels. The dust grains size distribution is represented by the mass fraction distribution (see text for details). Here only the five first concentric shells with radius [100,200,300,400,500] au are shown in upper panel with the corresponding average dust grain size distribution weighted by the gas density (blue line) in the lower panel. The dust grain size distributions in each cell defining the shell are shown (light grey lines) in the lower panel.

All Figures

Fig. 1 Snapshots taken from the RAMSES MHD models (ℳ0 and ℳ1) of a collapsing protostellar envelope (initial parameters can be found in the text) with dust grain size evolution computed by COALA. The maps show a 2000 au zoom in on the model when the central star has accreted 0.4 M⊙from the initial 2.5 M⊙core. The red cross represents the position of the star. The left column shows the total column gas density, Σ. The middle column shows the average grain size at the peak of the size distribution, apeak, weighted by the gas density along the line of sight (see the main text for further details). Similarly, the right column shows the average dust enrichment (ratio between the evolved dust-to-gas ratio, ϵ, and the initial dust-to-gas ratio, ϵ0) weighted by the gas density along the line of sight. In the top row, the square patterns near the sink are not numerical artifacts but result from the weighted average of the quantities along the line of sight.

In the text

Fig. 2 Snapshots and histograms of the dust enrichment at two given times of model ℳ1 with (top row) and without (bottom row) dust growth. The histograms are obtained with the dust enrichment calculated in each cell in a box with sides of 2000 au around the star. The two times are when M⋆ = 0.06 M⊙ (orange) and M⋆ = 0.4 M⊙ (blue). Dust growth and dynamics result in larger local variations in the dust enrichment, by a factor of greater than 2, in the envelope compared to the dust dynamics alone. The black line shows the histograms for the ℳ1 model with dust growth and 60 dust bins, highlighting the numerical convergence reached with 40-bin simulations. The red cross represents the position of the star. The variance value of ~0.01 of the histogram explains the low variation observed on the map for model ℳ1 when M⋆ = 0.4 M⊙.

In the text

	Fig. 3 Comparison of the computational times of the three simulations for model ℳ0: gas only (black), gas with 40 dust sizes (blue), and gas with 40 dust sizes and dust growth (orange). Using COALA in RAMSES only increases the global execution time by a factor of 1.7.
In the text

Fig. 4 Radial profiles of the dust grain size at the peak of the size distribution, apeak (see the main text for further details) and the grain-grain differential velocity of grains with size apeak from the turbulence models (see Appendix A.1), noted δvapeak, for the two snapshots shown in Fig. 1: models ℳ0 and ℳ1 at t0.4M⊙. Left panels: Radial profiles built by computing the mean apeak in concentric shells from the apeak map (for further details, see Fig. G.1 and the text). The mean dust grain size clearly increases with decreasing envelope radius, up to ~10 μm. The power law in model ℳ1 is higher than for model ℳ0. This shows that the grains grow faster in model ℳ1 since the gas turbulence tends to slow the collapse, leaving more time for grains to grow. Right panels: Dust grain populations (color-coded by apeak) at all envelope radii, plotted against the δvapeak they experience at this envelope location (in individual cells in the model). A vertical gray line represents the 100 au radius, approximating the ice line expected for such solar-type protostars. On each side of this line, orange and blue shaded areas represent the locations where the dust grains are expected to fragment (see text for further details).

In the text

Fig. 5 Time evolution of the fraction of the total dust mass in the envelope represented by each grain size bin. Initially, the grains from 0.19 to 0.27 μm represent ~20% of the total dust mass. After 100 kyr, these grains account for ~13% of the total dust mass. A part of the mass is represented by larger grains, such as grains with sizes of ~3–5 μm, which represent a steady-state value of ~1% of the total dust mass (white lines).

In the text

Fig. B.1 Benchmarks of the code COALA to solve the Smoluchowski coagulation equation (Eq. (1)) with the constant and the additive kernels. In the two plots, the y-axis are the dust grain mass distribution and the x-axis the grain masses in dimensionless unit. The mass range is chosen in order to have 19 orders of magnitude in mass similar to the mass range obtained from smin = 5 nm to smax = 1 nm. The black lines are the initial mass distribution. The blue lines are the exact solutions at a final time tfinal. The green circle are the numerical solutions for which the mass range is sampled in 40 mass bins. The numerical solutions follow accurately the exact solutions.

In the text

	Fig. C.1 Selected cube of 8000 au sides (grey square) in the simulations of around 32000 au side to perform the analysis in Fig. 4. Upper row: Total column gas density maps. Lower row: Grain size at the peak of the dust size distribution.
In the text

Fig. F.1 Comparison without dust growth (upper row) and with dust growth (lower row) for the model ℳ1 with M⋆ = 0.4M⊙. The red cross gives the position of the star. The left and middle columns show the total column gas density Σ and the average grain size at the peak of the dust size distribution apeak similarly to Fig. 1. For the model without growth, the apeak is the maximum size of the MRN distribution (~250 nm). The right column shows the evolution of the dust size distribution in a cell in the simulation.

In the text

Fig. G.1 Radial dependence of the dust grain size distribution, from the model ℳ1 with M⋆ = 0.4 M⊙. Concentric shells of varying radii from 100 to 4000 au (as shown in the upper panels) are used to compute the mean dust grain size distributions shown in the lower panels. The dust grains size distribution is represented by the mass fraction distribution (see text for details). Here only the five first concentric shells with radius [100,200,300,400,500] au are shown in upper panel with the corresponding average dust grain size distribution weighted by the gas density (blue line) in the lower panel. The dust grain size distributions in each cell defining the shell are shown (light grey lines) in the lower panel.

In the text