© The Authors 2025

1 Introduction

Planets around low-mass stars are formed from dust grains that consist of refractory materials and volatile ices. The grains initially reside in the interstellar medium, part of which transforms into a dense molecular cloud that contains a prestellar core (Tafalla 2018). After the protostar is born, the remainder of the core forms a protostellar envelope. Dust and gas fall from the envelope onto the protoplanetary disk (PPD). During this process, dust and gas are heated, up to a few hundred K. This can occur via two primary mechanisms. Firstly, heating is provided by a mild accretion shock at the centrifugal barrier near the outer rim of the PPD (Cassen & Moosman 1981; Sakai et al. 2014). Secondly, when the infall trajectory is near the outflow cavity, the grains are temporarily heated by protostellar radiation (Visser et al. 2009). Upon arrival on the disk, the gas and dust become shielded from protostellar and interstellar radiation and cool again. The disk material continues to experience heating events either through shocks in the disk or through flaring of the protostar (Armitage 2011).

The composition of interstellar ices can be reset during star formation before the ices are included into planets or comets (Jørgensen et al. 2020). The partial resetting of ice before it arrives in PPDs has been confirmed by astrochemical simulations (Visser et al. 2011; Drozdovskaya et al. 2014; Aota et al. 2015; Furuya et al. 2017) and in PPDs themselves (Podolak et al. 2011; Molyarova et al. 2018). These studies employed models that considered a single grain size. However, grains are expected to have a variety of sizes, shapes, and composition. Their size and material can cause differences in the dust grain temperatures T_d. A proven approach in numerical simulations is to discern the grains according to their size (Acharyya et al. 2011). For the modelling, the grain-size distribution can be represented by several discrete grain-size bins. These multi-grain models show that different temperatures for the grain-size bins affect the formation rate of the icy mantles and their eventual composition (Pauly & Garrod 2016).

The temperatures experienced by a gas-dust parcel infalling from the envelope towards the PPD can exceed the sublimation temperatures of the main ice species CO, CO₂, and H₂O. In this case, ice species first sublimate from the hottest grain population. Kalvāns & Freimanis (2024, hereafter KF24) found that observations of complex organic molecules (hot corinos) and carbon chains (warm carbon-chain chemistry) in low-mass protostars are best explained by an assortment of temperatures of the dust grains, rather than a single T_d for all grains.

The aim of this study is to explore how the rise and decrease in grain temperatures at typical infall conditions redistribute the ices between hotter and colder grain populations. In particular, we explore whether a heating-cooling cycle can induce an accumulation of ices in the coldest grain-size bin but leaves other grains ice poor. This situation briefly occurred during heating in the KF24 model. An asymmetric redistribution of ices like this can considerably affect the grain size and stickiness, which in turn affects the expected starting conditions for planet formation in PPDs. To fulfil the aim of this study, we applied the multigrain model of KF24 for a parcel of circumstellar matter infalling towards the inner region of a protostellar core. The parcel was heated before it cooled again, presumably during infall towards the inner region (PPD) of a protostellar core.

2 Methods

The study is based on the model of Kalvāns et al. (2024, hereinafter KKV24). The application of this model for protostellar envelope chemistry has illustrated the importance of the grain temperature distribution in circumstellar environments (KF24).

2.1 Chemical model

The recently updated multi-grain kinetic rate-equation astrochemical model ALCHEMIC-VENTA (originally introduced by Kalvāns 2015) now has an advanced treatment for molecular grain surface processes. We only include a brief description of the chemical aspects and refer to KKV24, where the chemical aspects are explained in detail and with the relevant references. Table 1 summarises the features of the code and the default values of the important physical and chemical parameters, as used in our Standard model.

A single major change was introduced in the code: The bulkice multi-layer approach for modelling the icy mantles on grains was removed. This was required for a physically correct description of rapid thermal desorption and re-adsorption of ice species onto grains. With a multi-layer approach, the rate of these processes could be artificially slowed down by the simulated transfer of molecules between the ice layers. Considering bulk ice does not substantially change the overall surface chemistry at elevated temperatures of the protostellar stage (Garrod 2013).

The surface chemistry starts with molecule adsorption onto the grain surface. Ice accumulates into five grain-size bins, and smaller grains have larger populations according to the MRN grain-size distribution, as shown in Table 2. All grains share a similar composite chemical make-up. The grains were assumed to be partially processed, as expected for a star-forming region. More precisely, the small grains were depleted onto the larger grains (Silsbee et al. 2020; Sipilä et al. 2020), whose abundance was increased by a factor of 1.25 to retain the MRN mass budget. Moreover, the surfaces of the refractory grains were assumed to be carbonaceous to account for the adsorption of polycyclic aromatic hydrocarbon molecules and polymerised organic residue that formed by the over-photoprocessing of ices during the previous star-forming activity (Hagen et al. 1979; Allamandola et al. 1988; Quirico et al. 2016). This assumed carbonaceous coating is not able to participate in strong inter-molecular hydrogen bonds with ice molecules, and this affects the chemical properties of bare grains (KKV24).

