Efficiency of non-thermal desorptions in cold-core conditions. Testing the sputtering of grain mantles induced by cosmic rays

Under cold conditions in dense cores, gas-phase molecules and atoms are depleted from the gas-phase to the surface of interstellar grains. Considering the time scales and physical conditions within these cores, a portion of these molecules has to be brought back into the gas-phase to explain their observation by milimeter telescopes. We tested the respective efficiencies of the different mechanisms commonly included in the models. We also tested the addition of sputtering of ice grain mantles via a collision with cosmic rays in the electronic stopping power regime. The ice sputtering induced by cosmic rays has been added to the Nautilus gas-grain model while the other processes were already present. Each of these processes were tested on a 1D physical structure determined by observations in TMC1 cold cores. The resulting 1D chemical structure was also compared to methanol gas-phase abundances observed in these cores. We found that all species are not sensitive in the same way to the non-thermal desorption mechanisms, and the sensitivity also depends on the physical conditions. Thus, it is mandatory to include all of them. Chemical desorption seems to be essential in reproducing the observations for H densities smaller than $4\times 10^4$~cm$^{-3}$, whereas sputtering is essential above this density. The models are, however, systematically below the observed methanol abundances. A more efficient chemical desorption and a more efficient sputtering could better reproduce the observations. In conclusion, the sputtering of ices by cosmic-rays collisions may be the most efficient desorption mechanism at high density (a few $10^4$~cm$^{-3}$ under the conditions studied here) in cold cores, whereas chemical desorption is still required at smaller densities. Additional works are needed on both mechanisms to assess their efficiency with respect to the main ice composition.


