© The Authors 2025

1 Introduction

The formation of a protostar and its surrounding protoplanetary disk is a complex process that involves numerous physical and chemical phenomena. One of the key components of these processes is the evolution of dust grains. These grains are formed in the cold and dense regions of the interstellar medium (ISM). They represent approximately 1% of the mass budget in the medium, and their size distribution is well modeled by a power law called the Mathis et al. (1977) distribution (MRN). They play a crucial role in the formation of the protostar and the disk, including their blackbody emission and opacity, which allow the cooling of the gas (Semenov et al. 2003; Omukai et al. 2005; Tsuribe & Omukai 2006). They also serve as the main formation sites for H₂ molecules, which heats up the medium (Gould & Salpeter 1963). Furthermore, their coupling with the magnetic field regulates its evolution via non-ideal magnetohydrodynamic (MHD) effects (Marchand et al. 2016; Tsukamoto & Okuzumi 2022). The size of the grains, and notably the presence of small grains, is also a key parameter as it can significantly impact the formation and shape of the disk (Zhao et al. 2016; Tsukamoto et al. 2023).

It is fundamental to understand the evolution of dust grains during the protostellar formation to model accurately all the physical processes that are related to them, including the dust opacity and the non-ideal MHD resistivities. The main effects that can impact the evolution of dust grains during the gravitational collapse of a molecular cloud are the interactions with the gas phase which remove (grain vaporization or sublimation) or add materials (grain growth or condensation), and notably through chemical reactions that occur on the surface of the grains (called chemisputtering), which leads ultimately to their total destruction when the temperature reaches a few thousand Kelvin. Lenzuni et al. (1995) performed an analysis of these effects in the case of a protostellar core contraction; they modeled the vaporization of three types of dust grains: carbonaceous, silicate, and aluminum oxide grains. The results of this study were then used to compute, among other things, the evolution of non-ideal MHD resistivities with temperature, and to build a resistivity table that has been used in recent numerical simulations (Marchand et al. 2016; Vaytet et al. 2018). Similarly, the modeling of the evolution of the dust opacity in numerical simulations is currently also made using an opacity table (Vaytet et al. 2013).

However, in this same study from Lenzuni et al. (1995), the authors highlight the fact that the dust-gas interactions could have timescales similar to the protostar formation timescale (the free-fall time). Thus, the dust vaporization, the dust opacity, and the non-ideal MHD resistivities cannot be solely functions of the temperature and density of the medium; instead, a dynamic approach has to be taken. One limitation of the study from Lenzuni et al. (1995) for the protostar formation is that they only considered a single evolution path for the temperature and the density of the medium. Since this evolution corresponds to a protostellar core contraction and not a gravitational collapse, their temperature and density increase much more slowly than can be encountered in a collapse simulation, as shown in Fig. 1. In young protostars, the continuous accretion and ejection and/or the rapid disk expansion can also cause a rapid increase and decrease in temperature and density, which could lead to different sublimation limits (see Tsukamoto et al. 2021; Morbidelli et al. 2024).

Thus, we propose to revisit the dust grain vaporization in the context of a protostar formation, with an approach similar to Lenzuni et al. (1995). In Sect. 2 we present the model and the method used to compute the evolution of the dust grains. In Sect. 3 we discuss the materials considered. In Sect. 4 we present the results of the computation of the dust grain evolution in a protostellar formation context. Finally, in Sect. 5 we discuss the implications of this study.

Fig. 1 Comparison of the evolutions of temperature (in green) and density (in red) used by (Lenzuni et al. 1995, solid line), with the evolution encounter in the central cell of a collapse simulation of a 1 solar mass cloud with the hydrodynamic code RAMSES (dashed line). The details of the simulation are given in Sect. 4.

2 Dust evolution modeling

To model the dust, we used the same set of assumptions as Lenzuni et al. (1995). A dust grain is defined as a collection of N monomers (which can be a molecule or an atom) of volume V_m and mass m_m , defining the volume of the grain as V = V_mN. We assumed that the grains are pure (i.e., each one is composed of only one kind of monomer), and we also limit our study to spherical grains, defining their radius a through the relation $V=V_{\text{m}}N=\frac{4}{3}\pi a^3$, and their surface area corresponding to 𝒜 = 4πa². Finally, we also supposed that the dust grains are always in thermal equilibrium with the gas phase, meaning that both share the same temperature T at each time.

Given this set of assumptions, we could then model the evolution of the number of monomers in one dust grain following the method of Gail & Sedlmayr (2013). We start from the equation $$\frac{\text{d}N}{ \text{d}t}=F_{\text{gr}}-F_{\text{vap}},$$(1)

where F_gr and F_vap are respectively the rates of addition and subtraction of monomers on the grain through reaction with the gas phase. It is more convenient to express them in terms of reaction rate per unit surface area J, so that F = 𝒜J. We can then use the relation between N and a to find the equation of evolution of the grain radius: $$\frac{\text{d}a}{ \text{d}t}=V_{\text{m}}\left(J_{\text{gr}}-J_{\text{vap}}\right).$$(2)

2.1 Grain growth

There are several ways to add monomers to a grain. It notably depends on whether the monomer is stable in the gas phase or not. The presence of other chemical species in the gas phase can modify the reaction path and then the growth rate.

We start with the simplest reaction that can add a monomer to the grain, which is $${\text{M}}_{N_{(\text{s})}}+{\text{M}}_{1{}_{(\text{g})}}\to {\text{M}}_{N+1{}_{(\text{s})},}$$(3)

where M_N(s) represents the grain with N monomers, and M_1(g) represents the monomers in the gas phase. The reaction rate of this process can be expressed as $$F_{\text{gr}}=\alpha f_{\text{col}},$$(4)

where f_col is the collision frequency between a monomHer and the grain, and α is the sticking coefficient, corresponding to the probability that a collision results in an absorption. The collision frequency in the case of a spherical grain can be expressed as $$f_{\text{col }}=\frac{1}{4}\mathcal{A}{v}_{\text{rel }}{n}_{\text{m}}=4\pi a^2\sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{m}}}}{n}_{\text{m}},$$(5)

where υ_rel is the mean relative thermal velocity between the dust grain and the monomers, which can be expressed as $\sqrt{\frac{8k_{\text{B}}T}{\pi m_{\text{m}}}}$¹, and n_m is the number density of monomers in the gas phase. We then find the reaction rate per unit surface area as $$J_{\text{gr}}=\alpha \sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{m}}}}{n}_{\text{m}}.$$(6)

We can then extend this case to a more general reaction that can add a monomer to the grain $${\text{M}}_{N{}_{(\text{s})}}+{\displaystyle \sum _{k=1}^{N_{\text{r}}}{v}_k}{\text{A}}_k\to {\text{M}}_{N+1{}_{(\text{s})}}+{\displaystyle \sum _{k=1}^{N_{\text{p}}}{\mu }_k}{\text{B}}_k,$$(7)

where the ν and µ parameters are the stoichiometric coefficients of the reaction and the A_k and B_k are the reactants and products, respectively. If this reaction has a limiting step in its kinetic path, we can estimate the growth rate by looking at the key reactant associated with this step. This reactant can be found through experimental results, but it is possible to get an idea by looking at the number densities n_k of each reactant and comparing the different n_k/ν_r,k ratios. If one is significantly smaller than the others, the species associated with it is a good candidate to be the key reactant. In this case, the growth rate per unit surface can be written similarly to our simple case and can be expressed as $$J_{\text{gr}}=\alpha \sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{key }}}}\frac{n_{\text{key }}}{v_{\text{key }}},$$(8)