The average temperature of dust depends on its composition, on the intensity of electromagnetic radiation, and on the grain size (Li & Draine 2001; Draine 2003; Cuppen et al. 2006; Li & Mann 2012). The latter two aspects were considered here. The temperature of the grain populations was determined by the grain radius, following the equation $$T_{d,a}={\left(\frac{a+b}{0.1}\right)}^y\times T_{d,0.1},$$(1)

where a is the radius of the refractory grain, b is the ice-mantle thickness (in µm), and T_d,0.1 is the temperature of 0.1 µm grains (Krügel 2007; Pauly & Garrod 2016). The exponent y is a variable. For simplicity, T_d,0.1 was assumed to be equal to T_gas (Sect. 2.2). T_gas has been shown to be similar to T_d of sub-micron grains by more complex models (Gavino et al. 2021).

The surface processes include binary reactions, photodissociation of the ice species, and six desorption mechanisms. These are sublimation, photodesorption by the interstellar UV radiation field (ISRF), and cosmic-ray (CR) induced photons, chemical desorption of the products of exothermic reactions, indirect chemical desorption by the frequent H+H→H₂ surface reaction, and CR-induced whole-grain heating. Out of these, sublimation primarily matters in the protostellar stage, while the other processes are mainly important for creating a feasible composition of the ice in the prestellar core that is inherited by the protostellar envelope. The rate of binary reactions and desorption is affected by the molecular desorption energy E_D, which varies depending on the surroundings of the molecule, whether it is a bare grain, strongly polar (H₂O), or hyper-volatile ices (CO).

As a dedicated chemical model, ALCHEMIC-VENTA has benefits for the calculation of the grain-size (ice-mantle thickness) evolution. It includes the UMIST RATE12 reaction net- work¹ and considers ice photoprocessing and desorption by several mechanisms, with rate coefficients that are updated according to recent experimental and theoretical results (KKV24). The most significant chemical aspect probably is the solid-phase conversion of CO into CO₂ and CH₃OH, that is, the chemical production of species with increasing sublimation temperatures T_subl, which is not easily reproduced by simplified functions (e.g., Rawlings & Williams 2021). This conversion is important because when the gas temperature T_gas rises, it allows repeated freeze-out of a major part of the C and O atoms into ice after the sublimation of the CO ice. Other chemical transformations of similar nature include N₂ → NH₃ and the CH₄ → C₂H₆→ longer carbon chains (KF24).

We studied events in a protostellar envelope that inherited its ices from the dark prestellar core. Table 3 lists the elemental abundances. Therefore, the initial conditions were set to be characteristic for these cores, with T_gas ≡ T_d,0.1 = 10 K, and the initial chemical composition was drawn from the prestellar stage model of KF24, using the elements of Table 3. Table 4 lists the initial abundance of the most important molecules (we used abundances relative to H (X/H), as used in the models, while in the results sections, the abundances are relative to H₂, n/n(H₂), as is often shown for astrochemical observations). The ISRF extinction was kept constant at A_V = 10 mag. The effects of the ISRF photons are still insignificant for the Standard model. The number density of H atoms was also kept constant, with n_H = 10⁸ cm⁻³ as the default value.

According to Eq. (1), T_d smoothly grows for all grains with increasing T_d,0.1, while the difference between T_d for grains in various size bins increases consistently. This may not always be the case as the incident spectrum of protostellar radiation is modified by absorption and scattering. Variations in the spectrum and intensity due to distance to the protostar and the amount of intervening material may change the relative efficiency for the heating of grains of different sizes. Furthermore, the optical properties of the grains also change when ices sublimate or refreeze onto their surfaces. These effects were not accounted for in our model and can affect the relative temperatures for grains of different sizes.

2.2 Model parameters

Visser et al. (2009); Drozdovskaya et al. (2014), and Furuya et al. (2017) indicated that for an infalling parcel of matter in a low-mass protostellar envelope, the rise and fall of the temperature occurs on timescales of about 20–100 kyr. We adopted the shorter timescale, divided equally between the heating (t_heat) and cooling (t_cool) episodes, as listed in Table 1. This is a conservative assessment because longer timescales would favour sublimation of ices and thus a more complete redistribution to colder grain-size bins. Therefore, t_heat and t_cool were not considered as variables in this study. The temperature increased from T_ini = 10 K to the peak gas temperature T_max over 10 kyr proportionally to the integration time squared t² (Garrod et al. 2008), while the decrease from T_max to the final temperature of T_end = 40 K followed a reverse curve. Because T_gas did not return to 10 K, the total integration time for the up-down temperature curve was 16.2 kyr, not 20 kyr. To study what happens after the heating episode, the simulations were continued for another 4 kyr with T_gas = T_end, up to a total time-span of 20 kyr.