Introduction
In recent years, our knowledge of the composition of the interstellar medium (ISM) and the process of star formation has entered a new area thanks to the powerful capabilities (in term of sensitivity, spatial resolutions, and new wavelength windows) of recent ground-based and space observatories (Spitzer, Herschel, ALMA, NOEMA, and the new EMIR receiver on the IRAM 30m). These new instruments have allowed for the detection of new molecules (in the gas-phase or in the ices), namely, either simple ones such as HCl + (De Luca et al. 2012) or complex ones such as (NH 2 ) 2 CO (Belloche et al. 2019). In addition to revealing the richness of the ISM chemistry, they have revealed the non-uniform distribution of their abundances from object to object or even within a certain type of object (with similar physical conditions). One of the key questions astrochemists want to answer is how these molecules in the gas-phase or in the ices form and whether we can reproduce their variability in abundance. The gas-phase composition of cold cores has been observed for a long time using ground-based telescopes in the mm and submm range. Spectral surveys have also revealed a diversity of relative abundances (Ohishi et al. 1992;Pratap et al. 1997;Dickens et al. 2000;Gratier et al. 2016). Current chemical models do a good job at reproducing most of the observed abundances (Wakelam et al. 2006), although there are still some molecules that have continue to present challenges, even though the chemistry has been considered thoroughly (e.g., CH 3 CCH and S 2 H, Hickson et al. 2016;Fuente et al. 2017). In recent years, however, the detection of complex organic molecules (such as HCOOCH 3 , CH 3 OCH 3 , and CH 3 CHO) in the cold environment of dark clouds has revived the interest for cold core chemistry (Bacmann et al. 2012;Vastel et al. 2014). Until this discovery, these molecules were believed to be formed only in the vicinity of hot protostars, where the high dust temperature would effectively promote their formation on the grains (Herbst & van Dishoeck 2009). Now both their formation at low temperature and the possibility to observe them in the cold gas-phase has again raised questions around their origin (Vasyunin & Herbst 2013;Balucani et al. 2015;Ruaud et al. 2015;Jiménez-Serra et al. 2016;Vasyunin et al. 2017) and determined their appearance at much earlier phases during the formation of stars. While there is a debate on the chemical path (gas-phase or surface) forming these molecules, it remains clear that some of them or their precursors need to desorb from the grains at low temperatures. Even for simple molecules commonly observed in the cold gas of dense cores, the physical conditions and the time scales produce a depletion of molecules onto the grains (see, e.g., Garrod et al. 2007). Astrochemical models simulating such environments need to invoke non-thermal desorption mechanisms to bring back into the gas-phase molecules that would have otherwise disappeared from the gas . Several mechanisms have been proposed and some of them have been studied in the laboratory, but their inclusion in astrochemical models is far from being perfect and their efficiency is far from being explicit. We cite a number of them here, but this is not an exhaustive list. The effect of cosmic rays on the ice chemistry has recently gained a lot of attention. For instance,  and  have proposed methods that include a cosmic-ray driven radiation surface chemistry into astrochemical models (see also Shingledecker et al. 2019Shingledecker et al. , 2020. Once a cosmic-ray particle collides with a grain, it produces radiolysis of species on the grains. Some products of these radiolysis can be supra-thermal and then very reactive, increasing the production of complex organic molecules. Such chemistry would certainly benefit from more experimental measurements. Kalvāns & Shmeld (2013) have developed a model in which the bulk of the ices (separated from the surface) is made of cavities in which the diffusion of the species is more efficient. The sizes and locations of these cavities are changed by the cosmic-ray impacts. Ivlev et al. (2015) has theoretically revisited the effect of impulsive spot heating induced by cosmic-ray collisions. These two last papers discuss a non-thermal desorption mechanism called "explosive desorption." In this process, when energy is brought to an ice, in which a fair amount of radicals would be frozen, this energy allows the sudden diffusion of these radicals. A chain of exothermic reactions between these radicals would then take place, provoking an explosive desorption of the surfaces. The idea was studied experimentally by D' Hendecourt et al. (1982) through the photolysis of ices. This process has been included in some models (see, e.g., Shen et al. 2004;Rawlings et al. 2013), however more experimental data on this process would be needed, as stated in these findings. Even the very exothermic formation of H 2 on the grains has been proposed to be responsible for non-thermal desorption localized around the location of the H 2 formation (Duley & Williams 1993;Willacy et al. 1994;Roberts et al. 2007). Such a process does not, however, appear to be efficient in the lab (Minissale et al. 2016;Wakelam et al. 2017b). The effect of fast grain rotation on the species desorption rates have been investigated by Hoang & Tung (2019). To explain some methanol gas-phase observations in a cold core, Harju et al. (2020) proposed the grain-grain collisions to be an efficient non-thermal desorption mechanism. In this paper, we test the respective efficiencies of the different mechanisms commonly included in the models (photodesorption, chemical desorption, and cosmic-ray-induced whole-grain heating) and tested the addition of sputtering of ice grain mantles via collision with cosmic rays in the electronic stopping power regime, leading to a localized thermal spike desorption (whose efficiency was recently measured). We describe in Section 2 the chemical and physical model used to study the efficiency of each of these processes. The model results are given in Section 3 and then compared to methanol observations in Section 4. In Section 5, we discuss some of the model assumptions and we present our conclusions in the last section.

General presentation
To study the non-thermal desorption processes, we used the Nautilus gas-grain code. Nautilus is widely detailed in Ruaud et al. (2016) and we only briefly describe it here. It is a numerical model that simulates the chemistry under interstellar conditions. For simplicity, only one grain size of 0.1µm is considered. The chemical processes are described as chemical reactions and the model computes the efficiency of each process based on a number of parameters. For bimolecular gas-phase reactions, for instance, the efficiency of each reaction is computed with a modified Arhenius temperature dependent law with three parameters. In some cases, these parameters are determined by experiments or theoretical calculations. In many other cases, the rate coefficients are guessed based on similarities with other systems. The considered gas-phase processes are bimolecular reactions (involving neutral-neutral and ion-neutral reactions), direct cosmicray ionization or dissociation, ionization or dissociation by UV photons, ionization or dissociation produced by photons induced by cosmic-ray interactions with the medium (Prasad & Tarafdar 1983), and electronic recombinations. Details on how the efficiency of each gas-phase process is computed are given in Wakelam et al. (2012). The chemical species in the gas-phase can be physisorbed on dust surfaces while colliding with the grains based on equations of Hasegawa et al. (1992). We used a sticking probability of 1, except for H and H 2 for which we use temperature dependent values (based on laboratory experiments from Matar et al. 2010;Chaabouni et al. 2012). For species on the surfaces, we make a distinction between the species in the most external layers (2 monolayers here) that we call surface species, and the species below these layers called mantle species. The species are adsorbed on the surface and become the mantle during the construction of the ices. Similarly, when the species desorb, only the species from the surface can desorb but surface species are gradually replaced by the mantle species. Once they are on the grains, these species can diffuse through thermal hopping or tunneling through diffusion barriers. For this particular model, we assumed that all species can diffuse through tunneling with an efficiency that depends on the mass of the species. The equations for diffusion through tunneling are based on Eq. 10 from Hasegawa et al. (1992). The thickness of the diffusion barrier is taken to be 2Å (Ásgeirsson et al. 2017). Both surface and mantle species are capable of diffusion but the process in the mantle is much less efficient. We assume a ratio of 0.4 for the surface species and 0.8 for the mantle species between the diffusion and binding energies. In addition, following Ghesquière et al. (2015), we assume that water drives the diffusion in the mantle, so that all mantle species with binding energy smaller than that of water are set to the water value. The surface reactions are assumed to proceed through the Langmuir-Hinshelwood mechanism. The probability of reaction is assumed to be one for exothermic and barrier less reactions. For reactions with barriers, the probability of reaction is computed taking into account the competition between diffusion and reaction as explained by Ruaud et al. (2016). This probability also depends on the efficiency of the reaction to tunnel through the chemical barrier and the reduced mass of the system. The width of the chemical barrier is taken to be 1Å. In addition to surface reactions between two adsorbed species, the model includes the ionization or dissociation of the surface and mantle species by UV photons and photons induced by cosmic-ray particles interactions with the gas. The values of the rate coefficient for these processes are the same as in the gas-phase because of a lack of laboratory or theoretical data. Kalvāns (2018) proposed applying a scaling factor to the photodissociation rates to use them for surface species. By comparing his model results with ice mantle observations, he determined that the best agreement was obtained for a scaling factor of 0.3. Such a result would, however, be strongly model-dependent. Thermal desorption is included in the model as well as a number of non-thermal desorption mechanisms that will be discussed below. In essence, the model described above is the same as in Ruaud et al. (2016), except for the diffusion through tunneling for all species. The gas and ice chemical networks are the same as in Wakelam et al. (2019).

Non-thermal desorption processes
The Nautilus model from Ruaud et al. (2016) includes three nonthermal desorption processes: photodesorption, chemical desorption, and cosmic-ray (CR) heating. We included in Nautilus a new process, which is the sputtering of the grains by the cosmicray particles (CR sputtering). Photo-desorption, CR heating, and chemical desorption are only allowed for surface species, while CR sputtering can occur for both surface and mantle species.

Photo-desorption
Photodesorption occurs when a single energetic UV photon absorbed close to the grain surface induces the desorption of molecules or radicals. Molecules at the surface can desorb immediately if the energy transferred by the UV photon overcome their binding energy. The energy can also be transferred to a neighboring atom or molecule, inducing an indirect desorption process. Molecules absorbing photons below the surface can also interact with upper layers of molecules or diffuse through the surface and desorb. Another possibility is called co-desorption, whereby the species below the surface takes with it the species above it. The process can in fact be even more complicated when including the photodissociation. For instance, Bertin et al. (2016) and Cruz-Diaz et al. (2016) showed experimentally that methanol would break and its fragments would desorb. Depending on the energy brought to the system, the photo-fragments could also stay close to the surface and recombine. For simplicity, we considered independently the dissociation (see Section 2) and the desorption. In other words, photo-desorption leaves the product intact but photo-dissociations on the surface can occur and then the products can photo-desorb. For the photo-desorption, we have used the simplified formalism for all species described in Ruaud et al. (2016): for the desorption rates induced by direct k des,UV and indirect k des,UV−CR UV photons (in s −1 ). The strength of the direct UV field is F UV = 1.0 × 10 8 photons cm −2 s −1 (Öberg et al. 2007). We adopt a secondary (cosmic ray induced) UV field of F UV = 10 3 photons cm −2 s −1 scaled to an ionisation rate of 10 −17 s −1 . This value is an adapted value considering the ∼3100 photons cm −2 s −1 for an ionisation rate of 3 × 10 −17 s −1 from Shen et al. (2004), whereas Prasad & Tarafdar (1983) gave about 2380 photons cm −2 s −1 for an ionisation rate of 3 × 10 −17 s −1 . In addition, S UV is the scaling factor for the UV radiation field and ζ the cosmic ray ionisation rate in s −1 . The yield (Y pd ) is assumed to be 10 −4 molecules per photons for all species (Andersson & van Dishoeck 2008). Furthermore, N site is the total number of surface sites on a grain while r dust is the radius of the grains. In our case, N site ∼ 1.2 × 10 6 and r dust = 0.1µm. The strength of the UV radiation field induced by the cosmic-rays is scaled with the local value of the CR ionization rate. The factor of 2 in the exponential of the photo-desorption rate by direct UV photons (Eq. 1) is taken from Roberge et al. (1991) and it takes into account the higher extinction of the grains in the UV wavelength range as compared to the visible (used to compute the Av). This is a mean value over the distribution of photons. A more robust approach would be to compute this scaling factor for each molecule of the ice by convolving the spectrum of photo-desorption for each molecule by the extinction curve of the grains. Such a calculation is, however, outside the scope of this paper, but should be considered in the future.

Grain heating by cosmic rays
The energy released by the collision between high velocity cosmic-ray particles and grains can produce a global heating of the grain or a localized one (Leger et al. 1985;Shen et al. 2004), which then cools down to its initial temperature. To include this process in Nautilus, we followed the simple formalism proposed by Hasegawa & Herbst (1993) for the whole grain heating: k des,cosmic = f CR × k des,therm (T peak ), with f CR being the ratio between the cooling time of the grain and the time between two collisions, while T peak is the peak temperature reached by the grain. The f CR and T peak parameters critically depend on the size of the grains as well as their nature and their coverage, as the cooling of the grains occurs mostly by molecular evaporation (Herbst & Cuppen 2006). In the current model, we kept the prescription by Hasegawa & Herbst (1993): f CR = 3.16 × 10 −19 and T peak = 70 K. This process is stochastic and our simple approach does not fully reproduce its complex nature. Kalvāns & Kalnin (2019, 2020a, for instance, studied this process in detail and proposed other ways to include it into astrochemical models. One limitation of our approach is that we used only one grain size, thus restricted to a representative grain size, in order to limit the number of free parameters. Herbst & Cuppen (2006) did a theoretical study of the efficiency of whole grain heating with the size and composition of the grains. They showed that for silicate grains, the peak temperature reached by small grains upon cosmic ray collision increases as the grain radius decreases. Iqbal & Wakelam (2018) showed that considering a grain size distribution in these complex astrochemical models increases strongly the importance of this process.

Chemical desorption
The energy released by exothermic surface reactions can produce partial evaporation of the products. We call this process chemical desorption. The exact description of the process is yet to be understood. As a result, several models have been proposed to include this process into kinetic models. Yamamoto et al. (2019) is the latest published model and compares its results with the two previously published ones: Garrod et al. (2007) and Minissale et al. (2016). The model proposed by Minissale et al. (2016) is the only one based on some experimental measurements but only for a few small systems (O + H, OH + H, and N + N). The three proposed models are all based on a number of unknown parameters. The model from Garrod et al. (2007) depends on the fraction of the energy released by the reaction that is transfered to the produced species. The model from Minissale et al. (2016) depends on the effective mass of the surface (which also defines the fraction of the energy retained by the product). The recent model from Yamamoto et al. (2019) is based on a different theory. It does not assume that the products retain some energy, which is then transformed to kinetic energy, but it assumes that the released energy heats the surface. So their model depends on the thermal diffusivity and specific heat of the surface. Yamamoto et al. (2019) compared the efficiency of the three models for bare grains and showed that this gave similar results. Minissale et al. (2016) showed experimentally that the chemical desorption would be very much less efficient on water ice surfaces than on bare grains. In the simulations presented here, we anticipate such concerns by stating that the grains are covered by water ices (with a significant fraction of CO 2 ). Since the model of Garrod does not explicitly depends on the nature of the surface (except by decreasing the a parameter to an unknown value) and that the model of Yamamoto does not give the parameters to be used for water ices, we included only the model from Minissale et al. (2016). Similarly to what we did in Wakelam et al. (2017b), the fraction of evaporation for singly produced species is computed following: with N the number of degree of freedom of the produced species (N = 3n with n the number of atoms in the produced species) and = (M−m) 2 (M+m) 2 the fraction of energy kept by the product with a mass m. M is the effective mass of the surface. We first compute the effective mass for bare grains of 120 amu and divide it by 10 to obtain the values for water ices. For the three measured systems, we used the measured values (f O+H = 30%, f OH+H = 25%, and f N+N = 50%). For channels producing more than one product, we set the chemical desorption efficiency to zero.

Sputtering by cosmic-rays
When CRs impact solids, many excited electron states are created. In less than a picosecond, these excited states relax in atomic motions and lead to a thermal spike (Toulemonde et al. 2000). As the result of the hot spot, some of the material is sputtered. The sputtering is scaling as the square of the stopping power, or the energy loss by thickness unit dE/dx (Dartois et al. 2015;Mejía et al. 2015). It sets heavy and low-energy CRs as the main contributor to this process.
We added this new type of process into Nautilus assuming that both surface and mantle species can desorb simultaneously. For simplicity, we assume that species desorb without breaking but this may not be the case all the time (as for photodesorption). The rate coefficient of these reactions are computed as follows: where ζ is the cosmic-ray ionisation rate of H 2 (in s −1 ), r dust the radius of a grain, and N site the number of sites on one grain. Y eff is the efficiency of desorption integrated over a cosmic-ray spectral distribution, which is a function of the number of layers of ices (n layers ) (Dartois et al. 2018). Y ∞ is the sputtering yield for thick ices and β, γ are two parameters associated to nature of the ice. For the physical structure adopted here, water ices remain the dominant ice component although CO 2 is also very abundant in the ice. This is due to the high dust temperature assumed (see discussion in Section 3.1). Sputtering on CO 2 ice is much more efficient than on water. We used the sputtering parameters for water ices: Y ∞ = 3.63 , β = 3.25, and γ = 0.57 (Dartois et al. 2018) and we test the values for CO 2 ices in Section 4. For the grain sputtering, in the experiments, the diameter of the craters made by the heavy (above C) cosmic-ray particles in the ices is on the order of a few nanometers. Thus, for grains with sizes above 10 nm, the measured yield should apply. The ice mantle thickness influence on the yield is already introduced in the formalism. For smaller sizes, the sputtering yield should be higher and capable of being summed up with the global heating process.

1D physical structure
To study the effect of the non-thermal desorption processes, we used the 1D cold core physical structure determined by Navarro-Almaida et al. (2020) from observations. These authors derived the density and gas temperature at several positions in the TMC1-C and TMC1-CP cores using CS and CS isotopologs and the RADEX radiative transfer code (van der Tak et al. 2007), together with the Markov Chain Monte Carlo (MCMC) approach (see, e.g., Rivière-Marichalar et al. 2019). The density structure of TMC1 was then obtained by fitting these densities to a Plummer-like analytical density profile, a widely used density profile for prestellar cores (Tafalla et al. 2002;Priestley et al. 2018). The visual extinction profiles of several positions in the cloud were obtained from the Herschel extinction maps in Navarro-Almaida et al. (2020). Assuming spherical symmetry and isotropic UV illumination, the UV shielding at every point inside the cloud is due to an extinction that was taken as half of that measured in the extinction maps. These values of extinction at each position were then interpolated using cubic spline interpolation. Finally, the dust and gas temperature profiles were determined by a cubic spline interpolation of the dust temperatures from the Herschel maps and the gas temperatures from the MCMC approach. The density, gas and dust temperatures, and visual extinctions are shown in Fig. 1 as a function of radius from the center of the core. The density starts around 5 × 10 3 cm −3 at 3.4 × 10 4 au, increases up to 6 × 10 4 cm −3 at 5000 au and then remains flat. The visual extinction is small outside (around 2) and increases up to 10 inside. The gas and dust temperatures decreases toward the inside. They are between 9 and 14.5 K, but the dust temperature is always slightly larger than the gas temperature. We added to this structure a radius dependent cosmic-ray ionization rate. Indeed, the CRs coming from the outside of the cloud are slowing down by the gas and following cosmic-ray ionization rate is decreasing with penetration (Padovani et al. 2009;Chabot 2016). Observations of cosmic-ray ionization rates using molecular ions, show a correlation with the column density (Neufeld & Wolfire 2017;Indriolo & McCall 2012). The H 2 cosmic-ray ionization rate (in s −1 ) was computed following the dashed red line fit presented in the lower panel of For low Av (≤ 0.5), we adopted a constant CR rate corresponding to a non attenuated CR flux of: We note that we are using a static 1D physical structure that does not evolve with time. This is clearly a simplification of the model, but we do not yet have a good constraint on the dynamical evolution of the studied cores. In addition, we are interested in comparing the efficiency of several chemical processes. Such a comparison would be more difficult to make if we added a time dependency on the physical conditions and moving structures.

Other model parameters
Starting from a mix of atoms (with abundances listed in Table 1 and apart from H, which is assumed to be in its molecular form), we ran the chemical model using the 1D physical structure and for a time span of 10 7 yr. The impact of the initial atomic hydrogen abundance is discussed in Section 5.1. The external incident Table 1. Initial abundances (with respect to the total proton density).

Species Abundance
He UV flux is assumed to be five times the Draine's field, as suggested by the observations . To independently study the effect of the non-thermal desorptions, we switched off all of them (no desorption model) and then turn them on one at a time. In the "no desorption model," the thermal desorption is still active, but we checked that in the conditions used here, it does not change the results if no desorption at all is assumed. The exceptions are H 2 and He, which have to be allowed to thermally desorb, otherwise most of the gas would end up depleted.

Model results
The four non-thermal processes presented in the previous section do not depend on the same quantities. The chemical desorption will depend on the abundance of the reactants and the efficiency of diffusion (which depends on the dust temperature).
The photo-desorption depends on both the visual extinction and the cosmic-ray ionisation rate, while the two last processes depend on the cosmic-ray ionization rate. In addition, all species will not be affected the same way. The effect for species essentially formed on the surface should be direct while the effect for species formed in the gas-phase is more complex as it can impact their gas-phase precursors. We separated several groups of molecules and we present the impact of the different non-thermal desorption processes. Figure 2 shows the abundance of the main constituents of the ices (H 2 O, CO, CO 2 , CH 3 OH, CH 4 , and NH 3 ) as a function of radius for an integration time of 6 × 10 5 yr, which is the typical age of an evolved pre-stellar core. We chose this time to emphasize the differences between the models. Indeed, at later times, the differences between the model results are much smaller and are negligible at 2 × 10 5 yr and earlier because the interactions with the grains are less efficient. The figures show that the abundance of the ice species increases toward the higher densities (smaller radii). All desorption mechanisms produce similar ice abundances at high density and similar to the model without desorption. Here, H 2 O and CO 2 dominate the ices for H densities larger than 10 4 cm −3 (Av = 4, inside 1500 au). Going outward, the CO 2 ice abundance drops and water clearly dominates until an Av of 4 (H density of about (5−6)×10 3 cm −3 ). The amount of water then depends on the non-thermal mechanism considered, chemical desorption being the most efficient to decrease it while sputtering is the least efficient. The large CO 2 abundance over CO reflects the grain temperature that is slightly above 10 K in the entire structure. The dust temperature used in these simulations comes from Herschel data (Malinen et al. 2012) and was derived by fitting the SED with gray-body emission (Navarro-Almaida et al. 2020). This procedure is known to overestimate (by about 1-2 K) the dust temperature in the center of the cores because of the contribution of the warmer grains in the external layers of the core surface that are located along the line of sight.

Main ice constituents
The effect of dust temperature is discussed in Section 5.2. At an Av lower than 3, hydrogenated species such as CH 4 , NH 3 , and CH 3 OH ices can be more abundant than CO 2 because there is more free hydrogen thanks to H 2 photo-dissociation.

Simple abundant species
In Fig. 3, we show the model results for simple gas-phase molecules often observed in cold cores. Except for water, these molecules are not directly formed on the grains but we can see that their gas-phase abundances are strongly sensitive to the nonthermal desorption processes because: 1) some precursors can be desorbed from the grains; and 2) after being formed in the gas-phase, they are depleted on the grains and non thermal desorption brings them back into the gas-phase. The CO gas-phase abundance is not strongly sensitive to the non-thermal desorptions. The other molecules (except for water) present a lower level of sensitivity in the outer parts, where the densities are smaller and, thus, the depletion is smaller. In the inner, denser regions, the larger gas-phase abundances are not produced by the same processes for all species. For CN, HCN, HC 3 N, and HCO + , cosmic-ray heating produces the largest abundances, then sputtering followed by chemical desorption and photo-desorption (the two last ones being equally efficient for CN). Here, CN, HCN (and HNC), and HC 3 N molecules are chemically linked. The CN molecule is mostly formed by N + CH, then can react with H + 3 to form HCNH + . While recombining with electrons, HCNH + produces HCN (and HNC); CN can also react with C 2 H 2 to form HC 3 N. The main effect of CR heating is to desorb CH 4 from the ices. Removing this process leads to results similar to the photo-desorption. Once in the gasphase, CH 4 participates to the ion-neutral chemistry by providing CH radicals. CH 4 for instance reacts with H + 3 to form CH + 5 , which recombines with electrons to produce CH 3 . CH 3 then reacts with atomic carbon to produce C 2 H 2 involved in the formation of HC 3 N. The increase of HCO + abundance is also due to the CH 4 CR heating desorption. One of the paths to form HCO + is through the reaction CO + CH + 5 . The gas-phase SO abundance is enhanced by both the chemical desorption and the sputtering in the inside. This molecule is formed by the neutral-neutral reaction S + OH and O + HS. In the case of chemical desorption, both the chemical desorption of HS and OH (during hydrogenation of S and O on the surfaces) are at the origin of the SO increase. We need to remove both of them to decrease the SO gas-phase abundance. In the case of sputtering, OH gas-phase abundance is highly increased. The CS gas-phase abundance is mostly enhanced by the chemical desorption and the CR heating. For the CR heating, it is again the desorption of CH 4 ice that is responsible for this increase.
In the CR heating model, gas-phase CS is mostly produced by three reactions: HCS + + e − , H + HCS, and H 3 CS + + e − . All three precursors, HCS, HCS + , and H 3 CS + , are formed by reactions between neutral or ionized atomic sulfur with CH 2 or directly CH 4 . Similarly to SO, it is the chemical desorption that increases most the gas-phase CS abundance and this is due to the larger HS abundance as it is formed by C + HS in the chemical desorption model. Figure 4 shows the abundance of CH 3 OH, CH 3 CHO, HCOOCH 3 , and CH 3 OCH 3 as a function of radius and for the same time as previously. Only three models are seen on the figure (except for CH 3 CHO) because the two other models (i.e., no desorption and CR heating) produce abundances in the gas-phase that are very small. These molecules are produced on the grains, or from precursors formed on the grains, and only efficiently evaporated by sputtering, chemical or photodesorption. Models with no desorption or only heating by CRs produces negligible gas-phase abundances of these species. The model with chemical desorption is the one producing the largest gas-phase abundance of HCOOCH 3 at all radii (although smaller than 10 −13 ). In our model, HCOOCH 3 is formed in the gas-phase by the reaction CH 3 OCH 2 + O → HCOOCH 3 + H. CH 3 OCH 2 is formed on the grain surfaces by the reaction s-C...CH 3 OH + s-H followed by chemical desorption of the produced species CH 3 OCH 2 . s-C...CH 3 OH is a Van Der Waals complex included in Nautilus by Ruaud et al. (2015).

