© O. Sipilä et al. 2020

1 Introduction

Interstellar dust grains play a large role in the chemical and physical evolution of molecular clouds that host starless and pre-stellar cores – the seeds of low-mass star formation. For example, the grains regulate the gas temperature in high-density regions, where collisions between grains and gaseous molecules are frequent (e.g., Goldsmith 2001); these grains also participate in shaping the chemical inventory in the clouds as new molecules are formed on their surface (see, e.g., Hama & Watanabe 2013 for a review).

The physical size of the dust grains varies from nanometers to some tenths of micrometers, or even larger depending on the environment. Various models of their size distribution and optical properties as a function of radius and bulk composition have been introduced in the literature (e.g., Mathis et al. 1977; Li & Draine 2001; Weingartner & Draine 2001; Jones et al. 2013). Yet in gas-grain chemical models used to predict molecular abundances in interstellar matter (ISM), it almost ubiquitously assumed that the grains are monodisperse; that is, they all have the same radii, which is usually taken to be around 0.1 μm (e.g., Hasegawa et al. 1992; Garrod & Herbst 2006; Taquet et al. 2012; Sipilä et al. 2015; Vasyunin et al. 2017). This assumption greatly simplifies the chemical modeling, but has negative consequences for the accuracy of the simulations. We highlight two issues here. First, grains are transiently heated to high temperatures by cosmic ray (CR) impacts which deposit energy into the grain; the amount of deposited energy depends, among other properties, on the physical size of the grain (Fano 1963; Léger et al. 1985). Nanometer-sized grains can be heated to hundreds of K, while sub-micrometer grains are only heated to a few tens of Kelvin (Herbst & Cuppen 2006). The transient heating allows efficient nonthermal desorption of molecules from the grain surface (Hasegawa & Herbst 1993). Second, the equilibrium temperature of a grain varies with its size. Considering that the efficiencies of chemical reactions on grain surfaces are strongly dependent on the grain temperature (e.g., Hasegawa et al. 1992), coupled with the fact that nonthermal desorption allows the surface material to interact with the gas phase, it is highly desirable to investigate the effect that size distribution may have on the chemical evolution of a star-forming object as a whole.

The effect of a grain size distribution on chemistry has been investigated in some previous studies, using models neglecting the effect of the size distribution on grain temperatures (Acharyya et al. 2011), including size-dependent equilibrium temperatures (Pauly & Garrod 2016; Ge et al. 2016) and grain-size dependence in CR heating efficiency (Zhao et al. 2018; Iqbal & Wakelam 2018). Of these studies, only the work of Iqbal & Wakelam (2018) considers size-dependent equilibrium temperatures and CR heating simultaneously, but their study was limited to a single set of physical conditions and lacked a treatment of size-dependence in grain cooling following a CR impact. Furthermore, the studies of Pauly & Garrod (2016) and Ge et al. (2016) considered a simple power-law dependence for the equilibrium dust temperature, or used temperature data appropriate for grains without ice mantles, respectively. Given that ice mantles already start to form on grains at low medium densities at which we can reasonably expect size-dependent variations in grain equilibrium temperatures, owing to the different response of the grain populations to heating by the interstellar radiation field (ISRF), there is a definite incentive to build a more complete model that ties the above together.

In this work, we present the first gas-grain chemical model that examines simultaneously the effects of grain-size-dependent CR desorption and grain equilibrium temperature on gas-phase and grain-surface abundances, applied to a variable set of physical conditions corresponding to those found in starless and pre-stellar cores. The heating and cooling timescales of the grains following CR impact events are scaled with the grain size, and the grain equilibrium temperatures are calculated with a radiative transfer code using dust models tailored for the present purposes. We track the ice abundances on each grain population explicitly to constrain the relative amounts of molecules as a function of grain size. We use the model to investigate ice abundances at the time of their formation in a low-density medium, and we also look into the relative ice abundances in high-density conditions preceding the protostellar stage.

The paper is structured as follows. Section 2 describes our model in detail. Section 3 presents the results of our model runs applied to single-point physical conditions and to the pre-stellar core L1544. We discuss the results in Sect. 4, where we also compare our work to the previous related studies available in the literature. Section 5 presents our conclusions. Appendices A–C present some additional results that complement the main text.

2 Model

In this paper, we investigate the effect of a grain size distribution on molecular abundances in physical conditions corresponding to those found in starless and pre-stellar cores. Specifically, we study how the CR desorption process and the grain equilibrium temperatureare affected by variations in grain size, and how these effects in turn influence gas-phase and grain-surface chemistry. This section presents the theoretical basis of our work and details the chemical and physical models that we employ.

2.1 Monodisperse grains versus grain size distribution

Gas-grain models of interstellar chemistry nearly universally assume monodisperse spherical dust grains; the grain radius is usually taken as 0.1 μm. In such a case, the grain abundance with respect to hydrogen can be calculated from (1)

where n_H is the total hydrogen number density (≈2n(H₂) inside dense clouds), m_H is the mass of the hydrogen atom, μ is the average particle mass, ρ_d is the grain material density, a_m is the monodisperse grain radius, and R_d is the dust-to-gas mass ratio. In this paper, we take μ = 1.4, ρ_d = 2.5 g cm⁻³, and R_d = 0.01.

If we instead consider a size distribution for the grains, a size scaling relation must be assumed. In this work, we adopt the grain size distribution from Mathis et al. (1977; hereafter MRN), where d n_d ∕da = Cn_Ha^−3.5. In this distribution,C is a normalization constant whose value can be derived under the assumption of a constant dust-to-gas mass ratio as follows: (2)

where a_max and a_min are the maximum and minimum radii considered for the distribution, respectively. The MRN distribution is continuous, but it can be easily discretized for inclusion in a chemical model by dividing the size interval into a number of grain size bins, thereby calculating the effective grain number density (3)

and the effective grain radius (4)

in each size bin i. Finally, the effective grain area in each bin is given by (5)

Presently we are interested in a comparison between models in which we assume either monodisperse grains or a grain size distribution. These are two inherently different cases. To facilitate the comparison, we set the monodisperse grain radius by requiring that the resulting total grain surface area () matches the total effective surface area of the size distribution (). In this paper we assume for the size distribution limits, leading to a monodisperse grain radius of a_m = 0.095μm based on the area constraint. The grains in the size distribution model are always assumed to be spherical.

Equalizing the total grain-surface area between the two models means that the (total) grain abundance is not the same in both cases (). An alternative approach would be to set the monodisperse grain radius so that the grain abundance equals the total abundance in the size distribution model. But we consider this approach inferior because it would lead to a discrepancy in the total grain surface area, which is of critical importance to the present study as adsorption rates depend on the grain area. We note that the dust-to-gas mass ratio is never varied from its fiducial value of 0.01, meaning that the size distribution and monodisperse grain models correspond to the same total dust mass.

2.2 Desorption induced by cosmic rays

