The Effects of Cosmic Rays on the Chemistry of Dense Cores

Cosmic rays are crucial for the chemistry of molecular clouds and their evolution. They provide essential ionizations, dissociations, heating and energy to the cold, dense cores. As cosmic rays pierce through the clouds they are attenuated and lose energy, which leads to a dependency on the column density of a system. The detailed effects these particles have on the central regions still needs to be fully understood. Here, we revisit how cosmic rays are treated in the UCLCHEM chemical modeling code by including both ionization rate and H2 dissociation rate dependencies alongside the production of cosmic ray induced excited species and we study in detail the effects of these treatments on the chemistry of pre-stellar cores. We find that these treatments can have significant effects on chemical abundances, up to several orders of magnitude, depending on physical conditions. The ionization dependency is the most significant treatment, influencing chemical abundances through increased presence of ionized species, grain desorptions and enhanced chemical reactions. Comparisons to chemical abundances derived from observations show the new treatments reproduce these observations better than the standard handling. It is clear that more advanced treatments of cosmic rays are essential to chemical models and that including this type of dependency provides more accurate chemical representations.


INTRODUCTION
Cosmic rays (CR) play a vital role in the chemistry of cold (10-30K), dense (> 10 2 cm −3 ) molecular clouds as they can pierce deep into them, unlike the interstellar UV radiation (for a review see Indriolo & McCall 2013). These high energy interstellar particles, primarily consisting of protons but can be heavier elements and electrons, have large energy ranges, up to ZeV energies (Blandford et al. 2014). Although the energies can be high, it is the lower energy CRs (≤ 1 TeV) that affect the dense interiors (Viti et al. 2013;Padovani et al. 2020). In these regions, CRs have a wide variety of effects, one of the most important being a producer of atomic hydrogen through the dissociation of H 2 (van der Werf et al. 1988;Montgomery et al. 1995;Li & Goldsmith 2003;Goldsmith & Li 2005;Padovani et al. 2018a). Other important effects are: being the dominant source of ionization; regulating the degree of coupling of the gas and the magnetic field; having an important role on the dynamics and the collapse timescale of collapsing clouds (e.g., Padovani et al. 2013Padovani et al. , 2014; providing heating and energy to dust grains (de Jong & Kamijo 1973;Kalvāns & Kalnin 2019;Sipilä et al. 2020Sipilä et al. , 2021Silsbee et al. 2021); producing internal UV photons (Prasad & Tarafdar 1983); may have a role on the charge distribution on dust grains (Ivlev et al. 2015); influencing disk growth (Kuffmeier et al. 2020); and affecting deuteration (Caselli et al. 2008). For example, each species ionised by a CR releases an electron. This secondary electron can cause further collisions, which in turn, depending on the energy, can induce more ionizations and heating ). If a secondary electron does not have enough energy to ionise a species, the species may become excited . Excited species produced arXiv:2206.11167v2 [astro-ph.GA] 5 Jul 2022 by CR bombardment have energy levels higher than their base counterparts, allowing these excited species to overcome some reaction barriers that would otherwise be difficult in cold environments. These species have been shown to drive more complex chemistry from reactions that can form interstellar complex organic molecules (Abplanalp et al. 2016).
Although CRs can pierce deep into the molecular clouds, they are still attenuated as they collide and lose energy. The denser the region is, the lower the CR ionization rate becomes (Padovani et al. 2018b). This leads to a dependency of the ionization rate on the density of a region, more precisely on the H 2 column density passed through by CRs.
As the Earth is shielded from the low-energy spectrum of CRs through solar modulation (see Potgieter 2013, for a review on solar modulation), measurements of the CR ionization rate taken from Earth are not indicative of measurements in the interstellar medium (ISM) and are in fact lower. Observations of molecules that are dependant or sensitive to the CR ionization rate (for example, H + 3 is produced from CR ionization of H 2 ) can be used as a tracer for the ionization rate (see Viti et al. 2013, for a review). The "typical" value for the CR ionization rate is often taken to be around the order of 10 −17 s −1 (e.g., Spitzer & Tomasko 1968;Solomon & Werner 1971;Herbst & Klemperer 1973;Li & Goldsmith 2003). It is necessary to note that while this may be known as the "typical" rate, observations shown environments with significantly higher rates. Diffuse clouds have been observed with ionization rates in the order of 10 −16 s −1 (Indriolo et al. 2007;Indriolo & McCall 2012) and rates of up to 10 −14 s −1 have been observed within the inner 300 parsecs of the Galactic centre (Oka et al. 2005;Le Petit, Franck et al. 2016). Recently both Voyager spacecrafts have passed beyond the heliopause, and have been observing lower energies of the CR spectrum (as low as 3 MeV for both nuclei and electrons) (Cummings et al. 2016;Stone et al. 2019). This data from the Voyager probes can be used to estimate a lower boundary for the ISM ionization rate (Ivlev et al. 2015;Padovani et al. 2018b). In fact, the local CR flux measured by the Voyager probes is thought to be unmodulated by the solar wind. However, the magnetometers on board the Voyager spacecrafts have not yet detected a change in the magnetic field direction, as would be expected if they had passed the heliopause (Gloeckler & Fisk 2015). Furthermore, the ionization rate using the fluxes from Voyager only give a lower limit to the observational estimates in nearby diffuse molecular clouds (e.g., Indriolo & McCall 2012).
The hydrogen chemistry of CRs is essential to the chemical evolution of a cloud. H + 3 is fundamental to the production of many polyatomic gas-phase molecules (Herbst & Klemperer 1973) and is formed through the CR ionization reaction: H 2 + CR → H + 2 +e − +CR and the subsequent reaction where H 2 ionization is the rate limiting reaction. In dense clouds H + 3 can then react via proton transfer with molecules such as CO (to form HCO + or HOC + ), O (forming OH + ), N 2 (forming HN + 2 ) and HD (forming H 2 D + ). See the review by Indriolo & McCall (2013) for a more in-depth summary.
CRs also dissociate molecular hydrogen in the ISM through the reaction: In high density regions this reaction is the only form of H 2 dissociation, as the UV photons for photodissociation cannot penetrate deep into the cloud. This reaction depends on the CR dissociation rate, which is often taken to be equal to the ionization rate. In chemical networks however, the rate is often lower than the "typical" value. In UMIST (McElroy et al. 2013) for example, the H 2 dissociation rate is 1.30×10 −18 s −1 . In Padovani et al. (2018a), it has been shown that the H 2 dissociation rate is higher than is often represented in chemical networks, is not a constant value and is not equal to the ionization rate. The rate is dependent on the secondary electrons produced from CR ionization and can be represented as a function of column density, similar to the CR ionization rate in Equation 1 (Padovani et al. 2018a).
As discussed (and as seen in similar works like Redaelli et al. (2021)) it is clear that CRs are extremely important to the chemistry of molecular clouds and their evolution, and hence it is essential that their effects are represented accurately within modern chemical models. This paper aims to improve the handling of CRs in gasgrain chemical models, by introducing both the CR ionization rate and the H 2 dissociation rate as functions of column density and to include the ability to produce excited species and their reactions on the grain. The chemical effects of these additions will be tested on models of collapsing cloud cores. These environments are crucial steps in the early stages of star formation and the effects of CRs on these objects where the gas density increases with time and changes the column density of the core, still need to be investigated. In section 2 we discuss the chemical modelling and detail the CR treatments we have included for this paper. In section 3 we describe the effects these treatments have on the chemical abundances of selected species and discuss the main processes involved in these changes, while in section 4 we summarize our findings.