Complex organic molecules observed in cold cores
The gas-phase abundances of CH 3 OH and CH 3 CHO are larger with the model including sputtering at high density while chemical desorption produces the largest abundances at radii larger than 5000 au. Finally, photo-desorption is the least efficient of the three models at high density but becomes much more efficient than sputtering at radii larger than 15000 au and more efficient than chemical desorption at radii larger than 27000 au. Methanol, CH 3 OH, is formed on the surfaces by successive hydrogenation of CO. The efficiency of the sputtering and photo-desorption are proportional to the abundance of CH 3 OH on the grains whereas chemical desorption depends on  2. Abundance of the main ice components as a function of radius (and visual extinction) for a time of 6 × 10 5 yr. The "no desorption" curve is almost the same as the "sputtering" one.
the abundance of the reactants. We note that in our model, the photo-desorption of methanol is not destructive while experiments by Bertin et al. (2016) showed that is should be partly destructive. In addition, sputtering of the entire mantle along with the surface is allowed, whereas chemical desorption and photo-desorption are only possible for species on the surfaces. Although the cosmic-ray ionization rate decreases inside the cloud, the flux of production of CH 3 OH is about ten times larger in the sputtering model at the inner point than with the chemical desorption model. Photo-desorption is more efficient at low visual extinction (i.e., in the outer part of the cloud). Although CH 3 CHO shows a similar sensitivity to the different desorptions, its formation path is different. In the chemical desorption model (at all radii), it is the gas-phase reaction O + C 2 H 5 that produces gas-phase CH 3 CHO and C 2 H 5 is produced mostly by H + C 3 H 7 , while C 3 H 7 is desorbed at the end of a long hydrogenation chain of surface reaction of C 3 . In the sputtering model, the O + C 2 H 5 reaction plays also a role (C 2 H 5 being evaporated by sputtering from the surfaces), especially in the outer parts, but the large increase of gas-phase CH 3 CHO inside is due to the dissociative recombination of CH 3 CHOH + , itself produced by the reaction CH 3 OCH 3 + H + . The CH 3 OCH 3 abundance at the highest density is indeed quite large, as seen in Fig. 4 for this model.
For CH 3 OCH 3 , it is the chemical desorption that is the least efficient whatever the radius while sputtering is the most efficient for radii inside 22000 au and photo-desorption outer this radius. CH 3 OCH 3 is formed on the surfaces via reactions such as s-H + s-CH 3 OCH 2 and s-CH 3 + s-CH 3 O. The surface abundance of this species is similar in all the models. The binding energy of the species results in a very small fraction of chemical desorption while the other mechanisms are more efficient.
The main production reactions discussed here were found by looking at the fluxes of the reactions producing these species. While doing various tests, we found that these COMs were particularly sensitive to the chemistry of Van Der Waals complexes from Ruaud et al. (2015) that we introduced into the network. By removing these processes, we obtain much fewer of these species both in the gas-phase and on the grains. This is particularly true for CH 3 OCH 3 , whose surface precursor CH 3 OCH 2 is formed by s-C...CH 3 OH + s-H → s-CH 3 OCH 2 . In addition, s-CH 3 CHO is also impacted because one of the channels producing this species on the surface involves C + s-CO → s-CCO  and C + s-H 2 CO → s-H 2 CCO. The species presented in the two previous (Sections 3.1 and 3.2) are not sensitive to this chemistry. We note that our model predictions for COMs cannot be compared to Jin & Garrod (2020) because we have used a different physical model, which is much less dense than theirs. Their chemical model includes new non-diffusive surface mechanisms to enhance the production of COMs and does not include the chemistry of Van Der Waals complexes from Ruaud et al. (2015). They found that the chemical desorption was the main efficient process to desorb COMs from the grains, mostly because of the H-abstraction from grain-surface COMs, followed by recombination, amplifying the chemical desorption. We include these re-actions as well and we find that sputtering could play a major role at the highest densities of our model.

