TY - JOUR
AU - Xi,, Hong-Wei
AB - ABSTRACT The temperatures of dust grains play important roles in the chemical evolution of molecular clouds. Unlike large grains, the temperature fluctuations of small grains induced by photons may be significant. Therefore, if the grain size distribution is included in astrochemical models, the temperatures of small dust grains may not be assumed to be constant. We simulate a full gas-grain reaction network with a set of dust grain radii using the classical MRN grain size distribution and include the temperature fluctuations of small dust grains. Monte Carlo method is used to simulate the real-time dust grain’s temperature fluctuations which is caused by the external low-energy photons and the internal cosmic ray induced secondary photons. The increase of dust grains radii as ice mantles accumulate on grain surfaces is also included in our models. We found that surface CO2 abundances in models with grain size distribution and temperature fluctuations are more than one order of magnitude larger than those with single grain size. Small amounts of terrestrial complex organic molecules (COMs) can also be formed on small grains due to the temperature spikes induced by external low-energy photons. However, cosmic ray induced secondary photons overheat small grains so that surface CO sublime and less radicals are formed on grains surfaces, thus the production of surface CO2 and COMs decreases by about one order of magnitude. The overheating of small grains can be offset by grain growth so that the formation of surface CO2 and COMs become more efficient. astrochemistry, ISM: abundances, ISM: molecules 1 INTRODUCTION It has long been recognized that interstellar dust grains play important roles for the formation of various molecules in various astronomical sources (Gould & Salpeter 1963; Charnley, Tielens & Millar 1992; Garrod & Herbst 2006). For instance, the efficient conversion from atomic hydrogen to molecular hydrogen was believed due to the interaction of atomic hydrogen and dust grains (Gould & Salpeter 1963; Hollenbach & Salpeter 1971; Biham et al. 2001). A well accepted H2 formation mechanism is that the hydrogen atoms first absorb on the grain surfaces, and then recombine to form H2 by the Langmuir–Hinshelwood mechanism. Many other species, including complex organic molecules (COMs), are believed to be formed by the same mechanism on grain surfaces (Garrod & Herbst 2006; Herbst & van Dishoeck 2009). Although dust grains have a sizedistribution in real astronomical sources (Mathis, Rumpl & Nordsieck 1977; Draine & Li 2001), most astrochemical models adopt a constant dust grain radius, which is |$0.1\,\mu \mathrm{ m}$|⁠, for the purpose of simplicity. Pauly & Garrod (2016) explained in detail why this approximation generally works reasonably well, but they also argued that there are limitations introduced by this approximation. Recently, the grain size distribution has been introduced into astrochemical modelling (Acharyya, Hassel & Herbst 2011; Ge, He & Li 2016; Pauly & Garrod 2016). It was found that surface chemistry on grains with a size distribution differs from that on grains with single radius 0.1μm. On the one hand, with the deposition of the gas phase species, the ice mantle on dust grains grows, thus the grain radius increase with time, which will increase the effective granular surface area and the rate of depletion of molecules (Acharyya et al. 2011). On the other hand, smaller grains will have more surface area because of their larger population. Pauly & Garrod (2016) showed that this effect can lead to the majority of ice mantle species reside on the smallest grains in the core collapse phase. Perhaps the most significant consequence of dust grain size distribution is the temperature difference of dust grains. Because the rate coefficients of surface reactions are strongly dependent on the dust grain temperatures, the ice mantle compositions are dependent on grain sizes. For instance, smaller dust grains are hotter than larger grains therefore species such as CO2 are more likely formed on smaller grains (Pauly & Garrod 2016). Moreover, smaller grains have more significant temperature fluctuations, which is called stochastic grain heating or single photon heating (Draine & Li 2001). Due to the small size, small grains have a smaller time-averaged vibrational energy (Draine & Li 2001), so that the grains whose radii are as small as |$0.005\,\mu \mathrm{ m}$| can be heated up to a few tens of kelvins by starlight photons, thus, the grain temperatures fluctuate between 5 and 50 K at |$A_v=0\,\mathrm{ mag}$| (Cuppen, Morata & Herbst 2006). To the best of our knowledge, so far astrochemical models that include stochastic heating of dust grains are only limited to the simplest molecular hydrogen formation models (Cuppen et al. 2006; Bron, Le Bourlot & Le Petit 2014). Stochastic heating may give us more new insights into the chemical evolution of ice mantle on dust grains. The temperature of the hottest grains in quite dark clouds in the models by Pauly & Garrod (2016) is only 15.6 K, while the transient temperature spikes induced by starlight can be more than 20 K. One particular interesting question is whether COMs, which were recently detected in cold prestellar cores (Öberg et al. 2010; Bacmann et al. 2012), can be formed on dust grains in cold cores because of the temperature spikes. The radicals that recombine to form COMs will not be mobile on grain surfaces until the temperatures of dust grains are more than 20 K (Garrod & Herbst 2006) therefore transient temperature spikes that are more than 20 K may help to form COMs in cold cores. In this paper, for the first time, we simulate a full gas-grain chemical network with physical conditions pertaining to cold dark cloud and include more explicit grain temperature fluctuations in our astrochemical models. The paper is organized as the following. The radiation inside cold cores will be introduced in Section 2 while the dust grain-size distribution and grain growth will be explained in Section 3. Heating and cooling of the grain are introduced in Section 4. The chemical models and simulation methods will be introduced in Section 5. We will show our results in Section 6. The comparison with observations and other models will be discussed in Section 7. Finally, in Section 8, we summarize our work and highlight the main conclusions. 2 PHOTON FLUX IN COLD CORES Following Cuppen et al. (2006), there are two types of photons based on their wavelength. Low-energy photons are those whose wavelength are between 250 nm and 1 cm, while high-energy photons are those whose wavelength are 91.2–250nm. Previous studies show that the average temperatures of dust grains are mainly controlled by external radiation field (Evans et al. 2001), and the heating of dust grains by internal radiation field was ignored (Zucconi, Walmsley & Galli 2001). However, for a cold dark molecular cloud with extinction |$A_v=10\,\mathrm{ mag}$|⁠, external high-energy photons can hardly penetrate into the inner part of the cloud, so the cosmic ray induced secondary photons are the dominant high-energy photons inside cold cores (Gredel et al. 1989). On the other hand, high-energy photons carry more energy, so they can cause more significant fluctuations. Because surface chemical reactions are very sensitive to dust grain temperature, we do consider the effect of high-energy photons heating in models as explained in the chemical model subsection. Dust grains are assumed to be silicate, which is used to choose the absorption coefficient data for the dust grains as in Cuppen et al. (2006). We follow Cuppen et al. (2006) to calculate the low-energy photon flux. The high-energy photon flux in our models are explained in detail here. We only consider the cosmic ray induced secondary photons for high-energy photons because they are the dominant ones. The spectrum of the cosmic ray induced secondary photons is complicated and consists of many lines (Shen et al. 2004). So, Shen et al. (2004) smoothed the spectrum between 850 to 1750 Å for the purpose of simplicity. Moreover, high-energy photons can do more than heating dust grains. High-energy photons can photodissociate ice mantle species or desorb surface species into gas phase so that the energy carried by high-energy photons is stored as chemical energy. Assuming one high-energy photon can only either photodissociate or desorb one surface species (Chang & Herbst 2014), the rate of high-energy photon flux that is stored as chemical energy is ∑ratei + ∑ratej, where i is for all photodissociation reactions and j is for all photodesorption reactions. For a dust grain with radius r, the rate of cosmic ray induced photons bombardment is QabsG0F0|$\pi$|r2, where G0 = 10−4 is the scaling factor for the cosmic ray induced photons (Shen et al. 2004), F0 = 108 cm−2 s−1 is the standard interstellar radiation field and Qabs is the wavelength dependent absorption coefficient for different sizes of silicate dust grains (Draine & Lee 1984; Draine 1985; Laor & Draine 1993), which is available online.1 So, a bare grain with radius r is bombarded by |$Q_{\rm{abs}}G_0F_0\mathrm{\pi}r^2 \Delta t$| high energy photons within the time period Δt. The rate of heating a dust grain with radius r by the cosmic ray induced photons is |$R_{\rm{pho}} = Q_{\rm{abs}}G_0F_0 {\mathrm \pi} r^2 - \sum{\rm{rate}_{\it{i}}} - \sum{\rm{rate}_{\it{j}}} $|. Because the spectrum of cosmic ray induced photons is complicated and we are only interested in the ice mantle composition change due to the temperature spikes induced by high-energy photons, we use the median value between 850 and 1750 Å, 1300 Å as the wavelength of cosmic ray induced photons in our models. In order to evaluate the impacts by the approximation, we also assume that the wavelengths of cosmic ray induced photons are uniformly distributed between 850 and 1750 Å in one model. 3 GRAIN SIZE DISTRIBUTION AND GRAIN GROWTH Dust grains in dense clouds may coagulate (Chokshi et al. 1993; Ysard et al. 2016). So, the distribution of dust grain size in dense clouds may differ from the classical MRN grain size distribution (Mathis et al. 1977), which is representative of the diffuse interstellar medium. However, the MRN size distribution has a simple analytical formula. Moreover, the grain size distribution in dense clouds is much more poorly known than in diffuse clouds. So, the MRN grain size distribution was still used by Acharyya et al. (2011) and Pauly & Garrod (2016) to study the effects of grain size distribution on the chemical evolution of dense clouds because simulation results based on the MRN size distribution can help as a first step to elucidate these effects. Therefore, we also assume that the grain-size distribution follows the classical MRN size distribution (Mathis et al. 1977), which is dn∝r−3.5dr. The largest and smallest radius are 0.25 and 0.005 μm, respectively. Following Acharyya et al. (2011) and Pauly & Garrod (2016), we initially divided the size distribution into five logarithmically equally spaced bins across the range of cross-section area. Then we calculate the mean cross-section area for each bin. Finally, the representative grain radius for each bin is calculated based on the mean cross-section area for that bin (Pauly & Garrod 2016). The 0th bin is the bin with |$r_{\mathrm{ max}0} =0.25\,\mu \mathrm{ m}$|⁠. The number of grains simulated for each bin followsthe relative abundances of grains according to the MRN dust distribution. The reason will be discussed in Section 5. We set the number of the grains in the 0th bin to be 1, so the number of dust grains in the ith bin, Ni is calculated as, \begin{eqnarray*} N_i=\frac{\int _{r_{\mathrm{ min}i}}^{r_{\mathrm{ maxi}}} r^{-3.5}\,\mathrm{ d}r}{\int _{r_{\mathrm{ min}0}}^{r_{\mathrm{ max}0}} r^{-3.5}\,\mathrm{ d}r}. \end{eqnarray*} (1) The initially calculated representative radii for the 0th through 4th bins are |$0.157\,\mu \mathrm{ m}$|⁠, |$0.0718\,\mu \mathrm{ m}$|⁠, |$0.0328\,\mu \mathrm{ m}$|⁠, |$0.0150\,\mu \mathrm{ m}$|⁠, and |$0.00687\,\mu \mathrm{ m}$|⁠, respectively. The numbers of dust grains in the 0th through 4th bins are found to be 1, 7, 49, 353, and 2500, respectively. Totally there are almost 3000 dust grains. The computational cost to simulate the chemical evolution of a chemical system with almost 3000 dust grains and the associated gas phase species is very expensive because we have limited CPUs. However, the number of total dust grains can be reduced based on the previous studies in order to reduce computational cost. Garrod & Pauly (2011) found surface chemistry is robust as long as the grain temperatures are not above 12 K or below 8 K. We found that the temperatures of grains in the 0th, 1st, and 2nd bins rarely fall out of the 8–12K range, if these grains are only heated by external low-energy photons. On the other hand, taking into account of the internal high-energy photons heating, the temperatures of grains in the 2nd bin can exceed 12 K while the temperatures of grains in the 0th and 1st bins are still within the range 8–12K. So the temperature fluctuations of the dust grains in the 0th, 1st, and 2nd bins can be ignored if we only consider the ambient radiation heating, however, the temperature fluctuations of the grains in the 2nd bins are significant if the heating by cosmic ray induced photons are considered. The surface chemistry on smaller dust grains may be different than that on much larger dust grains because the numbers of reactive species may be much less than one on the smaller grains, thus, the fluctuations of surface species numbers on smaller grains may change surface chemistry (Biham et al. 2001). However, it was found that the finite size effect is not important for a simple molecular hydrogen formation reaction network on dust grains with radius more than 0.05 μm (Biham et al. 2001) and becomes less important as the reaction network becomes more complicated (Chang & Herbst 2014). So, we can conclude that surface chemistry on the grains in the 0th and 1st bins are almost identical. Therefore, we can merge the 0th and 1st bins into a new 0th bin for a valid approximation. The mean cross-section area for the new 0th bin is calculated in order to calculate the representative grain radius for this bin (Pauly & Garrod 2016). The number of grains in the new 0th bin is reset to be one and the numbers of grains in other bins are calculated by equation (1). Table 1 summarizes the radii of the representative grains and the number of grains in eachbins. Table 1. Calculated grain-size distribution. Bin 0 1 2 3 Calculated radius/μm 0.086 997 0.032 83 0.015 01 0.006 87 Number 1 7 44 310 Bin 0 1 2 3 Calculated radius/μm 0.086 997 0.032 83 0.015 01 0.006 87 Number 1 7 44 310 Open in new tab Table 1. Calculated grain-size distribution. Bin 0 1 2 3 Calculated radius/μm 0.086 997 0.032 83 0.015 01 0.006 87 Number 1 7 44 310 Bin 0 1 2 3 Calculated radius/μm 0.086 997 0.032 83 0.015 01 0.006 87 Number 1 7 44 310 Open in new tab Grain growth due to the increase of ice mantle is also included in our models. With the gas phase species deposition, the grain radius will gradually become larger. The temperature fluctuation of dust grains will become less as ice mantle gradually becomes thicker. The depth of a monolayer of a grain is assumed as |$d_{\mathrm{ ML}}=1/\sqrt{A_s}$| (Pauly & Garrod 2016), where As is the site density of the grain. For a site density of 1.5 × 1015 cm−2, dML = 2.582 × 10−8 cm. 4 HEATING AND COOLING OF THE GRAIN Dust grains can also be directly heated by the bombardment of cosmic rays (Kalvans 2016). However, because the heating by cosmic rays is less frequent than that by photons, we focus on the temperature fluctuations induced by photons only in this work. The grain temperature fluctuations are caused by the absorption of low- and high-energy photons, which are treated as discrete random events in our simulations. The heat capacity is necessary to calculate the temperatures of grains that absorb photons. Following Cuppen et al. (2006), in which all the interstellar dust grains are assumed to be spherical silicate grains, the size-dependent heat capacity to convert photon energies to temperature spikes is, \begin{eqnarray*} c(T)=61.38 \,r^3\, T^3, \end{eqnarray*} (2) where r is the radius of a dust grain in μm while T is the temperature of a dust grain. The heat capacity is in eV K−1. Only radiative cooling of the grain is considered in our models in this work. The radiative cooling is treated with the usual continuous cooling approximation. We follow Cuppen et al. (2006) to calculate the temperatures of dust grains by radiative cooling. The sublimation of CO can also take away energy (Schutte & Greenberg 1991), however, CO sublimation cooling is not included in our models. We discuss the significance of CO sublimation cooling in Section 8. 5 CHEMICAL MODELS AND METHODS 5.1 Chemical models The chemical reaction network used in this work is from Hincelin et al. (2011), which is a full gas-grain reaction network. We slightly modify the reaction network as the following. First, the accretion of molecular hydrogen is absent in our reaction network in order to reduce the computational cost. Indeed, recent study shows that the absence of gas phase H2 accretion can only introduce negligible errors (Chang, Lu & Quan 2017). Secondly, we use the competition mechanism for the hydrogenation reactions that convert CO to methanol (Chang, Cuppen & Herbst 2007; Garrod & Pauly 2011; Chang & Herbst 2012). The hydrogenation of CO to form methanol has been studied extensively in laboratory experiments and the reaction barriers for the reactions H + CO and H + H2CO have been fitted using microscopic Monte Carlo method, which naturally includes the reaction–diffusion mechanism (Fuchs et al. 2009). The reaction barriers for H + CO and H + H2CO in Fuchs et al. (2009) are used in this work. We use a two-phase (gas phase and grain surface) model in which no distinction is made between the active layers and ice mantles for simplicity. The density is fixed to be, nH = 2 × 104 cm−3, while the gas temperature is fixed to be 10 K. The extinction is A|$v$| = 10 mag. The initial abundances used in this work are from Semenov et al. (2010) and areshown in Table 2. We keep the total dust to gas mass ratio to be 0.01. In total, there are NH = 7.11 × 1011 H nuclei in all cells of gas. Each grain is put in a cell of gas. The number of H nuclei in each cell is proportional to the cross-section area of the grain in the cell. The reason will be clear in the method section. The temperature fluctuation of dust grains with different sizes can be calculated with the methods introduced in Cuppen et al. (2006). To save CPU time, we neglect the temperature fluctuation of grains in the 0th bin and keep their temperature to be 10 K. Moreover, the temperature fluctuation of grains in the 1st bin induced by the low-energy photons are also ignored. In order to take into account the low-energy photons heating, the time-averaged temperature of grains in the 1st bin is calculated by a single run of temporal grain temperatures within a period of 1 yr, which is found to be 10.31 K. The temperature fluctuations of grains in the 1st bin induced by high-energy photons are included in models. All models are simulated for a period of 2 × 105yr which is the earliest time when gas phase species by model results agree best with observations in cold cores. Table 2. Initial abundances. Species Abundancea H2 0.5 He 9.00(−2) C+ 1.20(−4) O 2.56(−4) N 7.60(−5) P+ 2.00(−10) Na+ 2.00(−9) Cl+ 1.00(−9) Si+ 8.00(−9) S+ 8.00(−8) Mg+ 7.00(−9) Fe+ 3.00(−9) Species Abundancea H2 0.5 He 9.00(−2) C+ 1.20(−4) O 2.56(−4) N 7.60(−5) P+ 2.00(−10) Na+ 2.00(−9) Cl+ 1.00(−9) Si+ 8.00(−9) S+ 8.00(−8) Mg+ 7.00(−9) Fe+ 3.00(−9) aAbundance with respect to H, a(b) = a × 10b. Open in new tab Table 2. Initial abundances. Species Abundancea H2 0.5 He 9.00(−2) C+ 1.20(−4) O 2.56(−4) N 7.60(−5) P+ 2.00(−10) Na+ 2.00(−9) Cl+ 1.00(−9) Si+ 8.00(−9) S+ 8.00(−8) Mg+ 7.00(−9) Fe+ 3.00(−9) Species Abundancea H2 0.5 He 9.00(−2) C+ 1.20(−4) O 2.56(−4) N 7.60(−5) P+ 2.00(−10) Na+ 2.00(−9) Cl+ 1.00(−9) Si+ 8.00(−9) S+ 8.00(−8) Mg+ 7.00(−9) Fe+ 3.00(−9) aAbundance with respect to H, a(b) = a × 10b. Open in new tab The five models simulated in this work are summarized in Table 3. Model M1 includes the heating induced by low-energy photons only, while models M2 and M3 include both high-energy photons (cosmic ray induced UV photons) grain heating and low-energy photons (infrared wavelength photons) grain heating. The spectrum of the high-energy photons are approximated by a single energy (1300 Å) in model M2, while the wavelength of high-energy photons are assumed to be uniformly distributed between 850 and 1750 Å in model M3. Grain growth is not included in models M1, M2, and M3. The effect of grain growth is studied in model M4 in which dust grains are heated by both the high- and low-energy photons and grow in size as gas phase species accrete on grain surfaces. We also simulate a reference model, M5 in which stochastic heating and the growth of dust grains are absent. The radius of dust grains used in model M5 is the standard one, 0.1μm. The dust to gas mass ratio is also kept to be 0.01 in model M5 so that there are 6.16 × 1011 H nuclei in the single cell of gas. Table 3. Chemical models. Grain-size distribution Heating source Grain growth M1 Yes Low-energy photons No M2 Yes Low-energy photons and No high-energy photons (1300 Å) M3 Yes Low-energy photons and No high-energy photons (850–1750 Å) M4 Yes Low-energy photons and Yes high-energy photons (1300 Å) M5 No None No Grain-size distribution Heating source Grain growth M1 Yes Low-energy photons No M2 Yes Low-energy photons and No high-energy photons (1300 Å) M3 Yes Low-energy photons and No high-energy photons (850–1750 Å) M4 Yes Low-energy photons and Yes high-energy photons (1300 Å) M5 No None No Open in new tab Table 3. Chemical models. Grain-size distribution Heating source Grain growth M1 Yes Low-energy photons No M2 Yes Low-energy photons and No high-energy photons (1300 Å) M3 Yes Low-energy photons and No high-energy photons (850–1750 Å) M4 Yes Low-energy photons and Yes high-energy photons (1300 Å) M5 No None No Grain-size distribution Heating source Grain growth M1 Yes Low-energy photons No M2 Yes Low-energy photons and No high-energy photons (1300 Å) M3 Yes Low-energy photons and No high-energy photons (850–1750 Å) M4 Yes Low-energy photons and Yes high-energy photons (1300 Å) M5 No None No Open in new tab We define three terms to be used in the results section. The total population of a species, X, in all cells is noted as N(X) while the total population of the species X in all cells containing grains with the same radii ri is noted as |$\rm{N(X)}_{\it r_i}$|⁠, where i is 0, 1, 2, and 3. The radii r0, r1, r2, and r3 are 0.086997 μm, 0.03283 μm, 0.01501 μm, and 0.00687 μm, respectively. If the species X is a surface species, |$\rm{N(X)}_{\it r_i}$| is also the total population of the species X on all dust grains with radii ri. The total fractional abundance of a species, X, is defined as N(X)/NH. The fractional abundance of a species, X, in a type of cell is defined as |$\rm{N(X)}_{\it r_i}$|/(N|$\mathrm{ _H})_{r_i}$|⁠, where (N|$\mathrm{ _H})_{r_i}$| is the total population of H nuclei in all cells containing grains with the same radii ri. The type of cell is represented by the size of dust grains in the cell. Finally, the fractional abundance of a surface species, X, on grains with radii ri is defined as |$\rm{N(X)}_{\it r_i}$|/NH. 5.2 Methods The Gillespie algorithm is an ideal method to study stochastic chemical kinetics (Gillespie 1976). However, we found that the next reaction method (Gibson & Bruck 2000; Chang & Herbst 2012), which is a modified Gillespie algorithm, is more suitable to study the surface chemistry on stochastically heated grains. Therefore, the next reaction method is used in this work. The detailed explanation of the next reaction method can be found in Gibson & Bruck (2000) and Chang & Herbst (2012). We only briefly explain the method in the following. The absolute time when the next reaction occurs is used in the next reaction method. For gas phase reactions or surface reactions on dust grains with constant temperatures, the ith reaction will occur after time interval, τi = −ln X/Ri where X is a random number uniformly distributed within 0 and 1 while Ri is the reaction rate of the ith reaction. So the absolute time when the next ith reaction will occur is ti = t0i + τi , where t0i is the current time. If the reaction rate of the ith reaction changes from Ri to |$R_i^{^{\prime }}$| at time |$t_{0i}^{^{\prime }}$|⁠, where |$t_{0i}\lt t_{0i}^{^{\prime }} < t_i$|⁠, the absolute next reaction time for the ith reaction is updated as |$t_i^{^{\prime }}=t_{0i}^{^{\prime }} + \frac{(t_i - t_{0i}^{^{\prime }})R_i}{R_i^{^{\prime }}}$|⁠. The temperatures of dust grains always gradually cool down except when the grains are bombarded by star photons. Therefore, for the jth surface reaction on stochastically heated grains, the absolute next reaction time tj is calculated by the following equation, \begin{eqnarray*} -\ln Y=\int _{t_{0j}}^{t_{j}} R_j[T_d(t)]\, \mathrm{ d}t, \end{eqnarray*} (3) where Y is another random number uniformly distributed within 0 and 1, t0j is the current time, Rj is the reaction rate of the jth reaction and Td are the time dependent temperatures of dust grains that cool down by radiative cooling. Similarly, when the value of Rj changes because dust grains are heated by photons or other reactions change the abundances of the reactants of the jth reaction at time |$t_{0j}^{^{\prime }}$|⁠, the new absolute next reaction time |$t_j^{^{\prime }}$| is calculated as, \begin{eqnarray*} -\ln Y=\int _{t_{0j}}^{t_{0j}^{^{\prime }}} R_j[T_d(t)] \,\mathrm{ d}t + \int _{t_{0j}^{^{\prime }}}^{t_{j}^{^{\prime }}} R_j^{^{\prime }}[T_d(t)]\, \mathrm{ d}t \end{eqnarray*} (4) where |$R_j^{^{\prime }}$| is the new reaction rate of the jth reaction. Following Cuppen et al. (2006), we calculate the next arrival time of low-energy photons. The next arrival time of cosmic ray induced photons can be calculated in a manner similar to equation (3), using the rate of cosmic ray induced photons heating, Rpho discussed before. Thus, the absolute time when each event occurs can be found. We sort out the minimum absolute time and then execute the event, which is a chemical reaction or photon heating. If it is a chemical reaction, we change the population of reactants and products of the reaction and then update all the rates of reactions that the reactants and products of the reaction participate. If it is a photon heating event, we increase the temperature of dust grains and update all the rates of surface reactions whose reaction rates are dependent on surface temperatures. The process is repeated before the gas-grain system reaches the specified final time. In order to accelerate simulations, parallel computation is used. We can use one CPU to simulate the chemical evolution of surface ice and gas in one cell. In total we need almost 400 CPUs. However, we found that the computational costs to simulate the chemical evolution of different cells of gas and surface ice are different. For instance, the CPU time required by cells which contain grains in the 2nd bin is much more than that by the smallest cells which contain grains in the 3rd bin. In order to save CPUs, we use the multithread programming to simulate multiple cells of gas and ice in one CPU. One CPU can handle five of the smallest cells or one any other larger cell. Thus, the number of CPUs used is reduced. Totally 108 CPUs are used to perform simulations in this work. We assume the gas phase species are well mixed, so the number density of the same gas phase species in different cells of gas should be the same. Thus, the number of a gas phase species in different cells should be proportional to the volume of the cells. However, as the chemical system gradually evolves, the number density of gas phase species in different cells may not be the same because the chemical evolution of different cells may be different. The number of any gas phase species in each cell of gas is proportional to the surface area of the dust grain in each cell in order to ensure that the accretion of gas phase species does not introduce any difference of number density of gas phase species in different cells. However, chemical reactions other than accretion may introduce differences in the gas phase species number density. For instance, because of the temperature spikes on the smallest dust grains, surface H atoms are more likely to desorb from the smallest grains than the largest grains. Therefore, the number density of gas phase H in the cell of gas that contains the smallest grains may be higher than that in the cell of gas that contains the largest grain. Thus, we must mix gas phase species in different cells so that the number densities of each gas phase species in all cells are the same. Suppose we mix gas phase species in N1 cells that contain the smallest grains and N2 cells that contain the largest grains. The number density of gas phase species i after mixing is, \begin{eqnarray*} \bar{n}(i)=\frac{N_1 n(i)_1 + N_2n(i)_2}{N_1V_1 + N_2V_2}=\frac{\frac{N_1}{N_2} n(i)_1 + n(i)_2}{\frac{N_1}{N_2}V_1 + V_2}, \end{eqnarray*} (5) where V1 and V2 are the volumes of the cells containing the smallest and largest grains, respectively, while n(i)1 and n(i)2 are the population of gas phase species i in the cells containing the smallest and largest grains, respectively. Equation (5) shows that the ratio N1/N2 cannot be arbitrarily chosen in models because |$\bar{n}(i)$| is dependent on N1/N2. The relative abundances of grains according to the dust size distribution should be used to determine N1/N2. In order to keep the gas phase species well mixed, the message passing interface (MPI) scheme is used to exchange information between CPUs. The critical information is the number of each species in each cell. We redistribute the number of gas phase species in cells so that the number density of any species in different cells remain the same in simulations. We call the time interval mixing period later. On the other hand, the more information exchange between CPUs, the more CPU time is consumed for the information exchange by MPI. Therefore, the mixing period should be as long as possible. Different mixing periods are used in simulations to test the convergence of simulation results. The Taurus High Performance Computing system of Xinjiang Astronomical Observatory is used for the simulations in this work. It takes about one week to simulate models M1, M2, and M3 while the simulation of model M4 takes about two weeks., 6 RESULTS 6.1 Grain temperature fluctuations and grain growth Fig. 1 shows temperature fluctuations for different sizes of grains as a function of time in model M3, which includes both high- and low-energy photons heating. The wavelengths of high-energy photons are uniformly distributed between 850 and 1750 Å. We can see that high-energy photons heating events are much rarer than low-energy photons heating. For grains with radius 0.086 997 μm, we ignore temperature fluctuations because surface chemistry is robust within the range of temperature fluctuations as explained before. Similarly, for grains with radius 0.032 83 μm, we only consider the temperature fluctuations induced by the high-energy photons which can heat dust grains up to 13 K. For grains with radius 0.015 01μm and 0.00687 μm, both the high- and low-energy photons grain heatings are included. Low-energy photons can heat grains with radii 0.006 87 μm and 0.015 01 μm to more than 20 and 14 K, respectively. The temperature spikes induced by high-energy photons are much higher than that by low-energy photons. Grains with radii 0.006 87 μm and 0.015 01 μm can be heated by high-energy photons up to temperatures more than 35 and 20 K, respectively. Figure 1. Open in new tabDownload slide The temperatures as a function of time for different sizes of grains in model M3. Figure 1. Open in new tabDownload slide The temperatures as a function of time for different sizes of grains in model M3. The growth of grain radii is shown in Fig. 2, which corresponds to model M4. The radii of grains with initial radii 0.015 01 μm and 0.006 87 μm can increase up to about 0.02 μm and 0.01 μm, respectively, at the end of simulation. Figure 2. Open in new tabDownload slide The growth of grains with initial radii 0.015 01 μm and 0.006 87 μm in model M4. Figure 2. Open in new tabDownload slide The growth of grains with initial radii 0.015 01 μm and 0.006 87 μm in model M4. Fig. 3 shows the temperature fluctuations of dust grains at about 1.9 × 105yr in model M4, which includes the grain growth. We can see that the maximum temperatures decrease by a few kelvin for grains with initial radii 0.015 01 μm and 0.006 87 μm due to the grain growth. Figure 3. Open in new tabDownload slide The temperatures as a function of time for different sizes of grains in model M4 at about 1.9 × 105 yr. Figure 3. Open in new tabDownload slide The temperatures as a function of time for different sizes of grains in model M4 at about 1.9 × 105 yr. 6.2 Influence of mixing periods First, we test the convergence of simulation results by different mixing periods. Initially, the number of H2 in a cell is proportional to the volume of the cell. Moreover, because H2 is overwhelmingly more abundant than any other species, the formation or destruction of H2 by other species do not change much of the H2 abundances in each cell. Thus, H nuclei in each cell is proportional to the volume of the cell. Therefore, the fractional abundance of a species in a cell is proportional to the number density of the species. Thus, the fractional abundances in a cell instead of number density of gas phase species are used to test the convergence of simulation results. In order to find out the mixing period that is as long as possible, model M1 is simulated with mixing periods 103, 102, and 5 yr, respectively. We also simulate model M1 with an infinitely large mixing period so that no information is exchanged between different CPUs. Fig. 4 shows the fractional abundances of selected gas phase species in different types of cells as a function of time for different mixing periods. The different types of cells are represented by the size of dust grains in the cells. The influence of different length of mixing period can be seen from Fig. 4. With a mixing period that is infinitely large, because species such as H or N can easily desorb from grain surface when the temperatures of the smallest grains are well above 10 K, the fractional abundances of gas phase N in the cells with the smallest grains is higher than that with the largest dust grains whose temperatures are kept to be 10 K. For finite values of mixing periods which are no more than 103 yr, the fractional abundances of N in different cells are very close to each other. The fractional abundances of NH3 in different cells are almost identical as the mixing periods take finite values no more than 103yr. The difference of H abundances in different cells is not obvious, even the mixing period is infinitely large. Our simulation results shows that results do not change as much long as the mixing periods are not more than 103yr, which indicates that 103yr is ‘short’ enough so that the simulation results converge. So we always use a mixing period of 103yr in simulations. Figure 4. Open in new tabDownload slide Influence of mixing period on the fractional abundance of gas phase species in different cells. Figure 4. Open in new tabDownload slide Influence of mixing period on the fractional abundance of gas phase species in different cells. 6.3 Ice composition and surface COMs Fig. 5 shows the temporal evolution of the fractional ice composition (showing only major constituents) for different sizes of dust grains in model M1. The fractional ice composition is the percentage of that species with respect to the total number of major species in the ice mantle. We use the letter J to designate surface species hereafter. The major surface species are JH2O, JCO, JCO2, JCH4, JNH3, JH2CO, and JCH3OH. The fractional ice mantle compositions on the two largest grains are similar in model M1. Water ice is the most abundant granular species while surface CO is the second most abundant surface species. There is little JCO2 produced on these two largest grains in model M1 because the temperatures of these grains are constant and less than 12 K, so that the most efficient reaction to form JCO2, JCO + JOH → JCO2 + JH, does not happen easily. As grains become smaller, temperature fluctuations induced by photons are more significant, thus, the temperature of the smallest and the second smallest may exceed 12 K. We can see that the fractional ice composition changes significantly. First, the fractions of JCO2 on the two smallest grains increase, and even exceed the fraction of JCO at the end of the evolution on the smallest grains. Secondly, increasing dust grain temperature can increase the mobility of JH, thus, increases the rates at which JH atoms encounter other species. However, the reactions JH + JCO and JH + JH2CO become more difficult to happen as dust temperatures increase because we use competition mechanism for these two reactions. We can see that the fraction of JCH3OH increases as dust grains become smaller. The fraction of JCO decreases on the two smallest grains because JCO molecules are easier to be converted to other species due to the temperature spikes. We found that few JCO molecules can sublime regardless of the size of the grain in model M1. Finally, the fractional ice compositions of JH2O, JCH4, and JNH3 do not change much as grains become smaller. Figure 5. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in model M1. Figure 5. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in model M1. The temporal evolution of the fractional ice composition for different sizes of grains reported above is one realization because we run model M1 only once. Different realizations of the time evolution of the fractional ice composition in the same model might be different. Moreover, the fluctuation of the fractional ice composition in model M2 may be larger than that in the model M1 because of the larger temperature fluctuations in the model M2. In order to check the convergence of simulation results, we run the model M2 five times with different random seeds. We found that other than species with low abundances at early time, the fractional ice composition for different sizes of grains remain almost the same in the five different realizations. Therefore, the convergence of simulation results is achieved if we simulate each model just once. So, model results reported in this work are all from one realization of the temporal evolution of species abundances in each model. Influence of grain heating by cosmic ray induced photons can be clearly seen in Fig. 6. Although high-energy photons can heat the second largest grains and the second smallest grains to more than 12 and 20 K, respectively, the temperature spikes do not alter much the fractional ice composition on those grains. The fraction of JCO2 on the second smallest grains only slightly increases. The fractional ice composition on the smallest dust grains are strongly affected by the high-energy photons heating. First, we can see that high-energy photons grain heating does not help to produce JCO2 on the smallest grains. The fraction of JCO2 in the model M2 on the smallest grains is less than that in the model M1. Secondly, we can also see that as the fractional ice compositions evolve over time, the fraction of JCO dramatically decreases in model M2 because JCO can sublimate due to the high-temperature spikes induced by high-energy photons in model M2. That also explains why the extra high-energy heating cannot help the production of JCO2 on the smallest dust grains in model M2. Due to the decrease of the fraction of JCO in the model M2, the fraction of JH2CO also decreases in the model M2. The fraction of JCH4 on the smallest dust grains in the model M2 also quickly decreases as the fractional ice composition evolve over time because JC atoms can diffuse more quickly to react with surface species other than JH due to the temperature spike induced by high-energy photons, thus the number of surface carbon atoms that are hydrogenated to form JCH4 become less on the smallest grains due to high-energy photons heating. Figure 6. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in model M2. Figure 6. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in model M2. Fig. 7 shows the fractional ice composition (showing only major constituents) as a function of time in model M3. By comparing Figs 6 and 7, we can estimate the impact of the approximation that all high-energy photons have the same wavelength (1300 Å). The fractional ice composition as a function of time in the models M2 and M3 are similar. The only difference is that the fraction of JCH4 on the smallest grains in the model M2 is around a factor of 4 larger than that in the model M3. Dust grains in the model M3 can be heated to higher temperatures by high-energy photons, but some of the temperature spikes induced by high-energy photons in the model M3 are also lower than that in the model M2. Thus, the models M3 and M2 produce similar results. Figure 7. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in the model M3. Figure 7. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in the model M3. The influence of the grain growth on the fractional ice composition (showing only major constituents) can be seen in Fig. 8. The fractional ice compositions on the smallest grains are strongly affected by the grain growth, while the fractional ice compositions on larger grains are not much affected. Grain growth can decrease the temperature spikes of the overheated smallest grains, thus, offset the effect of high-energy photons heating. We can see that the fraction of JCO gradually increases after 104 yr, Because the radii of the smallest grains become larger, thus the temperature spikes induced by high-energy photon heating become smaller. Similarly, the fractions of JCO2 and JH2CO increase after 104 yr of evolution. Figure 8. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in the model M4. Figure 8. Open in new tabDownload slide The fractional ice composition (showing only major constituents) on different sizes of dust grains as a function of time in the model M4. Small amounts of terrestrial COMs can also be formed on dust grain surfaces in our models from M1 to M4. Fig. 9 shows the fractional abundances of JHCOOCH3 and JCH3OCH3 on grains with radius r as a function of time in those models. In order to form COMs on grain surfaces, the temperature spike of dust grains must be high enough so that radicals can diffuse. On the other hand, radicals are not easy to be formed on overheated grain surfaces, so temperature must be low enough so that radicals can be formed on grain surfaces. In the model M1, both JHCOOCH3 and JCH3OCH3 can only be formed on the smallest grains because radicals that recombine to form JHCOOCH3 and JCH3OCH3 cannot diffuse on larger grains. The smallest grains are not overheated in the model M1 so that sufficient surface radicals which recombine to form COMs are formed. The effect of high-energy photons heating on surface COMs formation is two-fold. First, in the model M2, small amounts of JHCOOCH3 and JCH3OCH3 molecules can also be formed on the second smallest grains due to the extra high-energy photons heating of grains. Secondly, the high-energy photons heating actually decreases the formation efficiency of JHCOOCH3 and JCH3OCH3 on the smallest grains because less JCH3O radicals are formed on the smallest grains in the model M2. The model M3 adopts a uniform distribution of the wavelength of high-energy photons, however, COM abundances in the model M3 are not much different than that in the model M2 in which all high-energy photons have the same energy. As the radii of grains increases, the temperature spikes become less. So, in the model M4, we can see that more JHCOOCH3 and JCH3OCH3 molecules are formed on the smallest grains, however, the formation of JHCOOCH3 and JCH3OCH3 on the second smallest grains becomes less efficient. Figure 9. Open in new tabDownload slide The fractional abundance for COMs on different sizes of dust grains as a function of time for different models. Figure 9. Open in new tabDownload slide The fractional abundance for COMs on different sizes of dust grains as a function of time for different models. The total fractional abundances of the major surface species and COMs in all models as a function of time are shown in Fig. 10. Since grain size distribution is not considered in model M5, the grain surface area in the model M5 is less than that in other models. Thus, less gas phase species accrete on dust grains in the model M5 than that in the other models. On the other hand, the total fractional abundances of surface species are also strongly affected by the temperature spikes. We can see that the total fractional abundances of JH2O and JCH4 are not much affected after a distribution of grain size, and stochastic heating of small grain are included in models. However, the total fractional abundances of JCO2 in the models M1–4 are more than one order of magnitude larger than that in the model M5 at 2 × 105 yr. Moreover, because the smallest grains are overheated in the models M2 and M3, overall JCO2 abundances in the models M1 and M4 are higher than that in the models M2 and M3. The total fractional abundances of JCH3OH in the models M1–4 are higher than that in the model M5 because the temperature spikes help to form methanol as explained before and the smaller grain surface area in the model M5. Moreover, the total fractional abundances of JCH3OH in the models M2 and M3 is lower than that in the models M1 and M4 because the smallest grains are overheated in the models M2 and M3 so that the formation of surface methanol in the models M2 and M3 is less efficient than that in the models M1 and M4. The terrestrial COMs, JHCOOCH3, and JCH3OCH3, can be formed in models from M1 to M4, but not in the M5 in which the dust grain temperatures are kept to be 10 K because its formation requires the surface diffusion ofradicals. Figure 10. Open in new tabDownload slide The total fractional abundances of major ice mantle species and surface COMs as a function of time. Figure 10. Open in new tabDownload slide The total fractional abundances of major ice mantle species and surface COMs as a function of time. 7 COMPARISON WITH OBSERVATIONS AND OTHER MODELS In this section, we compare our model results with the observational ice compositions towards background star Elias 16 and the model results from the previous studies. The selected models are MRN2 model in Acharyya et al. (2011) and 5G_T10_DIST in Pauly & Garrod (2016) because these two models adopt the same grain size distribution as that used in this work. The comparison is shown in the Table 4. Water ice abundances in all models are in good agreement with the observed values. The MRN2 model only includes a grain size distribution and grain growth while the temperatures of grains are fixed to be 10 K. Therefore, JCO2 abundances in the model MRN2 are more than one order of magnitude less than the observed values. The model 5G_T10_DIST includes the temperature distribution of dust grains so that the temperatures of small dust grains are more than 12 K, so JCO2 abundances in the 5G_T10_DIST is slightly more than the observed abundance. Our models M1, M2, and M3 take into account of the temperature fluctuations so that the temperature of small dust grains may exceed 12 K therefore significant amounts of JCO2 can also be produced. The abundances of JCO2 in the models M1 and M4 are close to the observed value while JCO2 abundances in the models M2 and M3 are about a factor of 4 smaller than the observed abundance. We can also see that abundances of JCO in our models M1–4 are larger than any other previous models and are slightly more abundant than the observed value. The observed abundances of JNH3 and JCH3OH towards background star Elias 16 have upper limits only while most of our model results are below this upper limit except that JCH3OH in the models M1 and M4 is slightly more abundant than the observed value. The abundances of JCH3OH in our models is about the same as that in the model 5G_T10_DIST, but less than that in the model MRN2. Our models from M1 to M4 produce less JNH3 than previous models do. Surface JCH4 and JH2CO have not been detected toward Elias 16. We can see that the abundance of JCH4 in the models M2 and M3 is the lowest among all models while the abundances of JH2CO in all our models from M1 to M4 are more than one order of magnitude higher than previous model results. Table 4. Comparison with observed ice compositions and previous model results. Species Elias 16 MRN2 MRN2 5G_T10_DIST M1 M2 M3 M4 time (yrs) 105 106 5 × 106 2 × 105 2 × 105 2 × 105 2 × 105 JH2O 6.4(−5) 1.85(−5) 7.0(−5) 1.59(−4) 7.72(−5) 9.43(−5) 9.40(−5) 8.24(−5) JCO 26 29 2(−3) 7.9 34 39 38 47 JCO2 20 1.4 0.27 36.2 24 5 6 15 JCH4 - 14 28 23.6 2 1 1 2 JCH3OH <3 4.4 16 2.9 5 1 1 4 JH2CO - 1.5(−4) 8.2(−8) 0.42 23 15 15 20 JNH3 <9 11 8.1 20.5 6 6 6 6 Species Elias 16 MRN2 MRN2 5G_T10_DIST M1 M2 M3 M4 time (yrs) 105 106 5 × 106 2 × 105 2 × 105 2 × 105 2 × 105 JH2O 6.4(−5) 1.85(−5) 7.0(−5) 1.59(−4) 7.72(−5) 9.43(−5) 9.40(−5) 8.24(−5) JCO 26 29 2(−3) 7.9 34 39 38 47 JCO2 20 1.4 0.27 36.2 24 5 6 15 JCH4 - 14 28 23.6 2 1 1 2 JCH3OH <3 4.4 16 2.9 5 1 1 4 JH2CO - 1.5(−4) 8.2(−8) 0.42 23 15 15 20 JNH3 <9 11 8.1 20.5 6 6 6 6 a(−b) means a × 10−b. Observational data are collated by Garrod, Wakelam & Herbst (2007). JH2O abundances are respect to H nuclei abundances, while the abundances of other species are given in the percentage of water ice abundances. The original JH2O abundances in the model MRN2 were respect to H2 and we have converted the abundances to that respect to H nuclei. Open in new tab Table 4. Comparison with observed ice compositions and previous model results. Species Elias 16 MRN2 MRN2 5G_T10_DIST M1 M2 M3 M4 time (yrs) 105 106 5 × 106 2 × 105 2 × 105 2 × 105 2 × 105 JH2O 6.4(−5) 1.85(−5) 7.0(−5) 1.59(−4) 7.72(−5) 9.43(−5) 9.40(−5) 8.24(−5) JCO 26 29 2(−3) 7.9 34 39 38 47 JCO2 20 1.4 0.27 36.2 24 5 6 15 JCH4 - 14 28 23.6 2 1 1 2 JCH3OH <3 4.4 16 2.9 5 1 1 4 JH2CO - 1.5(−4) 8.2(−8) 0.42 23 15 15 20 JNH3 <9 11 8.1 20.5 6 6 6 6 Species Elias 16 MRN2 MRN2 5G_T10_DIST M1 M2 M3 M4 time (yrs) 105 106 5 × 106 2 × 105 2 × 105 2 × 105 2 × 105 JH2O 6.4(−5) 1.85(−5) 7.0(−5) 1.59(−4) 7.72(−5) 9.43(−5) 9.40(−5) 8.24(−5) JCO 26 29 2(−3) 7.9 34 39 38 47 JCO2 20 1.4 0.27 36.2 24 5 6 15 JCH4 - 14 28 23.6 2 1 1 2 JCH3OH <3 4.4 16 2.9 5 1 1 4 JH2CO - 1.5(−4) 8.2(−8) 0.42 23 15 15 20 JNH3 <9 11 8.1 20.5 6 6 6 6 a(−b) means a × 10−b. Observational data are collated by Garrod, Wakelam & Herbst (2007). JH2O abundances are respect to H nuclei abundances, while the abundances of other species are given in the percentage of water ice abundances. The original JH2O abundances in the model MRN2 were respect to H2 and we have converted the abundances to that respect to H nuclei. Open in new tab Small amounts of COMs can be formed on the smaller grains in our models M1–4. The total fractional abundance of JHCOOCH3 and JCH3OCH3 can be as high as a few 10−8 and 10−9, respectively, in the model M4, which includes both stochastic heating of dust grains and grain growth. The formation of COMs on grains was not discussed in the previous models (Acharyya et al. 2011; Pauly & Garrod 2016). However, we do not expect COMs can be formed in the previous models because the temperatures of even the smallest grains are below 20 K, so radicals can hardly diffuse to formCOMs. 8 DISCUSSION AND SUMMARY We use the macroscopic Monte Carlo method to simulate gas-grain reaction networks under physical conditions pertaining to cold dark clouds. A distribution of grain sizes and stochastic heating of the smaller dust grains are included in the simulations. Five models are studied in this work. We simulate three models M1, M2, and M3 that do not include grain growth. Model M1 only considers dust grain heating by ambient low-energy photons while M2 and M3 also include the dust grain heating by cosmic ray induced UV photons inside molecular cloud. All cosmic ray induced UV photons have the same energy in the model M2 while the wavelength of the cosmic ray induced UV photons follows a uniform distribution in the model M3. Model M4 is simulated in order to find the effect of grain growth on the chemical evolution of ice mantles on stochastically heated grains. Finally, a reference model M5, in which there is no grain size distribution or stochastic heating of grains, is simulated for comparison. The fluctuations of the grain temperature dramatically alter the ice mantle compositions on the smaller dust grains. The abundances of JCO, JCO2, and COMs are more dependent on the temperature fluctuation of dust grains than any other species. The abundances of JCO2 increase significantly when stochastic heating is included in the models, which agree with previous studies which adopted size-dependent grain temperatures (Pauly & Garrod 2016). Moreover, our simulation results show that models that include low-energy photons only can already produce enough JCO2 that is consistent with the observational data toward Elias 16. Small amounts of terrestrial COMs can also be produced on the smallest grains due to the temperature spikes induced by low-energy photons. The influence of high-energy photons heating in models is two-fold. First, the smallest grains are overheated so that less JCO2 or COMs are formed. Moreover, JCO can easily sublime on the smallest grains. Secondly, the temperature spikes on the second largest grains are high enough so that more JCO2 and COMs are formed. The complicated spectrum of high-energy photons is approximated as either a single photon energy or a uniform distribution in our models. The choice of approximation does not have a large impact on the simulation results. Grain growth, which is included in model M4, is able to decrease the temperature spikes of the overheated dust grains as ice mantles gradually accumulate on dust grains. Thus, more COMs and JCO2 are formed on the smallest grains which are overheated by high-energy photons. Schutte & Greenberg (1991) argued that JCO sublimation cooling may be more important than radiative cooling if the temperatures of dust grains are above 26K. Their conclusion is based on the assumption that the whole grain surface is covered by volatile species (JCO). The overheated grains should cool down much more quickly if JCO sublimation cooling dominates over radiative cooling. However, we argue that JCO sublimation cooling is not likely to be important in our models because of the following reasons. The fraction of JCO on the smallest grains is less than 10−3 at the time 104 yr in model M2 as shown in Fig. 6. The total population of surface species on the smallest grains at the time 104 yr is a few monolayers. Therefore, the fraction of the surface covered by JCO on the smallest grains at the time 104 yr is less than 1 per cent in model M2. Similarly, we can estimate that the fraction is also less than 1 per cent at other time steps. Since the number of JCO desorbed from grain surfaces within fixed time interval is proportional to the population of JCO on grain surfaces, following Schutte & Greenberg (1991), we can estimate that the energy taken away by JCO sublimation on the smallest grains is less that 1 per cent of that by radiative cooling at 26 K. So radiative cooling still dominates over JCO sublimation cooling at 26 K because of the low coverage of JCO. Moreover, we can also estimate how much JCO sublimation cooling can decrease the temperature spike of the smallest bare grains in model M2. Following Schutte & Greenberg (1991), we assume the energy taken away by the sublimation of each JCO molecule is 960 K. Solving equation (2), the temperature of the smallest bare dust grain drops by ΔT ∼ 0.24 K for each JCO sublimation event at T = 26K. Equation (2) shows that the heat capacity of grains increases as T increases, thus ΔT should decrease at higher temperatures. On the other hand, we can estimate that each high-energy photon can only desorb a few JCO molecules from each of the smallest grains. We found that less than 80 000 JCO molecules desorb from each of the smallest dust grains during 2 × 105yr in model M2 while each of the smallest grains are bombarded by around 30 000 high-energy photons during the same time period. Since JCO can hardly sublime due to the temperature spike induced by low-energy photons, we can approximately estimate that each high-energy photon can only desorb less than three JCO molecules from each of the smallest grains. So JCO sublimation cooling can only lead to a temperature drop of < 0.72 K. Therefore, we can conclude that JCO sublimation cooling may not have a large impact in model M2. Moreover, the JCO sublimation cooling should be even less significant in the models that include grain growth because less JCO molecules sublime in these models. It is particularly interesting that COMs can be formed on the smallest grains when stochastic heating is included in models. The amount of COMs formed in models is small however. The fractional abundances of JHCOOCH3 and JCH3OCH3 are a few 10−8 and 10−9, respectively. The reason is that the abundances of radicals which can recombine are small on the smallest grains. Recent study shows that radicals can accumulate in ice mantles if we adopt a three phase model with photon penetration and bulk diffusion (Chang & Herbst 2014). Further research is necessary to study how COMs can be formed inside ice mantle via bulk diffusionmechanism. Traditionally, dust grains are assumed to be uniform in size with radii 0.1 μm, so the temperature fluctuations of dust grains are ignored in most astrochemical models. Our simulation results show that the temperature fluctuations of small grains whose radii are about 0.006 μm are large enough to alter the compositions of ice mantle. On the other hand, because of dust coagulation, the population of grains which are small enough to undergo significant temperature fluctuations may be less than that predicted by the MRN grain size distribution. Therefore, more study should be done to investigate the roles of small grain for the evolution of molecular clouds. We summarize our main results as the following: The temperatures of small grains inside cold molecular clouds are strongly affected by the ambient interstellar radiation and cosmic ray induced secondary photons. Considering the ambient interstellar radiation only, the temperature spikes are less than 16 K for the grains with radius 0.015 01 μm, while the temperature spikes can be more than 20 K for the smaller grains with radius 0.006 87 μm. The cosmic ray induced secondary photons can further increase the grain temperature fluctuation. The highest temperatures for grain with radii 0.015 01 μm and 0.006 87 μm can be more than 20 and 35 K, respectively. Dust grain temperature fluctuations can increase the production of JCO2. Overall, the formation efficiency of JCO2 on grain surfaces decreases if the cosmic ray induced secondary photons heating is included in models. The abundances of JCO on the smallest grains that are heated by cosmic ray induced secondary photons are much smaller than that on other sizes of grains because JCO can sublimate on these grains. Small amounts of terrestrial COMs are able to be formed on stochastically heated small grains because radicals that recombine to form COMs are mobile because of the temperature spike. Grain growth is able to decrease the temperature spikes of the overheated smallest dust grains, thus, the formation of JCO2 and COMs on the smallest grains becomes more efficient as the size of grains increases in models that include high-energy photons heating. The grain mantle chemistry is not much affected if a model adopts a uniform distribution of the wavelength of high-energy photons instead of the median wavelength. Footnotes 1 http://www.astro.princeton.edu/∼draine/dust/dust.diel.html ACKNOWLEDGEMENTS QC is a research fellow of the One-Hundred-Talent project of the Chinese Academy of Sciences. This work was funded by The National Natural Science foundation of China under grant 11673054. We thank Eric Herbst for critical reading of the manuscript and helpful discussions. The Taurus High Performance Computing system of Xinjiang Astronomical Observatory was used for the simulations. We thank our referee’s very constructive comments to improve the quality of the manuscript. REFERENCES Acharyya K. , Hassel G. E. , Herbst E. , 2011 , ApJ , 732 , 73 https://doi.org/10.1088/0004-637X/732/2/73 Crossref Search ADS Bacmann A. , Taquet V. , Faure A. , Kahane C. , Ceccarelli C. , 2012 , A&A , 541 , L12 Crossref Search ADS Biham O. , Furman I. , Pirronello V. , Vidali G. , 2001 , ApJ , 553 , 595 https://doi.org/10.1086/320975 Crossref Search ADS Bron E. , Le Bourlot J. , Le Petit F. , 2014 , A&A , 569 , A100 Crossref Search ADS Chang Q. , Cuppen H. M. , Herbst E. , 2007 , A&A , 469 , 973 Crossref Search ADS Chang Q. , Herbst E. , 2012 , ApJ , 759 , 147 https://doi.org/10.1088/0004-637X/759/2/147 Crossref Search ADS Chang Q. , Herbst E. , 2014 , ApJ , 787 , 135 https://doi.org/10.1088/0004-637X/787/2/135 Crossref Search ADS Chang Q. , Lu Y. , Quan D. , 2017 , ApJ , 851 , 68 https://doi.org/10.3847/1538-4357/aa99d9 Crossref Search ADS Charnley S. B. , Tielens A. G. G. M. , Millar T. J. , 1992 , ApJ , 399 , 71 https://doi.org/10.1086/186609 Crossref Search ADS Chokshi A. , Tielens A. G. G. M. , Hollenbach , D , 1993 , ApJ , 407 , 806 https://doi.org/10.1086/172562 Crossref Search ADS Cuppen H. M. , Morata O. , Herbst E. , 2006 , MNRAS , 367 , 1757 https://doi.org/10.1111/j.1365-2966.2006.10079.x Crossref Search ADS Draine B. T. , 1985 , ApJS , 57 , 587 https://doi.org/10.1086/191016 Crossref Search ADS Draine B. T. , Lee H. M. , 1984 , ApJ , 285 , 89 https://doi.org/10.1086/162480 Crossref Search ADS Draine B. T. , Li A. , 2001 , ApJ , 551 , 807 https://doi.org/10.1086/320227 Crossref Search ADS Evans N. J. II , Rawlings J. M. C. , Shirley Y. L. , Mundy L. G. , 2001 , ApJ , 557 , 193 https://doi.org/10.1086/321639 Crossref Search ADS Fuchs G. W. , Cuppen H. M. , Ioppolo S. , Romanzin C. , Bisschop S. E. , Andersson S. , van Dishoeck E. F. , Linnartz H. , 2009 , A&A , 505 , 629 Crossref Search ADS Garrod R. , Herbst E. , 2006 , A&A , 457 , 927 Crossref Search ADS Garrod R. T. , Pauly T. , 2011 , ApJ , 735 , 15 https://doi.org/10.1088/0004-637X/735/1/15 Crossref Search ADS Garrod R. T. , Wakelam V. , Herbst E. , 2007 , A&A , 467 , 1103 Crossref Search ADS Ge J. X. , He J. H. , Li A. , 2016 , MNRAS , 460 , L50 https://doi.org/10.1093/mnrasl/slw058 Crossref Search ADS Gibson M. A. , Bruck J. , 2000 , J. Phys. Chem. A , 104 , 1876 https://doi.org/10.1021/jp993732q Crossref Search ADS Gillespie D. T. , 1976 , J. Comput. Phys. , 22 , 403 https://doi.org/10.1016/0021-9991(76)90041-3 Crossref Search ADS Gould R. J. , Salpeter E. E. , 1963 , ApJ , 138 , 393 https://doi.org/10.1086/147654 Crossref Search ADS Gredel R. , Lepp S. , Dalgarno A. , Herbst E. , 1989 , ApJ , 347 , 289 https://doi.org/10.1086/168117 Crossref Search ADS Herbst E. , van Dishoeck E. F. , 2009 , ARA&A , 47 , 427 https://doi.org/10.1146/annurev-astro-082708-101654 Crossref Search ADS Hincelin U. , Wakelam V. , Hersant F. , Guilloteau S. , Loison J. C. , Honvault P. , Troe J. , 2011 , A&A , 530 , A61 Crossref Search ADS Hollenbach D. , Salpeter E. E. , 1971 , ApJ , 163 , 155 https://doi.org/10.1086/150754 Crossref Search ADS Kalvāns J. , 2016 , ApJS , 224 , 42 https://doi.org/10.3847/0067-0049/224/2/42 Crossref Search ADS Laor A. , Draine B. T. , 1993 , ApJ , 402 , 441 Crossref Search ADS Mathis J. S. , Rumpl W. , Nordsieck K. H. , 1977 , ApJ , 217 , 425 https://doi.org/10.1086/155591 Crossref Search ADS Öberg K. I. , Bottinelli S. , Jørgensen J. K. , van Dishoeck E. F. , 2010 , ApJ , 716 , 825 https://doi.org/10.1088/0004-637X/716/1/825 Crossref Search ADS Pauly T. , Garrod R. T. , 2016 , ApJ , 817 , 146 https://doi.org/10.3847/0004-637X/817/2/146 Crossref Search ADS Schutte W. A. , Greenberg J. M. , 1991 , A&A , 244 , 190 Semenov D. et al. , 2010 , A&A , 522 , A42 Crossref Search ADS Shen C. J. , Greenberg J. M. , Schutte W. A. , van Dishoeck E. F. , 2004 , A&A , 415 , 203 Crossref Search ADS Ysard N. , Köhler M. , Jones A. , Dartois E. , Godard M. , Gavilan L. , 2016 , A&A , 588 , A44 Crossref Search ADS Zucconi A. , Walmsley C. M. , Galli D. , 2001 , A&A , 376 , 650 Crossref Search ADS © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)
TI - Effect of stochastic grain heating on cold dense clouds chemistry
JO - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty1525
DA - 2018-09-21
UR - https://www.deepdyve.com/lp/oxford-university-press/effect-of-stochastic-grain-heating-on-cold-dense-clouds-chemistry-FPft3bk40s
SP - 2988
VL - 479
IS - 3
DP - DeepDyve
ER -