We assumed that the parcel did not enter dense, cold, and shielded inner parts of the PPD during the simulation. Thus, A_V and n_H did not rise significantly, and T_end remained relatively high at 40 K in Model Standard. Possible mixing with other infalling matter (nearby parcels) was not considered. The simplified conditions with fixed A_V, n_H, T_max, t_heat, and t_cool were not intended to reflect conditions in real infalling proto- stellar gas. Instead, the fixed conditions, roughly representative of protostellar envelopes, and the simple T_d dependence on the grain size allowed us to understand the transfer of volatile matter through gas between grain populations with different T_d without other complex background physics. In addition to the Standard model, described above and in Table 1, a number of other simulations were carried out. A physical parameter, n_H, T_max (and T_end), the irradiation intensity factor I, the time spent at an elevated temperature t_plat, y, or the initial abundance of ices, was modified in each of these models relative to Standard. Table 5 lists all the models and the changes in their parameter values. The considerations for the selection of changed parameters and their values are described below in Sect. 3, along with the modelling results.

3 Results

The gas-grain chemistry of the prestellar stage produces icy mantles that are inherited by the protostellar envelope. To study the sublimation of ices, it is convenient to classify the ice species into groups that share similar sublimation temperatures. These groups of species are primarily the hyper-volatile species N₂, CO, and CH₄ ices, which already sublimate at T_d > 22 K, CO₂ ice that sublimates at T_d > 55 K, and NH₃, CH₃OH, and H₂O ices, which start to sublimate when T_d exceeds about 90 K. The latter two groups include a variety of other compounds, such as C₂H₂, H₂S, and HCN, which sublimate together with CO₂, as well as HCOOH, H₂O₂, and HF, which sublimate together with or after H₂O ice. In what follows, the primary focus is on the total abundance of the species that constitute the groups and their chemical transformation into molecules that belong to another group.

The modelling of hundreds of species and thousands of reactions allowed us to simulate a series of chemical changes in the species and their phase transition that affect the grain-size evolution, such as the partial re-freeze of sublimated CO in the form of CO₂, the re-freeze of N₂ in the form of NH₃, and the re-freeze of CO₂ oxygen in the form of H₂O. At the end of the prestel- lar stage, all grain sizes had a b of about 64…82 ice molecule monolayers (MLs; Table 2). A few percent of C and O are in many organic species that contain several heavy atoms, each with abundances relative to H₂ (n/n(H₂)) below 10⁻⁶.

3.1 Standard results

We start with the results of the Standard model. These were further used as a reference for the comparison with other simulations. The ice sublimation started with the desorption of dinitrogen N₂. Half or 50% of the N₂ ice relative to its initial abundance was thermally desorbed at an integration time t = 2.4 kyr and a gas temperature T_des = 26 K. The sublimation was not momentary; 99% N₂ ice disappeared in the T_gas interval between 19 and 29 K (t = 1.5…2.7 kyr), while 80% N₂ ice loss was observed at 23…28 K (2.1…2.7 kyr). All ice species showed sublimation periods of several 10² yr, with a temperature rise to about 10 K. These slow sublimation processes occur because the grain temperatures are spread out and we did not use a single T_d value.

The 40 K value for T_end in the Standard model is too high for N₂ to re-freeze again; thus it does not have a value for its 50% re-adsorption temperature T_re. The ice molecule that experienced re-accumulation during the cooling stage in almost each simulation was carbon dioxide, CO₂. In the Standard model, the temperature T_re for 50% CO₂ re-adsorption (relative to its initial abundance) is T_re,_CO2 = 72K at t = 13.9 kyr. Re-accumulation of CO₂ ice from 10% to 90% of its initial abundance occured from 75 to 69 K (13.7… 14.1 kyr), while freeze-out from 1% to 99% of the initial abundance occured in the 80…68 K interval (13.4…14.2 kyr). These numbers show that ice re-accumulation also takes several hundred years. In the case of CO₂, the ice accumulation continued and reaches 120% relative to its initial abundance because it is produced via surface reactions. Similarly to N₂, carbon monoxide, CO, remained in the gas-phase at the end of model Standard. The abundance of CO slowly decreased as it was converted into CO₂ ice. This change is insignificant in model Standard (only 2% of CO), but becomes a major CO sink under irradiated conditions (Sect. 3.4). Table 6 lists the T_des and T_re values for all models, along with the initial n/n(H₂) for all the major ice species N₂, CO, CO₂, methanol, CH₃OH, ammonia, NH₃, and water ice, H₂O. As discussed below, T_des and T_re are key parameters in understanding what occurred in each simulation.

Figure 1 shows that the sublimation of ices occurs sequentially. It starts in the hottest grain-size bin. After they are sublimated from a hot grain-size bin, the former icy molecules quickly re-freeze onto grains whose T_d has not yet reached their thermal desorption temperature. Thus, the coolest grains accumulate most of the molecules that were sublimated from other grain-size bins, until the coolest grains are too hot and also sublimate their ices. This process is repeated for the hyper-volatile species CO, N₂, and CH₄ and for less volatile species, such as CO₂. In the case of H₂O, NH₃, and CH₃OH, this process occurs near T_max, and the sublimation from the largest, coldest grains is interrupted by the start of the cooling stage. Reactive species, such as hydrogen peroxide, H₂O₂, and formaldehyde, H₂CO, are initially produced via photoprocessing (of H₂O and CH₃OH ices, respectively) and are mostly destroyed at elevated temperatures even before full sublimation.