Hasegawa & Herbst (1993) derived an expression for the rate coefficient of a desorption induced by CRs (hereafter simply “CR desorption”) event, which for the desorption of a grain-surface species j takes the form (6)

where k_therm is the thermal evaporation rate coefficient at 70 K. This expression is based on the argument that energetic iron nuclei deposit enough energy when passing through a grain with radius 0.1 μm to heat it to 70 K, and that desorption effectively takes place at 70 K; the molecular evaporation cools the grain back to 10 K. The f(70 K) efficiency term is the grain “duty cycle” at 70 K, defined as the ratio of the grain cooling timescale (R_cool ~ 10⁻⁵ s) to the time interval between successive CR heating events for 0.1 μm grains, which Hasegawa & Herbst (1993) calculated to be R_heat = 3.16 × 10¹³ s based on the iron CR flux given by Léger et al. (1985). The heating event interval scales as a⁻² (Léger et al. 1985), while the cooling timescale is independent of grain radius; the latter scales with the temperature as exp (−E_b∕T), where E_b is the binding energy of a typical volatile grain-surface species, which Hasegawa & Herbst (1993) fixed to 1200 K. We thus obtain the following for the grain-size-dependent duty cycle: (7)

where T_max is the radius-dependent transient maximum temperature that the grain reaches upon a CR impact, f(0.1 μm, 70 K) = 3.16 × 10⁻¹⁹ is the duty cycle for 0.1 μm grains (Hasegawa & Herbst 1993), and the last two terms represent the scaling of R_cool and R_heat, respectively. The quantity T_max is derived from the grain heat capacity (Léger et al. 1985). Small grains are transiently heated to temperatures in excess of 150 K (the upper limit in Léger et al. 1985); we used the extrapolated heat capacity from Zhao et al. (2018) to calculate T_max for the smallest grains in the distribution.

We divided the grain size distribution into five logaritmically spaced size bins and derived the effective grain radius and abundance in each bin using the formulae presented in Sect. 2.1. The transient maximum temperatures, successive heating intervals, and duty cycles were then calculated based on the effective grain radii. The results of these calculations are collected in Table 1, in which we also list for reference the CR desorption rate coefficient of CO, taking 1150 K for the CO binding energy; CR desorption rate coefficients increase with grain radius. The duty cycle f(a, T_max) varies by more than ten orders of magnitude over the size distribution and obtains a very low value for the smallest grains, but the CR desorption rate coefficient varies only by approximately two orders of magnitude over the distribution; this is because the low duty cycle values are compensated by the high transient maximum temperatures for the smallest grains. The number of grain size bins has an effect on our results, but only for the relative ice abundances in the various grain size bins and not for the total ice abundances; we discuss this issue in more detail in Appendix A.

2.3 Size-dependent equilibrium grain temperature

The grain size affects the equilibrium temperature (T_dust) of the grains. Chemical reactions on grain surfaces are very sensitive to variations in T_dust owing to two factors: first, the thermal diffusion rate on the surface scales as , where is the binding energy of species j; and second,barrier-mediated surface reactions are scaled analogously by the inverse of T_dust. Therefore, allowing the equilibrium grain temperature to vary over the size distribution affects the reactivity on the different grain populations. For this reason we explicitly tracked the chemical abundances on the grains in each size bin; the chemical model thus contains five versions of all surface species in our current setup of five grain size bins.

In several previous works concerning models of the structure of starless and pre-stellar cores (e.g., Sipilä et al. 2017), we used opacity data from Ossenkopf & Henning (1994) to derive T_dust profiles as a function of density. In the present work, we cannot employ the same approach because the Ossenkopf & Henning (1994) tables are derived over a broad range of grain sizes (from 5 × 10⁻³ to 0.25 μm), and we are considering a setup with multiple grain size bins. We therefore generated new sets of absorption and scattering cross sections for the five size bins shown in Table 1, and over a broad distribution from 0.03 to 0.3 μm, using the SIGMA code (Lefèvre et al. 2020¹; Woitke et al. 2016; Min et al. 2005). We assumed that the grains consist of a carbon/silicate core in a 50/50 volume ratio (matching the grain material density of ρ_d = 2.5 g cm⁻³ in the chemical code) and are coated with a water ice layer with a volume fraction half that of the core. The refractive index data is taken from Ossenkopf & Henning (1994) (included in SIGMA) such that this setup corresponds roughly to the “thin ices” model of Ossenkopf & Henning (1994). Radiative transfer calculations to determine T_dust profiles were carried out using the CRT code (Juvela 2005); for these, we assumed that the core is surrounded by an external layer corresponding to a visual extinction of A_V = 1 mag, and we took the ISRF spectrum from Black (1994). The results of these calculations are presented in Sect. 3.2.

2.4 Additional model details

The results presented in this paper pertain to four individual model cases. These are collected in Table 2. Model M1represents the standard case of monodisperse grains, which is almost universally adopted in gas-grain models of ISM chemistry. In model M2, we instead consider a grain size distribution, assuming that T_max and the duty cycle are constant for all grain sizes. This model is included specifically to illustrate the effect of ignoring thescaling of the CR desorption rate coefficients with grain radius and temperature. In this case, T_max,m = 73.86 K and f(0.095 μm, T_max,m) = 1.16 × 10⁻¹⁹ are, respectively, the temperature and duty cycle obtained for the monodisperse grains in model M1, so that models M1 and M2 are comparable. Model M3 represents the fully radius-dependent distribution model as detailed above (cf. Table 1). Model M4 introduces into M3 the grain-size-variable T_dust calculated using CRT.

We used the gas-grain chemical code discussed in Sipilä et al. (2019; see also references therein), updated to accommodate for the grain size distribution. In this paper we employed point-like physical models to demonstrate the effects on the chemistry due to different grain setups. For simplicity we adopted a constant temperature T_gas = T_dust = 10 K in these models; (model M4 is not considered), and a constant visual extinction A_V = 10 mag. We varied the density of the medium between n(H₂) = 10⁵ and 10⁷ cm⁻³.

In Sipilä et al. (2019), we presented the results of an extensive benchmarking study attempting to reproduce the observed distribution of ammonia in the pre-stellar core L1544, for which we used the physical model by Keto & Caselli (2010). We construct a similar model in this work for the purposes of studying the size distribution effects in the context of a real core that has a nonconstant density and temperature structure. We derive time-dependent chemical abundance profiles in L1544 with dust temperature profiles designed specifically for the present work (see Sect. 3.2) in connection with the density and gas temperature structure from Keto & Caselli (2010). We divide the so-obtained physical source model into concentric shells and calculated chemical evolution separately in each shell. This process of deriving time-dependent abundance profiles is described in more detail in Sipilä et al. (2019); we employ here the same initial chemical abundances (Table 3) and the same chemical networks as in that paper². An A_V profile, important for the photoreactions in the chemical model, was constructed assuming an external visual extinction of (matching the radiative calculations). We stress that the L1544 core model is adopted just so we can explore the influence of the grain size distribution on chemical evolution using a core model with a high central density; in this work we do not attempt to reproduce any observations toward L1544, and our present aims could be fulfilled using any generic core model, such as a Bonnor-Ebert sphere.