Considering abundances
Methanol is an interesting species as its formation is almost entirely on the grains as can be seen by the very small gas-phase abundance obtained with the model without any non-thermal desorption (Fig. 4). The different non-thermal desorption processes studied in the previous section are potentially efficient at the same time with respective efficiencies although for methanol sputtering, chemical desorption, and maybe photo-desorption dominate. We have run an additional model in which we included all non-thermal desorption mechanisms with the same 1D physical structure. In Fig. 5, we show the gas-phase abundance of methanol computed by the model for different times (between 10 5 and 10 6 yr) as a function of density (rather than radius as shown in the previous section). We show the model results when each of the important non-thermal desorption mechanism is added or all of them. In the same figure, we also plot observed methanol abundances. After 10 6 yr, the gas-phase CH 3 OH abundance evolves less quickly. The observed methanol abundances have been derived from observations of four rotational transitions of methanol toward TMC-1C, TMC-1(CP), TMC-1(NH3) at 3mm as a part of the "Gas phase Elemental abundances in Molecular Clouds" (GEMS) IRAM 30m large program . The observed lines have been modeled within the software CASSIS 1 (developed by IRAP-UPS/CNRS, Vastel et al. 2015) with the RADEX 2 non-LTE radiative transfer code (van der Tak et al. 2007), using the Markov Chain Monte Carlo method (MCMC), more details will be presented in an upcoming paper (Spezzano et al. in prep). Observed column densities at each cloud position, with the H 2 column densities, are given in Table A.1.
The modeled gas-phase methanol is high at high density at 10 5 yr and decreases with time. At low density, on the contrary, the abundance is small and increases with time. The observations show a rather flat abundance of methanol with density, with a small increase toward the small densities. With respect to the observations, the model with all desorption mechanisms and times between 3.6 × 10 5 and 10 6 yr best reproduces the observations. Around 3.6×10 5 yr, the observations are best reproduced at high density thanks to the cosmic-ray sputtering while at smaller densities, the observations are best reproduced for later times thanks 1 http://cassis.irap.omp.eu 2 http://www.strw.leidenuniv.nl/∼moldata/radex.html to the chemical desorption. At 3.6 × 10 5 yr, the modeled and observed abundances are in agreement within a factor of 10. An underproduction of gas-phase methanol at small densities ( < 2 × 10 4 cm −3 ) could indicate a more efficient chemical desorption under these conditions as it would be if the water ice coverage of the grains were small. Experiments from Minissale et al. (2016) showed a much more efficient chemical desorption for bare grains and the efficiency would decrease as the grains are covered with water. Navarro-Almaida et al. (2020) found that the observed H 2 S gas-phase abundance in the same regions with densities smaller than 2 × 10 4 cm −3 could only be reproduced by the chemical model assuming the high values of chemical desorption for bare grains. They suggested that this density would represent a change in the chemical composition of the surface of the grains. Our results on methanol goes in the same direction. Using the current prescription of Minissale et al. (2016) for water ice surfaces, the fraction of produced CH 3 OH that desorb in the gas-phase is 0.06%. We changed this fraction in the model and found that the observations at low density can be reproduced with an efficiency ten times larger (0.6%). The experiments of Minissale et al. (2016) could provide upper limits on the CH 3 OH chemical desorption on water ice or bare surfaces of 8%, which is much greater than what we would need. Alternatively, laboratory experiments on CR sputtering have shown to be more efficient with CO or CO 2 ices, rather than pure water ices, and our model indeed predict an ice CO 2 abundance as high as water and even larger under some conditions. In summary, to reproduce the flat CH 3 OH abundance as a function of density for a single time, we need to change the efficiency of the chemical desorption or the CR sputtering with the radius, which would be consistent with a change in the ice composition. There are too many uncertain parameters to be able to constrain quantitatively the efficiency of these processes this way but to illustrate the model sensitivity, we have run three additional models here representing extreme cases. The first one uses Minissale's bare grain prescription for the chemical desorption Article number, page 9 of 19  and the second model uses the sputtering parameters for CO 2 ices (Y ∞ = 21.9, β = 56.3, and γ = 0.6, Dartois et al. 2021). For the third model, we have considered all processes and the prescriptions for bare grains and CO 2 ices. These three models are shown in Fig. 6 for two different times (3.6×10 5 yr and 1×10 6 yr) and noted as "chem des high", "sputtering high", and "all desorptions high". In the same figure, we also report the previous models. The sputtering from an ice mostly composed of CO 2 is much more efficient than from water ice producing at high density a CH 3 OH gas-phase abundance almost ten times larger. The enhanced chemical desorption produces larger CH 3 OH gasphase abundance than the enhanced sputtering at all densities. At a small density (< 2 × 10 4 cm −3 ), a high efficient chemical desorption, as would be expected to occur on bare grains, seems necessary to reproduce the high observed gas-phase methanol abundance. At larger densities (> 2 × 10 4 cm −3 ), however, the enhanced chemical desorption seems overly efficient, while the enhanced sputtering appears to be a major asset.
The modeling results are strongly dependent on the physical conditions. The external UV field for instance was set to 5 in Draine units based on observational constraints. Decreasing this value increases the methanol abundance. The 1D physical model used here does not cover the points observed at the highest densities (n H = 8 × 10 4 cm −3 and higher). The fact that sputtering produces larger methanol abundance than chemical desorption occurs at higher density with these models than the less efficients ones (above n H = 6 × 10 4 cm −3 instead of 4 × 10 4 cm −3 ). We checked that this tendency seen on the curves is confirmed by running higher density models.