where n_key , m_key , and νkey are respectively the number density of the key reactant, its mass, and its associated stoichiometric coefficient.

2.2 Grain vaporization

The removal of monomers from the grain can be caused either by ejection of a monomer from the surface of the grain by thermal agitation (thermal vaporization or free vaporization) or by a reaction with the gas phase (chemisputtering). In both cases these processes can be described as the reverse of reactions (3) and (7). A general reaction describing the evolution of the number of monomers in the grain can be written as $${\text{M}}_{N(\text{s})}+{\displaystyle \sum _{k=1}^{N_{\text{r}}}{v}_k}{\text{A}}_k\rightleftharpoons {\text{M}}_{N+1(\text{s})}+{\displaystyle \sum _{k=1}^{N_{\text{p}}}{\mu }_k}{\text{B}}_k.$$(9)

The presence of species Bk on the right-hand side indicates whether this is free vaporization (no B_k) or chemisputtering (at least one B_k). There are two ways to estimate the rate of vaporization. We can use the detailed balance principle, but in the case of chemisputtering, we can also use a similar approach to the growth model with a key species that is responsible for the chemisputtering. The second approach is more accurate, but requires experimental data to determine the reaction probability.

2.2.1 Detailed balance principle

This principle states that at equilibrium, the growth rate and the vaporization rate are equal for each reaction occurring at the surface of the grain. This principle can be expressed in terms of rate per unit surface area as $$J_{\text{vap }}=J_{\text{gr},\text{\hspace{0.17em}\hspace{0.17em}eq}}=\alpha \sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{key }}}}\frac{{\stackrel{\circ }{n}}_{\text{key }}}{v_{\text{key }}},$$(10)

where ${\stackrel{\circ }{n}}_{\text{key}}$ is the number density of the key reactant (defined in Sect. 2.1) when the thermodynamical equilibrium is reached. The total rate of evolution of the grain can then be written as $$\begin{array}{l}{J}_{\text{tot }}=J_{\text{gr }}-J_{\text{vap }}=\alpha \sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{key }}}}\frac{n_{\text{key }}}{v_{\text{key }}}\left(1-\frac{{\stackrel{\circ }{n}}_{\text{key}}}{n_{\text{key}}}\right)\\ =\alpha \frac{p_{\text{key }}}{v_{\text{key }}\sqrt{2\pi m_{\text{key }}{k}_{\text{B}}T}}\left(1-\frac{{\stackrel{\circ }{p}}_{\text{key }}}{p_{\text{key }}}\right),\end{array}$$(11)

where p_key = n_keyk_B T is the partial pressure of the key reactant.

All we need now is a method for estimating the number density and pressure at equilibrium of the key reactant. We assume that we have a state out of equilibrium, with the partial pressure of the key species denoted p_key and the partial pressure of all the reactants and products denoted $p_{{\text{A}}_k}$ (with k ≠ key) and $p_{{\text{B}}_k}$ . We denote by ${\stackrel{\circ }{p}}_{\text{key }}$ the partial pressure of the key reactant that is needed to reach equilibrium, given the partial pressure of all the other reactants and products. To compute ${\stackrel{\circ }{p}}_{\text{key }}$, we use thermochemical considerations: at equilibrium between the gas phase and the grain, reaction (9) obeys the law of mass action which states that we have the following relation between the different partial pressures of the different species: $${\displaystyle \prod _{k=1}^{N_{\text{p}}}{p}_{{\text{B}}_k}^{{\mu }_k}}={\stackrel{\circ }{p}}_{\text{key}}{}^{v_{\text{key}}}{e}^{-{\text{Δ}}_{\text{r}}G/RT}{\displaystyle \prod _{k=1,k\ne \text{key}}^{N_{\text{r}}}{p}_{{\text{A}}_k}^{v_k}}.$$(12)

Here Δ_rG is the Gibbs energy (or free enthalpy) per mole of reaction (9), and can be expressed as $${\text{Δ}}_{\text{r}}G={\text{Δ}}_{\text{f}}G\left({\text{M}}_{1(\text{s})}\right)+{\displaystyle \sum _{k=1}^{N_{\text{p}}}{\mu }_k}{\text{Δ}}_{\text{f}}G\left({\text{B}}_k\right)-{\displaystyle \sum _{k=1}^{N_{\text{r}}}{v}_k}{\text{Δ}}_{\text{f}}G\left({\text{A}}_k\right),$$(13)

where Δ_f G is the Gibbs energy of formation per mole and M_1(s) represents a monomer on the dust grain.

From the relation (12), we just need to isolate ${\stackrel{\circ }{p}}_{\text{key }}$ and divide by the actual partial pressure of the key species p_key to obtain $$\frac{{\stackrel{\circ }{p}}_{\text{key }}}{p_{\text{key }}}={\left(\frac{{\displaystyle {\prod }_{k=1}^{N_{\text{p}}}{p}_{{\text{B}}_k}^{{\mu }_k}}}{e^{-{\text{Δ}}_{\text{r}}G/RT}{\displaystyle {\prod }_{k=1}^{N_{\text{r}}}{p}_{{\text{A}}_k}^{v_k}}}\right)}^{1/v_{\text{key}}},$$(14)

which gives us the ratio which appears in Eq. (11). This ratio can be interpreted as the inverse of the pseudoactivity of the monomer on the grain (Gail & Sedlmayr 2013), which is often also called the supersaturation ratio. Thus, just by knowing all the partial pressures of the different species involved in the reaction, Eq. (14) allows us to compute the total rate of evolution of the grain (11).

2.2.2 Kinetic approach

The detailed balance principle approach of the vaporization works well when the current state is not far from equilibrium. In situations where vaporization completely dominates the growth in reaction (9), the thermodynamic approach is limited as it does not take into account kinetic considerations for the vaporization of the grain (the α parameter is defined for the condensation kinetic process) that can reduce the effective vaporization rate, notably in the case of chemisputtering. In this situation we model the vaporization rate per unit surface area as Lenzuni et al. (1995) did, with the expression $$J_{\text{vap }}=Y_{\text{key }}\sqrt{\frac{k_{\text{B}}T}{2\pi m_{\text{key ,vap }}}}\frac{n_{\text{key,vap }}}{{\mu }_{\text{key,vap }}},$$(15)

where n_key,vap and µ_key,vap are associated with the key species involved in the vaporization process, and the Y_key is the yield coefficient (or the reaction probability, similar to the sticking coefficient) of the key species, which needs to be determined experimentally. This expression is very similar to expression (10), but instead focuses on the presence of species B_key, which is responsible for the chemisputtering of the grain. This approach is more accurate than the detailed balance principle, but we need to have experimental data to determine the yield coefficient.

For this work, we used Eq. (15) to compute the total rate of vaporization of the grain if experimental data were available, otherwise we used the detailed balance principle (10).

2.3 Numerical method