3 Results

3.1 Effect of grain size distribution and size-dependent CR desorption on molecular abundances

We proceed with a step-by-step analysis to uncover the effect that the grain size distribution has on molecular abundances in dense core conditions. Our starting point is the comparison of models M1 and M2; the twoare equivalent except for the introduction of a distribution of grain sizes in the latter model. Figure 1 shows the abundances of some selected gas-phase species as a function of time at three H₂ densities. The plot demonstrates clearly that simply introducing a distribution of grain sizes, and keeping the CR desorption scheme unaltered, has no significant effect on the abundances of the various species. This conclusion holds for all species in our model and thus is not limited to the species that we selected for Fig. 1, although the move to a distribution of grain sizes naturally means that grain-surface species are spread out over several populations of grains (see below).

The effect of the grain size distribution becomes evident as we introduce the more realistic case of grain-size-dependent transient maximum temperatures and duty cycles (model M3), although the magnitude of the effect varies with each species. Fig. 2 shows a comparison of models M1 and M3. The size distribution model predicts an increase in abundances for the various species until t ~ 10⁵ yr, compared tothe monodisperse grains model, with an increasing effect toward higher densities. At longer timescales the situation is – in general – reversed. These differences can be understood by careful examination of grain-surface molecular distributions and their effect on the gas phase. The surface distributions can be highly nontrivial when the grain size is allowed to vary, depending on the species and on the time.

We show in Fig. 3 the abundances of CO, N₂, and H₂O ice as functions of time. Let us consider N₂ in the lowest-density model (n(H₂) = 10⁵ cm⁻³) at t = 10⁴ yr. At this time, gas-phase N₂ is formed mainly through the N + CN→C + N₂ reaction in both cases (M1 and M3). However, the contribution of CR-desorbed icy N₂ to the gas-phase N₂ abundance is only less than 1% in model M1, but about 10% in model M3. This accounts for the slightly higher gas-phase N₂ abundance in the latter model. Formation of N₂ in the ice competes against hydrogenation and desorption of atomic N. Hydrogenation is the most efficient destruction process of atomic N on all grains, but its relative efficiency diminishes toward larger grains because of the associated increase of CR desorption rate coefficients (Table 1). As the grain size increases, atomic H is desorbed more effectively, boosting N₂ production through the N + N association reaction. However, for the largest grains, CR desorption of atomic N is efficient enough to start inhibiting N₂ formation. Therefore, there exists a “sweet spot” for grain-surface N₂ production, which in our model occurs for grains slightly larger than 0.1 micron. This is clearly evident in Fig. 3, where grain-surface N₂ is in model M3 mainly on the grains in the fourth size bin (a_eff = 0.146 μm) at t = 10⁴ yr. The CR desorption rate coefficients in model M1 are smaller than those on the largest grains in model M3, and so the total effect of CR desorption is greater in the latter model at early times. As the medium density is increased, CR desorption events occur more often, and gas-phase N₂ formation is dominated by CR desorption in both models M1 and M3; the N₂ ice abundance is highest on the largest grains only for a period of a few hundred years at n(H₂) = 10⁷ cm⁻³.

At later times, between t ~ 10⁵ yr and a few × 10⁶ yr in the model with n(H₂) = 10⁵ cm⁻³, model M3 predicts a lower gas-phase N₂ abundance (Fig. 2), even though N₂ ice is more abundant than in model M1. The difference is again directly related to grain chemistry. The main formation pathway for gas-phase N₂ at late times is N₂ H⁺ + e⁻→N₂ + H, and its main destruction pathway is , constituting a loop between N₂ and N₂H⁺. This applies to both models M1 and M3, and in fact the abundance of is, at late times, somewhat higher in model M3 than in M1. An inspection of Fig. 3 reveals that N₂ ice is abundant only on the smallest grains at t > 10⁵ yr because of the high CR desorption efficiency on large grains, and so the net CR desorption rates are at late times higher in model M1 than in model M3.

The distribution of CO ice is straightforward: the larger the grain, the lower the CO ice abundance. CO is formed mainly in the gas phase, and the majority of CO ice is material deposited from the gas phase through adsorption (unlike N₂ ice which is formed in situ at early times). Adsorption rates are largest on small grains with the greatest surface area, and hence the CO ice abundance simply decreases with increasing grain radius. After t ~ 10⁵ yr when CR desorption becomes important, larger abundance variations between grains of different radii develop owing to the increase of the CR desorption efficiency toward larger grains (Table 1). This effect does not impact CO as greatly as N₂, owing to itshigher binding energy. The total CO ice abundance profile displays a clear peak at n(H₂) = 10⁵ cm⁻³, while remaining flat throughout the time evolution at n(H₂) = 10⁷ cm⁻³. The efficiency of CO ice processing is tied to the abundance of gas-phase atomic H, which is roughly inversely proportional to the density of the medium. At n(H₂) = 10⁵ cm⁻³, the total CO ice abundance decreases past ~10⁵ yr as a result of hydrogenation processes that convert CO into formaldehyde and further into methanol. The production of methane is also effective. At higher densities, fewer H atoms are available for hydrogenation, and the decline in CO abundance occurs only to a rather limited degree; the abundances of formaldehyde and methanol ice are correspondingly much lower than at n(H₂) = 10⁵ cm⁻³.

Much like that of CO, the distribution of H₂O ice is in these physical conditions a simple decreasing function of the grain size label. However, because of the high binding energy of water ice (5700 K in our model), no significant time-dependence is exhibited in the abundance ratios between ices on the different grain populations.

Figure 2 shows that the gas-phase ammonia abundance is slightly boosted by the introduction of a grain size distribution at high density, and does not experience a similar drop with time with respect to model M1 as the CO and N₂ abundances do. In general, how the abundance of a given species responds to alterations in the grain model is dependent on several factors, such as whether that species is formed mostly in the gas phase or on grains and the volatility of the species. (We return to this issue in Sect. 4.2.) Nevertheless it is clear already from the very small sample of species represented in the preceding figures that switching from monodisperse grains to a size distribution has a marked effect on the overall chemistry, especially at high density. However, we note that in this work we started from initially atomic gas (except for hydrogen and deuterium which are initially in H₂ and HD) to better bring out the density dependence in the effect that the different grain models have on the chemistry. This of course represents an extreme case, and in reality we expect the gas density to increase gradually as the pre-stellar core forms and ultimately collapses, which may mean that differences between models M1 and M3 are not as great as those shown in Fig. 2.

Fig. 1 Abundances of selected gas-phase species, indicated in the figure, as a function of time for three medium densities. The solid lines correspond to model M1, while the dashed lines correspond to model M2.