Figure 2 shows that the sequential sublimation from the five grain-size bins produces a characteristic seesaw pattern for the gas-phase abundances of the less volatile molecules CH₃OH, NH₃, and H₂O in the heating stage of the Standard model (KF24). Table 2 lists the temperatures of the five grainsize bins when T_gas = T_max. As the temperature drops, only H₂O, NH₃, CH₃OH, and CO₂ re-freeze before T_end = 40 K. The re-formation of icy mantles exclusively occurs onto the largest, coldest grains. For these, the ice thickness b_0.232 reaches 363 MLs at the end of the simulation (t = 20 kyr), while it is below 1 ML for all the other grain-size bins.

The abundances of the re-accumulated ices are slightly changed compared to their initial values, as listed in Table 6: 106% for H₂O, 102% for CH₃OH, 105% for NH₃, and 120% for CO₂ ice. The increase arises from CO and N₂, which are converted into less volatile molecules that can be stored on icy grains. Similar solid-phase final abundances were also obtained for other models, unless noted otherwise.

3.2 Gas density

The hydrogen atom number density n_H has an assumed default value of 10⁸ cm⁻³ in Model Standard. Practically all of the H atoms are in H₂ molecules. The value of n_H can reach about 10¹² cm⁻³ during the transfer from the envelope to the disk (Visser et al. 2011). To investigate how these higher densities affect the transfer of ices between grain-size bins, we also considered n_H = 10¹⁰ cm⁻³ in Model Dens10 and 10¹² cm⁻³ in model Dens12. In addition, observations of complex organic molecules and warm carbon-chains in protostellar envelopes (Cazaux et al. 2003; Sakai et al. 2008) indicate that significant heating also occurs in less dense regions. For these, we used Models Dens6 and Dens7 with n_H of 10⁶ and 10⁷ cm⁻³, respectively.

In the higher-density Models Dens10 and Dens12, which are relevant for PPD conditions, the T_d differences for sub- µm grains are expected to be smaller than those derived with y = −1/6 (Podolak et al. 2011). To reflect a strong gas-grain thermal coupling, we assumed that y = −1/20 in these models. With this approach, the T_d difference for the smallest and largest grains at T_gas = 10 K is below 0.5 K (lower than that calculated with the sophisticated model of Gavino et al. 2021) and 16 K at T_gas = 150 K. The smaller differences in T_d may enable re-adsorption of gas-phase volatiles onto the surfaces of other grain-size bins, not just the coldest-grain population, because the other grain-size bins could also have sufficiently low T_d.

With a higher n_H, the molecules more often stick to the grains, which allows a faster growth of the icy mantles. In this way, higher n_H shifts the ice sublimation to higher temperatures. With y assumed to be −1/20, the lowest T_d is now higher relative to T_gas, which induces sublimation at lower T_gas. Table 6 shows that volatiles sublimate at higher T_gas in Model Dens10, where even H₂O ice is sublimated, unlike in the Standard model. In Model Dens12 the rate of molecule adsorption onto grains, which is proportional to n_H², is so high that no significant degradation of the H₂O ice abundance is observed. In other words, at a higher n_H, a higher T_d is needed for the thermal desorption rate to be higher than the rate of molecules that stick to grains. Moreover, some re-freeze near the end of the simulation at T_end = 40 K can also be experienced by the hyper-volatile molecules. This is shown by the re-accumulation of CO in Model Dens12, where the end abundance of CO ice reaches 2% of its initial abundance. In the two high-density models, a full migration of the ices to the coldest, largest 0.232 µm grains occurs, and b₀.₂₃₂ reaches 365 and 370MLs for Dens10 and Dens12, respectively.

In Models Dens6 and Dens7, the number density of gas and dust is lower by factors of 100 and 10 relative to the Standard model, respectively. The rate of molecules and atoms that stick to the grains is correspondingly lower by factors of 10 000 and 100. As a result, Fig. 2 shows that the coldest-grain population is unable to adsorb the volatiles before other grain-size bins become cold enough to grow their own icy mantles. In Model Dens7 most of the ice mass is still located on the large a = 0.232 µm grains with an ice thickness b_0.232 = 279 MLs. The next largest bin with a = 0.146 µm has b_0.46 = 96 MLs, the 0.092 µm grains have b_0.092 = 14 MLs, while for the two smallest and hottest grain-size bins (0.058 and 0.037 µm) b is below 1 ML. For Model Dens6, the ice mass is evenly distributed between all grain-size bins with an ice thickness between b_0.232 = 91 MLs and b_0.037 = 34 MLs. Table 6 shows that the sublimation of ices and especially their re-accumulation occurs at lower temperatures for the low-density models.

Table 6 shows that there is a major difference for T_re between the Standard model on one hand and Models Dens6, Dens7, and Dens10 on the other hand. The reason for the low T_re values is that these models desorb all their ice at T_max, exposing the bare-grain surface. The bare surfaces were assumed to be carbonaceous in our model, and did not form hydrogen bonds (KKV24). This means that the H₂O, NH₃, and CH₃OH surface molecules sublimate above T_d ≈ 80 K. A significant adsorption of these molecules is prevented until the Td of the coldest grain population allows accumulating ≈1 ML of ice. After this point, the ice accumulation is quick. An extreme case here is CH₃OH, which only re-freezes to 48% of the initial abundance because the gas in Model Dens6 is destroyed. On the other hand, most CH₃OH ice does not sublimate in Model Dens12 and can therefore reach an abundance of 144% relative to the initial value. The surplus abundance is reached during the first few thousand years because CO is hydrogenated on cold grains.