Considering column densities
In order to compare the model to the observations in a different way, we reconstructed the theoretical column densities predicted by our model. To do so, we assume a spherical geometry and integrated the column density along each line of sight (see Navarro-Almaida et al. 2021, for more details on the method). We then obtained the column densities as a function of distance from the centre of the clouds shown in Fig 7. The observed radii have been computed for a source distance of 140 pc. We have four observed positions and three cores. The theoretical column densities decreases rapidly as we move away from the centre of the clouds. The conclusions on the most efficient processes as a function of radius (or density)  . Column densities of gas-phase methanol as a function of distance (in au) to the centre of the source. The different graphs represent different times. The lines in each box represent the models in which one of the non-thermal desorption processes has been added or all of them. The points are the observed column densities as described in the text and listed in Table A.1. discussed in the previous section stand here. Our model still underestimate the observed column densities at the larger radii (smaller densities). At early times (1-2×10 5 yr), the inner column densities are better reproduced while the outer column densities (at 16800 au) are strongly underestimated. At later times, the inner predicted column density is smaller but the outer ones are larger. The predicted column density profile is then flatter. The overall observed column densities (which are rather constant across the clouds) are better reproduced for the times around 3.6 − 6 × 10 5 yr. In that case, the predicted column densities are below the observed ones by less than a factor of 10. From the observational point of view, the column densities are what is measured (with various assumptions for the radiative transfer analysis). The observed abundances used in the previous section were obtained by dividing these column densities with H 2 column densities computed from Herschel dust observations at the same positions and with a similar beam size. When computing the observed abundances with respect to H 2 , the main hypothesis is that the methanol lines and the dust observations probe the same column of material. When we compared the observed and predicted abundances, we used the local densities measured at each position (and deduced from the molecular excitation). Since our model, at the positions of interests, has the 1D structure playing a minor role (the Av is high enough to prevent the direct UV photons to play a major role), the comparison does not depend on the assumed 1D structure, which may not represent the three sources altogether. When comparing the column densities, we are much more dependent on the assumed density profile.