Fig. 2 As described Fig. 1, but in this figure the dashed lines correspond to model M3.

Fig. 3 Distributions of CO (top row), N2 (middle row), and H2O (bottom row) ice as functions of time at three different medium densities. The green solid lines correspond to model M1. The red andblue lines correspond to model M3; the red lines represent the abundances on each grain population, while the blue line represents the sum over the individual grain populations. The green and blue lines overlap almost perfectly for H2O ice. The effective grain radius increases with bin label as indicated in the top right panel.

3.2 Grain-size-dependent equilibrium dust temperature and molecular abundances in L1544

Next, we examine the differences in the predictions of models M1, M3, and M4 when applied to the pre-stellar core L1544. Figure 4 presents T_dust profiles obtained using CRT, when applying the dust model discussed in Sect. 2.3 to the L1544 density structure. We also show the T_dust profile of Keto & Caselli (2010) for reference. The T_dust profile calculated over the full grain size distribution is used in models M1 and M3. In the case of model M4, the equilibrium T_dust increases with grain size, which holds throughout the range of densities in the core. The difference between the temperatures of the smallest and largest grains (ΔT) is the largest between n(H₂) = 10⁴ cm⁻³ and n(H₂) = 10⁵ cm⁻³, where ΔT ~ 1 K, and in general ΔT is close to 1 K in a broad density range from the lowest densities up to ~10⁵ cm⁻³. The low endof this range represents the medium densities where ices form in the first place, while at the high end molecular freeze-out becomes important. The T_dust profile calculated over the entire distribution is largely similar to those of the individual grain populations, but presents a somewhat different slope in the n(H₂) = 10⁴ cm⁻³ to n(H₂) = 10⁶ cm⁻³ density range. We note that the dust temperature profiles calculated here are clearly different from that obtained by Keto & Caselli (2010). The new dust model constructed for the present work is not meant to reproduce observed temperatures or (submillimeter) intensities toward L1544, but is merely designed to investigate the effect of variations in equilibrium T_dust on ice abundances.

A gas temperature profile must be supplied to carry out chemical calculations using the L1544 density structure. The gas temperatureaffects grain-surface chemistry through its influence on adsorption rates, which depend on the thermal speed of the molecules, and it naturally impacts the synthesis of molecules in the gas phase. Our aim is to investigate the effect that the various approaches to the treatment of grain heating have on surface chemistry, and thus for models M1 to M4 to be maximally comparable,for simplicity, we adopt the gas temperature profile of Keto & Caselli (2010) in all cases. It would be possible to calculate new gas temperature profiles for the present purposes using, for example, the hydrodynamics code presented inSipilä & Caselli (2018) and Sipilä et al. (2019), but because the gas temperature also depends on the dust temperature through the dust-gas collisional coupling, the gas temperature profiles in models M3 and M4 would not be equal. Moreover, carrying out this calculation is not straightforward for model M4 where there are several dust temperature components to consider. We thus retain the Keto & Caselli (2010) gas temperature profile for L1544 so that the effect of varying dust temperature can be more easily distinguished between the models.

Figure 5 shows selected molecular abundance profiles in the L1544 model at different durations of chemical evolution starting from the (mostly) atomic initial state. The gas-phase N₂ and CO abundances follow the trends shown in Fig. 2: N₂ is strongly depleted at high densities (n(H₂) ≳ 10⁵ cm⁻³ in models M3 and M4 from ~10⁵ yr onward), while the CO abundance takes a somewhat longer time to drop below the model M1 value. There is no significant difference between models M3 and M4 for the gas-phase species.

Deviation from the results shown in Fig. 3 is evident for N₂ ice; the abundance is a factor of 2–3 lower in models M3 and M4 than in model M1. Water ice presents similar behavior. This effect is tied to the response of grain chemistry to the temperature (~6 versus 10 K in Fig. 3). Decreasing the temperature allows atomic hydrogen to subsist longer on the grain surfaces, increasing hydrogenation rates. There is a consequent overall boost to the abundances of, for example, nitrogen hydride ices when the temperature is decreased from 10 to 6 K. The difference between models M1 and the two size distribution models is due to grain-size-dependent CR desorption efficiency. In models M3 and M4, in which the CR desorption rates depend on the grain size, the (average) CR desorption rate of atomic hydrogen is lower than in model M1, leading to higher nitrogen hydride ice abundances and hence a lower N₂ ice abundance in the size distribution models. The gas-phase abundance of ammonia decreases with the more efficient destruction of N₂ ice in models M3 and M4 (Fig. 5). The response of the models to changes in T_dust is highly nontrivial and highlights the need for precise dust temperature modeling as opposed to assuming that the temperature is constant (typically 10 K) in the shielded high-density regions inside pre-stellar cores.

Figure 6 presents a breakdown of the abundances of CO, N₂, and H₂ O ice as functions of medium density at t = 10⁴ yr and t = 10⁵ yr in models M3 and M4. We show the relative ice abundance in each grain size bin, that is, the abundance in each bin divided by the total abundance (sum over all bins). We also show in one panel the relative abundances of each grain population for reference. The plot clearly indicates that even though there is no significant difference in the total ice abundances predicted by the two models (Fig. 5), there are large variations in the distributions over the grain sizes depending on time and the medium density.

The distribution of ices over the grain populations tends to be a straightforward function of the grain size in model M3 over the entire density range, that is, the ice abundances decrease with grain radius. We highlight two notable features: First, for N₂ at early times we recover the same trend evident in Fig. 3 where the ice is most abundant in bin 4; we note, however, that this trend does not appear at n(H₂) ≲ 10⁴ cm⁻³); and, second, the relative water ice abundance is heavily dependent on the medium density, in stark contrast to the results shown in Fig. 3. This is caused by the grain temperature gradient across the core – in the center where T_dust is low, CR desorption dominates over hydrogenation processes on large grains. Hence the water ice content is heavily concentrated on small grains, whereas the abundances of carbon-containing molecules tend to be relatively higher on large grains.

The abundance profiles are much more complicated at the low medium densities in model M4, where the differences in T_dust between the size bins strongly affects the efficiencies of chemical processes on the grain surfaces. Most of the N₂ and CO ice tends to be incorporated in the ice in size bins 3 and 4, which is in stark contrast to the distribution in model M3. Water ice, being the end product of hydrogenation processes and very strongly bound to the grain surface, does not exhibit similar complicated size-dependent abundance variations. The predictions of models M3 and M4 are in near-perfect agreement at densities above n(H₂) ~ 10⁵ cm⁻³, where the temperature difference between the smallest and the largest grains is small.

