Effect of grain size distribution and size-dependent grain heating on molecular abundances in starless and pre-stellar cores

We present a new gas-grain chemical model to constrain the effect of grain size distribution on molecular abundances in starless and pre-stellar cores. We introduce grain-size dependence simultaneously for cosmic-ray (CR)-induced desorption efficiency and for grain equilibrium temperatures. We keep explicit track of ice abundances on a set of grain populations. We find that the size-dependent CR desorption efficiency affects ice abundances in a highly non-trivial way that depends on the molecule. Species that originate in the gas phase follow a simple pattern where the ice abundance is highest on the smallest grains (the most abundant in the distribution). Some molecules, such as HCN, are instead concentrated on large grains throughout the time evolution, while others (like $\rm N_2$) are initially concentrated on large grains, but at late times on small grains, due to grain-size-dependent competition between desorption and hydrogenation. Most of the water ice is on small grains at high medium density ($n({\rm H_2}) \gtrsim 10^6 \, \rm cm^{-3}$), where the water ice fraction, with respect to total water ice reservoir, can be as low as $\sim 10^{-3}$ on large (>0.1 $\mu$m) grains. Allowing the grain equilibrium temperature to vary with grain size induces strong variations in relative ice abundances in low-density conditions where the interstellar radiation field and in particular its ultraviolet component are not attenuated. Our study implies consequences not only for the initial formation of ices preceding the starless core stage, but also for the relative ice abundances on the grain populations going into the protostellar stage. In particular, if the smallest grains can lose their mantles due to grain-grain collisions as the core is collapsing, the ice composition in the beginning of the protostellar stage could be very different to that in the pre-collapse phase.


Introduction
Interstellar dust grains play a large role in the chemical and physical evolution of molecular clouds that host starless and prestellar 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), and they also participate in the shaping of the chemical inventory in the clouds as new molecules are formed on their surface (see for example 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 for 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 the interstellar matter (ISM), it almost ubiquitously assumed that the grains are monodisperse, that is, they all have the same radius, 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 it has negative consequences for the accuracy of the simulations. We highlight two issues here. First, grains are transiently heated to high temperatures by cos-mic 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 K (Herbst & Cuppen 2006). The transient heating allows efficient non-thermal 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 non-thermal desorption allows the surface material to interact with the gas phase, it is highly desirable to investigate the effect that a size distribution may have on the chemical evolution of a star-forming object as a whole.
Indeed 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 including grain-size-dependence in CR heating efficiency (Zhao et al. 2018;Iqbal & Wakelam 2018). Of these studies, the work of Iqbal & Wakelam (2018) is the only one that 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 A&A proofs: manuscript no. sizeDistribution 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 start to form on grains already at low medium densities where one 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 prestellar 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 in order 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 Section 4, where we also compare our work to the previous related studies available in the literature. Section 5 presents our conclusions. Appendices A to C present some additional results that complement the main text.

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 temperature are 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.

Monodisperse grains versus grain size distribution
Gas-grain models of interstellar chemistry nearly universally assume monodisperse spherical dust grains, with the grain radius usually taken as 0.1 µm. In such a case, the grain abundance with respect to hydrogen can be calculated from where n H is the total hydrogen number density (≈ 2n(H 2 ) 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 −3 , and R d = 0.01. If one considers instead 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 dn d /da = Cn H a −3.5 . Here C is a normalization constant whose value can be derived under the assumption of a constant dust-to-gas mass ratio: where a max and a min are respectively the maximum and minimum radii considered for the distribution. The MRN distribution is continuous, but it can be easily discretized for inclusion in a chemical model by dividing the size interval [a min , a max ] into a number of grain size bins, calculating then the effective grain number density and the effective grain radius in each size bin i. Finally, the effective grain area in each bin is given by Presently we are interested in a comparison between models where 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 (X d,m πa 2 m ) matches the total effective surface area of the size distribution (Σ i (n d σ/n H ) i eff ). In this paper we assume [a min , a max ] = 0.03 µm, 0.3 µm 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 (X d,m Σ i X i d ). 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 as 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. Hasegawa & Herbst (1993) derived an expression for the rate coefficient of a CR-induced desorption (hereafter simply "CR desorption") event, which for the desorption of a grain-surface species j takes the form

CR-induced desorption
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 −5 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 13 s based on the iron CR flux given by Léger et al. (1985). The heating event interval scales as a −2 (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 for the grain-size-dependent duty cycle: 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 −19 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. T max is derived from the grain heat capacity (Léger et al. 1985). Small grains will be transiently heated to temperatures in excess of 150 K (the upper limit in Léger et al. 1985); here we use 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, where we also show 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.

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: 1) The thermal diffusion rate on the surface scales as is the binding energy of species j, and 2) Barrier-mediated surface reactions scale analogously by the inverse of T dust . Therefore, allowing the equilibrium grain temperature to vary over the size distribution will affect the reactivity on the different grain populations. For this reason we track explicitly 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 prestellar cores (e.g., Sipilä et al. 2017), we have 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 −3 µm to 0.25 µm), and we are presently considering a setup with multiple grain size bins. We therefore generated new sets of absorption and scattering cross sections for the five size bins displayed in Table 1, and over a broad distribution from 0.03 to 0.3 µm, using the SIGMA code (Lefèvre et al. (subm.) 1 ; 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 −3 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) so 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.

