{'fraction_non_alphanumeric': 0.03218893904162677, 'fraction_numerical': 0.021656864086080335, 'mean_word_length': 4.287091719898759, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Atlas of Cosmic ray-induced astrochemistry

Cosmic-rays are the primary initiators of interstellar chemistry, and getting a better understanding of the varying impact they have on the chemistry of interstellar clouds throughout the Milky Way will not only expand our understanding of interstellar medium chemistry in our own galaxy, but also aid in extra-galactic studies. This work uses the ALCHEMIC astrochemical modeling code to perform numerical simulations of chemistry for a range of ionization rates. We study the impact of variations in the cosmic ray ionization rate on molecular abundances under idealized conditions, given by constant temperatures and a fixed density of 1e4 cm-3. As part of this study we examine whether observations of molecular abundances can be used to infer the cosmic ray ionization rate in such a simplified case. We find that intense cosmic-ray ionisation results in molecules, in particular the large and complex ones, being largely dissociated, and the medium becoming increasingly atomic. Individual species have limitations in their use as probes of the cosmic ray ionization rate. At early time (<1 Myrs) ions such as N2H+ and HOC+ make the best probes, while at later times, neutral species such as HNCO and SO stand out, in particular due to their large abundance variations. It is however by combining species into pairs that we find the best probes. Molecular ions such as N2H+ combined with different neutral species can provide probe candidates that outmatch individual species, in particular N2H+/C4H, N2H+/C2H, HOC+/O and HOC+/HNCO. These still have limitations to their functional range, but are more functional as probes than individual species previously used.


INTRODUCTION
Astrochemistry is a complex web of processes that are hard to disentangle, as physical parameters affect the chemistry of the interstellar medium (ISM) in different ways. The cosmic-ray (CR) ionization rate, ζ CR , is the important ingredient that leads the way for ion-neutral reactions by ionizing neutral atoms and molecules. In particular, deep inside molecular clouds CR particles are the sole source of ionization, as the thick material shields the inner regions from ionizing UV photons (Indriolo et al. 2009). In extreme environments CR ionization can drive the formation of species but also destroy them as molecules are ionized and quickly dissociated due to a large amount of free electrons.
CRs are predominantly charged particles, with protons being the dominating component, and smaller contributions from electrons, positrons, and bare nuclei of helium and other heavy elements (see, e.g. Beringer et al. 2012, for a thorough review).
But the origin of CRs is still not fully understood. It is impossible to directly pinpoint the point of origin as the charged CR particles 1 are deflected and redirected by the Galactic magnetic field.
Additionally, disentangling the direct CR contribution is complicated due to the interactions CR particles suffer on their way to Earth, and their interactions with X-rays.
CRs can be studied with a variety of methods. Astrochemical methods, as discussed here, provide one opportunity. Specifically, these methods are sensitive to the ionization rate caused by CR interactions. Other methods are summarized by, e.g. Grenier et al. (2015). It is, for example, possible to directly detect the protons, electrons, and other constituents of the CR sea in particle detectors.
Unfortunately, the trajectories of CR particles with energies below a few GeV are heavily influenced by the solar wind. One must correct for this so-called solar modulation if one wishes to study these lower-energy CRs via direct detection. This possibly excludes CR data from Voyager 1, which potentially has left the heliosphere. Alternatively, one can explore CRs via indirect methods, i.e., by studying the electromagnetic radiation CRs produce as they propagate through interstellar space.
Synchrotron radiation from CRs moving in magnetic fields is one example. CRs also produce γ rays 1 In the following we use the terms cosmic rays (CRs) and cosmic-ray particles synonymously. when they interact with ISM particles: the collisions produce pions that decay in a variety of ways.
The γ ray observations allow us to study CRs down to energies ∼ 1 GeV (e.g. Yang et al. 2014).
CRs with energies 0.1 GeV are the primary agent responsible for the ionization in the dense sectors of the ISM that are shielded well from UV radiation (e.g., Fig. 14 of Padovani et al. 2009).
Note that the properties of these CRs are poorly constrained by many of the aforementioned experiments: direct detection experiments and investigations of γ rays can only provide reliable statements for energies ≥ 1 GeV. This means that astrochemical studies that are sensitive to the ionization rate can deliver meaningful constraints on low-energy CRs, but it also means that CR-influenced astrochemistry is badly constrained by other observations, such as those of gamma-rays.
The importance of these CR ionization events for initiating interstellar chemistry means that observations of molecules become useful probes for ζ CR . ISM ion-molecule chemistry begins with the formation of H + 3 , which is directly linked to the CR ionization rate as H + 3 can then interact with neutral atoms and small molecules to form larger molecules, such as HCO + , OH + and H 2 O + , which are thus indirectly linked to the CR ionization rate. For an overview of the ion chemistry in the interstellar medium, we refer the reader to Snow & Bierbaum (2008).
As such, observations of ions and other species are essentially indirect probes of the cosmic-ray ionization through its effects on the interstellar chemistry. By observing and interpreting the abundances of molecules in the ISM we can therefore constrain ζ CR . This can, for example, be done using statistical methods, such as χ 2 minimization for fitting parameters. In the case of a few well-understood molecules with a simple chemistry, it is even possible to estimate the CR ionization rate directly from observed molecular abundances (see, e.g. Caselli et al. 1998;Neufeld et al. 2010;Indriolo et al. 2012Indriolo et al. , 2015.
But first we need to get a better understanding of CR-driven chemistry in the ISM by exploring how and which molecules are affected by the CR ionization rate the most. The aim of this paper is therefore to seek better probes that can improve future surveys of the CR ionization rate by studying the distinctive effect that it has on a selection of species often found in large surveys. To develop a first overview here we focus on a -rather idealized -situation in the parameter space of cloud structure. We consider molecules in a medium with constant gas and dust temperature at a density of 10 4 cm −3 that develop out of some initial state assumed to be the same for all calculations. Reality is likely to be more intricate. However, our approach reduces the problem to a tangible one, and the calculations are certainly suited to revealing the astrochemical complexity encountered even under the most simple conditions.
In Section 2 we discuss previous observational measurements. In Section 3 we present our model and a comparative model of derived measurement techniques. The results of the models are separated into discussions of the general trends of the increasing ζ CR in Section 4, and single species and species pairs as effective probes in Section 5 and 6 respectively. We discuss our conclusions in Section 7.

OBSERVATIONAL MEASUREMENTS
The rate of CR ionization events varies depending on the environment and its physical properties.
In Figure 1 we plot CR ionization rates and H 2 column densities derived from observations toward diffuse clouds, low-mass cores, and protostellar envelopes. We also include the models derived by Padovani et al. (2009), who considered two different determinations of the steady-state local interstellar spectrum (LIS) for CR protons by Webber (1998, "minimum";W98) and Moskalenko et al. (2002, "maximum";M02), and two electron LISs by Strong et al. (2000), C00 for their "conventional" model C and E00 for their steeper ("SE") model. We combine the two proton and electron LIS results in the four lines in Figure 1.
Following Herbst & Klemperer (1973), a CR ionization rate of 3 × 10 −17 s −1 is often adopted for dense cores in the solar neighborhood. One should, however, note that quite a range of values for the CR ionization rate has been reported in the literature (e.g. Caselli et al. 1998;Williams et al. 1998;Maret & Bergin 2007). In particular, exploration of the large sample studied by Caselli et al. (1998) yields rates in the range 10 −18 to 10 −16 s −1 . Such estimates are subject to a range of uncertaintiesof which some are examined in this paper. Still, Caselli et al. claim Rachford et al. (2002), and Indriolo et al. (2007Indriolo et al. ( , 2015; low-mass cores are from Caselli et al. (1998), Williams et al. (1998), Alves et al. (2001), and Maret & Bergin (2007); massive protostellar envelopes are from de Boisanger et al. (1996), van der Tak et al.
(2000), van der Tak & van Dishoeck (2000), Doty et al. (2002), and Shaw et al. (2008). The red squares are observations of low-mass cores by Caselli et al. (1998), the gray boxes are observations of diffuse clouds toward Sagittarius B2 by Indriolo et al. (2015), and the gray box shows the range of values derived from the sample of low-mass cores by Williams et al. (1998). The solid lines are from fits derived by Padovani et al. (2009) from proton (Webber 1998, W98), (Moskalenko et al. 2002, M02) and electron (Strong et al. 2000, C00, E00) LIS data. This is an updated version of Figure 15 in Padovani et al. (2009). by an actual variation in the CR ionization rate. If so, a variation of ζ CR by two orders of magnitude in the solar neighborhood is certainly worthy of future investigations.
Many previous studies have sought to quantify the relationship between the CR ionization rate and measurables such as column densities of different ions. In Figure 1 we plot the results of observations for determining the CR ionization rate in different objects. Atomic and molecular diffuse clouds themselves lie in the range 10 −16 − 10 −15 s −1 , while dense molecular clouds and low-mass cores lie in the lower range 5×10 −18 − 1 × 10 −15 s −1 .

MODEL AND MEASUREMENT TECHNIQUES
One of the goals of this study is to present our adopted model and then compare it to equations derived for determining the CR ionization rate from observational quantities.

Model
In this study we have utilized the gas-grain chemical model "ALCHEMIC" developed by Semenov et al. (2010), who also presented a detailed description of the code and its performance. Below we give a brief explanation of the code. A UV photodesorption yield for surface species of 10 −3 is assumed (see, e.g. Öberg et al. 2009a,b). The self-shielding of H 2 from photodissociation is calculated by Equation 37 with Draine & Bertoldi (1996), assuming a total extinction of 0.5 mag with a total column density 1.1 × 10 21 cm −2 (Güver &Özel 2009). The shielding of CO by dust grains, H 2 , and its self-shielding, are calculated using the precomputed table from Lee et al. (1996). We use the chemical network from Albertsson et al. (2013), which was updated for the ortho-para chemistry of H 2 , H + 2 and H + 3 by Albertsson et al. (2014). For initial abundances we adopt the "low metal" initial abundances from Graedel et al. (1982), meaning that the medium begins as half molecular  . Gas and dust temperature as a function of ζ CR , based on heating and cooling rates calculating using the procedure described in Goldsmith (2001). Table 1. Selection of observationally relevant molecules of primary concern in this paper. CH 3 CCH is renamed to C 3 H 4 in the adopted chemical network. In Table 1 we have compiled a list of observationally relevant molecules that we seek to investigate as possible probes. This selection is motivated by a number of studies that are broadly representative of observational investigations of molecules at millimeter wavelengths. Specifically, the line surveys by Foster et al. (2011Foster et al. ( , 2013, Frau et al. (2012), Jackson et al. (2013), and Watanabe et al. (2014) give a broad overview of species that are easily detected in molecular clouds. We include additional molecules that are relevant for particular astrophysical reasons. These species are taken from the following studies: Goldsmith (2001) and Clark et al. (2013), who investigated cooling via spectral lines; Caselli et al. (1998), who introduced DCO + as a tool for CR determination; Armijos-Abendaño Based on Figure 1 we will investigate CR ionization rates in the range 5×10 −18 −1×10 −15 s −1 distributed in 25 logarithmically separated steps. Because CRs are heat sources, we include a routine to calculate the CR heating and cooling rates using the procedure described in Goldsmith (2001) to adjust the temperature based on the typical temperature of the modeled object. For these calculations we adopt a dust temperature of 10 K, consider negligible gas-dust coupling, and calculate gas temperature from balancing the heating and cooling rates based on the depletion factor at a density of 1×10 4 cm −3 . This gives us a improved estimate on the temperature structure of our object.  CR ionization rate. We adopt = 0.07 ± 0.04 (Indriolo et al. 2012). The gray dashed line represents a perfect agreement between the CR ionization rate relations and the actual modeled CR ionization rate.

Molecules
So far, one of the most practiced methods for constraining ζ CR is using ions and other species as indirect probes by interpreting their abundances based on our understanding of the interstellar chemistry. In Figure 3 we plot the calculated versus the modeled CR ionization rate. Equations are based on studies by Caselli et al. (1998) for dense clouds (top panel), and Neufeld et al. (2010) and Indriolo et al. (2015) for diffuse clouds (bottom panel). The dotted line represents a perfect agreement between the equations and our model. The relation by Caselli et al. (1998) fits our model results quite well, within a factor 2-5. We also see changes between 1 and 10 Myr, but the chemistry does not change significantly after 10 Myr.
The equations from Neufeld et al. (2010) and Indriolo et al. (2015) have a much larger disagreement with the models for diffuse clouds. But if we exclude the lower end of the CR ionization rate, 10 −16 s −1 as CR ionization rates are higher than those in dense clouds, they agree within a factor of 2-5.
It appears that there is not much variation already after 1 Myr for diffuse clouds. Caselli et al. (1998) determined the ionization degree and cosmic ionization rate for 24 dark cloud cores based on observations by Butner et al. (1995). Employing simple steady-state models for deriving HCO + /CO ratios (see, e.g. Wootten et al. 1979;Guelin et al. 1982), they derived the following relation: where δ is the reaction rate for H + 3 (and H 2 D + ) due to reactions with neutral species such as CO and O. Caselli et al. (1998) is the depletion factor, i.e. the abundance fraction of gaseous CO over total (gas phase plus ice mantle) molecules.
[X] denotes the absolute abundance of species X. The rates k i of other reactions included are as follows: For the comparison we adopted a TMC 1-like molecular cloud model with n H = 10 4 cm −3 , A V = 10 magnitudes, n(e − ) calculated from the models, and compared our model to Equation 3.
Modeling of dense clouds is easier, as we had high extinctions that block out UV photons and other additional processes, and we are left with CR photons as a dominant process. But some of the results presented here can also be transferred to diffuse clouds, and for comparison we also briefly discuss CR where k j are the reaction rates for the following reactions: and is an efficiency factor introduced by Neufeld et al. (2010) that accounts for the fact that not every CR ionization event results in the net production of an OH + molecule (e.g. = 0.07 ± 0.04; Indriolo et al. 2012). For this comparison we assume a similar diffuse cloud model to Indriolo et al. (2015) with n H = 35 cm −3 , T = 100 K, n(e − ) = 5.25 × 10 −3 cm −3 , = 0.07 ± 0.03 and A V = 1 magnitude.
Because of the intense conditions, these environments are often predominantly atomic. Hence, we initiate our models with atomic abundances, including hydrogen.
It is clear that the equation overestimates ζ CR . But this overestimation decreases with time and at higher modeled ζ CR it even underestimates the rate. It shows the problem with assuming steady-state chemistry. Hollenbach et al. (2012) argued, for example, that steady-state is achieved on a timescale ∼ 10 9 yr/n(H), corresponding to ∼ 30 Myr at the density of 35 cm −3 adopted here. It is conceivable that gas in the diffuse ISM is subject to relatively constant conditions over long periods of time. In that case estimates following the Indriolo et al. (2015) formalism do yield correct results. Still, the caveat of steady-state chemistry should be kept in mind. UV photons are contributing significantly to the production of electrons in diffuse clouds, and this is taken into account in our model, but not the equations.
While the two relations reproduce our model results quite well, there are still discrepancies. Our aim is to outline the use of molecular probes of ζ CR that are now routinely accessible with advanced facilities such as the Atacama Pathfinder Experiment (APEX) or the Atacama Large Millimeter/submillimeter Array (ALMA). In the next sections we will review and discuss the chemical species that are the most promising as probes of the CR ionization rate over the full ionization range and chemical ages in dense clouds. The identified list of species might also help with finding possible probes in circumstances where the standard probes are not available.

GENERAL TRENDS
We begin by looking at the overall trends in the chemistry to get an idea about the effects an increase in ζ CR has on the chemistry of molecular clouds. Because our selection of species is small, we have calculated these trends based on the full chemical network, consisting of approximately 1300 species. In Figure 4 we look at the ratio of atomic versus molecular abundances (  The atomic versus molecular abundance ratio in the left panel is calculated as the ratio between the mass-weighted sums of the two groups: Up to 10 kyr the ratio is fairly flat and not affected by the increasing ζ CR . It is not until later that the effects starts to show. At 100 kyr the material has become overall more molecular, with the small exception for the highest ζ CR ∼ 10 −15 s −1 for which molecular dissociation has already begun dominating the chemistry. At 1 Myr the trend from 100 kyr continues, and from here on the molecular dissociation stepwise begins dominating with time, which results in a monotonic increase in the ratio with ζ CR . As such, for old objects 1 Myr, this could be a very powerful probe. It requires a detailed examination of the abundances of atoms versus molecules, which we leave for a future study. The effects of molecular dissociation are much clearer when we consider the abundance ratio between ions and neutral species, shown in the right panel of Figure 4. Charged species can carry either a positive (+) or negative (-) charge, and we calculate the mass-weighted ratio as The distribution clearly shows that ions are sensitive to an increase in ζ CR , as we see a shallow but clear trend of the ratio increasing with ζ CR . But it is the large increase at 1 Myr for ζ CR 3 × 10 −16 s −1 that is the strongest sign of the molecular dissociation, where the number of atomic ions is greatly increased. Aside from that, we see that the ratios decrease overall with time, meaning that the material is, with time, slowly becoming more neutral: part of the atomic material is reforming neutral molecules. Though the trend increases for old objects at the highest ζ CR ∼ 10 −15 s −1 and the trends are monotonic at every time step, the variance is very small. Again, this should be investigated by future studies.
We will also look at the effects of an increase in both density and time. Because temperature is directly linked to ζ CR through the calculated cooling and heating rates, we do not investigate this parameter separately.
In Figure 5 we plot abundances for a variety of commonly observed molecules as a function of time for three different densities (rows; 10 4 , 10 6 and 10 8 cm −3 ) and two different CR ionization rates (columns; 10 −17 s −1 and 10 −15 s −1 ).
Increasing the density has the strongest effect on the chemical time scale. In Figure 5 we see the abundances increasing, in some cases monotonically, until they reach a peak value. This peak moves to earlier times with increasing density. After this point, the chemistry becomes more complex. Compared to previously only being dominated by ion-neutral reactions, surface chemistry and neutral-neutral reactions begin to play a significant role. An increase in small ions, in particular atomic ions, also affects the chemical equilibrium values.  With time, often after approximately 1 or 10 Myr, depending on the ionization rate, the system reaches an equilibrium state. The equilibrium state of the species' abundances largely depends on the adsorption of species onto grain surfaces, which becomes more efficient with higher densities and ζ CR . The effect from an increase in ζ CR is indirect, as it also increases the temperature, and thus the amount of material desorbed into the gas phase for some species, but also increases the amount of ions, which has significant effects on the chemistry.
In comparison, ζ CR can have very significant influence on the chemistry, in particular on the steadystate abundances. But density has a similar effect between low and high ζ CR . That means that density only competes with time to effect the chemistry, and its effect is much less significant. We conclude that density is not as important as time for the chemical evolution, and in order to simplify the analysis, we exclude it from the study by keeping it constant at n H = 10 4 cm −3 .

INDIVIDUAL TRACERS
Here, we look at the overall trends in abundances for the species listed in Table 1. Species must fulfill a number of conditions to be useful as probes of the CR ionization rate. First, molecules must be sufficiently abundant to be observable; in our case we consider a relative abundance of 10 −15 to be the limit of observability. Second, ideal probes must show a monotonous increase or decrease of their abundances versus ζ CR . Third, these abundance variations as a function of ζ CR are ideally large so that differences can easily be observed. While the species that we consider in Table 1 have been proven to be observable, and hence fulfill the first criterion, we will take into account the two other criteria.

Probe score
The usefulness of an individual or pair of species will be discretized to a probe score, ranging between 0 and 1, based on the above criteria. The score is defined by the multiplication of two factors; the fractional range in ζ CR where the species can be used as a probe, and the logarithm of the maximum abundance variations. By multiplying the two factors, which we have named as the f unctional and the variation factors, respectively, we calculate the probe score.
The functional factor will be calculated based on the most dominant trend between positive and negative trends. In order to do this, we need to determine trends in the data and as such, need to define what constitutes a significant variation. All CR ionization events will not result in the production of the target species, represented by for the case of OH + in Section 3.2. We consider a conservative effectiveness of 50%. This means that for an increase in ζ CR by one order of magnitude, abundance should increase by at least a factor of 5, and anything less signifies a plateau. After plateaus have been identified we use Theil-Sen regression to determine the remaining trends in the distribution. This is a robust regression method from the python package scikit-learn that ignores outliers.
A plateau limits the functional range of a probe, but a repeated switch between positive and negative trends will also have a great impact on how useful a probe is. If it switches once, we get a bow-shaped distribution, and as long as we can clearly distinguish between the positive and negative sides of the distribution, the species is still useful as a probe. However, if the trend changes more than once, it severely reduces how useful the species is as a probe, and in these cases we therefore halve the probe score.
The variation factor in the probe score is defined as the logarithm of the ratio between the maximum and the minimum abundance, R = log 10 (max[X]/min[X]), for species X. Abundance variations will be restricted, as variations beyond a certain factor will not have a significant impact. In this study we consider this limit to be a factor of 100, based on our observational experience. We still want to be able to separate and rank species that have extremely large abundance variation. For this reason the variation factor will be linear with ζ CR up to a factor 0.1 below the maximum allowed variance.
After that, we adopt a relation ∝ x/(1 + x), where x is the logarithm of the abundance variation.
Finally, the variation factor is normalized.

Survey of observable species
We start with an overview of our selection of species. In Figure 6 the abundances are shown separately for 0.01, 0.1, and 1 Myr (top to bottom) in so-called "violin plots". Abundances are shown relative to their median values and for comparison the mean abundance is indicated with a green line.
The median values of each species are printed at the bottom of each panel. The gray shapes show the probability density function, a non-parametric estimate of the probability distribution, determining the probability of specific abundances to be found at certain abundance ranges.
The probe scores of individual species are shown at the bottom of each panel. They are color-coded based on the four quartiles of the distribution of calculated probe scores for all species in Table 1    The species in Figure 6 are arranged such that ions are shown first, followed by neutral molecules, both with an increasing number of atoms. Though ions are strongly connected to ζ CR , their probe scores are telling a different story, with many examples of having low probe scores. No ions show a persistently high score. This is due to a mix of small abundance variations and sensitivity to the increasing electron density that follows the increase in ζ CR . N 2 H + and HOC + both stand out as great probes at times ¡1 Myr, and to some lesser extent also HCO + , but at 1 Myr C + is the only useful probe.
Amongst the neutrals, the trends vary a lot, and we find good probe scores for all different types and sizes of species. Amongst the atoms, C stands out as an excellent probe at 1 Myr, just like C + , but at earlier times its probe scores also suffer from small abundance variations. Amongst diatomic species, SO is the sole species that shows promise, and in fact at all three time steps. The triatomic species all have below average scores, but amongst the more complex species we find several probes, and amongst these, HNCO and C 3 H 4 are the most promising.

Selected tracers
As such we have identified eight promising probe candidates that we will take a closer look at by dissecting their chemistry and evaluating their fit as probes of ζ CR . In Figure 7 we illustrate variations in abundances for these selected species and in the following section individually discuss their chemistry. In Table 2 we summarize the regions where the individual species are useful as probes for ζ CR .
We compare our model results with the observations of the starless core L1498 by Tafalla et al. (2006), in order to ascertain the accuracy of our numerical models. Tafalla et al. simultaneously observed N 2 H + , HCO + and SO from our list of species. We have marked the observed abundances in Figure 7 with ×, adopting ζ CR = 3 × 10 −16 s −1 derived by Caselli et al. (1998). The comparison points to a chemical age of ∼10 kyr or 1 Myr. These very different ages illustrate the uncertainties involved when only using a small number of species. The latter is closer to the derived lower age limit by Maret et al. (2013) of a few hundred-thousand years. Discrepancies are likely due to Maret et al. adopting the density structure from Tafalla et al. (2004) instead of a single density point as in this paper.
Atomic ionized carbon (C + ) does not increase strongly enough at 10 kyr, and at high ζ CR it ends with a plateau. At 100 kyr molecular dissociation starts to affect the abundances, resulting in a wavy pattern during this transitional period. It is not until ∼ 1 Myr that this ion becomes a promising probe. But once again, due to molecular dissociation, its abundance variation is approximately five orders of magnitude.
N 2 H + is linked to ζ CR through its formation from H + 3 reacting with N 2 . At late chemical ages (≈ 1 Myr) it slowly reaches a peak abundance at ζ CR ∼ 10 −16 s −1 , after which the abundance quickly drops several orders of magnitude due to molecular dissociation. This limits its usefulness as a probe at late times, which we see reflected in its probing score, while at earlier times it is a great probe with monotonic and large increasing trends.
HCO + is a commonly used probe for ζ CR . But as we see in the second row of Figure 7, its isomer HOC + proves to be a much more consistent probe. At earlier times < 1 Myr, both HCO + and HOC + show a strong increase with ζ CR . It is at 1 Myr that the two ions differ. While HCO + suffers greatly from molecular dissociation, and at higher ζ CR there is also an increase in o-H 2 , which converts HCO + to HOC + . As such, HCO + has a plateau until ζ CR ∼ 3 × 10 −16 s −1 when the dissociation quickly drops abundances. Meanwhile HOC + benefits from a small but monotonic increase in abundance and only suffers from a plateau at ζ CR ∼ 10 −15 s −1 . As such, HOC + is clearly a more capable ζ CR probe.
Note that the abundance of HOC + is lower than HCO + by a factor of 10-100 depending on the time and ζ CR . This, of course, makes observations of HOC + more challenging than observations of HCO + .
Atomic carbon (C) shows a wide range of behaviors over the three time steps. At 10 kyr it shows a slow decrease in abundance, so slow that it is considered a plateau up until ζ CR ≈ 10 −16 s −1 . The decrease in the C abundance is much more significant at 100 kyr, dropping approximately three orders of magnitude until ζ CR 10 −16 s −1 , after which molecular dissociation slowly begins to dominate and increases the C abundance. This results in a much higher probe score, but because of the arc-shape, it is a functional probe over a very limited range in ζ CR . At 1 Myr the molecular dissociation dominates the chemistry of atomic carbon, resulting in a monotonic positive trend. Compared to C + , the two are great probes at 1 Myr, but atomic carbon is also a useful probe at 100 kyr.
SO has a clear positive trend at 10 kyr, but at 100 kyr it suffers from a plateau at ζ CR 10 −16 s −1 , again due to molecular dissociation. At 1 Myr this trend also spreads to lower ζ CR , but the abundance is monotonic aside from a brief plateau at low ζ CR 10 −17 s −1 . As such, SO is showing promise at all time steps, aside from a restricted functional range at 100 kyr.
HNCO is a very peculiar case. At 10 kyr it is slightly restricted by a small plateau, but at 100 kyr it is restricted to only being useful up to ζ CR ∼ 3 × 10 −17 s −1 . The peculiarity arises at 1 Myr, when its abundance drops by almost eight orders of magnitude. This feature alone makes it interesting as a probe, in particular for more distant or diffuse targets. Reactions with H + and He + split the molecule into smaller parts, which are not so easy to reform. Considering that uncertainties in observational measurements increase for distant targets, HNCO shows great promise due to its large variations for for chemical ages 1 Myr. The large variations means it is not as sensitive to these errors as other probes and we can get more precise estimates of ζ CR .
The final example is C 3 H 4 , another molecule primarily formed on the surface from C 3 being hydrogenated. At 10 kyr and 100 kyr, we see a sharp increase in the abundance at ζ CR ∼ 5 × 10 −17 s −1 that is most prominent at 10 kyr. This is caused by the sudden increased inflows of carbon and hydrogen atoms on the dust surfaces and by the dissociation of H 2 and longer carbon chains C m for m> 2. Despite this bump, C 3 H 4 shows an almost monotonic increase at 10 kyr, while at 100 kyr it is much more restricted at ζ CR 10 −16 s −1 due to reactions with molecular ions such as H 3 O + , H + 3 and HCO + . Finally, at 1 Myr the molecule shows fluctuations in its abundance, similar to a sinus curve, which is due to destruction from reacting with molecular ions at low ζ CR , while at ζ CR 5 × 10 −17 s −1 it begins to reform into C 4 H 2 instead thanks to the strong increase in atomic carbon. Compared to HNCO, C 3 H 4 reforms much more efficiently. This is because it reacts with ions produce molecular ions such as C 3 H + m , with m ≥ 2, which can quickly reform into C 3 H 4 . However, in this case it is a drawback, keeping its abundance steady and the molecule not useful as a probe for older objects 1 Myr.

Species
Probe area 10 kyr 100 kyr 1 Myr Table 2. Listing of the regions where selected individual species show the most promise as probes of ζ CR .
The color-coding specifies how restricted the functional range of the probe is: none (green), small (yellow), moderate (orange), and severe (red). The symbols specify the nature of the abundance variation: ticks ( ) indicate only one monotonic trend, tildes (∼) signify some level of restriction, and bars (|) denote a symmetric arc-shaped distribution centered near ∼ 10 −16 s −1 .

SPECIES PAIRS
If we look at pairs of species as probes, it is possible that they can remove restrictions or troublesome behaviors that individual species might have. Because forming pairs based on the selected species in Table 1 will create a large number of pairs, we will focus our discussion here on the species pairs with probe scores from the top quartile (scores ≥ 0.66) in at least two of the time steps. These are listed in Table 3, sorted by their mean score. From the probe scores we can immediately see that the species pairs perform much better than than the individual species. Many of these pairs have their lowest score at 100 kyr, but there are also many examples of species that have their lowest probe score at 1 Myr. Only a few exceptions have problems at 10 kyr: NH 2 D/HC 3 N, SO/C 4 H and C + /O, but they perform much better at the other time steps.
From the more than 200 species pairs generated, only 44 have at least 2 green probe scores amongst the 3 time steps. Not surprisingly, the majority of these pairs include an ion (22/44). The most commonly listed ion is N 2 H + (12/44), which follows from it being one of the most promising individual probes. After this, HCO + and HOC + (4/44 and 5/44, respectively) are significantly less common.
Amongst the neutrals, HNCO (12/44) is most common, mostly owing to its large abundance variation.  Table 3. Listing of species pairs in this study that show, in at least two out of the three time steps, a high probing score corresponding to the top quartile of the total distribution. The species pairs are sorted by their mean probing score. The green scores represent the top quartiles with scores ≥ 0.66, the yellow scores represent the second quartiles with scores ≥ 0.53, the orange scores represent the third quartile with scores ≥ 0.31, and the red scores represent any non-zero values.
It is closely followed by C (8/44) and HC 3 N (6/44). With the exception of HC 3 N, these belongs to the best performing individual species that we identified in Section 5.3.
As an individual tracer, the mean score of HC 3 N is pulled down because it is not affected by the increasing ζ CR at 10 kyr. Otherwise HC 3 N performs well with green scores at other time steps. When paired with other species that are strongly affected by ζ CR at 10 kyr, we get a much more promising probe.
Looking at the size of the species pairs, it appears that small differences of 1 − 2 atoms result in the best probes. Looking at differences in atomic weight between the two species, there is a grouping around 13−17 in atomic mass, identified as the difference between one heavier element (C, N, O), and possibly also a lighter element (H, D). These small differences suggest that the formation process for either partner of a species pair should have a length difference of one, or at most two, steps.
However, the two individual species in the pairs often are formed from different formation pathways.
The pairs N 2 H + /C 4 H and N 2 H + /C 2 H, are excellent examples of this, where one species belongs to nitrogen chemistry, and the other to carbon chemistry. In short: similarity is preferential, but they should not be too similar.
In Figure 8 we plot the abundances of the top 10 scoring species pairs from The remaining pairs are limited by arc-shaped distributions or extensive plateaus at 100 kyr, and either have a monotonic trend, or small limitations by plateaus or an arc-shaped distribution skewed Species Probe area 10 kyr 100 kyr 1 Myr Table 4. Listing of the regions of the most promising species pairs based on their average probe scores.
The color-coding specifies how restricted the functional range of the probe is: no (green), small (yellow), moderate (orange), and severe (red). The symbols specify the nature of the abundance variation: ticks ( ) indicate only one monotonic trend, tildes (∼) signify some level of restriction, and bars (|) represent an arc-shaped distribution.
toward the high end of the ζ CR range. Even with these limitations, they are capable as probes, as long as one knows enough about the targets to avoid the problematic time scales and areas in the ζ CR range.
We can clearly see that there are promising probes amongst these species pairs, and in Table 4 we list a summary of how well the most promising species pairs perform. Given that our species selection in Table 1 is largely driven by observations, most of the abundance ratios listed in Table 4 can be observed using (sub)millimeter-wavelength telescopes. Species like C 4 H, C 3 H 2 , and HOC + might require significant integration times, but they are generally observable. The situation is different for ratios involving atoms, though. Observations of C at 492 GHz from the ground are only possible in exceptionally high and dry sites.

CONCLUSIONS
We have looked at the potential of individual and pairs of species as probes for ζ CR in dense clouds.
The study is limited to species commonly observed in large-scale surveys and for ζ CR = 5 × 10 −18 − 1×10 −15 s −1 at the time steps 0.01, 0.1, and 1.0 Myr. Our study is a pilot intended to be expanded to include a larger sample of species to prepare for the new generation of telescopes and future surveys.
Ions are strongly linked to an increase in ζ CR , but molecular ions suffer from molecular dissociation at high ζ CR 10 −16 s −1 and 1 Myr. This means that primarily atoms are increasing in abundance beyond this point, and as we look at later time steps, the molecular dissociation is dominating at even lower ζ CR . Atomic ions continuously increase with ζ CR , and even more so when the molecular dissociation is dominating. Because the abundance of some atomic ions, in particular H + , C + , He + , as well as those of some abundant molecular ions (e.g. H + 3 , HCO + , H 3 O + ), are very reactive, the atomic/molecular as well as ion/neutral ratios will affect the chemistry, and the former in particular could possibly also be used as a probe.
The probe candidates show a range of behaviors depending on their chemistry and in many cases it is the complexity of the species that determines how sensitive they will be to ζ CR . It controls the abundances of atomic and molecular ions as well as electrons, which are the driving components of the destructive pathways. Large, complex molecules with long formation pathways will decrease or even dilute the effects of the steadily increasing ζ CR , while abundances are low because reforming them takes a long time.
In order to infer how great an individual or pair of species can probe ζ CR , we identify different trends in the abundance variations and use them to calculate a probe score. It is based on the multiplication of two factors: one relating to the f unctional range of the probe, as in how much of the ζ CR range it can be used as a probe, and one factor relating to the magnitude of abundance variation that the species experience over the studied ζ CR range. We find that all probe candidates investigated in this study have limitations. That is, none of the individual or pairs of species are perfect probes for the complete investigated ζ CR range and at all three time steps, but a few come very close with only small limitations.
Amongst individual species, the ions N 2 H + , HCO + and HOC + are the most capable probes. HNCO is a particularly interesting molecule, as it shows an incredibly large range of abundance variation at 1 Myr of almost eight orders of magnitude. This sensitivity is caused by its formation being slow and largely associated with surface chemistry. This may make it a useful probe for more distant targets, as it will not be as sensitive to large errors that can be expected for distant objects.
Species pairs show more promise. Most of them show a monotonic trend with an increasing ζ CR at 10 kyr and 1 Myr, while many suffer from restrictions during the transitional period at 100 kyr, where many trends are reversed as molecular dissociation begins to affect the chemistry.
Amongst the pairs, there are a few that stand out. N 2 H + /C 4 H only shows restrictions at 100 kyr for high ζ CR 3 × 10 −16 s −1 and at 1 Myr for low ζ CR 8 × 10 −18 s −1 . N 2 H + /C 2 H on the other hand has no restrictions at 1 Myr, but is more limited at 100 kyr to ζ CR 10 −16 s −1 . N 2 H + /C 2 H, however, on the other hand has no restriction at 1 Myr, but a larger restriction at 100 kyr due to the peak of the arc being shifted to lower ζ CR .
Many species pairs are restricted or not very useful as probes at 100 kyr due to low variations or the trends not being monotonic. But there is an equal number of species pairs that have their lowest probe score at 1 Myr. However, the highest scoring species pairs are almost exclusively belong to the former type.
A number of pairs with HNCO as a partner have the very particular property that their abundance varies by almost eight orders of magnitude at 1 Myr. This can be particularly useful in studies of distant targets, where larger abundance variations could improve the accuracy of the estimated ζ CR .
We have summarized the functional ranges of individual species and species pairs, so that future surveys and observations can target the right probes for the right environments. While this study focused on dense clouds, at least some of these probes are likely also applicable in diffuse clouds and other astrochemical environments. This can be investigated in a future study, and the current one could also be extended to include more species to prepare for the next generation of observational facilities.