A notable result, shown in Fig. 6, is that the relative ice abundances do not follow the relative abundances of the grains themselves. The general trend is that the smallest grains are under-represented in this regard; that is, the relative ice abundance is lower than the corresponding grain abundance. In the other size bins, the trend is reversed, except for the largest grains at high density. The consequences of introducing a grain-size-variable equilibrium T_dust on the relative ice abundances are very interesting, because the effect is most pronounced in the density regime where ices start to form, in regions (partly) shielded from external UV radiation. The initial chemical abundances in our present model, in which the gas is initially (mostly) atomic, are not appropriate for the inner, high-density regions of the core model. However, our results do imply potential repercussions to the ice composition in the inner regions of pre-stellar cores as well because the gas density increases through (quasi-static) contraction starting from lower-density gas. We test the effect of different initial conditions in Sect. 4.1, in which we also explore the role of the external radiation field on chemical evolution because of its effect on the equilibrium dust temperatures. Finally, we note that the thickness of the ice mantle that accumulates on the surface of a grain depends on the grain radius as well. We include for completeness some additional discussion on this topic in Appendix B.

Fig. 4 Equilibrium dust temperature profiles as a function of density, calculated for the five distinct grain populations (black lines; model M4) or over the entire distribution (green; models M1 and M3). The effective grain radius increases with bin label as indicated in the plot. The red line shows the Tdust profile from Keto & Caselli (2010). The inset shows a zoom-in of the n(H2) = 104 cm−3 to n(H2) = 105 cm−3 density regime, excluding the result for the full distribution.

Fig. 5 Abundances of selected gas-phase species (top row) and grain-surface species (labeled with asterisks; bottom row) as functions of medium density in the L1544 model. The abundances are displayed at three distinct time steps, ranging from 104 yr (left) to 106 yr (right). The different line styles correspond to model M1 (solid lines), model M3 (dashed lines), and model M4 (dotted lines). The ice abundances represent sums over all populations in models M3 and M4.

4 Discussion

4.1 Effect of initial conditions

We already pointed out above that the initial chemical composition adopted in this paper is not appropriate in the center of the core where the medium density is high (n(H₂) > 10⁶ cm⁻³). Starting instead from an initially (mostly) molecular gas composition would influence the timescales related to the development of ice chemistry, both for the total abundances and for the relative abundances on the various populations of grains. In a realistic scenario, the core structure is expected to develop gradually from a cloud with a lower average density, and chemical evolution takes place throughout the evolution of the physical structure. To emulate the chemical evolution during the core formation process, we constructed a simple two-stage model, in which the chemistry is first let to evolve over some duration of time in a low-density environment roughly corresponding to molecular clouds, after which we extracted the abundances of all species and used them as initial abundances for the core model. A similar process has been adopted in other works in the literature, for example, to investigate the chemical properties of the starless core L1521E (Nagy et al. 2019). For the molecular cloud properties, we took the values at the edge of the L1544 core model: n(H₂) = 2.66 × 10² cm⁻³, visual extinction A_V = 1 mag, T_gas = 18.7 K, and . We considered an array of grain temperatures, that is, we employed our model M4 for this test. We ran the molecular cloud model for t = 10⁶ yr (stage 1), and then ran the core model for an additional t = 10⁶ yr (stage 2).

Figure 7 shows the results of these calculations for the standard model M4 and its variant where chemical evolution starts from the abundances obtained at the end of stage 1. Given enough time for chemicalevolution in the core phase (using the L1544 physical structure), the results in the two cases are virtually identical at low densities. Abundance differences of a factor of a few appear at densities above n(H₂) ~ 10⁵ cm⁻³. Using molecular initial abundances does not have a clearly predictable effect on the results; switching the initial abundances causes some species to increase in abundance (e.g., water and ammonia), while others decrease in abundance (e.g., CO, N₂). The results of the comparison are also affected by how long the gas is let to evolve in the molecular cloud. We have shown in Sipilä et al. (2019) that running such a low-density model only for t = 10⁵ yr does not have an appreciable effect on the chemistry in the core stage. The discussion is therefore affected by uncertainties involved in the timescales of various processes, such as the duration of core formation and the lifetime in the pre-stellar phase. Abundance variations of a factor of a few across (a part of) the core, while seemingly small, may have a significant effect on observable line emission depending on critical densities of the lines (Sipilä & Caselli 2018).

Another property that has a crucial influence on our results in the outer core is the strength of the ISRF. If the field strength is increased, in particular in the UV, there are implications not only to the absolute dust temperatures but also to the relative dust temperatures in the size bins. We investigated this issue by running another variant of the M4 model in which the core is subjected to the unattenuated ISRF (that is, ). Figure 8 shows the T_dust profiles obtained when omitting the visual extinction external to the core. In this case a turnover point appears at a density of a few times 10³ cm⁻³; at densities lower than this, the grain equilibrium temperatures decrease with grain size. This is because the smallest grains absorb UV wavelengths efficiently (we return to this point in Appendix C). For densities above a few times n(H₂) ~ 10⁴ cm⁻³, where UV photons do not penetrate, the smallest grains are the coolest and the dust temperatures in the models with or 1 mag are very close to each other.

The modification in T_dust naturally affects molecular abundances as well. To quantify this issue, we ran model M4, setting the dust temperature profiles to those shown in Fig. 8, but retaining the assumption in the chemical model. This is an obvious inconsistency between the radiative transfer and chemical models, but it allows us to quantify more properly the effect of modifying T_dust alone, as moving from to 0 mag would also translate to a significant increase in the efficiency of photoreactions, thereby affecting the chemical evolution in the gas phase and in the ice. Figure 9 shows the relative ice abundances in the models with or 0 mag, concentrating on the low-density outer part of the core. Once again the response of the ice to changes in the grain properties depends on the molecule. For CO in the model with , the relative ice abundance distribution in each size bin is straightforward and, notably, CO ice is most abundant in bin 1 even though that grain population has the highest equilibrium temperature at low medium densities. N₂ ice however shows again its great sensitivity to the physical conditions as far as its formation is concerned, while water ice does not respond strongly to the changes in T_dust.

We reiterate that in these tests we neglected the effect of the one magnitude decrease in visual extinction on photochemistry, and that we would expect the (ice) chemistry to change drastically if was also decreased in the chemical model. Furthermore, it must also be noted that Fig. 9 shows the relative ice abundances in each grain size bin and the total ice abundances are very low in these conditions. In addition, the high dust temperaturein the outskirts of the core means that thermal desorption of relatively weakly bound species (such as atomic H and N) is comparable to, or overpowers, CR desorption. Nevertheless, our results clearly show that the ice composition can be strongly affected by the introduction of a grain size distribution, further underlining the demand for a detailed treatment of chemical evolution throughout the formation of a core, as opposed to assuming a static physical structure. Even such a detailed model is, unfortunately, affected by uncertainties from many sources, such as the initial elemental abundances.

Fig. 6 Relative abundances in each grain size bin as functions of density in the L1544 model at t = 104 yr (top) and at t = 105 yr (bottom). The red lines correspond to model M3, while the blue lines correspond to model M4. The green lines show the relative abundance of the grains in each bin, shown only in the bottom left panel for clarity (but note that the y-axis scaling is different in each subplot). The effective grain radius increases with bin label.