Additional model details
The results presented in this paper pertain to four individual model cases. These are collected in Table 2. Model M1 represents the standard case of monodisperse grains, which is al- Grain size distribution with radii and abundances given in Table 1, but T max and duty cycle set to constant values for all grain sizes (see main text) M3 Grain size distribution with parameters given in Table 1  M4 As M3, but equilibrium T dust depends on grain size most universally adopted in gas-grain models of ISM chemistry.
In model M2, we consider instead 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 the scaling of the CR desorption rate coefficients with grain radius and temperature. Here, T max,m = 73.86 K and f (0.095 µm, T max,m ) = 1.16 × 10 −19 are respectively the temperature and duty cycle that one obtains 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 use 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 employ point-like physical models to demonstrate the effects on the chemistry due to different grain setups. For simplicity we adopt in these models a constant temperature T gas = T dust = 10 K (that is, model M4 is not considered), and a constant visual extinction A V = 10 mag. The density of the medium is varied between n(H 2 ) = 10 5 cm −3 and 10 7 cm −3 .
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 here for the purposes of studying the size distribution effects in the context of a real core that has a non-constant density and temperature structure. We derive timedependent chemical abundance profiles in L1544 by using dust temperature profiles designed specifically for the present work (see Sect. 3.2) in connection with the density and gas temper- Table 3. Initial abundances (with respect to n H ) used in the chemical modeling.
Species Abundance H 2 5.00 × 10 −1 He 9.00 × 10 −2 C + 7.30 × 10 −5 N 5.30 × 10 −5 O 1.76 × 10 −4 S + 8.00 × 10 −8 Si + 8.00 × 10 −9 Na + 2.00 × 10 −9 Mg + 7.00 × 10 −9 Fe + 3.00 × 10 −9 P + 2.00 × 10 −10 Cl + 1.00 × 10 −9 ature structure from Keto & Caselli (2010). We divide the soobtained physical source model into concentric shells and calculate 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 2 . An A V profile, important for the photoreactions in the chemical model, was constructed assuming an ex- ternal visual extinction of A ext V = 1 mag (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 -here 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.

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 two are equivalent except of 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 2 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 5 yr, compared to the 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 non-trivial 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 2 , and H 2 O ice as functions of time. Let us consider N 2 in the lowest-density model (n(H 2 ) = 10 5 cm −3 ) at t = 10 4 yr. At this time, gas-phase N 2 is formed mainly through the N + CN −→ C + N 2 reaction in both cases (M1 and M3). However, the contribution of CR-desorbed icy N 2 to the gas-phase N 2 abundance is only less than 1% in model M1, but about 10% in model M3. This accounts for the slightly higher gas-phase N 2 abundance in the latter model. N 2 formation 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 2 production through the N + N association reaction. However, for the largest grains, CR desorption of atomic N is efficient enough to start inhibiting N 2 formation. Therefore, there exists a "sweet spot" for grain-surface N 2 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 2 is in model M3 mainly on the grains in the fourth size bin (a eff = 0.146 µm) at t = 10 4 yr. 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 larger in the latter model at early times. As the medium density is increased, CR desorption events occur more often, and gas-phase N 2 formation is dominated by CR desorption in both models M1 and M3; the N 2 ice abundance is highest on the largest grains only for a period of a few hundred years at n(H 2 ) = 10 7 cm −3 .
At later times, between t ∼ 10 5 yr and a few × 10 6 yr in the model with n(H 2 ) = 10 5 cm −3 , model M3 predicts a lower gasphase N 2 abundance (Fig. 2), even though N 2 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 2 at late times is N 2 H + + e − −→ N 2 + H, and its main destruction pathway is H + 3 + N 2 −→ N 2 H + + H 2 , constituting a loop between N 2 and N 2 H + . This applies to both models M1 and M3, and in fact the abundance of H + 3 is, at late times, somewhat higher in model M3 than in M1. Inspection of Fig. 3 reveals that N 2 ice is abundant only on the smallest grains at t > 10 5 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 A&A proofs: manuscript no. sizeDistribution the gas phase, and the majority of CO ice is material deposited from the gas phase through adsorption (unlike N 2 ice which is formed in situ at early times). Adsorption rates are largest on small grains which have the greatest surface area, and hence the CO ice abundance simply decreases with increasing grain radius. After t ∼ 10 5 yr when CR desorption becomes important, larger abundance variations between grains of different radius develop due to the increase of the CR desorption efficiency toward larger grains (Table 1). CO is not affected by this effect as greatly as N 2 , owing to its higher binding energy. The total CO ice abundance profile displays a clear peak at n(H 2 ) = 10 5 cm −3 , while remaining flat throughout the time evolution at n(H 2 ) = 10 7 cm −3 . The efficiency of CO ice processing is tied to the abundance of gasphase atomic H, which is roughly inversely proportional to the medium density. At n(H 2 ) = 10 5 cm −3 , the total CO ice abundance decreases past ∼ 10 5 yr due to hydrogenation processes that convert CO into formaldehyde and further into methanol. The production of methane is also effective. At higher densities, less 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 2 ) = 10 5 cm −3 .
Much like that of CO, the distribution of H 2 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 2 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, on the volatility of the species, etc. (We return to this issue in Sect. 4.2.) Nevertheless it is clear already from the very small sample of species displayed 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 here we started from initially atomic gas (except for hydrogen and deuterium which are initially in H 2 and HD) in order 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 one expects the  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 T dust profile from Keto & Caselli (2010). The inset shows a zoom-in of the n(H 2 ) = 10 4 cm −3 to n(H 2 ) = 10 5 cm −3 density regime, excluding the result for the full distribution.
gas density to increase gradually as the pre-stellar core forms and ultimately collapses, which may mean that differences between models M1 and M3 will not be as great as those displayed in Fig. 2. 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 2 ) = 10 4 cm −3 and n(H 2 ) = 10 5 cm −3 , 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 5 cm −3 . The low end of 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 2 ) = 10 4 cm −3 to n(H 2 ) = 10 6 cm −3 density range. We note that the dust temperature profiles calculated here are clearly different from the one 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 in order to carry out chemical calculations using the L1544 density structure. The gas temperature affects 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 different approaches to the treatment of grain heating have on surface chemistry, and thus for models M1 to M4 to be maximally comparable, we adopt for simplicity 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 in Sipilä & 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 2 and CO abundances follow the trends displayed in Fig. 2: N 2 is strongly depleted at high densities (n(H 2 ) 10 5 cm −3 in models M3 and M4 from ∼ 10 5 yr onwards), 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 displayed in Fig. 3 is evident for N 2 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 K 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 where 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 2 ice abundance in the size distribution models. The gas-phase abundance of ammonia decreases with the more efficient destruction of N 2 ice in models M3 and M4 (Fig. 5). The response of the models to changes in T dust is highly non-trivial 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 2 , and H 2 O ice as functions of medium density at t = 10 4 yr and t = 10 5 yr in models M3 and M4. Here, 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 displays clearly 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: 1) For N 2 at Article number, page 7 of 16 A&A proofs: manuscript no. sizeDistribution early times we recover the same trend evident in Fig. 3 where the ice is most abundant in bin 4 (note however that this trend does not appear at n(H 2 ) 10 4 cm −3 ), and 2) 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 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 2 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 2 ) ∼ 10 5 cm −3 , where the temperature difference between the smallest and the largest grains is small.
A notable result displayed 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 underrepresented 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, as 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, where we also explore the role of the external radiation field on chemical evolution due to 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. istic scenario, the core structure is expected to develop gradually from a cloud with a lower average density, and chemical evolution will take 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, where 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 extract the abundances of all species and use 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 ) = 2.66 × 10 2 cm −3 , visual extinction A V = 1 mag, T gas = 18.7 K, and T dust = [13.3, 13.5, 13.7, 13.9, 13.9] K. 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 6 yr (stage 1), and then ran the core model for an additional t = 10 6 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 chemical evolution 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 2 ) ∼ 10 5 cm −3 . 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, ammonia), while others decrease in abundance (e.g., CO, N 2 ). 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 for t = 10 5 yr only does not have an appreciable effect on the chemistry in the core stage. The discussion is therefore affected by uncertainties involved in the time-scales 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 where the core is subjected to the unattenuated ISRF (that is, A ext V = 0 mag). 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 3 cm −3 ; 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 2 ) ∼ 10 4 cm −3 where UV photons do not penetrate, the smallest grains are the coolest, and the dust temperatures in the models with A ext V = 0 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 A ext V = 1 mag 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 A ext V = 1 mag to 0 mag would also translate to a significant increase in the efficiency of photoreactions, affecting the chemical evolution in the gas phase and in the ice. Figure 9 shows the relative ice abundances in the models with A ext V = 1 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 A ext V = 0 mag, the relative ice abundance distribution in each size bin is straightforward and, notably, CO ice is most abundant in bin 1 despite that grain population having the highest equilibrium temperature at low medium densities. N 2 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 A ext V was decreased also in the chemical model. Furthermore, it must also be noted that: 1) Fig. 9 shows the relative ice abundances in each grain size bin, and that the total ice abundances are very low in these conditions, and 2) The high dust temperature in 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.