Effect of the initial atomic abundance
Starting with no atomic initial abundance is a very common simplification of cold core chemical modeling. Molecular hydrogen is the first molecule to be formed in the interstellar medium. Its formation involves several micro-physical processes on top of interstellar grains, going from low energetic processes (sticking, diffusion and Langmuir-Hinshelwood mechanism) to more energetic ones (Eley-Rideal and "hot atom"). The efficiency of the H 2 formation depends strongly on the nature of the surfaces and their shape (see Wakelam et al. 2017a, for a review of on the H 2 formation in the ISM). In astrochemical models of cold cores with large networks, namely, models that are not dedicated to the formation of H 2 itself (such as Biham et al. 2001;Cuppen et al. 2010;Cazaux et al. 2011;Bron et al. 2014), only the Langmuir-Hinshelwood mechanism is considered for the formation of H 2 , which may underestimate the formation efficiency of H 2 at moderate density. From an observational point of view, H 2 can only directly be observed in wavelengths (UV, near-IR, and mid-IR) and is not accessible in dense cold regions. In such regions, we can have some estimates of the H 2 abundance, and so the remaining atomic hydrogen, by measuring the residual atomic hydrogen fraction via HI Narrow Self-Absorption (HINSA) observations. Using such a method, Li & Goldsmith (2003) and Goldsmith & Li (2005) estimated the ratio between H and H 2 between 10 −4 and a few 10 −3 . Using a H/H 2 ratio of 10 −3 would lead to a H abundance of 5×10 −4 with respect to the total proton density. We ran our models for various values of initial hydrogen (5 × 10 −4 , 10 −3 , and 10 −2 ). Initial atomic hydrogen abundances of 5 × 10 −4 and 10 −3 produce minor differences as compared to the results presented in the previous sections. Starting with a higher abundance of 10 −2 (which may not be realistic) would strongly affect the methanol abundance (much less than the other species discussed in this paper) at early times and for densities larger than a few 10 3 cm −3 . At a H density of 2 × 10 4 cm −3 , the methanol gas-phase abundance could be increased by more than one order of magnitude with respect to the one computed by our "all desorption" model at 10 5 yr. This difference fades with time and disappears at 10 6 yr. For the present physical structure, an initial H/H 2 ratio different from one could help increase the modeled abundances for intermediate spatial points but only if the initial H abundance is high and only for early times (see Fig. B.1, to be compared to Fig. 5).