Now that we have solved Eq. (2), which allows us to compute the evolution of the grain radius, we need to compute this evolution for a given evolution of temperature T(t) and density ρ(t) for the gas phase surrounding the grain. As stated before, we make the assumption that the grain is always at thermal equilibrium with the gas phase. The numerical method is summarized in Fig. 2.

First of all, the computation of the radius evolution rate J_tot needs as input the composition of the gas phase, or at least the quantity of each species involved in the reactions with the grain. To this end, we use the code FASTCHEM2 (Stock et al. 2022), which allows us to compute the chemical composition of the gas phase at equilibrium, given a temperature, pressure, and the relative abundance of each chemical element. This means that we make the assumption that the gas phase is always at chemical equilibrium, which can lead to some inaccuracies (see Sect. 5 for more details).

For this calculation, the only requirement is to set the abundance of each element in the gas phase at the beginning of the computation, allowing us to determine the composition of the gas phase with FASTCHEM2. For this we need to set a global elemental abundance for the gas–dust mixture, and set the initial dust quantity. Then, we just have to remove the elements composing the dust from the global elemental abundance to get the composition of the gas phase. Then, to integrate Eq. (2), we need an integration scheme. We use the Euler method, but better methods such as Runge-Kutta could be used.

The choice of the integration time step Δt is crucial. To do this we need to consider the different timescales that exist in our problem: the temperature evolution timescale Δt_T, and the radius evolution timescale Δt_a. The density evolution timescale τ_ρ has the same order of magnitude as Δt_T for the different gas evolutions we consider in this work, and is therefore not considered. To evaluate these timescales, we need to define a maximum temperature step ΔT and a maximum radius step Δa. It is then possible to define the timescales as $$\begin{array}{l}\text{Δ}{t}_T=\frac{\text{Δ}T}{\frac{ \text{d}T}{ \text{d}t}}\\ \text{Δ}{t}_a=\frac{\text{Δ}a}{\frac{ \text{d}a}{ \text{d}t}}=\frac{\text{Δ}a}{V_{\text{m}}{J}_{\text{tot}}}.\end{array}$$(16)

Several cases can occur:

1. Δt_T ⪅ Δt_a: the grain evolves at the same speed as or more slowly than the temperature. The time step is then set as the smallest value between the two timescales, and we perform a single integration step.

2. Δt_T ≫ Δt_a: the evolution of the grain is very fast compared to the evolution of the temperature; thus, we can consider that it reaches its equilibrium at a constant temperature. As we want to limit computation time, we first increase the time by Δt_T to change the temperature and density, such that the grain is set out of equilibrium. We then set Δt = Δt_a, or a fraction of it, and perform multiple steps until the grain radius reaches its equilibrium with the gas phase a_eq(T). To know when the equilibrium is reached, we compare the derivative of the radius with a parameter ϵ, and we stop the integration when the derivative is smaller than ϵ .

If the grain radius changes during the time step, it also means that it has modified the composition of the gas phase by removing or adding the materials composing the monomer. The number of monomers in the dust grains that change during the time step Δt can be computed from the radius variation and is equal to $$\text{Δ}N=\frac{1}{V_{\text{m}}}\frac{4}{3}\pi \left(a{(t+\text{Δ}t)}^3-a{(t)}^3\right).$$(17)

From this we can remove or add from the gas phase the chemical elements forming the monomers, and then recompute the new chemical composition of the gas with FASTCHEM2.

It should be noted that the right-hand side of Eq. (2) does not depend on the radius of the grain. It is then possible to evolve multiple grains with different radii at the same time, such that we can mimic the evolution of a grain size distribution.

Fig. 2 Schematization of the numerical method used to compute the evolution of a grain radius. When the radius a evolves, we use Eq. (2) to compute the radius time derivative. The computation of the gas composition is done with FASTCHEM2 (Stock et al. 2022).

3 Dust materials

To apply this model to a protostellar formation problem, we need to set the composition of the dust that is present during the process. We choose a solar composition for the chemical elemental abundances, such that the three main species of grain that are present are carbonaceous grains, silicates (specifically olivine here), and aluminum oxides. For the initial dust quantity, we follow Lenzuni et al. (1995), assuming that most of the Si, Mg, and Al are locked in the grains, such that their elemental abundances set the dust quantities of silicates and aluminum oxides, while for the carbonaceous grains, it is 70% of the carbon is locked in the dust. Some of their properties are summarized in Table 1. For all the reactions that appear in this section, the values for the Gibbs energies of reaction are taken from Chase (1998) and are given in Appendix A. We describe below how each of these species interacts with the gas phase.

3.1 Carbon

The monomers composing the carbonaceous grains are mainly carbon atoms. The first reaction with the gas phase that can then occur at the surface of this kind of grain is the adsorption of carbon atoms alone or in small groups of i atoms (called i- mers, with i typically equal to 2 or 3) in the gas phase, and their subsequent vaporization. This corresponds to the reaction $${\text{C}}_{N(\text{s})}+\frac{1}{i}{\text{C}}_{i(\text{g})}\rightleftharpoons {\text{C}}_{N+1(\text{s})}.$$(18)

This reaction can be treated with the thermodynamic approach described in Sect. 2.2.2, but it is also possible to use a kinetic approach, as described by Grassi et al. (2017) or Stahler et al. (1981), who considered the probability that a monomer can escape the grain using the Debye frequency, which takes into account the structure of the grain. Our calculations show that both approaches are equivalent, so we keep the thermodynamic approach for this study. Finally, the sticking coefficient for this reaction is about 0.3 as shown in Table 3 of Lenzuni et al. (1995).

Carbonaceous grains can also be subject to chemisputtering by three chemical species in the gas phase: atomic hydrogen, molecular hydrogen, and water molecules. The chemisputtering by atomic hydrogen takes the global form $${\text{C}}_{N+m(\text{s})}+n{\text{H}}_{(\text{g})}\to {\text{C}}_{N(\text{s})}+{\text{C}}_m{\text{H}}_{n(\text{g})}.$$(19)

This kind of chemisputtering has been extensively studied, including studies by Barlow & Silk (1977) and Draine (1979). The global process is summarized in the work of Lenzuni et al. (1995). Since these reactions are dominated by the vaporization process, we use Eq. (15) to compute the global radius evolution through hydrogen sputtering, and we use Eq. (24) of Lenzuni et al. (1995) for the yield coefficient: $$Y_{\text{H}}=Y_{\max }\exp \left(-{\left[\frac{\Vert T-T_{\max }\Vert }{\sigma }\right]}^{1.75}\right).$$(20)

Here Y_max = 3.15 × 10⁻², T_max = 625 K, and σ = 1.35 K.

The H₂O molecules can also be a source of chemisputtering, as shown by Stahler et al. (1981). We again use the same yield coefficient as in Lenzuni et al. (1995) with their Eq. (27), which is $$Y_{{\text{H}}_2\text{O}}=Y_0\exp \left(-\frac{T_0}{T}\right),$$(21)

with Y₀ = 6.85 × 10⁵ and T₀ = 29 000 K.

Finally, the chemisputtering by molecular hydrogen H₂ is also possible. This chemisputtering is not taken into account by Lenzuni et al. (1995) because of the lack of data. In our case we use the survey of Krakowski & Olander (1970), which compiles various measurements of the yield coefficient $Y_{{\text{H}}_2}$ as a function of the temperature. The results are shown in Fig. 3.