General trends in molecular abundances
In the above, we showed for N 2 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 2 CO and CH 3 OH all present ice abundance profiles where the abundance decreases monotonously with increasing grain radius. However, HCN and HNC are both most abun- bin 1 bin 2 bin 3 bin 4 bin 5 Fig. 8. As Fig. 4, but the core is subjected to the unattenuated ISRF. dant 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 3 H 2 is a mixed case as it presents an enhanced abundance in bin 4 over that in bin 3, even though the abundance in the other size bins decreases monotonously with increasing grain size. As explained above, these effects are tied to attributes such as binding energies and volatilities, and also the gas-phase/grain-surface origin of the species. No correlation to the chemical relationship between pairs of species exists. For example, H 2 CO and CH 3 OH (both produced mainly by the hydrogenation of CO) behave similarly to CO, but CO 2 presents a distribution over the grain populations that is similar to that of H 2 O. This is in qualitative agreement with the connection between H 2 O ice and CO 2 ice that has been established observationally (see for example Boogert et al. 2015). However, we note that the fraction of CO 2 ice that our model predicts (with respect to the most abundant ice constituent H 2 O) is very low, implying that a revision of CO 2 formation in our model is needed in order 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 have 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., subm.) 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 time-scales of molecular freeze-out and core collapse. If the collapse proceeds faster than the required time for the readsorption 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 where no grain-grain collisions occur. Also, the ice content would be different than in the case of no collision (or otherwise)-induced desorption, owing to the fact that 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 4 cm −3 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.