MODELLING
The chemical code selected for this paper is UCLCHEM (Holdship et al. 2017). UCLCHEM is a time dependent gas-grain chemical code, written in modern Fortran. UCLCHEM is an open sourced chemical code, freely available for use and modification. It is diverse in use due to its modular nature. Specific environments (shocks, cores, collapses) each have their own physics module. UCLCHEM uses separate gas and grain networks. The default gas phase network used is the UMIST RATE12 1 network, described in McElroy et al. (2013) and is used for this paper. The grain network used is described in section 2.1.3 below. For more detailed information on UCLCHEM see Holdship et al. (2017)  In Padovani et al. (2018b) a polynomial fit was developed to express the dependency of the CR ionization rate on column density. We have implemented such a fit into UCLCHEM.
where k is an integer from 0 to 9, c k is the fitting coefficient and N is the column density.
Equation 1 is used to calculate the ionization rate at each time step and the calculated rate is used in all chemical reactions that involve the CR ionization rate. Table A1 gives two sets of fitting coefficients. One, labelled as model L, describes the trend of the ionization rate as a function of the column density obtained by using the Voyager data; the other, labelled as model 1 www.udfa.net H, represents the average value of the ionization rate in diffuse clouds (Neufeld et al. 2010;Shaw et al. 2008;Indriolo & McCall 2012;Neufeld & Wolfire 2017 Padovani et al. (2018a) evaluated the H 2 dissociation rate based on the same CR interstellar spectra used to compute the ionization rate. In the following, the dissociation rate is parameterised as a function of the column density by using Equation 1 with the coefficients listed in Table A1. Similarly to the CR ionization treatment, the dissociation rate can now be calculated at each time step of the model. However, this handling of the dissociation rate can only be activated if the CR ionization dependency is also activated and is automatically set to the same model (L or H) as the CR ionization dependency. Figure 1 shows how the CR ionization rate and H 2 dissociation rate differ from the 'standard' handling of UCLCHEM (i.e. the fixed, user-defined value) under increasing density. This particular example shows all three models at a ×1 ionization factor. UCLCHEM handles the ionization rate in multiples of 1.3×10 −17 s −1 , so an ionization factor of ×1 will correspond to an ionization rate of 1.3×10 −17 s −1 and a H 2 dissociation rate of 1.3×10 −18 s −1 for the standard handling. UCLCHEM calculates the column density by multiplying the size of the cloud (in parsecs) by the total hydrogen density (cm −3 ).

Excited Species
UCLCHEM uses a user defined grain network, separate from the gas phase network. The default grain network that is provided with UCLCHEM handles some basic CR and photon interactions, freeze out reactions and chemical desorption and diffusion reactions. This network was used for this paper, with the additions of excited species production and reactions due to CRs: these excited species are added using the principles described in detail in  and used in .
The underlying principles are that CR bombardments of a solid species generally have one of the following outcomes: Where A is the target species, B and C are dissociated products and * represents an excited species. The reaction rates for these interactions are defined in  and follow the formula: Where G Rn is the radiochemical yield for the reaction pathway Rn, (Rn being R1 to R4 above), S e is the electronic stopping cross section. φ ST is the integrated Spitzer-Tomasko CR flux (from 0.3 MeV to 100 GeV) and has a value of 8.6 cm −2 s −1 . ζ is a CR ionization rate scaling factor the CR flux.
Once produced, an excited species will either react with another solid species, or relax back to the ground energy state. The excited species reaction proceeds at the rate: Where f br is the branching ratio, v 0 is the vibrational frequency of the species, N site is the number of physisorption sites on the grain and n dust is the dust density.

Modelling pre-stellar cores
The effects of these additions will be studied in the cases of pre-stellar cores. Pre-stellar cores represent the early stages of low-mass star formation and have densities in the range of ∼10 4 -10 7 cm −3 , depending among other things, on their evolutionary stage. The models will be set to mimic these regions; in each case, the model will start at an initial density of 10 2 cm −3 . At ∼10 6 years, the models will collapse in free fall to a specific final density. To cover the density range, the models have four possible final densities: 10 4 , 10 5 , 10 6 or 10 7 cm −3 . After collapsing, the models are set to run with static conditions until a final time of 10 8 years is reached, in order to investigate the chemical evolution over time. To determine the influence that temperature and radiation field may have with the CR ionization dependency, each will be varied independently (see Table 1 for values). Note: UCLCHEM assumes 1 Habing to be the Galactic radiation field strength. The chemical species that will be analysed are: H 2 O, CS, NH 3(grain) , N 2 H + , NH 3 , CO (grain) , HCO + , H 2 O (grain) and CO 2(grain) . These species are important as some act as tracers of the gasphase of pre-stellar cores, their regions and physical conditions (CS, N 2 H + ) (Lee et al. 1999), some are key species for grain chemistry and the chemical complexity (NH 3 , NH 3(grain) ) (Rodgers & Charnley 2001) and others are some of the most abundant species found in these regions (H 2 O (grain) ,CO (grain) ,CO 2(grain) )(Öberg et al. 2011).
The CR ionization dependency, H 2 dissociation dependency and excited species production and reactions will all be tested individually as well as combined. To test the effects of the CR ionization dependency, under each condition, models will be run with the CR ionization dependence turned off and compared to the same conditions with the L and the H model dependencies activated. The H 2 dissociation rate dependency can only be activated when the L or H ionization model is also selected. To test the influence of the H 2 dissociation rate, L and H models with the dissociation dependency disabled will be compared to the same model with the dissociation dependency activated. The excited species can be activated without the CR or H 2 dissociation rate dependencies. As such their effects will be examined independently of the other additions and then combined together. Table 1 shows the summary of the parameters investigated, their values and their descriptions. In total, a grid of 280 models were run.

RESULTS
When discussing the influence of the CR ionization rate dependency on the chemistry of pre-stellar cores, our simulated core evolution is split into three phases. The pre-collapse phase covers the period of up to ∼ 10 6 years and represents the period of time leading up to the beginning of the cloud collapse (abundances are examined at a time of 10 5 years). This phase of the model has a gas density of 10 2 cm −3 , until the collapse phase where density begins to increase. The cloud collapse phase represents the time in which the cloud is undergoing collapse in free fall and occurs between ∼10 6 -6×10 6 years, depending on the final density. The density here is increasing over this period from the initial density to the selected final density of the model. The postcollapse represents the period of static density, after the cloud collapses to the designated final density. The post-collapse phase will always have a constant density equal to the selected final density parameter of the individual model. In order to assess trends in our simulations, we set a lower limit for an "observable" fractional abundances of 10 −13 ; below this fractional abundance, changes across the parameter space will be considered irrelevant. Additionally, any changes in abundances that are below a factor of 3 will not be discussed, as these differences are not likely observable. Figure 2 shows the effects of adding the ionization rate density dependency on chemical abundances for a final density of 10 4 cm −3 , a temperature of 10 K, a radiation field strength of 1 Habing and an un-adjusted ionization rate factor of ×1. Table A2 summarises the abundance trends. In the pre-collapse phase, the addition of the ionization rate dependency results in reduced abundances for all our selected species. The same trend is seen for both the L and H models with the H model having enhanced effects (i.e. larger reductions in abundances, up to 3 orders of magnitude, see Table A2). We note that during this phase of the pre-stellar core evolution, the CR ionization rate is indeed higher by almost a factor of ∼ 20 for the H model compared to the value used for the standard model, which proves to be very destructive in the early stages.

Density Dependent ionization Rates
This destruction comes from both the increased presence of ionised species and increased grain desorptions. H + and He + are two ions that play important roles in gaseous destruction. For example, with CS, the main destruction route during the precollapse under the standard model comes from its photodissociation into S and C with some contribution from photoionization and reactions with ionised species. The increased ionization rates for the L and the H models result in increased abundances for the H + and the He + ions. The H + ion plays a part in the CS destruction through the route: The importance of this reaction is enhanced with the L model and is significantly more dominant for the H model. The H model also has higher destruction contributions from the He + ion via the routes: He + + CS → S + + C + He He + + CS → S + C + + He A second example of this is H 2 O, where the destruction in the standard model has some contributions from C + and H 2 O reactions but is mostly dominated by H 2 O photodissociation into OH and H. Similarly to CS, the increased abundances of H + at higher ionization rates drove the destruction of H 2 O through the reaction: While gaseous species are seeing reduced abundances due to ions, solid species are seeing these reductions via CR induced desorptions. For example, H 2 O (grain) in the standard model is primarily destroyed by CR induced UV desorptions. With the increased ionization rates of the L model, direct CR desorption also begins to take place. The even higher ionization rates of the H model result in both the CR induced UV desorption and the direct CR desorptions becoming more efficient, resulting in the reduced abundances seen.  During the post-collapse, species only see notable abundance changes with the H model. Gas phase H 2 O and CS have increased abundances while the solid phase CO (grain) and CO 2(grain) show decreased abundances. In this case the solid phase species tend to undergo more destruction, with CO (grain) and CO 2(grain) decreasing in abundances by over 3 orders of magnitude and the gaseous species see increases up to a factor of 20. While these large decreases in abundances are mainly caused by CR induced desorptions, there are also contributions from reduced formation rates. CO freezeout to CO (grain) is the dominant formation route during the cloud collapse. Under the H model this formation method is significantly inhibited during the collapse phase, reducing the amount of freeze-out taking place. CO 2(grain) is also affected by the CO freeze-out inhibition. The primary formation route for CO 2(grain) comes from the diffusion of CO (grain) and OH (grain) . Reduced abundances for both of these species with the H model inhibit the amount of CO 2(grain) formation. The combination of less formation and more desorption result in these significant decreases seen with the H model. These increased desorptions can also influence the gas phase species. After the collapse, the primary H + destruction route for H 2 O is no longer efficient; this fact coupled with the increased desorptions of H 2 O (grain) , result in the increased postcollapse abundances seen. Other species, like CS, are not as reliant on desorptions. After the collapse the primary formation route for CS comes from the photodissociation of H 2 CS into CS and H 2 . The H 2 CS molecule also shows increased abundances for the L model and significantly so for the H model, which in turn leads to more efficient photodissociation.
When the final density is increased, this reduces the effects of the ionization dependency addition. CO (grain) and CO 2(grain) are prime examples of this effect, as the large decreases in abundances seen at 10 4 cm −3 are no longer present at higher densities (abundance changes are now under an order of magnitude). The reduced ionization rates as density increases is the main cause of this feature through reduced CR desorptions. Also, under these conditions CO freeze-out during the collapse is not inhibited with the increased ionization rates. This, along with the lower desorptions, lead to the reduced effects of the ionization dependency.
3.1.1. Effects due to temperature variations Figure 3 shows the L and H models with an increased initial temperature of 20 K, while Table A3 summarises the effects on the L and H models with initial temperatures of 20 K and 30 K. During the pre-collapse phase, when the temperature is increased, only CS and HCO + show notable changes. Both still show a reduction in abundance with the ionization dependency, but the reduction is lower at higher temperatures (CS is reduced by over 2 orders of magnitude and HCO + is reduced by up to a factor 6 for the H model). At higher temperatures only H 2 O, CS and HCO + have abundances above the set limit of 10 −13 . For the post-collapse, increasing the temperature to 20 K results in larger abundance changes from the H model for solid and gas phase NH 3 . However, CS, CO (grain) and CO 2(grain) instead see less of a change than at 10 K (significant in the cases of CO (grain) and CO 2(grain) , where the large reduction in abundance at 10 K H model is no longer seen). At 30 K only CS and NH 3 have their abundance change by over a factor of 3 (both by a factor of ∼8) and in both cases these changes are less than they are at 20 K for the L model. Increasing the density at these temperatures has a similar effect as at 10 K (i.e. reduced changes as density increases).
Gas phase NH 3 abundances are strongly influenced by the solid phase NH 3 abundances. Gas phase formation comes completely from grain desorptions. NH 3(grain) however, at 10 K under all models, is formed via H (grain) and NH 2(grain) diffusion. As temperature increases this reaction becomes less dominant, particularly during the pre-collapse phase and with the H model at 30 K. In this case there is no formation from the diffusion reaction, which results in NH 3(grain) having abundances below the set limit.
CO 2(grain) , at higher temperatures, relies less on CO (grain) diffusion. During the post-collapse the primary formation rate comes from the diffusion reaction: The reduced abundance changes for CO 2(grain) here are a result of the H 2 CO (grain) abundances. Under the L model, H 2 CO (grain) has higher abundances, resulting in more diffusion, which is balanced out by the increased desorption of CO 2(grain) , leading to little change from the standard model. The H model on the other hand, still sees a reduction in abundance, which is a result of the increased desorptions, but the reduction is not as severe as at 10 K due to the H model having increased O (grain) abundances for more diffusion.
3.1.2. Effects due to variations in the radiation field strength Figure 4 shows the differences between the standard model and the models where the new treatment of the cosmic ionization rate is included, when the radiation field is increased by a factor of 10. Table A3 summarizes the results for enhancing the radiation field by a factor of 10 and 100. In the pre-collapse phase, enhancing the radiation factor by a factor of 10 reduces the abundance changes produced by the ionization dependency, which are further reduced when the radiation field is increased by a factor of 100 (abundance changes are up to 1 order of magnitude, see Table A3). During the post-collapse phase, in general, increasing the radiation field enhances the effects of the ionization rate dependency. At a radiation field of 100 Habing only H 2 O, CS, NH 3 , and HCO + are above the 10 −13 threshold and have increased abundances compared to the standard model. H 2 O, NH 3 and HCO + changes are enhanced by ∼1, ∼2 and ∼1 orders of magnitude, respectively. The large abundance increase for NH 3 comes from desorption from the grains. Under the standard and the L model, NH 3(grain) abundances are below the set limit. This is not the case for the H model. The increased grain abundance here is due to the diffusion of H (grain) and NH 2(grain) . Under these conditions, both of these species have significantly higher abundances with the H model than the L or stan-dard. This increases NH 3(grain) formation which then can desorb into the gas phase. For HCO + , under the increased radiation field there are two main formation routes, photoionization of HCO and the H + 3 reaction: Under the L and H models, both CO and H + 3 have much higher abundances, leading to the increased production of HCO + . When the density is increased above 10 4 cm −3 , there are no significant differences between the effects of the ionization dependency at the standard radiation field strength and the effects at increased strengths. Figure 5 shows the effects of increasing the ionization rate by a factor of 10 on the chemical abundances with the ionization rate dependency. Table A3 summarizes the results of increasing the rate by factors of 10 and 100. Increasing the initial ionization rate in this manner proves to be very destructive both with and without the ionization dependency, particularly in the pre-collapse period. Many species have abundances below the set limit with this increased ionization rate factors. Those that are visible show much larger abundance reductions (up to 2 orders of magnitude, see Table A3) than at an un-adjusted rate. During the post-collapse phase, increasing the ionization rate by a factor of 10 and 100, increases the influence of the ionization dependency. Species that had no notable changes under the standard rate now show reduced abundances. CO (grain) and CO 2(grain) again are particularly affected, showing reductions of up to 7 orders of magnitude. These effects are more enhanced with the ×100 initial ionization rate. As with the other conditions, cores of higher densities show reduced effects from the ionization dependency, even with the enhanced initial CR ionization rate. These trends originate for the same reasons as for the models with a standard initial CR ionization rate but are more pronounced (e.g. increased ions and desorption). Table A4 summarises the only cases where the H 2 dissociation rate has any effect. In short, under an enhanced radiation field of 100 times the galactic one, an increase of abundance was seen for the solid species (by a factor of ∼4), at a low density (10 4 cm −3 ) and only for the L model. The other notable effect is seen at an increased initial ionization rate of ×10, for CO (grain) , also at 10 4 cm −3 with the L model, where a significant increase of abundance (∼3 orders of magnitude) is seen. This large increase in abundance can be traced to the diffusion of H (grain) and CO (grain) into HCO (grain) . Under these conditions this is the dominant destruction route for CO (grain) . Under the H 2 dissociation rate this reaction pathway is severely inhibited, reducing the destruction of CO (grain) during and after the collapse. Abundances of H (grain) here are also lower for the H 2 dissociation dependency model, which may explain the inhibition. Table A4 summarises the only conditions where the inclusion of the excited species had a notable effect on the abundances of the selected species. At an increased temperature of 30 K, N 2 H + shows a reduced abundance by a factor of 4 at a density of 10 5 cm −3 . Increasing the initial ionization rate by a factor of 10 reduces CO 2(grain) abundances at a density of 10 4 cm −3 . The most significant effects come from increasing the ionization by a factor of 100. While the higher ionization rates provide more excitations, the increased destruction of the species is not only from their excitation and subsequent reactions. CO (grain) for example is also affected heavily by the H (grain) and CO (grain) diffusion reaction. The addition of the excited species also produce higher abundances of H (grain) , which increased the amount of CO (grain) diffusion.

Density Dependency with Excited Species
In this section models with both the ionization rate and dissociation rate dependencies activated are compared with and without the inclusion of excited species. Table A5 show the effects under "standard" conditions (10 K initial temperature, radiation field strength of 1 Habing and an initial ionization factor of ×1). CO (grain) and CO 2(grain) are the main species affected, and have reduced abundances when the excited species are included. These reduced abundances are caused via the same destructive methods as discussed in the previous subsection.

Effects due to parameter variations
When the excited species are included with the CR ionization and H 2 dissociation in the chemical models, the effects of varying the temperature and radiation field strength are reduced, while the effects of varying the ionization factor are increased. As such these effects are quickly summarised here.
Increasing the initial temperature and increasing the radiation field strength both inhibit the effects of including the excited species. Under a higher temperature, including the excited species only has an effect on N 2 H + (this is at 30 K and a density of 10 7 cm −3 with the H model only) where the species shows an increased abundance. Increasing the radiation field strength only has an inhibition effect on lower densities. At 10 5 cm −3 and above, there are no differences between 1, 10 and 100 Habing.
Including the excited species with increased initial ionization rates of 10 and 100 times standard, has a greater effect than at the ×1 value. The abundance changes are both larger and seen for more species. At a 10 times standard initial ionization rate, effects are only seen with the H model at 10 5 cm −3 and 10 6 cm −3 . In these conditions, most species see reduced abundances (several orders of magnitude for CS, CO (grain) , HCO + and CO 2(grain) ). Under an increased rate of 100 times standard, the species see larger reductions, under the same densities but with the L model instead. Similar reductions are also seen here at 10 7 cm −3 with both the L and H models.

Comparison to observations
In this section, we qualitatively compare our models to a set of observations from the Cyanopolyyne peak ("CP" or Core D; Hirahara et al. (1992)) of the molecular cloud TMC-1 which is thought to currently undergoing rapid core formation (Choi et al. 2017). This core was chosen as it has been well studied and its density is expected to be around 10 4 cm −3 with a temperature of about 10 K. The density of 10 4 cm −3 is a good candidate for this study, as the effects of the ionization rate dependency are greater due to the low density.
Both von Procházka & Millar (2020) and Fuente et al. (2019) report chemical abundances of the CP region. Fuente et al. (2019) reports molecular abundances derived from observations using the IRAM 30m (3mm and 2mm) and the Yebes 40m telescopes. Depending on the setup, the IRAM 30m has a spatial resolution of ∼29" and the Yebes 40m has a HPBW of 42" or 84". On the other hand, von Procházka & Millar (2020) reports a collection of upper and lower abundance limits obtained from ∼20 other studies (see Table 2 in von Procházka & Millar (2020)), with the emphasis that focusing on upper and lower limits somewhat mitigates the errors in observations and modelling. Figure 6 shows the data reported in Fuente et al. (2019), compared to the UCLCHEM models at 10 K and a final density of 10 4 cm −3 with and without the inclusion of the CR ionization rate dependency during the post-collapse phase. The "standard" UCLCHEM model for CS and HCO + tends to display abundances that are near the lower limit, while the L is closer to the central values. The H model over predicts the abundance for CS and under predicts the abundance for HCO + . With N 2 H + all models under predict the abundance, with the L model being the closest and the H model under predicting the most. Figure 7 shows data reported in von Procházka & Millar (2020), also compared to the UCLCHEM models at 10 K and a final density of 10 4 cm −3 with and without the inclusion of the CR ionization rate dependency during the post-collapse phase. In the cases of HCO + , N 2 H + and NH 3 , for all models the post-collapse abundances are within in the upper and lower limits of the observations. In von Procházka & Millar (2020) only an upper limit is noted for H 2 O. In this case the H model exceeds the stated limit while the L and basic model do not.
In the case of this region in TMC-1 and the compared species, the L model of the ionization dependency appears to perform the best of the three models. The H model over or under predicts the abundances on several occasions, suggesting this upper limit for the ionization rate may in fact be too high.
One caveat that must be addressed is how UCLCHEM handles grains. This version of UCLCHEM considers a grain to be a single layer (i.e. no distinction between grain surface and bulk). It is therefore necessary to speculate on the effects a multi-layer grain approach may have on these results. Species in the bulk are somewhat shielded from CR impacts and the subsequent desorptions. As CR desorptions are critical to abundance changes seen for species like CO (grain) and CO 2(grain) , it is likely the abundance changes seen with the inclusion of the CR ionization dependency will be less significant. Excited species in the bulk are also more protected. Desorptions from excitations and excited reactions would be reduced with a greater emphasis on relaxations, again reducing the effects we see in our models.

SUMMARY
In this paper we improve the treatment of cosmic rays in the gas-grain time dependent chemical code UCLCHEM by including the dependency of the cosmic ray ionization and H 2 dissociation rates on the column density of the gas, as well as the excited species due to the cosmic rays on the grains. We then evaluate the effects of these additions on the chemistry of pre-stellar

Species Conditions Behaviour
Ionization rate dependency H2O 10 6 cm −3 , H model Abundance increase larger than other densities. N2H + 10 5 cm −3 , H model Only density to see a notable change (increased abundance). NH3 10 6 cm −3 , H model Only density to see a notable change (increased abundance). HCO + 10 6 cm −3 , H model Abundance increase larger than other densities.
Ionization rate dependency with parameter variations N2H + 20 K, 10 6 cm −3 , H model Abundance change greater than other densities. 30 K, 10 6 cm −3 , H model Only condition to see abundance change at 30 K. Abundance not as reduced as at 20 K.
NH3 (grain) 30 K, 10 5 cm −3 , H model Only condition to see abundance change at 30 K. Abundance further reduced than at 20 K.
CO2 (grain) 30 K, 10 5 cm −3 , H model Only condition to see abundance change at 30 K. Abundance further reduced than at 20 K.
CS 100 times radiation field, 10 4 cm −3 , H model Abundance is not as reduced as at 10 times radiation field.
H2O (grain) 100 times initial ionization rate, 10 5 cm −3 , H model Reduction in abundance is greater than at 10 4 cm −3 . Table 2. Table showing the post-collapse species, conditions and the behaviour that do not follow the general trends of the ionization dependency.  cores. It is evident that the cosmic ray ionization rate dependency on the column density of the core is the most influential of the treatments, with the inclusion of excited species on the grains playing roles only under specific conditions. Our conclusions can be summarized as follows: • In the low densities of the precollapse phase (∼10 2 cm −3 ) the ionization rate dependency is very destructive due to CR induced desorptions and the production of chemically important ions.
• After the core collapses, the inclusion of the dependency of the CR ionization rate on the column density of the core leads to increased grain desorptions, which decreases solid species abundances (and subsequently increases gaseous species abundances), species like H 2 O, CO (grain) and CO 2(grain) are particularly affected by this. Other gaseous species, like CS show increased abundances from dissociations of larger molecules like H 2 CS.
• Changing the physical parameters of the cloud alters the impact of the new treatments in a nontrivial manner. Higher densities have lower ionization rates with the dependency, reducing the abundance changes seen for all species. Increasing the temperature also has a similar effect for CS, CO (grain) and CO 2(grain) (increased formation rates balance out the destruction from the CRs), while NH 3 shows lower abundances due to less NH 3(grain) formation and subsequent desorption. Increasing the radiation field strength enhances the effects of the ionization dependency, which occurs as a result of grain and gas formation routes. NH 3 for example sees increased abundances with the H model due to higher NH 3(grain) grain formation and desorption, while HCO + sees larger changes due to formation in the gas via H + 3 and CO.
• The H 2 dissociation rate dependency and the inclusion of excited species only affect the chemistry of some of the investigated species under specific conditions. The H 2 dissociation dependency increases abundances of some solid species for the L model under two conditions, a ×10 ionization rate and a ×100 radiation field strength. CO (grain) sees these abundance increase due to the inhibition of its primary destruction route, from a reduced abundance of the H (grain) reactant. The excited species, reduce solid abundances at higher ionization rates, particularly with CO 2(grain) and CO (grain) . While the excitations and subsequent reactions reduced the solid abundances, destruction also comes from reactions with H (grain) (which also sees higher abundances under these conditions).
• Chemical models with and without the ionization dependency were compared to molecular abundances in the TMC-1 cyanopolyyne peak from Fuente et al. (2019)

ACKNOWLEDGEMENTS
We would like to thank Jonathan Holdship for his contribution to the chemical modelling, particularly with UCLCHEM's analysis tools. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No 811312 for the project "Astro-Chemical Origins" (ACO). This work is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme MOPPEX 833460.