Fig. 1 Evolution of n/n(H2) for the major ice molecules on the surfaces of the five grain-size bins in the Standard model. The thicker lines show larger grains.

Fig. 2 Effects of nH on the ice sublimation and re-freeze with Models Dens6, Dens7, Standard, Dens10, and Dens12 from top to bottom. Left: ice thickness (black lines) changes with Td (red lines in the background). The thicker lines show larger grains. Right: n/n(H2) for ice (solid lines) and gas-phase (dotted lines) species.

Fig. 3 Effects of changes related to Tmax on the ice sublimation and re-freeze with models Tmax50 (top) and Tmax100 (middle) with Tmax 50 and 100 K, respectively. The bottom plots show the results from Model Plateau with tplat = 10 kyr. Left: ice thickness (black lines) and grain temperature (red; the thicker lines show larger grains). Right: n/n(H2) of major species in gas (dotted) and ice (solid lines). This is to be compared with Model Standard in Fig. 2, where Tmax = 150 K and tplat = 0. The scale changes for the secondary y-axis (Td) in the left plots.

3.3 Maximum temperature and its duration

Different trajectories from the envelope to the disk can have different T_max. Our default value of 150 K is among the highest (Drozdovskaya et al. 2014; Furuya et al. 2017) and was chosen because water ice is sublimated from most grain populations at this temperature, which allows the transfer of most of the ice mass between the grain-size bins. Lower maximum gas temperatures of 50 and 100 K were considered in models Tmax50 and Tmax100, respectively, at the same t_heat and t_cool.

Figure 3 shows that only the CO and N₂ ices sublimate in Model Tmax50. Both molecules also re-freeze onto the coldest grain-size bin surface before T_end = 20 K is reached. The ice thickness for the 0.232 µm grains changes from 64 MLs initially to 276 MLs at the end of the simulation. In Model Tmax100, CO₂ ice also sublimates and re-accumulates, while N₂ ice does not re-freeze at T_end = 30 K, and finally, b_0.232 = 251 MLs.

The Tmax50 and Tmax100 models differ in another way. The heating rate is lower by factors of 3 and 2, respectively, relative to the Standard model. This is because T_max is lower, while t_heat remains constant. Table 6 shows that this change has no substantial effect on T_des for the N₂ and CO ices, while for CO₂ ice, T_des is lower by 5 K.

A third model that considered changes in the heating timescales is Plateau, where we assumed that the parcel of infalling matter spends time t_plat = 10 kyr at T_max = 150 K. This simple change allowed us to investigate the redistribution of ices when all major volatiles were sublimated from all grain-size bins. This scenario never occurs in Model Standard. The first 10 kyr of Plateau are a copy of the same period in Model Standard. However, when the T_max plateau starts, H₂O ice is depleted below the 1 ML levels already during the first 500 years. The only other major change that occurs during t_plat is the destruction of CH₃OH in the gas (e.g., Millar et al. 1991). Upon re-freeze, the abundance of CH₃OH ice is restored to only 31% of its initial abundance (i.e. T_re,CH3OH is never reached). This decline of methanol is faster than that in the model of Garrod (2013), for example, because the CR ionization rate ζ_CR in our models is higher and because we used a different gas-phase chemical network.

The re-accumulation of the H₂O, NH₃, and CO₂ ices is more complete than that of CH₃ OH in Model Plateau, and similarly to the case of Standard, it exclusively occurs on the large, cold grains, where it reaches final abundances similar to those of Model Standard. The final b_0.232 is 362MLs. Figure 3 shows that ice re-accumulation occurs about 2250 yr after the end of the T_max plateau stage, when T_d for the large grains is only 84 K (corresponding to the re-accumulation gas temperature T_re,H2O = 96 K). This T_d is significantly below the sublimation temperature of H₂O ice (T_d ≈ 100 K). The ice accumulation occurs late because hydrogen bonds between ice species and bare-grain surfaces are lacking, similarly to models Dens6 and Dens7 (Sect. 3.2).

3.4 Irradiation

With I we regulated the intensity of ISRF and CRs. The value of I was assumed to be unity for the Standard model, where A_V = 10 mag and ζ_CR = 4.1 × 10¹⁷ s⁻¹ (Ivlev et al. 2015, KKV24). A detailed investigation of changes in radiation for infalling envelope material (see e.g. Drozdovskaya et al. 2014) is beyond the scope of this study. Our aim here was to determine the overall trend for the ice evolution under elevated irradiation. Therefore, highly irradiated conditions were represented with a single model that differed from Standard with I set to 100. This value shows significant chemical effects in the solid and gas phases, but it is not excessively destructive for the 10³..10⁴ yr long chemical timescales considered in the simulations. In the case of ISRF, I = 100 was achieved by lowering A_V from 10 to 5, while a 100-fold increase in CRs is expected to mimic some effects of energetic particles from the protostar. The same increase was also attributed to CR-induced photons.