Comparison to previous works & 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 desorption in their model (Zhao et al. 2018;Iqbal & Wakelam 2018). Here 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), and we have 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) did. 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 grainsize-dependence in CR desorption rate coefficients. They 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 Article number, page 11 of 16 A&A proofs: manuscript no. sizeDistribution from those presented here (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) did, 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 various 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 op-tical properties of the grains. Here 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 evolution where no significant ice formation has yet taken place. Conversely, at high densities the ices will 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.

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 cosmic ray (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 slower 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 kept explicit track of 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 displaying the crucial parameters as a function of grain radius. We also considered a model where 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 non-trivial consequences for the ice content. Some molecules are most abundant on the smallest grains (e.g., CO, water, CO 2 ), 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 2 however is formed at early times almost exclusively on grain surfaces where the competition between desorption and hydrogenation influences the N 2 formation efficiency as a function of grain size. As a consequence, the grain-surface N 2 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 on large grains. Some other molecules, such as HCN, can be most abundant on large grains over a very long period of time (even > 10 6 yr).
Allowing the equilibrium grain temperature to vary with grain radius introduces a temperature gradient of about one K between the smallest and the 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 where 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 temperature on the ice composition at low medium density, where ices begin to form, implies that the ice content at high density could be influenced as well 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 different size 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 so in a future work we will investigate the time-evolution of A&A proofs: manuscript no. sizeDistribution ices in a collapsing core, and will also explore the evaporation of ice mantles on small grains due to grain-grain collisions.
O. Sipilä, B. Zhao, and P. Caselli: Effect of grain size distribution on molecular abundances A&A proofs: manuscript no. sizeDistribution 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 6 yr to grow just two MLs of ice at n(H 2 ) = 10 3 cm −3 . Second, in the monodisperse grains model (M1), the ice thickness does not evolve much after t = 10 5 yr, and peaks around n(H 2 ) = 10 5 cm −3 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 example N 2 ice is at times more abundant on large grains than on small grains, is not due to the ices accumulating at varying rates.