3.2 Silicates

The monomer composing the silicate grains is the group Mg2SiO₄ (forsterite). Several studies have reported experimental data for the vaporization of forsterite in an H₂ atmosphere (Hashimoto 1990; Nagahara & Ozawa 1994, 1996; Tsuchiyama et al. 1998; Kuroda & Hashimoto 2002). Gail & Sedlmayr (1999) highlighted the study of Nagahara & Ozawa (1996) and discussed it as a chemical sputtering reaction whereby H₂ reacts with forsterite and enhances the vaporization compared with a vacuum vaporization where there is no reactive gas. This dichotomy in reaction mechanisms led these authors to describe the vaporization of olivine with two mechanisms that are summarized by two reactions: the free vaporization $${\text{Mg}}_2{\text{SiO}}_{4(\text{s})}\rightleftharpoons 2\text{Mg}+\text{SiO}+\frac{1}{2}{\text{O}}_2,$$(22)

and the chemisputtering by H₂ $${\text{Mg}}_2{\text{SiO}}_{4(\text{s})}+3{\text{H}}_2\rightleftharpoons 2\text{Mg}+\text{SiO}+3{\text{H}}_2\text{O}.$$(23)

For both reactions we use Eq. (11) to compute the vaporization and growth rate. The sticking coefficient α is in practice close to 0.1 for SiO, as shown in Shornikov et al. (1998); Fedkin et al. (2006), where values between 0.05 and 0.15 have been found.

Fig. 3 Yield coefficient of the chemisputtering by molecular hydrogen as a function of the temperature from the survey of Krakowski & Olander (1970).

3.3 Aluminum oxides

The monomer composing the aluminum oxide grains is the group Al₂O₃ . In the presence of H₂, there are three main reactions that can occur, shown to be predominant by thermodynamic calculations with the thermochemical software and the database FactSage (Bale et al. 2016): $${\text{Al}}_2{\text{O}}_{3(\text{s})}+3{\text{H}}_2\rightleftharpoons 2\text{Al}+3{\text{H}}_2\text{O,}$$(24) $${\text{Al}}_2{\text{O}}_{3(\text{s})}+2{\text{H}}_2\rightleftharpoons {\text{Al}}_2\text{O}+2{\text{H}}_2\text{O,}$$(25) $${\text{Al}}_2{\text{O}}_{3(\text{s})}+{\text{H}}_2\rightleftharpoons 2\text{AlO}+2{\text{H}}_2\text{O}.$$(26)

They are all chemisputtering reactions by H₂ . Similarly to the case of silicate, we use Eq. (11). We can also consider the free vaporization of aluminum oxide, which is described by the reaction $${\text{Al}}_2{\text{O}}_{3(\text{s})}\rightleftharpoons 2\text{Al}+\frac{3}{2}{\text{O}}_2.$$(27)

The values for the sticking coefficients range between 0.2 and 0.3, according to Burns (1966).

Fig. 4 Color map showing the ratio between the chemisputtering vaporization and the free vaporization flux as a function of temperature and density of the gas, with a saturation at 10−5 and 105. Red corresponds to the area where the chemisputtering vaporization flux dominates, and blue corresponds to the area where the free vaporization flux dominates. The yellow dotted line represents the barotropic law (see Eq. (30)), and the green dotted line represents a typical disk model from Andrews et al. (2009, see our Sect. 5.2 and Eqs. (41) and (42)). The cyan solid line represents the limit where the vaporization flux is equal to the growth flux. The vaporization dominates for higher temperatures (above the line) and the growth dominates for lower temperatures (below the line). The carbon grains do not have this line as vaporization always dominates for the range of temperatures and density of this plot. The cyan dashed line represents the sublimation limit prescription of the silicates used in Isella & Natta (2005). Its expression is given in Eq. (31). Finally, the orange solid line represents the limit where the vaporization timescale is equal to the free-fall timescale. For higher temperatures (above the line), the vaporization timescale is smaller than the free-fall timescale, and for lower temperatures (below the line), the vaporization timescale is larger than the free-fall timescale. The two other orange lines represent the limit where the vaporization timescale is equal to the free-fall timescale, but multiplied by a factor of 100 (dot-dashed) and 0.01 (dashed). The jump in the orange lines on the middle panel is due to a change in the key reactant in reaction (23), being Mg for lower densities and SiO for higher densities. The vaporization timescales are computed with a grain radius of 5 × 10−7 cm in Eq. (28).

3.4 Timescales

It is possible to compute the characteristic timescale of variation of the grain radius for each of these reactions as a function of temperature. Via Eq. (2), we can compute the characteristic time of destruction of a grain of radius a₀ by a given reaction as $${\tau }_{\text{grain }}=\frac{a_0}{V_{\text{m}}{J}_{\text{tot}}}.$$(28)

This allows us to have a first idea of the importance of each reaction in the evolution of the grain. In the context of gravitational collapse, the free-fall time is the characteristic time of evolution of the temperature and density of the gas, and is given by $${\tau }_{\text{ff}}=\sqrt{\frac{3\pi }{32G\rho }}.$$(29)

In the context of gravitational collapse, the barotropic law from Machida et al. (2006) gives a good idea of the evolution of the temperature as a function of the density. It is defined as $$T=T_0\sqrt{1+{\left(\frac{\rho }{{\rho }_1}\right)}^{2{𝑔}_1}}{\left(1+\frac{\rho }{{\rho }_2}\right)}^{{𝑔}_2}{\left(1+\frac{\rho }{{\rho }_3}\right)}^{{𝑔}_3},$$(30)

with the critical densities ρ₁ = 3.9 × 10⁻¹³ g cm⁻³, ρ₂ = 3.9 × 10⁻⁸ g cm⁻³, and ρ₃ = 3.9 × 10⁻³ g cm⁻³, and the coefficients ɡ₁ = 0.4, ɡ₂ = −0.3, and ɡ₃ = 0.56667.

From this we can obtain a good idea of the sublimation limit of each grain material by examining Fig. 4. We note that, in the context of gravitational collapse (dotted black line), chemisputtering dominates vaporization for all three types of grain.

Then, there are two conditions that determine if a dust grain is sublimated for a given temperature and density: the vaporization flux must be greater than the growth flux, such that vaporization dominates growth (the vaporization-growth limit), and the vaporization timescale must be smaller than the dynamical timescale of the gas phase (the dynamical limit), so that the grain has time to evaporate at this temperature. In the context of increasing temperature, the sublimation limit for each material in Fig. 4 corresponds to the higher curve between the orange and black solid lines, representing the dynamical limit and the vaporization-growth limit, respectively.

These considerations demonstrate that in the context of gravitational collapse (black dotted line), the vaporization of silicate and aluminum oxide grains is determined by the vaporizationgrowth limit. We see that our model reproduces quite well the silicate sublimation limit used in Isella & Natta (2005, the blue dashed line), which is based on observational data from Pollack et al. (1994) and gives the sublimation temperature of silicate grains to be equal to $$T_{\text{sub }}=2000{\left(\frac{{\rho }_{\text{gas}}}{1 \text{g} {\text{cm}}^{-3}}\right)}^{1.95\times {10}^{-2}} \text{K}.$$(31)