Figure 4 and Table 6 show that evolution of the amount of ice and its thickness at the end of the simulation in Model Irrad100 is generally similar, although the underlying chemical processes differ significantly. The abundance of the ice molecules is reduced by three processes: thermal desorption, photodesorption, and ice photoprocessing. The latter results in fragments of photodissociated solid-phase species that combine into new molecules, part of which is desorbed. Thus, the loss of ices in this model can no longer be described as simply ‘sublimation’. Under the strong irradiation, the chemistry of major elements and species changes, as summarised below.

• Carbon and oxygen. Initially, CO ice is rapidly destroyed by photoprocesses and is already depleted at 18 K, when CO sublimation does not yet occur. H₂O ice and CO gas are combined into CO₂ and CH₃ OH ices, which both sublimate and are mostly destroyed in the gas. During the cooling stage, only CO₂ ice is restored, and at the end time of 20 kyr, the abundance of CO₂ has become comparable to that of H₂O ice, with n/n(H₂) of 1.7 × 10⁻⁴ and 2.4 × 10⁻⁴, respectively.

• Organic molecules. The synthesis of CH₃OH via hydrogenation of CO becomes ineffective when T_d > 24 K because atomic H, produced by H₂O photodissociation or adsorbed from the gas-phase, sublimates too rapidly. CH₃OH is fully desorbed and then destroyed by the gasphase chemistry, and does not re-form onto grains during the second half of the simulation. The gas-phase chemistry produces hydrocarbons, such as ethane C₂H₆, and those with desublimation temperatures above 40 K accumulate onto the grain surfaces.

• Nitrogen. As N₂ sublimates, N is slowly accumulated into the much less volatile NH₃ ice (Kalvāns & Kalnin 2019). During the cooling stage and especially towards the end of the simulation, other ice species, such as HCN and HC₃N, become abundant, although most N atoms still go into NH₃ ice.

During the simulation end phase with constant T_end = 40 K, the ice thickness continues to rise steadily and reaches b_0.232 = 356 MLs for the large grains and <25 MLs for other grain-size bins. This continued grain growth occurs because the major gasphase species CO and N₂ are constantly directly destroyed by radiation or ions, such as He⁺, which is produced by CRs. The ionization induced by ISRF is less significant because of molecular (self-) shielding effects. Some of the resulting atoms and atomic ions recombine into less volatile molecules that stick to the grain surfaces, where they continue to increase the ice thickness.

3.5 Grain temperature trend

We adopted a default value of y = −1/6 in the Standard model, that is, small grains with a large total surface area were hotter than large grains. This trend is supported by Draine (2003); Cuppen et al. (2006), and Krügel (2007). The opposite trend (e.g., Heese et al. 2017) with y = 1/6 was considered with Model Tpow1/6. KF24 concluded that the least likely scenario is that all grains have similar temperatures because T_d is affected not only by the grain size, but even more so by their material (see Sect. 1 and Woitke et al. 2024).

Figure 4 shows that the main difference in the results of Models Standard and Tpow1/6 is the early and rapid ice sublimation, followed by a late and rapid re-accumulation in Tpow1/6. This is shown quantitatively with the sublimation temperatures in Table 6. All ice is sublimated at 150 K. The first ice layer reappears on the small grains when their T_d drops to 81 K (at T_gas = 95 K). Similarly to Models Dens6, Dens7, and Plateau, the ice re-accumulation is delayed by the assumed non-existence of hydrogen bonds on the bare-grain surface. After the first ice ML has formed, the small grain-size bin adsorbs all the hydrogen-bonding less volatile ices within 10 yr.

The accumulation of all the hydrogen-bonded species increases the radius of the a = 0.037 µm grains to a + b = 0.067 µm, which is higher than the second-smallest grain-size bin with a = 0.058 µm. With y > 0 and size growth, the T_d of the small grains suddenly rises from 81 to 88 K before it continues to decline, as programmed in the code (Sect. 2.2). The rise delays the accumulation of the next volatile, CO₂, which explains why T_re,CO2 = 63 K in Model Tpow1/6 is lower than that for Model Standard, where it is 72 K. (We recall that T_des and T_re are gas temperatures that are indicative of the assortment for T_d of the five grain-size bins.) Similarly to Model Standard, the CO and N₂ ices do not re-accumulate in Model Tpow1/6.

While the temperature calculations for grains with different sizes and compositions indicate a significant spread in T_d (see Sect. 2.1), we also explored Model Tpow-1/20 with a narrow T_d range with y taken to be −1/20, similarly to Models Dens10 and Dens12. These similar grain temperatures may occur when all grains happen to share a similar composition and size. A narrow T_d spread reduces the possibilities for a redistribution of the ice mass between the grain-size bins (Sect. 3.2). Figure 4 shows that a full ice removal occurs in Model Tpow-1/20. The T_des and T_re values are lower than those of Model Standard because of the higher minimum T_d for large grains. During the cooling stage, the recovery of H₂O, NH₃, and CH₃OH ices occurs, and the additional delay is induced by the bare-grain surface, which does not participate in hydrogen bonds and thus has lower thermal desorption temperatures for these species (Sect. 3.2).