Effect of the dust temperature
For the results presented here, we used the dust temperature obtained from Herschel observations (following Fuente et al. 2019;Navarro-Almaida et al. 2020). The obtained 1D dust profile is always above 10 K at all radii (see Fig. 1) and above the gas temperature. The surface chemistry can be very dependent on this parameter. As such, our predicted abundance of ice CO is low while the CO 2 ice contains a large fraction of the oxygen (as well as water). This has consequences for the predicted ice and gas methanol. We redid all our simulations with a dust temperature equal to that of the gas, which mostly changes in the high-density regions (≥ 2 × 10 4 cm −3 ). In Figs. C.1 to C.3, we show the abundance of the same species as in Section 3 (Figs. 2 to 4) computed with a dust temperature equal to the gas one at each radius of the 1D structure. There are a number of differences between the two sets of models at high density. As such, the following results are focused on these regions. Thus, we refer to the standard model as the one in which the dust temperature is determined by Herschel observations (and shown in the remainder of the paper). The most obvious difference is the larger CO ice abundance, which propagates to the CH 3 OH ice (and H 2 CO ice not shown in the figures). As a consequence, the CO 2 ice is predicted much smaller going inward the cold core. The smaller grain temperature decreases the efficiency of the CO 2 ice formation (which possesses a barrier) and so more CO is available for H 2 CO and CH 3 OH. At 6 × 10 5 yr, the inner CH 3 OH ice abundance is two times larger with the smaller grain temperature (whatever the model presented on the figures as they all give the same abundance). Another consequence of the colder grains is a higher H 2 O ice abundance by a factor of two, which is produced by a higher H 2 abundance on the grains (forming water by reacting with OH on the surface). The higher H 2 abundance in the ices is due to less efficient thermal desorption. Indeed, in our models, H 2 and He are both allowed to thermally desorb and considering their low binding energies, this process is efficient in our standard model, whereas it is negligible for H 2 in this lower dust temperature model. The CH 4 and NH 3 ices are not changed. The lower dust temperature produces some significant changes on the simple molecules commonly observed in the gas-phase (comparing Fig. 3 Ruaud et al. 2015). Overall, a higher grain coverage of CO produces a slight increase of the level of oxygen in the gas-phase, increasing the gas-phase abundances of SO, H 2 O and HCO + . The higher abundance of HCN, in the "photodesorption" and "no desorption models," is also due to a greater abundance of oxygen in the gas-phase but less directly. In our standard model, part of the HCN is destroyed by Si + (to form SiNC + ). In the colder dust model, the increase of OH abundance in the gas-phase reduces the Si + abundance by reacting with it to form SiO + . To find these complex chemical effects, we have tested the model each time by switching off and on the processes. The methanol-increased abundance on the grains propagates to the gas-phase (comparing Fig. 4 to Fig. C.3). The other three COMs discussed in this paper (CH 3 CHO, CH 3 OCH 3 , and HCOOCH 3 ) are not significantly affected by the smaller grain temperature. The overall effects of the various non-thermal desorption processes previously discussed stand. The results presented in section 4 are not significantly changed either. The methanol gas-phase abundance is just slightly more abundant at high density (see Fig. C.4).