Fig. 7 As shown in Fig. 5, but in this figure indicating results from model M4 (solid lines) and model M4 with low-density cloud initial conditions (dotted lines). The time runs from the beginning of stage 2 (see main text).

Fig. 8 As shown in Fig. 4, but in this figure the core is subjected to the unattenuated ISRF.

4.2 General trends in molecular abundances

In the above, we showed for N₂ and CO ice that the abundance distributions on different grain populations in model M3 are straightforward at late times: abundances decrease as the grain radius increases. This finding does not extend to all species. Fig. 10 shows the distributions of selected additional molecules in model M3. H₂ CO and CH₃OH all present ice abundance profiles in which the abundance decreases monotonously with increasing grain radius. However, HCN and HNC are both most abundant in bin 4 virtually throughout the time evolution, and in the case of HNC, the fifth size bin dominates over the smallest three over a long time period of time. c − C₃H₂ is a mixed case as it presents an enhanced abundance in bin 4 over that in bin 3, even though the abundance in the othersize bins decreases monotonously with increasing grain size. As explained above, these effects are tied to attributes such as binding energies and volatilities as well as the gas-phase/grain-surface origin of the species. No correlation to the chemical relationship between pairs of species exists. For example, H₂ CO and CH₃OH (both produced mainly by the hydrogenation of CO) behave similarly to CO, but CO₂ presents a distribution over the grain populations that is similar to that of H₂O. This is in qualitative agreement with the connection between H₂O ice and CO₂ ice that has been established observationally (see for example Boogert et al. 2015). However, we note that the fraction of CO₂ ice that our model predicts (with respect to the most abundant ice constituent H₂O) is very low, implying that a revision of CO₂ formation inour model is needed to draw firm conclusions on this molecule.

Our results imply interesting consequences for the initial chemical composition of forming protostellar systems. This statement can be divided into two categories. First, introducing a distribution of grain sizes influences both gas-phase and ice abundances in high-density regions deep inside pre-stellar cores (also affecting the strength of the gas-dust collisional coupling; Ivlev et al. 2019). This alone is a clear indication of the need to consider a size distribution as opposed to simply assuming that all grains are the same size. Second, the relative chemical composition of the ice on grains of different sizes is such that at high medium densities most of the molecular material is on the smallest grains. It is likely that the smallest grains are rapidly removed from the distribution owing to impacts with larger grains at high densities regulated by ambipolar diffusion (Silsbee et al. 2020) as the collapse progresses. The heat generated in such collision events may lead to the small grains shedding their mantles as they are heated to higher temperatures than the large grains. The efficiency of this process depends on the relative velocity of the colliding grains, and it may also be that other mechanisms in addition to ambipolar diffusion would be needed to enhance the desorption efficiency on small grains enough to impact gas-phase abundances. Regardless, the critical question is then related to the timescales of molecular freeze-out and core collapse. If the collapse proceeds faster than the required time for the re-adsorption of previously evaporated molecules onto the (large) grains, the total amount of ice going into the prototellar phase would be smaller than in the case in which no grain-grain collisions occur. Also, the ice content would be different than in the case of desorption that is not induced by collisions (or otherwise), because some species can be most abundant on large grains as opposed to small grains. This effect could conceivably account for the observed presence of gas-phase ammonia in the center of L1544 (Crapsi et al. 2007) and, at least partly, for the observed abundance differentiation in a variety of molecules around its center (Spezzano et al. 2017). We will investigate this important issue in a future work.

The introduction of the grain size distribution affects the ionization fraction of the cloud, because the high amount of small grains leads to an overabundance of negatively ionized grains as compared to models considering monodisperse grains. We find that the size distribution models predict, for medium densities above a few × 10⁴ cm⁻³ where grain processes become important, a factor of a few less free electrons than model M1 does, depending on the time. A more quantitative analysis is, however, challenging because of the way we have defined the size distribution: we fixed the total grain area to correspond to that in the monodisperse grains model, which means that the overall number density of grains is higher in the size distribution model.

Fig. 9 Relative abundances in each grain size bin as functions of density in the L1544 model at t = 104 yr (top) and at t = 105 yr (bottom) for model M4. The blue lines correspond to , while the redlines correspond to . The effective grain radius increases with bin label. The x-axis range differs from that in Fig. 6. The dotted black line indicates the approximate position of the turnover point of the dust temperatures (Fig. 8).

Fig. 10 Distributions of selected molecules in the ice (indicated on top of each panel) as functions of time, in model M3 at n(H2) = 105 cm−3 and T = 10 K. The red lines represent the abundances on each grain population, while the blue line represents the sum over the individual grain populations. The effective grain radius increases with bin label as indicated in the plot.

4.3 Comparison to previous works and restrictions of the present model

There are, to our knowledge, two previous works in the literature that have discussed the effect of a grain size distribution on molecular abundances in physical conditions similar to those in the present paper, including grain-size dependence for CR desorptionin their model (Zhao et al. 2018; Iqbal & Wakelam 2018). We briefly compare our present approach to these two papers. We note that other authors have studied the chemical effects arising from a grain size distribution without modifying the CR desorption efficiency as a function of grain radius (e.g., Acharyya et al. 2011; Pauly & Garrod 2016).

The CR desorption rate coefficients calculated in the present work are derived analogously to the approach of Zhao et al. (2018). We verified that, for the same minimum and maximum grain size limits and discretization of the distribution,we obtain the same duty cycles and transient maximum grain temperatures as Zhao et al. (2018). However, there is a significant difference in the chemical models as they did not consider grain-surface chemical reactions at all, but allowed molecules to adsorp onto and desorb from grains. Therefore our present results are not comparable to theirs, and we do not, for example, recover the result from Zhao et al. (2018) that the abundances of nitrogen carriers in the gas phase could be enhanced by the introduction of a grain size distribution. We did confirm through benchmarking calculations that if we disregard grain-surface reactions and use the same chemical network as Zhao et al. (2018), we obtain similar overall results; some differences related to input parameters remain.

Iqbal & Wakelam (2018) studied chemical abundances in physical conditions corresponding to cold cores by employing a gas-grain chemical model that includes a treatment of grain-size dependence in CR desorption rate coefficients. These authors did not however account for the changes in the grain cooling timescale induced by the variation in the transient maximum grain temperature, and consequently the duty cycles and hence the CR desorption rate coefficients that they derived are very different from those presented in this work (and in Zhao et al. 2018). Also, they did not fix the total grain surface area between the monodisperse and size distribution models as we do in this work, making our respective approaches incompatible. Nevertheless, we obtain broadly similar results in some respects as Iqbal & Wakelam (2018), for example, that gas-phase abundances in models with and without grain size distribution typically agree within an order of magnitude.