Fig. 4 Effects of a 100-fold increase in irradiation with Model Irrad100 (top) and the reverse grain temperature trend with Model Tpow1/6 (bottom) on ice loss and re-accumulation. Left: ice thickness (black lines) changes with Td (red). The thicker lines show larger grains. Right: ice (solid lines) and gas-phase (dotted lines) abundances of the major ice species relative to H2. This is to be compared with Model Standard in Fig. 2.

3.6 Distribution of the ices on the grains

Interstellar ice may be unevenly distributed among grains, and some grain-size bins may have thicker ice at the expense of others (Iqbal & Wakelam 2018; Silsbee et al. 2021). In order to determine whether this phenomenon affects the distribution of ice among the grain-size bins after a heating-cooling cycle, we considered two models in which all ice was either located on the surface of the smallest 0.037 µm grains (Model Ice_small) or of the largest 0.232 µm grains (Model Ice_large). In terms of the ice distribution, these two models can be viewed as fictional illustrating cases.

Figure 5 shows that the extreme initial distribution of ices on the surface of a single grain-size bin in Models Ice_small and Ice_large does not change the end result, the ice accumulation onto the coldest grains, that we obtained with the other models. The phase transition in Model Ice_large is similar to that in a single-grain model. All ices start, sublimate, and re-adsorb onto the surfaces of the coldest grain-size bin. No appreciable ice layer ever appears on the surface of other grainsize bins, where the ice thickness does not exceed 0.4 MLs at any point in time. The overall evolution of the chemical abundances closely resembles that of Model Standard, except for the lack of the seesaw pattern in the gas phase (Sect. 3.1). Table 6 shows that the T_des values are higher by a few degrees than those of the Standard model because no previous sublimation from warmer grains was possible in Model Ice_large. After T_max, all information on previous ice distribution is lost, and Models Standard and Ice_large show very similar results in the cooling stage.

The initial distribution of ices in Model Ice_large resembles the final distribution in Model Standard. Moreover, the final distribution of the ices of Model Ice_large is similar to the initial distribution of this same model. In other words, Model Ice_large can be viewed as the next heating-cooling cycle of a gas parcel after the cycle of Model Standard. This second cycle again results in ice that is exclusively accumulated on the surface of the coldest grain-size bin. This statement excludes the hypervolatile species that do not accumulate because T_end = 40 K prevents their adsorption. Nevertheless, Model Tmax50 indicates that these species re-accumulate in a similar manner, but onto the coldest grains.

Model Ice_small displays a more complex evolution than Ice_large, which nevertheless arrives at the same end result. The hyper-volatile species and CO₂ sublimate from the small 0.037 µm and re-freeze onto the surface of the larger, colder grains before they sublimate fully. This is similar to Model Standard, where a major part of the ice also is on the 0.037 µm grains. This process is not repeated for H₂O, NH₃, and CH₃OH because the larger and colder grains have bare, carbonaceous surfaces by the definition of Model Ice_small. This means that E_D for these species is significantly lower, and they cannot accumulate. Therefore, full ice sublimation occurs in Model Ice_small already when H₂O and the associated ices disappear from the 0.037 µm grains at T_gas = 114 K and T_0.037 = 134 K. This results in a delayed re-accumulation of the ices, as discussed in Sect. 3.2. The early desorption and late adsorption, compared to Model Standard, can be seen quantitatively with the T_des and T_re values in Table 6.

Fig. 5 Results of Model Ice_small, where all ice is on the surface of a = 0.037 µm grains (top) and Model Ice_large with ice on the 0.232 µm grains (bottom). Left: ice thickness (black lines) changes with Td (red). The thicker lines show larger grains. Right: ice (solid lines) and gas-phase (dotted lines) abundances of the major ice species relative to H2. This is to be compared with Model Standard in Fig. 2.

3.7 Summary

The results of Model Standard and of most of the other models show that, first, during the heating period, the molecules that were sublimated from warmer grains re-settle into the coldest grain-size bin before full sublimation occurs. Second, during the cooling period, ice irreversibly and exclusively re-accumulates into the coldest bin. This redistribution of ice into the coldest grain-size bin is valid for gas densities from 10⁷ cm⁻³, which is characteristic of dense, heated envelopes to 10¹² cm⁻³ in PPD interior. The large 0.232 µm grains have 12% of the entire bare- grain surface area in the model. For the redistribution to occur, the majority of the heated grains must exceed the sublimation temperature of the ice species. The hyper-volatile ices N₂ and CO disappear when T_d exceeds ≈34 K, and they are likely redistributed for all infall trajectories. For CO₂ ice, this temperature is >60 K, while for less volatile ices, such as H₂O, it is ≈124K. These T_d correspond to the temperature of coldest available grains, which means that T_gas will always be higher.