For carbon grains, we see that the dynamical limit is the relevant factor for the sublimation limit. We also observe that if we change the value of the dynamical timescale by a factor of 0.01 and 100, it significantly changes the sublimation limit of the carbon grains, highlighting the crucial role of the dynamical evolution of the gas in determining the evolution of carbon grains, while it only slightly affects the dynamical limit of silicate and aluminum oxide grains.

To understand why the carbon material is very dependent on the dynamical evolution of the gas, we can examine the effects of each reaction that can destroy a carbon grain. In Fig. 5 we compare the free-fall time and different characteristic times of destruction of a carbon grain with a radius of 5 × 10⁻⁷ cm (the smallest grains in our size distribution) as a function of temperature, assuming a barotropic law. We observe that the chemisputtering by H atoms has a range of temperatures where the characteristic time has approximately the same order of magnitude as the free-fall time, between 900 and 1100 K. This has the consequence of reducing the steepness of the total characteristic vaporization time, resulting in significant changes in the crossing temperature with the different characteristic timescales. This suggests that the dynamical approach is necessary to accurately model the evolution of carbon grains.

Fig. 5 Evolution of the characteristic time of destruction of a carbon grain with a radius of 5 × 10−7 cm as a function of temperature for each reaction. The free-fall time is also shown, assuming a barotropic law for temperature-density relation (Machida et al. 2006).

4 Dust evolution computation

Now that we have defined the content of our dust grains, we can try to apply our model to different temperature and density evolutions T(t) and ρ(t), which we call trajectories from now on.

4.1 Initial setup

Before computing the evolution of the dust grains, we need to set the initial conditions for the quantity of dust grains and the size distribution. For the quantity of each material, we use the values of the dust-to-gas ratios in Table 1 to compute their total density from the initial gas density of the trajectory.

For the size distribution, we use the power-law distribution from Mathis et al. (1977), also known as the MRN distribution. According to the MRN distribution, between a_min = 5 × 10⁻⁷ cm and a_max = 2.5 × 10⁻⁵ cm, the grain size distribution is nonzero and can be expressed as $$\text{d}n(a)=C{a}^{\lambda }\text{d}a,$$(32)

where dn is the number density of grains that have a radius contained between a and a + da, λ is the power-law index and is close to −3.5 according to Mathis et al. (1977), and C is a normalization constant. This constant can be computed from the total dust density ρ_dust as $${\rho }_{\text{dust }}={\displaystyle {\int }_{a_{\min }}^{a_{\max }}{\rho }_{\text{grain }}}\frac{4}{3}\pi a^{3}C{a}^{\lambda }\text{d}a.$$(33)

The computation of this integral gives the expression of C, $$C=\frac{{\rho }_{\text{dust }}}{{\rho }_{\text{grain }}}\frac{3(\lambda +4)}{4\pi a_{\min }^{{}_{{}^{\lambda +4}}}\left({\xi }^{\lambda +4}-1\right)},$$(34)

where ξ = a_max/a_min. To compute the evolution of the grain size distribution, and also the non-ideal MHD resistivities and the opacities, it is useful to work directly with a discretized version of the MRN distribution. To do this, we choose a number of bins N_bin and we define the radius limit between the bin i and i + 1 (where i is between 1 and N_bin) as $a_{\text{lim},i}=a_{\min }{\xi }^{i/N_{\text{bin}}}$, so that they are logarithmically spaced. The number of grains in bin i is then given by $$n_i=\frac{C}{\lambda +1}\left(a_{\lim ,i}^{\lambda +1}-a_{\lim ,i-1}^{\lambda +1}\right).$$(35)

We now choose that all grains in each bin have the same radius a_i. The choice of this radius a_i depends on whether we want to conserve the total mass of the grains or their total surface area. However, if we take a sufficiently high number of bins, the difference between the two choices becomes negligible, and we can simply choose the geometric mean of the bin, so that $a_i=\sqrt{a_{\lim ,i}{a}_{\lim ,i-1}}$. In our case, we set N_bin = 100, such that the error in the total dust mass and the total surface area of the MRN distribution is less than 0.04%. This process results in the grain repartition represented by the blue dot in Fig. 6. We can now apply the method described in Sect. 4.1 for the evolution of each radius in the distribution.

Fig. 6 Dust size distribution for aluminum oxide grains at the beginning of the computation (in blue) and at the end of the computation (in red) for the trajectory from Bhandare et al. (2024) shown in Fig. 8.

4.2 Dust evolution

We can now compute the evolution of the MRN dust distribution for a given trajectory, which consists of a temperature T(t) and density ρ(t) evolution. We show the results of the computation for three different trajectories. The first one is the trajectory from Lenzuni et al. (1995), which is shown in Fig. 1.

The second one is a trajectory that follows the evolution of the central computational cell of a collapse simulation performed with the AMR hydrodynamic code RAMSES (Teyssier 2002). The numerical and physical setup is similar to the one used in Ahmad et al. (2023): we start from a one solar mass cloud with a ratio of thermal to gravitational energy of α = 0.25, without any rotation or magnetic field. The evolution in temperature and density of the central cell is shown in Fig. 1.

The third is the trajectory of a one-micron grain extracted from a simulation of the collapse of a solar mass molecular cloud performed by Bhandare et al. (2024), shown in Fig. 8. The grain trajectories selected here are the ones that experience an ejection from the center of the system (see Fig. 9), causing a fast decrease in temperature and density. The evolution of the dust-to-gas ratio for each grain species and for each trajectory is shown in Fig. 7.

First of all, we see that the results for the Lenzuni et al. (1995) trajectory (left panel in Fig. 7) are in good agreement with what was obtained in their article. The only difference is the temperature of full destruction of aluminum oxide, which is around 1610 K in our computation, whereas it is around 1720 K for Lenzuni et al. (1995). This difference probably arises from the difference in the modeling of the vaporization process. In our case we use a set of three reactions (24)–(26), while in their case they computed the chemical equilibrium of ${\text{Al}}_2{\text{O}}_{3_{\left(\text{s}\right)}}$ between the chemical species Al, AlOH, Al₂O, AlH, H₂, and H₂O.

The main difference between Lenzuni’s trajectory and the trajectory from the collapse simulation is the speed of evolution of temperature and density (see Fig. 1), where the latter evolves much faster. This difference directly impacts the results of the computation, particularly for the carbonaceous grains. Since the vaporization reactions for this type of grain have a timescale close to the evolution of temperature, as shown in Fig. 5, a change in the speed of temperature variation modifies the temperature of full destruction of the carbonaceous grains. For Lenzuni’s trajectory it is around 1030 K, while it becomes around 1120 K for the collapse simulation. There is no real impact on the other two species of grains since the timescales for their reactions are much shorter than the evolution of temperature.