Conclusions
In this work, we used the full gas-grain model Nautilus and included a new non-thermal desorption process, which is the sputtering of ices by cosmic-ray particles, experimentally studied in the laboratory by Dartois et al. (2018). We tested its efficiency with respect to the non-thermal desorption mechanisms commonly included in astrochemical models: chemical desorption, grain heating by cosmic-rays, and photo-desorption.
We also tested our model prediction with methanol observations in the TMC-1C, TMC-1(CP), and TMC-1(NH3) cold cores. We focused the discussions on the main ice components, simple molecules usually observed in cold cores (CO, CN, CS, SO, HCN, HC3N, and HCO + ), and complex organic molecules (COMs such as CH 3 OH, CH 3 CHO, CH 3 OCH 3 , and HCOOCH 3 ). Our conclusions are as follows: -We found that all species are not sensitive in the same way to the non-thermal desorption mechanisms and the sensitivity also depends on the physical conditions. It is then mandatory to include all of them. Grains heating by cosmic rays appears to impact significantly some of the species formed in the gasphase (such as CN, HCN, HC 3 N, and HCO + ) through the desorption of CH 4 ices at high density. For molecules formed on the grains (such as H 2 O, CH 3 OH, CH 3 OCH 3 ), the sputtering of the ices induced by cosmic-ray collision dominates the desorption at high density while photo-desorption dominates at low density (although resulting the gas-phase is much less high as the ice abundances are smaller). -A direct comparison of the gas-phase methanol in TMC-1C, TMC-1(CP), and TMC-1(NH3) cold cores with our model predictions shows that our models can reproduce the observations within a factor of 10 at 3.6 × 10 5 yr for all densities. Chemical desorption seem essential to reproduce the observations for H densities smaller than 4×10 4 cm −3 , while sputtering is essential above this density. -The models are, however, systematically below the observed abundances. Considering a more efficient chemical desorption, as measured on bare grains by Minissale et al. (2016), appears to head in the right direction for low densities but also seems a bit extreme. Using a higher efficiency for the sputtering, as measured in CO 2 ices, increases the methanol abundance in the gas at high density. -Comparing CH 3 OH observed and modeled column densities, rather than abundances, leads to similar conclusions except that the inner column densities are fully reproduced at early times (1-2×10 5 yr); however, in that case, the most external column densities are strongly underestimated by our models.
The best agreement at all radii is obtained for times between 3.6×10 5 and 6×10 5 yr. In that case, the predicted column densities is flatter with the radius but below the observed column density by less than a factor of 10. -Considering that a small fraction of the hydrogen is not in H 2 at the beginning of the chemical calculation only impacts the methanol abundance if this fraction is as high as H/H 2 = 2 × 10 −2 , which would be unrealistic. -In our simulations, we used the dust temperature observed with Herschel, which is a little above the gas temperature determined based on the molecular excitations. Setting the dust temperature equal to the gas one decreases the dust temperature in the inner regions of the cores by a few Kelvin. Such changes will affect the ice reservoir by producing less CO 2 , increasing the CO and CH 3 OH abundances. This higher CH 3 OH ice abundance propagates in the gas-phase by increasing the gas-phase abundance by a factor of a few. The change in CO ice coverage slightly impacts the oxygen quantity in the gas phase because, in our model, the atomic oxygen physisorbed on CO ice partly releases O and OH in the gas phase either by interaction with direct and secondary UV photons and hydrogenation.
In conclusion, the sputtering of ices by cosmic-ray collisions may be the most efficient desorption mechanism at high density (a few 10 4 cm −3 under the conditions studied here) in cold cores, while chemical desorption is still needed at smaller densities. Additional studies are needed on both chemical desorption and CR sputtering to assess their efficiency with respect to the main ice composition (especially in the presence of mixtures).  ) in the case of the model starting with initial abundances of zero, 10 −3 , and 10 −2 . Two times are shown: 3.6 × 10 5 and 6 × 10 5 yr. The points are the observed abundances as described in the text. In all cases, all the desorption mechanisms were included.  Abundance of the main ice components as a function of radius (and visual extinction) for a time of 6 × 10 5 yr, and the model in which the dust temperature is set equal to the gas one. The "no desorption" curve is almost the same as the "sputtering" one.