Under conditions of high irradiation, CO and H₂O ices are converted into CO₂. Therefore, more of the ice mass can sublimate and re-freeze at medium T_d of ≈60 K (up to 73 K in Model Dens12). Irradiation does not change the general pattern of ice redistribution onto the coldest grains. This is also true for models with altered grain-size T_d distributions, for instance, for Model Tpow1/6, where the largest grains have the highest T_d and Model Tpow-1/20 with a narrow T_d distribution. Unlike Model Standard, full ice sublimation occurs in both models, and the redistribution occurs when the infalling parcel of matter cools again.

A final aspect is that the re-freeze of the H₂O-like ices occurs at much lower T_d in cases of full ice sublimation. This is because the bare-grain surface was assumed to be carbonaceous and does not form hydrogen bonds in our model. Therefore, T_re, which is necessary for re-forming the first full ice ML, is lower by about 30…40 K than T_des. In a real protostellar envelope, the surfaces of the refractory grains can also consist of silicates or other materials. Along with the T_d differences between grain materials (KF24), this means that the grain surface material is another variable that can promote accumulation of ice onto a minority of grains with specific properties.

4 Conclusions

To summarise the results, the selective accumulation of ice onto the surface of the coldest-grain population during a heatingcooling cycle is a strong rule. It cannot be easily changed by variations in the physical or chemical conditions of matter falling onto a PPD from a protostellar envelope. This means that most grains enter the PPD ice poor, while a small minority are ice rich. This uneven ice distribution can affect the grain coagulation, which is more likely and at higher collision speeds, when the grains are covered with an ice layer (Gundlach & Blum 2015; Nietiadi et al. 2017, 2020). As far as we know, the specific case of collisions between ice-covered and bare refractory grains has not yet been investigated. The relative grain velocities can also be significantly affected by their size differences (Ormel et al. 2009; Silsbee et al. 2020), which can be exaggerated or reduced by the selective accumulation of ice.

In prestellar cores, grains grow by accumulating icy mantles (Whittet et al. 2001), while rapid further grain growth by coagulation is observed already in protostellar envelopes (Miotello et al. 2014; Testi et al. 2014). Coagulation is affected by the dual nature of icy and bare grains that are present in the same place and at the same time. Heating events affecting PPDs, such as the event observed by Burns et al. (2023), can also be expected to irreversibly produce dual grain populations (bare and ice rich), which affects the behaviour of grains in the disk, as described by Podolak et al. (2011). The spread-out snowlines in PPDs, predicted by Gavino et al. (2021), likely do not occur, as the ices always sublimate from and freeze out onto the same population of the coldest grains as in Model Ice_large.

The full sublimation of ice and its subsequent reaccumulation will result in fully stratified icy mantles with species of increasing volatility towards the surface. This aspect was not simulated in this model. This ice structure facilitates the rapid sublimation of ice components with no need to diffuse through layers of less volatile ices. Moreover, the segregation into relatively pure ices may hamper the synthesis of some complex molecules, including organics. For example, one route for CH₃OH production is by photoprocessing a H₂O-CO icy mixture, which cannot be expected to be active in stratified ices.

Acknowledgements

This research is funded by the Latvian Science Council grant “Desorption of icy molecules in the interstellar medium (DIMD)”, project No. lzp-2021/1-0076. I am grateful for the support of Ventspils City Council. This research has made use of NASA’s Astrophysics Data System. I thank the anonymous referee for efforts that significantly improved the paper.

References

All Tables

All Figures

	Fig. 1 Evolution of n/n(H2) for the major ice molecules on the surfaces of the five grain-size bins in the Standard model. The thicker lines show larger grains.
In the text

	Fig. 2 Effects of nH on the ice sublimation and re-freeze with Models Dens6, Dens7, Standard, Dens10, and Dens12 from top to bottom. Left: ice thickness (black lines) changes with Td (red lines in the background). The thicker lines show larger grains. Right: n/n(H2) for ice (solid lines) and gas-phase (dotted lines) species.
In the text

Fig. 3 Effects of changes related to Tmax on the ice sublimation and re-freeze with models Tmax50 (top) and Tmax100 (middle) with Tmax 50 and 100 K, respectively. The bottom plots show the results from Model Plateau with tplat = 10 kyr. Left: ice thickness (black lines) and grain temperature (red; the thicker lines show larger grains). Right: n/n(H2) of major species in gas (dotted) and ice (solid lines). This is to be compared with Model Standard in Fig. 2, where Tmax = 150 K and tplat = 0. The scale changes for the secondary y-axis (Td) in the left plots.

In the text

Fig. 4 Effects of a 100-fold increase in irradiation with Model Irrad100 (top) and the reverse grain temperature trend with Model Tpow1/6 (bottom) on ice loss and re-accumulation. Left: ice thickness (black lines) changes with Td (red). The thicker lines show larger grains. Right: ice (solid lines) and gas-phase (dotted lines) abundances of the major ice species relative to H2. This is to be compared with Model Standard in Fig. 2.

In the text

Fig. 5 Results of Model Ice_small, where all ice is on the surface of a = 0.037 µm grains (top) and Model Ice_large with ice on the 0.232 µm grains (bottom). Left: ice thickness (black lines) changes with Td (red). The thicker lines show larger grains. Right: ice (solid lines) and gas-phase (dotted lines) abundances of the major ice species relative to H2. This is to be compared with Model Standard in Fig. 2.

In the text