In the case of the trajectory from the simulation of Bhandare et al. (2024) shown in Fig. 8, there are two regimes. It starts with the temperature and density increasing at a rate similar to the collapse simulation, but at a certain point, the growth rate suddenly increases until it reaches a maximum, and then it decreases at the same rate (caused by the ejection of the grains from the center of the collapse). The extremely fast increase in temperature implies a shift of the temperature of full destruction of the carbonaceous grains to around 1310 K (higher than the vaporization limit of silicates, which is around 1290 K). The subsequent decrease in temperature before destroying all the aluminum oxide grains leads to a sudden growth of the latter grains, recovering all their mass when the temperature is sufficiently low. However, as the carbonaceous grains and silicates have already been completely evaporated, their population remains at zero, resulting in a drastic change in the dust quantity for temperatures lower than 1300 K compared to before the temperature spike. This result is due to the absence of nucleation in our model, and the impossibility of recondensation of carbonaceous and silicate materials on the remaining aluminum oxide grains, which is discussed in Sect. 5.1. The recondensation of the aluminum oxide grains is also visible in the size distribution shown in Fig. 6. As the recondensation can only happen on the grains that are still present in the gas, we obtain at the end of the trajectories a size distribution that is devoid of its smallest grains.

For all the other trajectories shown in Fig. 8, the results are similar and mainly depend on the maximum temperature reached. If this temperature is below 1300 K, all types of grains recover their initial mass. If it is between 1300 and 1600 K, the results are similar to what was just discussed. However, if the temperature exceeds 1600 K, the grains are completely destroyed and there is no recondensation when the temperature decreases.

Fig. 7 Dust-to-gas ratio as a function of temperature for three different trajectories. For panel a we use the trajectory from Lenzuni et al. (1995), shown in Fig. 1. For panel b we use a trajectory from a collapse simulation, also shown in Fig. 1. For panel c we use a trajectory from the collapse simulation of Bhandare et al. (2024), shown in Fig. 8. The arrows in the right panel indicate the time evolution of the system.

Fig. 8 Evolutions of the temperature in panel a and the density in panel b of the gas surrounding one-micron grains in the collapse simulation of Bhandare et al. (2024). The 55 trajectories selected from the simulation are those where the grains experience a maximum temperature between 1000 and 2000 K. Each trajectory starts at a temperature of 600 K. The red solid line represents the trajectory chosen in this article to compute the dust evolution. The grain position evolution is shown in Fig. 9.

Fig. 9 Evolution of the grain positions associated with the temperature and density evolution from panels a and b of Fig. 8.

Fig. 10 Evolutions of the three non-ideal MHD resistivities as a function of temperature for the trajectory from Lenzuni et al. (1995) shown in Fig. 1. On the left is shown the ambipolar diffusion resistivity ηad, in the middle is the Hall effect resistivity ηH, and on the right is the Ohmic resistivity ηΩ. In each plot the resistivity computed through the dust evolution computation is shown as a solid line, and the resistivity computed through the Marchand et al. (2016) vaporization modeling is shown as a dashed line.

4.3 Non-ideal MHD resistivities

The three non-ideal MHD resistivities that appear in the nonideal MHD equations (Marchand et al. 2016) can be computed from a given discretized grain size distribution. For this computation, we use a routine using the analytical derivation of the grain ionization from Marchand et al. (2021). This routine needs as input a discretized grain size distribution, the temperature T , the numerical density of particles in the gas phase n_gas, the ionization rate x_i , and the magnetic field magnitude B. The temperature is given directly by the trajectory we used to compute the dust vaporization, and the numerical density of particles in the gas phase is given by the gas density and the mean molecular weight of the gas μ_gas, so that n_gas = ρ/ (μ_gasm_H), with μ_gas = 2.34 for a solar composition, and m_H is the mass of hydrogen. For the ionization rate and the magnetic field magnitude we use the same prescription as Marchand et al. (2016), we set a constant ionization rate x_i = 10⁻¹⁷ s⁻¹, and we follow Li et al. (2011) such that the magnetic field scales as $$B=1.43\times {10}^{-7}\sqrt{\frac{n_{\text{gas}}}{1 {\text{cm}}^{-3}}}\text{\hspace{0.17em}G}.$$(36)

The goal of this section is to compare the results of our vaporization computation with those of the vaporization modeling by Marchand et al. (2016). In this article it is assumed that the grain sizes always follow an MRN distribution, and it is assumed that the grain quantity is a function only of the temperature: each type of grain evaporates linearly within a given temperature range. These ranges were chosen based on the results of Lenzuni et al. (1995): carbon evaporates between 750 and 1100 K, silicates between 1200 and 1300 K, and aluminum oxides between 1600 and 1700 K. As our computation gives slightly different results for Lenzuni’s trajectory, we chose to modify the vaporization temperature ranges so that the strict comparison of the two resistivity models is meaningful. We set the vaporization temperature ranges to be between 800 and 950 K for the carbon grains, between 1250 and 1300 K for the silicate grains, and between 1530 and 1620 K for the aluminum oxide grains. With this choice, we observe in Fig. 10 that the resistivities computed from our vaporization computation for Lenzuni’s trajectory are almost identical to the resistivities computed with the simple vaporization model described above.

However, if we maintain the same vaporization model and consider the Bhandare et al. (2024) simulation trajectory, we obersve some differences in our computations (see Fig. 11). First of all, the late vaporization of the carbonaceous grains generates higher resistivities for temperatures between 950 and 1300 K. When the temperature decreases below 1300 K, the resistivities are much lower (by a factor of 100) than the simple vaporization model. In the Marchand et al. (2016) vaporization model, all the grains recover their mass when the temperature decreases sufficiently, which does not occur in our computation since the materials cannot recondense if no more grain is present. Additionally, the resistivities between 1300 and 1550 K just after the spike in temperatures are also different by a factor of 10 between the model and the computation, while the aluminum oxide grains recovered their initial mass (see Fig. 7). The reason for this is shown in Fig. 6: the spike in temperature removes the smallest grains. Even if the total mass of the aluminum oxide grains is recovered, this causes the grain number to be reduced, thereby decreasing the resistivities.

Fig. 11 Same as Fig. 10, but for the trajectory from Bhandare et al. (2024) shown in Fig. 8. The arrows indicate the time evolution of the system.

Fig. 12 Planck and Rosseland mean opacities as a function of temperature for three different trajectories. For panel a we use the trajectory from Lenzuni et al. (1995), shown in Fig. 1. For panel b we use a trajectory from a collapse simulation, also shown in Fig. 1. For panel c we use a trajectory from the collapse simulation of Bhandare et al. (2024) shown in Fig. 8.

4.4 Opacities

It is also possible to compute the opacities of a given dust distribution. For this purpose, we use the library DSHARP (Birnstiel et al. 2018). By taking a grain radius a and the complex refractive index $\underset{\_}{\text{n}}=n+ik$ of the material forming the grain, it can compute for a given frequency f the opacity per unit of mass κ(f, a). For a mixture of N_species grain species with N_radii different radii, the total opacity per unit of mass is given by $${\kappa }_{\text{tot}}(f)={\displaystyle \sum _i^{N_{\text{species}}}}\frac{{{\displaystyle \sum }}_j^{N_{\text{radii}}}{n}_{i,j}{a}_j^3{\kappa }_i\left(f,a_j\right)}{{{\displaystyle \sum }}_j^{N_{\text{radi}}}{n}_{i,j}{a}_j^3},$$(37)

with n_i,j the number density of grains of species i and radius j.