We adopted in this work the grain size distribution from Mathis et al. (1977), in which the grain abundance has a simple power-law dependence on the grain size. Switching to another size distribution model where the relative amounts of small and large grains are different from the present case would naturally affect the ice compositions on the different grain populations. However, as long as the total grain surface area was kept constant with respect to our fiducial monodisperse grains model, we would expect no significant effect on gas-phase abundances due to such a switch. We do not expect our conclusions to change in a fundamental way if the size-dependence of the distribution is varied, unless there is a strong overabundance of large grains due to grain growth.

Another property that affects our results, especially pertaining to the relative ice abundances, is the adopted optical properties of the grains. In this work we constructed custom grain models mimicking the (unprocessed) thin ices case of Ossenkopf & Henning (1994). The new models are assumed to apply at all densities in the model cloud; that is, we do not modify the grain optical properties as a function of chemical evolution. The differences in grain equilibrium temperature over the size distribution are the result of density-dependent attenuation of external radiation. Obviously, the assumption of thin ices on the grains, as far as optical properties are concerned, is inappropriate at low medium densities early into the cloud evolutionwhere no significant ice formation has yet taken place. Conversely, at high densities the ices quickly grow thicker than what is assumed in our opacity model. It is possible with our current model setup to modify the dust opacities time-dependently based on the ice thickness (separately in each grain size bin), simulating the real-time effect of evolving ice layers on the optical properties of the grains, as a function of position in the cloud. Such an undertaking is however beyond the scope of the present paper and will be presented in a future work.

5 Conclusions

We investigated the effect of a grain size distribution on molecular abundances in physical conditions corresponding to starless and pre-stellar cores, using a gas-grain chemical model in which the CR desorption efficiency depends on the size of the grain. Small grains are heated to high temperatures by impinging CRs, but they also cool fast. Large grains cool more slowly and are heated by only a few tens of K by CRs. Consequently, the CR desorption rate coefficients are higher on large grains than on small grains. We divided the MRN size distribution into five distinct grain size bins and explicitly tracked the ice abundances in each bin as a function of time. The size distribution can be implemented in a chemical model in a straightforward manner, and we provided a table showing the crucial parameters as a function of grain radius. We also considered a model in which the grain equilibrium temperature is allowed to vary as a function of grain radius. For this, we determined the grain temperatures with a radiative transfer code using dust models custom built for the present purposes.

We found that the grain-size-variable CR desorption efficiency has highly nontrivial consequences for the ice content. Some molecules are most abundant on the smallest grains (e.g., CO, water, and CO₂), while others can be most abundant on larger grains (e.g., HCN and HNC) – although not necessarily on the largest grains in the distribution. This effect is tied to the main formation mechanism of a given molecule. CO, for example, is mainly produced in the gas phase and its grain-surface abundance is a simple decreasing function of the grain size. N₂ however is formed at early times almost exclusively on grain surfaces where the competition between desorption and hydrogenation influences the N₂ formation efficiency as a function of grain size. As a consequence, the grain-surface N₂ abundance is at early times highest on the larger grains in the distribution, but is later concentrated on small grains because of the high CR desorption efficiency onlarge grains. Some other molecules, such as HCN, can be most abundant on large grains over a very long period of time (even > 10⁶ yr).

Allowing the equilibrium grain temperature to vary with grain radius introduces a temperature gradient of about one K between the smallest and largest grains, depending on the physical conditions. In regions where external UV radiation penetrates, the smallest grains are the warmest of the distribution, but the situation is reversed in the inner, well-shielded regions inside starless and pre-stellar cores. The temperature gradient has a strong effect on the relative abundances of the ices on the different grain populations at low medium densities. However, allowing or disallowing variations in the equilibrium grain temperature does not have a large effect on the total abundances of grain-surface molecules or on gas-phase abundances; the practical difference is only in the relative abundances on the grain populations.

In the present paper we only considered cases in which there is no evolution of the physical conditions from low density to high density through gravitational contraction. The strong influence of the variations in equilibrium grain temperatureon the ice composition at low medium density, where ices begin to form, implies that the ice content at high densitycould also be influenced as the parental cloud contracts and a core forms. Another fundamental finding in the present work is that the distribution of molecules on grains of various sizes is not uniform. We propose that grain-grain collisions in infalling gas (possibly boosted by ambipolar diffusion) that result in grain coagulation and the (partial) evaporation of the ice mantles on small grains, could contribute to the observed chemical differentiation toward the centers of starless and pre-stellar cores, for example, in L1544 (Spezzano et al. 2017), L183 (Swade 1989; Lattanzi et al. 2020), and TMC-1 (Pratap et al. 1997). These findings imply great consequences for the chemical composition in star-forming clouds, and therefore in a future work we will investigate the time evolution of ices in a collapsing core, and will also explore the evaporation of ice mantles on small grains due to grain-grain collisions.

Acknowledgements

The authors thank the referee, Dr. Tyler Pauly, for constructive comments on the manuscript that helped to improve the presentation. The support by the Max Planck Society is gratefully acknowledged. O.S. thanks Mika Juvela and Mika Saajasto for helpful discussions related to the grain models.

Appendix A Variations in the number of grain size bins

The choice of the number of discrete grain sizes naturally affects ice abundances. When the number of grain size bins is increased, the materialadsorbed from the gas phase is spread out over an increasing number of grains. However, as long as the limits of the size distribution (minimum and maximum grain radius) are kept fixed along with the other critical parameters such as the dust-to-gas mass ratio, the chosen number of size bins does not affect the total grain surface area, and hence should not have an effect on our results at least when the bin number is greater than a few.

We compare in Fig. A.1 the ice abundances in our fiducial model M3 (five grain size bins) against variants of model M3 where the grain size bin number is either three, eight, or ten. In the model with the fewest size bins, some very small deviations from the fiducial model are found at early times, but none are apparent after t ≳ 10⁴ yr. Increasing the number of size bins decreases the abundances on the individual grain populations, but the total abundance (sum over all grain populations) is unaffected by the choice of the number of grain sizes. We also recover similar chemical behavior for the ices as in the fiducial case, for example, that the N₂ ice abundance is at early times dominated by the larger grains in the distribution. The abundances of gas-phase molecules are similarly unaffected by the number of grain size bins. These results justify our choice of five size bins as being enough to obtain accurate results, while simultaneously limiting the required computational time, which increases nonlinearly with each added grain size bin.

Fig. A.1 Distributions of CO (top row), N2 (middle row), and H2O (bottom row) ice as functions of time at a fixed medium density of n(H2) = 106 cm−3. Left, middle, and right panels: results for model M3 using 3, 8, or 10 grain size bins, respectively. The blue solid lines indicate the total ice abundance, while the black crosses indicate the corresponding total abundance in our fiducial M3 model. The red, orange, and green lines represent the breakdown of the total ice abundance in the grain size bins as shown in the top right panel. This breakdown is omitted for the fiducial M3 model for clarity.