There are several ways to compute a mean opacity over the frequency. The two most common are the Planck and Rosseland mean opacities, which appear in the equations for the radiative transfer and are then important for the modeling of the protostar formation. The Planck opacity is defined as $${\kappa }_{\text{P}}=\frac{{{\displaystyle \int }}^{\text{}}\kappa (f)B(f,T)\text{d\hspace{0.17em}}f}{{{\displaystyle \int }}^{\text{}}B(f,T)\text{d\hspace{0.17em}}f},$$(38)

where B (f, T) is the Planck function. The Rosseland opacity is defined as $$\frac{1}{{\kappa }_{\text{R}}}=\frac{{{\displaystyle \int }}^{\text{}}\frac{1}{\kappa (f)}\frac{\partial B(f,T)}{\partial T}\text{d\hspace{0.17em}}f}{{{\displaystyle \int }}^{\text{}}\frac{\partial B(f,T)}{\partial T}\text{d\hspace{0.17em}}f}.$$(39)

We can track the evolution of the Planck and Rosseland mean opacities for each trajectory. For the values of the complex refractive index, we use the values from Zubko et al. (1996) for carbonaceous grains, Draine (2003) for silicate grains, and Eriksson et al. (1981) and Kischkat et al. (2012) for aluminum oxide grains. The results are shown in Fig. 12. We observe that each dust vaporization phase reduces the total opacity of the medium. There are not many differences between the Lenzuni et al. (1995) trajectory and the collapse simulation trajectory, except for a shift in the destruction of carbonaceous grains to higher temperatures, which also leads to a slight shift in the opacities. For the Bhandare et al. (2024) simulation trajectory, the simultaneous vaporization of carbonaceous and silicate grains results in a significant drop in the opacities around 1300 K. The recondensation of aluminum oxide grains as the temperature starts to decrease leads to a slight increase in the opacities, as the radius of surviving grains increases.

5 Discussion

5.1 Limitations

The main limitations that we considered for this study are as follows:

1. We limited the study to pure grain, implying that recondensation of a given material is not possible on other kinds of grains. It could also change the vaporization process since a given material could be protected from chemisputtering or free vaporization by another material covering it, one that evaporates at higher temperatures.

2. We did not take into account nucleation, meaning that we did not create any new grains. This could potentially allow fully destroyed species to start to recondense.

3. We supposed that the gas phase is in chemical equilibrium, which is probably not the case. It is notably argued in Lenzuni et al. (1995) that the production of CO molecules from CH_n molecules, product of the chemisputtering of the carbonaceous grains, is too slow to happen during the grain vaporization process. This kinetic limitation allows the carbon atoms to recondense more easily on the grains, shifting upward the temperature of full destruction of carbonaceous grains. This could be corrected in the future by computing the chemical evolution of the gas phase, with the code ULCHEM for instance (Holdship et al. 2017).

4. The model used to compute non-ideal MHD resistivities (Marchand et al. 2021) is different from the model developed in Marchand et al. (2016), which is used to compute the resistivity table used in numerical simulations (Vaytet et al. 2018). The main reasons for these differences are that graingrain collisions and thermionic emission (significant from 800 K) are not taken into account in the Marchand et al. (2021) model, and the resistivity table from Marchand et al. (2016) is limited to grains with only one charge for its computation. These three effects imply that, in disk regions, the table gives resistivity values that are lower than in the model used in this paper. However, these differences do not change the global results and conclusions of Sect. 4.3.

5. We do not know how the quantity of dust materials that has been ejected and experienced partial sublimation compares with the remaining dust in the disk. As the ejected dust goes back into the outer regions of the disk (Tsukamoto et al. 2021), this could potentially change the MHD resistivities if the dust reflux from the ejection is not negligible compared to the remaining dust in the envelope.

5.2 Disk application

Another context where we could apply our grain evolutionary model could be the protoplanetary disks. The main differences with the protostellar collapse is the characteristic timescale of evolution of the gas temperature and density: it is not the freefall time but the Keplerian time, which corresponds to the time to perform a full rotation around the central star. The Keplerian time is given by the Kepler’s third law and is expressed as $${\tau }_{\text{K}}=\sqrt{\frac{4{\pi }^2{r}^3}{G{M}_{\star }}},$$(40)

with M_★ the mass of the central star and r the distance to the central star. We can perform the same timescale analysis as in Sect. 3.4 for the protostellar disk. We use for that a typical disk model orbiting a one solar mass central star from Andrews et al. (2009) which is in good agreement with observations, and gives the following expression for temperature and middle density of the disk $$T=280{\left(\frac{r}{1\text{\hspace{0.17em}AU}}\right)}^{-1}\text{\hspace{0.17em}\hspace{0.17em}K},$$(41) $$\rho =6\times {10}^{-10}{\left(\frac{r}{1\text{\hspace{0.17em}\hspace{0.17em}AU}}\right)}^{-9/4}\text{g\hspace{0.17em}\hspace{0.17em}}{\text{cm}}^{-3}.$$(42)

The temperature-density implicit relation is shown in Fig. 4 (dotted green line). We see that the curve is quite close to the barotropic law for temperature lower than 1600 K, such that the sublimation limit for silicate and aluminum oxide grains should remain the same as in our collapse application. However, the carbonaceous grains are affected by the dynamical evolution of the gas, and the Keplerian time is much longer than the free-fall time so we could expect a difference in the sublimation limit.

It is possible again to compute the timescale of destruction of a carbon grain of radius a₀ for each kind of reaction. The results for the smallest grain of the MRN distribution are shown in Fig. 13, with a comparison with the Keplerian time. The crossing of the total timescale (red dotted line) with the Keplerian time curve gives us an idea of the radius at which carbon grains should be destroyed. The chemisputtering by H₂ and H₂O give us a threshold radius of around 8 × 10⁻² AU.

The chemisputtering by H₂ has a radius range of around 1 × 10⁻¹ AU, where it is the main reaction destroying the grains, and it is a factor of 10 longer than the Keplerian time. This has the consequence of reducing the steepness of the total vaporization timescale. As a result, there is a larger range of radii around 1 × 10⁻¹ AU where the total vaporization timescale is greater than the Keplerian time, but only by a factor of 10–100. This means that in a few revolutions, the carbon grains can be greatly affected by the H-sputtering when they are in this range of radii from the central stars. This suggests that a dynamical approach to the vaporization of carbon grains in the protoplanetary disk may be necessary as the H-sputtering prevents us from determining a specific radius of destruction.

Fig. 13 Evolution of the characteristic time of destruction of a carbon grain with a radius of 5 × 10−7 cm as a function of the radius for each reaction, for the Andrews et al. (2009) disk model. The Keplerian time is also shown for comparison.

6 Conclusion

In this article we presented a model for computing the evolution of the dust grains in a protostellar collapse. On the one hand, we showed that the dynamical evolution of the gas is crucial for the evolution of the carbonaceous grains, notably because the chemisputtering by H atoms has a range of temperatures where the characteristic time has the same order of magnitude as the free-fall time. On the other hand, the silicate and aluminum oxide grains sublimate at a relatively fixed temperature, but the evolution of their size distribution is crucial in order to compute the non-ideal MHD resistivities and the opacities. We also showed that the vaporization of the grains in complex trajectories with decreasing temperature can have a significant impact on the nonideal MHD resistivities and the opacities of the medium, which are not adequately captured by a simple temperature-dependent vaporization model. Additionally, we discussed the limitations of our model and the potential applications of our model to protoplanetary disks. Despite these limitations, the dynamical approach gives results that are quite different from the current one used in numerical simulations, notably in the case of fast temperature and density variation. This suggests that in the context of protostellar collapse or of dust grain ejection from outflows, the dynamical approach is necessary.

A potential follow-up to this work could be to include the chemical evolution of the gas phase, to take into account the nucleation of new grains, and to apply the model to protoplan- etary disks. In this way, we could implement the vaporization of the grains in hydrodynamical codes to have a more realistic model of dust evolution in the protostellar collapse. The works of Lombart et al. (2024) on grain coagulation and fragmentation could be extended to take into account the chemisputtering of the grains, and be coupled with various hydrodynamical codes, such as RAMSES (Teyssier 2002) and IDEFIX (Lesur et al. 2023).

Acknowledgements

We thank the anonymous referee for their useful comments that have improved the quality of this paper. This work has received funding from the French Agence Nationale de la Recherche (ANR) through the projects PROMETHEE (ANR-22-CE31-0020), and DISKBUILD (ANR-20-CE490006). We thank Asmita Bhandare for providing us with the data of the collapse simulation from their Bhandare et al. (2024) paper. We finally thank Sébastien Charnoz for the enriching discussions we had.

Appendix A Thermochemical data

We list here the values of all the parameter needed to compute the vaporization of the grains. To compute the Gibbs energy of formation ∆_fG of the different chemical species involved in the model, we use the standard enthalpy of formation ∆_f⁰H and the standard entropy S⁰ of the species. The Gibbs energy of formation at a given temperature T is then given by $${\text{Δ}}_{\text{f}}G={\text{Δ}}_{\text{f}}{}^{0}H-T{S}^0.$$(A.1)

The data for the species involved in reactions (22) to (27) are shown in Table A.1. The values for the sticking coefficients α used in the model are also shown in Table A.2 for information, but their precise values have little impact on the model results as long as the order of magnitude chosen is of the good order (here, between 0.1-0.3). Finally, for the other gas species that are not directly involved in the vaporization process, but which still appear in the gas equilibrium computation by FASTCHEM2 (Stock et al. 2022), we use the data table available (file logK.dat) in the FASTCHEM GitHub repository², which is build using mainly the data from Chase (1998).

References

All Tables

All Figures

	Fig. 1 Comparison of the evolutions of temperature (in green) and density (in red) used by (Lenzuni et al. 1995, solid line), with the evolution encounter in the central cell of a collapse simulation of a 1 solar mass cloud with the hydrodynamic code RAMSES (dashed line). The details of the simulation are given in Sect. 4.
In the text

	Fig. 2 Schematization of the numerical method used to compute the evolution of a grain radius. When the radius a evolves, we use Eq. (2) to compute the radius time derivative. The computation of the gas composition is done with FASTCHEM2 (Stock et al. 2022).
In the text

	Fig. 3 Yield coefficient of the chemisputtering by molecular hydrogen as a function of the temperature from the survey of Krakowski & Olander (1970).
In the text

Fig. 4 Color map showing the ratio between the chemisputtering vaporization and the free vaporization flux as a function of temperature and density of the gas, with a saturation at 10−5 and 105. Red corresponds to the area where the chemisputtering vaporization flux dominates, and blue corresponds to the area where the free vaporization flux dominates. The yellow dotted line represents the barotropic law (see Eq. (30)), and the green dotted line represents a typical disk model from Andrews et al. (2009, see our Sect. 5.2 and Eqs. (41) and (42)). The cyan solid line represents the limit where the vaporization flux is equal to the growth flux. The vaporization dominates for higher temperatures (above the line) and the growth dominates for lower temperatures (below the line). The carbon grains do not have this line as vaporization always dominates for the range of temperatures and density of this plot. The cyan dashed line represents the sublimation limit prescription of the silicates used in Isella & Natta (2005). Its expression is given in Eq. (31). Finally, the orange solid line represents the limit where the vaporization timescale is equal to the free-fall timescale. For higher temperatures (above the line), the vaporization timescale is smaller than the free-fall timescale, and for lower temperatures (below the line), the vaporization timescale is larger than the free-fall timescale. The two other orange lines represent the limit where the vaporization timescale is equal to the free-fall timescale, but multiplied by a factor of 100 (dot-dashed) and 0.01 (dashed). The jump in the orange lines on the middle panel is due to a change in the key reactant in reaction (23), being Mg for lower densities and SiO for higher densities. The vaporization timescales are computed with a grain radius of 5 × 10−7 cm in Eq. (28).

In the text

	Fig. 5 Evolution of the characteristic time of destruction of a carbon grain with a radius of 5 × 10−7 cm as a function of temperature for each reaction. The free-fall time is also shown, assuming a barotropic law for temperature-density relation (Machida et al. 2006).
In the text

	Fig. 6 Dust size distribution for aluminum oxide grains at the beginning of the computation (in blue) and at the end of the computation (in red) for the trajectory from Bhandare et al. (2024) shown in Fig. 8.
In the text

Fig. 7 Dust-to-gas ratio as a function of temperature for three different trajectories. For panel a we use the trajectory from Lenzuni et al. (1995), shown in Fig. 1. For panel b we use a trajectory from a collapse simulation, also shown in Fig. 1. For panel c we use a trajectory from the collapse simulation of Bhandare et al. (2024), shown in Fig. 8. The arrows in the right panel indicate the time evolution of the system.

In the text

Fig. 8 Evolutions of the temperature in panel a and the density in panel b of the gas surrounding one-micron grains in the collapse simulation of Bhandare et al. (2024). The 55 trajectories selected from the simulation are those where the grains experience a maximum temperature between 1000 and 2000 K. Each trajectory starts at a temperature of 600 K. The red solid line represents the trajectory chosen in this article to compute the dust evolution. The grain position evolution is shown in Fig. 9.

In the text

	Fig. 9 Evolution of the grain positions associated with the temperature and density evolution from panels a and b of Fig. 8.
In the text

Fig. 10 Evolutions of the three non-ideal MHD resistivities as a function of temperature for the trajectory from Lenzuni et al. (1995) shown in Fig. 1. On the left is shown the ambipolar diffusion resistivity ηad, in the middle is the Hall effect resistivity ηH, and on the right is the Ohmic resistivity ηΩ. In each plot the resistivity computed through the dust evolution computation is shown as a solid line, and the resistivity computed through the Marchand et al. (2016) vaporization modeling is shown as a dashed line.

In the text

	Fig. 11 Same as Fig. 10, but for the trajectory from Bhandare et al. (2024) shown in Fig. 8. The arrows indicate the time evolution of the system.
In the text

	Fig. 12 Planck and Rosseland mean opacities as a function of temperature for three different trajectories. For panel a we use the trajectory from Lenzuni et al. (1995), shown in Fig. 1. For panel b we use a trajectory from a collapse simulation, also shown in Fig. 1. For panel c we use a trajectory from the collapse simulation of Bhandare et al. (2024) shown in Fig. 8.
In the text

	Fig. 13 Evolution of the characteristic time of destruction of a carbon grain with a radius of 5 × 10−7 cm as a function of the radius for each reaction, for the Andrews et al. (2009) disk model. The Keplerian time is also shown for comparison.
In the text