Appendix B Ice thickness as a function of grain size

In the main text, we emphasize several times that the ice abundances are the largest on the smallest grains, which is sensible given that they represent the largest portion of the total grain surface. Because the number of binding sites on the grain surface decreases with grain radius, we expect an accompanying decrease in the ice thickness. Figure B.1 shows the ice thickness in monolayers (ML) for models M1, M3, and M4, as a function of density in the L1544 core model. We highlight three features of the plot. First, the ice accumulates very slowly at low density; for example, it takes ~ 10⁶ yr to grow just two MLs of ice at n(H₂) = 10³ cm⁻³. Second, in the monodisperse grains model (M1), the ice thickness does not evolve much after t = 10⁵ yr, and peaks around n(H₂) = 10⁵ cm⁻³ because the CR desorption efficiency increases with density, being especially important at the highest densities. In models M3 and M4, a similar peak is evident for all but the smallest grains, for which the ice thickness keeps increasing with density owing to the inefficiency of CR desorption on small grains. For the other grain sizes, the peak shifts to lower densities as the grain size increases, in accordance with the increase in CR desorption efficiency. There is no significant difference between the ice thickness in models M3 and M4. Third, the ice thickness is a simple decreasing function of grain size at all densities. This means that the high variability of the relative ice content seen, for example, in Fig. 6, where for instance N₂ ice is at times more abundant on large grains than on small grains, is not due to the ices accumulating at varying rates.

Fig. B.1 Ice thickness in monolayers (MLs) as a function of density in L1544 at three different times. The green, red, and blue lines correspond to models M1, M3, and M4, respectively. The grain size increases with bin label.

Appendix C Grain-size-dependent optical depth

In Sect. 4.1 we tested the effect of extinction on the size-dependent equilibrium dust temperatures. Here, we complement that discussion by showing, in Fig. C.1, the optical depth per hydrogen column density for the six grain models constructed for this paper via SIGMA. The plot serves to explain the density-dependence in the equilibrium grain temperatures shown in Figs. 4 and 8 in the main text. In our fiducial model where , the UV part of the ISRF spectrum is already strongly attenuated at the core edge, and the grain temperatures are dictated by the optical properties of the grains at longer wavelengths. The optical depth in the near- to mid-infrared wavelength range is a simple decreasing function of grain size. The equilibrium grain temperatures are thus highest for the largest grains when there is appreciable extinction external to the core. However, when this extinction is removed from the model, the grain temperatures in the outer core are dominated by the optical properties in the UV part of the spectrum, and consequently the smallest grains are then the warmest.

Fig. C.1 Optical depth per hydrogen column density for the six grain models constructed for this paper, that is, the five models corresponding to the size bins and the one calculated over the entire size distribution. The grain size increases with bin label. Right-hand panel: a zoom-in of the short-wavelength part of the spectrum.

References

All Tables

All Figures

	Fig. 1 Abundances of selected gas-phase species, indicated in the figure, as a function of time for three medium densities. The solid lines correspond to model M1, while the dashed lines correspond to model M2.
In the text

	Fig. 2 As described Fig. 1, but in this figure the dashed lines correspond to model M3.
In the text

Fig. 3 Distributions of CO (top row), N2 (middle row), and H2O (bottom row) ice as functions of time at three different medium densities. The green solid lines correspond to model M1. The red andblue lines correspond to model M3; the red lines represent the abundances on each grain population, while the blue line represents the sum over the individual grain populations. The green and blue lines overlap almost perfectly for H2O ice. The effective grain radius increases with bin label as indicated in the top right panel.

In the text

Fig. 4 Equilibrium dust temperature profiles as a function of density, calculated for the five distinct grain populations (black lines; model M4) or over the entire distribution (green; models M1 and M3). The effective grain radius increases with bin label as indicated in the plot. The red line shows the Tdust profile from Keto & Caselli (2010). The inset shows a zoom-in of the n(H2) = 104 cm−3 to n(H2) = 105 cm−3 density regime, excluding the result for the full distribution.

In the text

Fig. 5 Abundances of selected gas-phase species (top row) and grain-surface species (labeled with asterisks; bottom row) as functions of medium density in the L1544 model. The abundances are displayed at three distinct time steps, ranging from 104 yr (left) to 106 yr (right). The different line styles correspond to model M1 (solid lines), model M3 (dashed lines), and model M4 (dotted lines). The ice abundances represent sums over all populations in models M3 and M4.

In the text

Fig. 6 Relative abundances in each grain size bin as functions of density in the L1544 model at t = 104 yr (top) and at t = 105 yr (bottom). The red lines correspond to model M3, while the blue lines correspond to model M4. The green lines show the relative abundance of the grains in each bin, shown only in the bottom left panel for clarity (but note that the y-axis scaling is different in each subplot). The effective grain radius increases with bin label.

In the text

	Fig. 7 As shown in Fig. 5, but in this figure indicating results from model M4 (solid lines) and model M4 with low-density cloud initial conditions (dotted lines). The time runs from the beginning of stage 2 (see main text).
In the text

	Fig. 8 As shown in Fig. 4, but in this figure the core is subjected to the unattenuated ISRF.
In the text

Fig. 9 Relative abundances in each grain size bin as functions of density in the L1544 model at t = 104 yr (top) and at t = 105 yr (bottom) for model M4. The blue lines correspond to , while the redlines correspond to . The effective grain radius increases with bin label. The x-axis range differs from that in Fig. 6. The dotted black line indicates the approximate position of the turnover point of the dust temperatures (Fig. 8).

In the text

	Fig. 10 Distributions of selected molecules in the ice (indicated on top of each panel) as functions of time, in model M3 at n(H2) = 105 cm−3 and T = 10 K. The red lines represent the abundances on each grain population, while the blue line represents the sum over the individual grain populations. The effective grain radius increases with bin label as indicated in the plot.
In the text

Fig. A.1 Distributions of CO (top row), N2 (middle row), and H2O (bottom row) ice as functions of time at a fixed medium density of n(H2) = 106 cm−3. Left, middle, and right panels: results for model M3 using 3, 8, or 10 grain size bins, respectively. The blue solid lines indicate the total ice abundance, while the black crosses indicate the corresponding total abundance in our fiducial M3 model. The red, orange, and green lines represent the breakdown of the total ice abundance in the grain size bins as shown in the top right panel. This breakdown is omitted for the fiducial M3 model for clarity.

In the text

	Fig. B.1 Ice thickness in monolayers (MLs) as a function of density in L1544 at three different times. The green, red, and blue lines correspond to models M1, M3, and M4, respectively. The grain size increases with bin label.
In the text

	Fig. C.1 Optical depth per hydrogen column density for the six grain models constructed for this paper, that is, the five models corresponding to the size bins and the one calculated over the entire size distribution. The grain size increases with bin label. Right-hand panel: a zoom-in of the short-wavelength part of the spectrum.
In the text