Abundance ratios of OH/CO and HCO+/CO as probes of the cosmic ray ionization rate in diffuse clouds

The cosmic-ray ionization rate (CRIR, $\zeta_2$) is one of the key parameters controlling the formation and destruction of various molecules in molecular clouds. However, the current most commonly used CRIR tracers, such as H$_3^+$, OH$^+$, and H$_2$O$^+$, are hard to detect and require the presence of background massive stars for absorption measurements. In this work, we propose an alternative method to infer the CRIR in diffuse clouds using the abundance ratios of OH/CO and HCO$^+$/CO. We have analyzed the response of chemical abundances of CO, OH, and HCO$^+$ on various environmental parameters of the interstellar medium in diffuse clouds and found that their abundances are proportional to $\zeta_2$. Our analytic expressions give an excellent calculation of the abundance of OH for $\zeta_2$ $\leq$10$^{-15}$ s$^{-1}$, which are potentially useful for modelling chemistry in hydrodynamical simulations. The abundances of OH and HCO$^+$ were found to monotonically decrease with increasing density, while the CO abundance shows the opposite trend. With high-sensitivity absorption transitions of both CO (1--0) and (2--1) lines from ALMA, we have derived the H$_2$ number densities ($n_{\rm H_2}$) toward 4 line-of-sights (LOSs); assuming a kinetic temperature of $T_{\rm k}=50\,{\rm K}$, we find a range of (0.14$\pm$0.03--1.2$\pm$0.1)$\times$10$^2$ cm$^{-3}$}. By comparing the observed and modelled HCO$^+$/CO ratios, we find that $\zeta_2$ in our diffuse gas sample is in the { range of $1.0_{-1.0}^{+14.8}$ $\times$10$^{-16}- 2.5_{-2.4}^{+1.4}$ $\times$10$^{-15}$ s$^{-1}$. This is $\sim$2 times higher than the average value measured at higher extinction, supporting an attenuation of CRs as suggested by theoretical models.

While ζ 2 cannot be easily observed directly, the use of tracers is favored instead. H + 3 is one of the most commonly used tracers of the CRIR. It is produced by the CR ionization of H 2 and it is destroyed through reactions with abundant neutral species (e.g., CO, O) and electrons Geballe et al. 1999;Dalgarno 2006;Indriolo & McCall 2012). The CRIR can be derived once the gas temperature, volume density (n H ), and column densities (e.g. H 2 , CO) are known. Other molecules that are directly relevant to the H + 3 chemistry are considered as potential probes of the CRIR, such as HCO + and DCO + (Guelin et al. 1982;van der Tak & van Dishoeck 2000;Caselli et al. 1998), OH + (Indriolo et al. 2015; Bacalla et al. 2019), and H 2 O + Neufeld & Wolfire 2017;Bialy et al. 2019). Measuring ζ 2 with ions such as H + 3 , OH + , and H 2 O + , requires the presence of bright massive stars in the background, which is not very common. Furthermore, the deuterium species can only be detected in high extinction regions due to its relatively low abundance. Although the above methodology can provide reasonable measurements, it is difficult and somehow impossible for general use.
However, oxygen-bearing molecules have -in principlethe potential to constrain the CRIR; most of the formation of oxygen-bearing species (e.g., OH, HCO + ) starts from the hydrogenation of ionized oxygen (O + ). Since the ionization potential of atomic oxygen (13.62 eV) is very close to that of atomic hydrogen (13.6 eV), the majority of oxygen in the cold neutral medium (CNM) is in its atomic form. Thus, in such FUV-shielded regions, O + and oxygen-bearing species are the indirect products of CRs.
OH has long been proposed to be an alternative tracer of molecular gas due to its fairly constant abundance in the ISM (Liszt & Lucas 1998, 2002Li et al. 2018).
The thermal emission and absorp-tion lines of OH 18 cm have been detected extensively toward the CNM, which is extended to the outskirts of the molecular clouds and where the CO emission is faint or undetectable (Turner 1979;Magnani et al. 1988;Wannier et al. 1993;Cotten et al. 2012;Li et al. 2018;Busch et al. 2021). HCO + , in the traditional sense, is a dense gas tracer in the ISM. However, recent observations found that HCO + is ubiquitous in diffuse and translucent clouds (Pety et al. 2017;Luo et al. 2020), especially with absorption measurements against strong continuum sources (e.g., quasars, H ii regions).
In this paper, we combine absorption observations of HCO + , and absorption and emission observations of CO to calculate their column densities along each LOS. We attempt to find the chemical connection between these oxygen-bearing molecules (CO, OH, and HCO + ) that are ubiquitous in diffuse clouds. We investigate the variance of chemical abundances of oxygenbearing molecules under different environmental parameters (e.g., gas volume density, FUV intensity, CRIR), especially the potential connection of their molecular abundance ratios as a probe of CRIR.
This paper is organized as follows. The observations and archival data used in this work are presented in Section 2. We present the results of column densities and constraints on the gas density in Section 3. In Section 4, we perform photodissociation region (PDR) modelling of CO, OH, and HCO + under various environmental parameters. We discuss the chemistry of CO, OH, and HCO + in the diffuse cloud and derive the chemical abundances in chemical equilibrium in Section 5. We derive the abundance ratios of OH/CO and HCO + /CO and constrain the CRIR by combining the observed HCO +to-CO abundance ratio and chemical models in the lowdensity (n H ∼ 10 2 cm −3 ) diffuse cloud in Section 6. The main results and conclusions are summarized in Section 7.
2. OBSERVATIONS AND DATA 2.1. HCO + (1-0) and CO (1-0) The HCO + (1-0) and CO (1-0) absorption observations toward 13 strong continuum sources were carried out during Apr. 2016 to May. 2016 with ALMA (project ID: 2015.1.00503.S, PI: L. Bronfman). The calibration of the raw data was performed using the Common Astronomy Software Applications (McMullin et al. 2007). Self-calibration was performed toward 2 sources (3C454.3 and 3C120) to increase the signal-to-noise ratio and eliminate the spectral contamination from bandpass calibrators. For each line-of-sight (LOS), we decompose the absorption spectra (τ ν ) into different Gaussian components to derive the column density at each velocity component. A detailed description of observations and data reduction can be found in Luo et al. (2020).
The HCO + (1-0) integrated optical depth (in units of km s −1 ) and CO (1-0) emission toward another 15 sources are taken from Table 1 in Liszt & Gerin (2023). The original data of HCO + absorption was presented in , Liszt & Lucas (2000), and Liszt & Gerin (2018). We exclude sources in which the brightness temperature of CO is bright (T mb ≥ 5 K) since we only focus on diffuse LOSs.

Reddening E(B-V) data
The E(B-V) data is taken from Green et al. (2019), in which the values were derived by combining stellar photometry from Pan-STARRS 1 and 2MASS, and parallaxes from Gaia. We use the E(B-V) values to derive the total gas column density along the LOS.

Constraints on the gas densities
For those components where both CO (1-0) and (2-1) transitions are available, we use the non-LTE radiative transfer code radex (van der Tak et al. 2007) to constrain the gas densities along the LOSs. We account for H 2 as the main colliding partner of CO. We perform several models covering H 2 volume densities of 10 1 ≤ n H2 ≤ 10 4 cm −3 , CO column densities of 10 11 ≤ N CO ≤ 10 16 cm −2 and kinetic temperatures of 10 1 ≤ T k ≤ 10 2 K. We then find the optimum solutions by maximizing the likelihood function using the Markov chain Monte Carlo (MCMC) method, which is encoded in emcee (Foreman-Mackey et al. 2013). The likelihood function is defined as: where the τ i obs and σ i obs are the optical depth and uncertainty of the observed i-th transition, respectively. τ i model is the modelled optical depth by radex.
The representative optimum results of n H2 under different T k are shown in the left panel of Fig. 1. As T k increases by a factor of 10, the resultant n H2 decreases by a factor of 2∼5. However, the column densities are less influenced by T k and vary by <4% (for CO) and <0.1% (for HCO + ) (with respect to the values at T k =50 K, as shown in the middle and right panels of Fig. 1.) Table 1 summarises the derived n H2 and n H of each velocity component at T k = 10, 50, and 100 K. The derived n H2 is in the range of (0.11±0.01-4.4±0.4)×10 2 cm −3 .

Calculation of column densities
For the rest of the components where only CO (1-0) and HCO + (1-0) are available, we consider an excitation temperature of T ex = 4 K for CO and T ex = 2.73 K for HCO + (e.g. Godard et al. 2010;Luo et al. 2020).
The column densities are calculated with (Mangum & Shirley 2015): where |µ lu | 2 is the dipole matrix element, Q rot is the rotational partition function, g u is the degeneracy of the upper energy level, and E u is the energy of the upper energy level. For each transition, the |µ lu | 2 , E u , and the rest frequency ν are taken from the CDMS database (Müller et al. 2001(Müller et al. , 2005. For linear molecules, the partition function is given by (McDowell 1987): where B 0 is the rigid rotor rotation constant. The calculated column densities of CO and HCO + are summarised in Table 2.

PHOTODISSOCIATION REGION MODELLING
Variations in chemical abundances can be used as a diagnostic tool for estimating the environmental parameters of the ISM. To better understand the response of the abundances of the molecular species we examine (CO, OH, HCO + ) and which are ubiquitously detected in diffuse and translucent clouds, we perform chemical simulations under a range of different environmental parameters (e.g. varying densities, FUV intensities, and cosmic-ray ionization rates). For these calculations, we use the publicly available code 3d-pdr 2 (Bisbas et al. 2012). Symbol "*" represents high confidence fit from 3d-pdr models (see Section 6.3).  Table 3. The gas phase elemental abundances relative to hydrogen in the PDR models.
Elements Abundance Elements Abundance In our modelling, we use a suite of one-dimensional uniform slabs with densities of n H = 10 m cm −3 where m = 2.0, 2.5, 3.0, 3.5, interacting with four FUV intensities χ/χ 0 =1,10,30,100 (where χ is normalized to the spectrum of Draine 1978), and four CRIR ζ 2 = (5, 10, 50, 100)×10 −17 s −1 . The FUV radiation field is considered as plane-parallel that impinges from one direction. The diffuse component of radiation is not accounted for. The maximum value of A V in our simulations is 20 mag. The gas-phase element abundances relative to hydrogen are shown in Table 3. We set all of the element carbon in the ionized phase with an abundance of 1.4×10 −4 and 60% of the total hydrogen in the molecular phase as the initial conditions in our PDR models (see Röllig et al. 2007, for further details). Throughout the text, the diffuse cloud is referred to as molecular gas with n H ≤500 cm −3 and A V ≤1 mag, the translucent cloud is referred to 500≤ n H ≤5000 cm −3 and A V ≤5 mag. Figure 2 shows our simulation results for a fixed ζ 2 = 10 −16 s −1 , in which the relative abundance of each species to H 2 is defined as the ratio of the corresponding column densities along the LOS. The variance of CO abundance is shown in the first row of Fig. 2. CO is mostly photodissociated, especially at lower extinction regions (e.g., A V <1 mag). As a consequence, the abundance of CO at low A V decreases by two orders of magnitude (from ∼ 10 −7 to ∼ 10 −9 ) as the FUV intensity increases from χ/χ 0 =1 to 100. The abundance of CO increases with increasing A V (e.g., 0.6≤ A V ≤5 mag) due to FUV shielding by dust and also self-shielding. However, the abundance of CO slightly decreases by a factor of a few from low to intermediate extinctions (A V ≤0.6 mag), which is due to the decrease of the abundance of its precursors (see Section 5). For a fixed FUV intensity, the abundance of CO increases with increasing density. The divergence of CO abundance at different densities is smaller (within a factor of a few, from 5×10 −8 to 5×10 −7 for χ/χ 0 = 1) at low extinction (A V <1 mag) but larger, over an order of magnitude, at high A V . At high-density or low FUV intensity models, the abundance of CO reaches the canonical value (e.g., ∼10 −4 as in the solar neighborhood) at lower extinction (A V = 2 mag for χ/χ 0 = 1 and n H = 10 3.5 cm −3 ) than that of low-density or high FUV intensity models (A V = 5 mag for χ/χ 0 = 10 and n H = 10 3.5 cm −3 ).
On the contrary, the abundance of OH has a different pattern to that of CO, as shown in the second row of Fig. 2. At low A V (e.g., 2 mag), the abundance of OH increases with the increasing FUV intensity, especially for high-density simulations (e.g., n H = 10 3.5 cm −3 ). It increases from ∼10 −9 to 10 −8 as the FUV intensity increases from χ/χ 0 =1 to 100. The abundance of OH decreases with the increasing A V when A V <2 mag. The slope of the decreasing trend of OH abundance with A V becomes larger for lower densities and higher FUV intensities. As can be seen from Fig. 2, the abundance of OH decreases by a factor of 2 for the high-density (n H = 10 3.5 cm −3 ) and low FUV intensity (χ/χ 0 =1) model, while it decreases by two orders of magnitude for the low-density (n H = 10 2 cm −3 ) and high FUV intensity (χ/χ 0 =100) model.
The knee point where the OH abundance trend turns, is approximately at an A V of 2 mag, after which it monotonically decreases with the increase of density below the A V value, while it converges above the A V value (except for n H = 10 3.5 cm −3 ). Particularly, the divergence of OH abundance at low FUV intensities is larger than that of at high FUV intensity. With the density increase from 10 2 cm −3 to 10 3 cm −3 , the abundance of OH decreases by nearly two orders of magnitude at χ/χ 0 =1, while it decreases only less than one order of magnitude at χ/χ 0 =100. Furthermore, the abundance of HCO + follows a similar trend to that of OH at low FUV intensity; HCO + decreases with the increase of density and A V (A V <2 mag). While this trend is maintained for the highdensity (e.g., n H >10 3 cm −3 ) and low FUV intensity models (e.g., χ/χ 0 <100), the abundance of HCO + declines at low extinction for low-density and high FUV intensity models.
At high FUV intensities, HCO + does not form as efficiently as OH does. The difference between HCO + and OH is most likely due to the interruption of the formation of HCO + through the OH channel at high FUV intensities (see Section 5). Figure 3 shows our simulations with a constant FUV intensity (χ/χ 0 =1). The abundances of CO, OH, and HCO + have a similar trend as ζ 2 varies: they increase with increasing CRIR at low extinction (e.g., A V < 1 mag for CO at n H =10 3 cm −3 ), while they have an inverse relation at high extinctions. The abundance of CO increases by an order of magnitude as the CRIR increase from 10 −16 to 10 −15 s −1 at low ex- . The chemical abundances of CO, OH, and HCO + as a function of AV using 3d-pdr for four densities (nH = 10 2 , 10 2.5 , 10 3 , and 10 3.5 cm −3 ) and four FUV intensities (χ/χ0 = 1, 10, 30, 100) for a fixed CRIR of ζ2 = 10 −16 s −1 .
tinctions for low-density models (e.g., n H ≤10 3 cm −3 ). The A V value of the "inverse point" where the dependence of CO abundance on the CRIR turns opposite, is higher (A V ∼5 mag) for low-density models and lower (A V ∼0.3 mag) for high-density models. For the n H = 10 3.5 cm −3 models, an order of magnitude increase in CRIR (from ζ 2 = 10 −16 s −1 to 10 −15 s −1 ) affects the CO abundance from a factor of ∼ 2 at low A V to 10 2 at higher A V .
As an important ISM environmental parameter, the CRIR has a high impact on the abundance of OH. In particular, it is almost proportional to ζ 2 at low A V . As seen in the second row of Fig. 3, the abundance of OH increases by an order of magnitude as ζ 2 increases from 10 −16 to 10 −15 s −1 at low A V (e.g., A V <1 mag) for all four densities explored. The A V value of the "inverse point" is larger than that of CO, which decreases with increasing density. At high extinction (A V >2 mag) and high-densities (n H >10 3 cm −3 ), the abundance of OH decreases with the increasing ζ 2 .
The dependence of HCO + abundance on the CRIR is similar to that of OH. With the CRIR increase by an order of magnitude, the abundance of HCO + increases by an order of magnitude at low A V for n H ≤10 3 cm −3 . However, at low-density models (n H = 10 2 cm −3 ), an increase of CRIR from 5×10 −16 s −1 to 10 −15 s −1 only results in a slight increase (≤30%) on the abundance of HCO + , meaning that the abundance of HCO + would hit a maximum value in this case (see Section 5.2 for more discussion). At high-densities (n H ∼10 3.5 cm −3 ), the increase of CRIR by an order of magnitude increases HCO + abundance by a factor of ∼5. The A V value of the "inverse point" is the same as that of OH.

The abundance of OH in diffuse cloud
In diffuse clouds, the formation of OH can be traced back from two channels that form OH + . The first one starts from the charge transfer reaction: with O + forming OH + through reaction: In addition to the reactions R1 and R2, atomic oxygen can directly react with H + 3 to form OH + (the second channel): Once OH + has formed, it can hydrogenate to form the precursors of OH, namely H 2 O + and H 3 O + : in which reaction R4 is the main destruction process for OH + . OH is formed through electron recombination reactions: and From reaction R11, H 2 O also contributes to the formation of OH by photodissociation: Reaction R12 is also the main destruction path of H 2 O. The minor destruction path of H 2 O is through C + in diffuse clouds: The destruction of OH in diffuse clouds includes photodissociation: reaction with C + : and reaction with H + : The OH abundance (x(OH)=n(OH)/n H ) in chemical equilibrium can be solved when the formation and destruction processes reach a balance (see detailed derivation in Appendix A). In diffuse clouds, photodissociation dominates the destruction of OH, and reactions R15 and R16 play a minor role. Since the abundance of C + is an order of magnitude higher than H + when the ζ 2 is not high (e.g., ζ 2 10 −15 s −1 , Le Petit et al. 2016b), we ignore here the contribution from R17 (the last term in equation A17). Thus, the abundance of OH can be written 3 as: where φ is the term in the square brackets of equation 4 and θ is in the form of equation A15. Since the term δ is small (∼0.1) in average diffuse ISM conditions (e.g., T =30 K, A V = 1 mag), φ ≈1. To a good approximation, equation 4 can be simplified as: . Since there is an anti-correlation between x(OH) and n H as density increases, the abundance of OH decreases monotonically ( similar to that shown in Fig. 2). On the other hand, if the UV radiation field increase or the extinction decreases, the gas temperature will increase. Therefore, all reaction rates in equation 5 will decrease and, therefore, the abundance of OH will also increase.
As can be seen from equation 5, the abundance of OH is proportional to ζ 2 in the diffuse cloud. Figure 4 shows the resultant abundance of OH (A V = 1 mag, n H = 10 2 cm −2 ) from both equations 5 as well as the which are not included in our analysis 4 . Equation 5 represents a good approximation ( with <60% uncertainty) for ζ 2 <10 −15 s −1 when compared with 3d-pdr models. However, in the PDR simulations, the abundance of OH reaches the maximum at ζ 2 ∼ 2×10 −15 s −1 . For higher ζ 2 , the OH abundance decreases. This is because the destruction of OH at high CRIR is dominated by H + (Reaction R17) rather than by photodissociation and C + (Reactions R14∼R16), which has been ignored in deriving equation 5. Similar work using isothermal simulations has been reported by Bisbas et al. (2017, their Fig. 10), in which the OH abundance peaks at ζ 2 /n H ≈ 2×10 −17 cm 3 s −1 .
Though there are limitations, we still highlight the usage of equation 5 at ζ 2 <10 −15 s −1 instead of running chemical networks, especially when coupled with hydrodynamical simulations, as this dramatically reduces the overall computational expense.

The abundance of HCO + in diffuse cloud
The formation of HCO + starts from ion-neutral reactions: 4 The reaction rate of O and H 2 is <<10 −30 cm 3 s −1 at a gas temperature below 100 K, which can be ignored. However, once the temperature increases above 200 K, the formation rate increases to 10 −20 cm 3 s −1 , which is comparable to that in equation A16.
where CO + is the result of reaction R15. Reaction R18 is very efficient as almost all CO + forms HCO +5 . The destruction of HCO + is always dominated by free electrons: Thus, in chemical equilibrium, the HCO + abundance (x (HCO + )) can be written as: By substituting equations A1-A17 to the above and relating x (H 2 O) with x (OH), we obtain: where is: For low values of ζ 2 (e.g., ζ 2 10 −15 s −1 ), the production of electrons results primarily from the ionization of atomic carbon. The electron abundance can be, thus, approximated with the C + abundance, e.g. x (e) = x (C + ) (Goldsmith 2001). If the ionization degree of the cloud is large or the gas density is low, x(O)/x(C + ) ≈ 2, and the abundance of CH is an order of magnitude lower than that of OH. Thus, the contribution from the last term of equation 7 is small. Equation 7 can, then, be simplified as: For higher ζ 2 (e.g., ζ 2 10 −15 s −1 ), a significant fraction of electrons are produced by H + (Le Petit et al. 2016b). In this case, equation 7 should be applied.
However, we should still keep in mind that at low extinctions and high FUV intensities, the destruction of CO + by H and free electrons becomes more efficient than that of H 2 (Reaction R18). The formation of HCO + is, therefore, interrupted. The response of HCO + abundance does not follow a similar trend to that of OH at low A V and high FUV intensities (as shown in Fig. 2). In such a case, equation 9 is not valid. Overall, equation 9 shows that x(HCO + ) is proportional to x(OH) in diffuse clouds. High-sensitivity absorption observations found that the integrated optical depth of HCO + has a tight linear relation with that of OH (N (HCO + )/N (OH) = 0.03, . Considering a diffuse cloud with gas temperature ranging from 50 to 100 K, at an extinction range of 0.5 − 1.5 mag, the resultant HCO + /OH ratio is nearly constant (∼3.5×10 −3 , see Fig. 5). Note that while this value is an order of magnitude lower than the above observations, the abundance of OH measured in diffuse gas has large uncertainty between different LOSs (see Section 5.4 for more discussion).
Additionally, most of the HCO + observations of  have high optical depths (τ >1). This means that the formation and destruction of HCO + may be different from what we considered here. Furthermore, the CH abundance in translucent clouds is likely to be comparable to that of OH (∼4×10 −8 , Liszt & Lucas 2002;Sheffer et al. 2008), meaning that the HCO + /OH ratio could be underestimated by equation 9.
Finally, in the case of a high CO abundance, the latter becomes the precursor of HCO + through reaction: In this case, the abundance of HCO + may be underestimated.

The abundance of CO in diffuse cloud
There are three major formation paths of CO: Reactions R16, R21, and The destruction of CO in the diffuse and translucent clouds is dominated by photodissociation: and by the interaction with He + CO + He + → C + + O + He.
The CO abundance (x (CO)) in chemical equilibrium can be written as: η CO denotes the destruction rate of CO by He + , which is proportional to the CRIR. If we substitute x(HCO + ) as x(OH) using equation 7, we obtain: The CO abundance is proportional to the gas density. In low-density diffuse gas where photodissociation dominates the destruction of CO, the OH channels dominate the formation of CO (Luo et al. 2023). Equation 11 can be simplified as : Since x (OH) is proportional to the CRIR, x (CO) increases with the increasing CRIR.

Comparison between the observed molecular abundances and model predictions
Radio emission line observations toward high latitude diffuse clouds (0.4 A V 1.1 mag) have reported an abundance of OH between 1.6×10 −7 and 4×10 −6 (Magnani et al. 1988). Absorption measurements toward quasars at radio wavelengths suggest a relative abundance ratio to total H column density (N OH /N H ) of 2.5 − 5 × 10 −8 at A V of 1 mag (Crutcher 1979;. Considering that the molecular fraction is 10%-20% at such LOSs , the abundance of OH with respect to H 2 is in the range of ∼10 −7 to 10 −6 . Recent radio observations by Tang et al. (2021) covering a broad range of A V (0.2−60 mag) suggest that the abundance of OH is higher (∼10 −6 ) at low A V and lower (∼10 −7 ) at higher A V . The high spatial resolution (∼0.12 pc) observations by  in the Taurus boundary also found a decreasing trend of OH abundance with the A V .
Our models show that the abundance of OH strongly depends on the density distribution and the ISM environmental parameters (see Section 4, Fig. 2 and 3) and can vary by more than two orders of magnitude (∼ 10 −9 − 10 −7 at the same A V ). As can be seen from equation 5, x(OH) is anti-correlated with k pd (R22) and n H . This is consistent with the OH survey in nearby molecular clouds of Tang et al. (2021). Considering the aforementioned large scatter of OH abundance in diffuse clouds, we adopt an OH abundance in the range of 10 −7 to 10 −6 as the "typical value". Thus, for any given FUV intensity (1≤ χ/χ 0 ≤100) and density (10 2 < n H ≤ 10 3.5 cm −3 ), the model underestimates the abundance of OH at an A V 0.2−2 mag if a ζ 2 = 10 −16 s −1 is adopted (Fig. 2). As seen from Fig. 3, the CRIR should be no less than 10 −16 s −1 if we are to reproduce the observed abundance of OH in diffuse clouds.
Due to the sub-thermal excitation of HCO + transitions in low-density gas (Godard et al. 2010;Luo et al. 2020), its abundance can vary over an order of magnitude without precise measurement of both the excitation temperature and optical depth (Luo et al. 2023). The abundance of HCO + in diffuse clouds has been measured frequently through absorption observations against strong continuum sources (e.g., quasars, H ii regions). Despite the different methods in obtaining the column density of H 2 , the abundance of HCO + is fairly constant in diffuse gas ((1.7-3.1)×10 −9 , Liszt & Gerin 2016;Gerin et al. 2019;Luo et al. 2020). At a CRIR of 10 −16 s −1 , our models underestimate the abundance of HCO + in diffuse gas in all density and FUV intensity ranges explored (Fig. 2). Therefore, to reproduce the observed abundance of HCO + , the gas density should be approximately in the range of 10 2 ≤ n H ≤ 10 2.5 cm −3 and the CRIR should be ζ 2 > 10 −16 s −1 (Fig. 3). The inferred density is consistent with the inferred density in Section 3.1 and those by radio and UV absorptions (e.g., 80−160 cm −3 , Goldsmith 2013; Liszt & Gerin 2016;Luo et al. 2020).
In diffuse clouds, most of the carbon is in the form of C + or C 0 due to insufficient shielding from UV photons. Similar to HCO + , the CO low-J transitions are usu-ally sub-thermally excited in the low-density diffuse gas (Goldsmith et al. 2008;Luo et al. 2020). The abundance of CO in diffuse clouds can vary by two orders of magnitude in different environments (e.g., 2.6×10 −8 -2×10 −5 from UV absorption measurements).
Considering all the above, we find that the abundances of OH, HCO + , and CO in diffuse clouds suggest an ISM environment with low gas densities (n H 10 2.5 cm −3 ) and high CRIRs (ζ 2 10 −16 s −1 ). A more quantitative analysis is discussed below ( §6.3).

THE ABUNDANCE RATIOS AND CRIR IN DIFFUSE CLOUDS
As can be seen from equations 5 and 9, the abundances of OH and HCO + are proportional to the CRIR in diffuse gas, which implies they can be used to constrain CRIR with a given density. However, molecular hydrogen does not emit radiation that can be observed from radio telescopes due to the lack of permanent dipole moment, the measurement of an accurate H 2 column density in the CNM, as well as its exact abundance, is very hard. Instead, it is possible to put constraints on the CRIR by combining the observed molecular column densities and their ratios with chemical models.

The abundance ratio of OH/CO
The abundance ratio of OH/CO can be obtained from equation 11: The second term in the square brackets can be safely ignored in the diffuse gas 7 : .
(14) Figure 6 shows the predicted abundance ratio of OH/CO with the column density of CO (N CO ) from 1D slab model simulations (n H = 10 2 cm −3 ). The abundance ratio of OH/CO decreases as the increasing N CO . The abundance ratio of OH/CO will increase with the increasing FUV intensity by a factor of a few at low N CO , while it remains approximately constant at high N CO . This is because, at a low FUV intensity (e.g., χ/χ 0 ∼ 1), the gas temperature is significantly lower than that of a moderate FUV intensity (e.g., χ/χ 0 ∼ 10). The decreasing gas temperature would increase the reaction rates of R15, R16, and R19, leading to a lower OH/CO ratio. However, at high N CO , CO molecules exist mainly in well-shielded regions, in which the gas temperature does not significantly increase even when high external FUV intensities exist.
As shown in Fig. 6, the abundance ratio of OH/CO monotonically increases with the increasing CRIR at a given N CO . Since high CRIR heats the gas, the reaction rates of R15, R16, and R19 will decrease.
At the meantime, η CO will increase as the increasing CRIR, leading to a higher OH/CO ratio.
The results shown in Fig. 6 indicate that once the gas density can be constrained from line ratios (e.g., C i/CO, CO rotational line ladder), the CRIR can be inferred from OH and CO observations without knowledge of their exact abundance relative to H 2 . This is potentially useful in high spatial resolution radio spectral line observations, especially when the column densities of H 2 cannot be determined.
However, the main challenge of generalizing the use of OH and CO as a probe of CRIR is that the excitation temperature of OH is usually within a few K above the Galactic synchrotron background (Li et al. 2018). This will lead to an order of magnitude uncertainty in N OH through emission lines if T ex varies between 0.1-1.0 K above T bg (note that N OH ∝ T ex /(T ex -T bg )). Absorption measurements are less suffered from the above issues (N OH ∝ T ex ), while the optical depth of OH is usually between a few to tens of 10 −3 (more than an order of magnitude lower than that of HCO + ). This affects the feasibility of detecting the absorption of OH in a relatively short integration time even with the 7 The term k pd (R24)+η CO is comparable to k R20 + k R23 at A V = 1 mag, while 1/n(H) is apparently a few orders of magnitude higher than x (CH) x (CO) x (O), thus, the second term can be ignored.
JVLA. With the proposed capability of high-sensitivity radio telescopes such as FAST, and SKA in the future (McClure-Griffiths et al. 2015), it is possible to constrain the CRIR through OH and CO in both the Milky Way and external galaxies.
6.2. The abundance ratio of HCO + /CO Contrary to OH, HCO + can be easily detected in diffuse clouds through absorption measurements against strong continuum sources (e.g., a few to tens of minutes with ALMA), and has less uncertainty in calculating its column density.
Following the same way, if we replace equation 9 with equation 14, the abundance ratio of HCO + /CO can be written as: (15) Figure 7 shows the predicted abundance ratio of HCO + /CO with N CO from 1D slab model simulations (n H = 10 2 cm −3 ), overlaid with the observed values from absorption measurements. Red dots are the observed values, and different line curves represent PDR models under different CRIRs. As seen from Fig. 7, the majority (22 out of the total 26) of the observational values are within 10 −16 ζ 2 10 −15 s −1 .
At a given N CO , the abundance ratio of HCO + /CO will first increase with increasing CRIR. This is similar to the variance of OH/CO as we have explained in Section 6.1. With the increasing photodissociation rate (R24) and the gas temperatures, reaction rates R15, R16, R19, and R23 decrease while the reaction rates of electron recombination (R21) increase. Thus, the increasing trend of HCO + /CO is slowing down because of the increasing destruction rate of HCO + by free electrons.
Consequently, there is a plateau of the HCO + /CO ratio (as well as the abundance of HCO + ), where the increasing CRIR no longer increases the HCO + /CO ratio.
However, the HCO + /CO ratio decreases at an even higher CRIR. This is because the most important precursors of HCO + (R18), such as CO + , are destroyed by atomic hydrogen and free electrons at high temperature, reducing the formation efficiency of HCO + . The abundance ratio of HCO + /CO reaches a maximum when the CRIR is approximately 10 −15 s −1 for n H = 10 2 cm −3 .
−18 to −14 with a step size of 0.07 dex and log 10 (n H /cm −3 ) = 1 -3.5 with a step size of 0.1 dex, interacting with three FUV intensities (χ/χ 0 = 1, 5, 10). We vary ζ 2 and n H to find the optimum model by minimizing the reduced χ 2 function: where f i obs and σ i obs are the observed molecular column densities and their corresponding uncertainties of the ith species, respectively. f i model is the column density from PDR models and N = 1 is the degree of freedom.
We show two representative χ 2 red distributions from our fitting in Figure 8. The higher values of χ 2 red are colored with bright blue and the lower colored with dark blue. White color indicates values beyond the maximum in this colour-bar. Red contours represent the boundary of the best-fit parameter where the deviation between modeled values and observations is 1σ. We treat a fitting result with a minimal value of χ 2 red 1 (upper panel of Figure 8) as a "high confidence" fit, and a fitting result with a minimal value of χ 2 red 1 (lower panel of Figure 8) as a "low confidence" fit.
As seen from Figure 8, the best-fit values drift to both higher ζ 2 and n H as the increasing FUV intensity. However, at high FUV intensity, the best-fit value of n H from PDR models is much higher than the allowable density range by MCMC runs (Table 1. On the other hand, increasing χ/χ 0 would increase T k . In that case, the inferred gas density should be lower (as seen in Figure 1). Thus, it is not likely that FUV intensity is high for our sources.
For T d = 15 K, the derived FUV intensity is χ/χ 0 = 1.4. A variation on T d of ±2 K would only result in a difference on χ/χ 0 by a factor of 2. In the following, we only consider the fitting results at χ/χ 0 = 1. The fitting results of all sources are summarised in Table 2, in which high confidence fit parameters are labelled with a "*" at the end of each row. The ζ 2 in our sample (high confidence results) is in the range of The n H is in the range of (0.1-3.2)×10 2 cm −3 , which is consistent with that obtained from MCMC runs in Section 3.1.
In Fig. 9, we plot the Luo et al. (2023) ζ 2 toward IC 348 measured from HCO + and CO measurements, and the Indriolo & McCall (2012) ζ 2 from H + 3 measurements. For comparison, we also include the Padovani et al. (2022) CR attenuation models L , H , and with low-energy spectral slope = −1.2. To plot ζ 2 as a function of E(B-V), we convert N H2 to E(B-V) by assuming that the atomic fraction f atomic follows the power law: (Luo et al. 2023), and the f atomic = 0.64 8 at A V ≤1 mag 9 . The average value of CRIR in our samples (red dots) is (8.0±6.4)×10 −16 s −1 . This is ∼ 2 times higher than the value measured with H + 3 toward nearby massive stars (ζ 2 = 3.5 +5.3 −3.0 ×10 −16 s −1 ) and toward IC 348 (4.7±1.5 ×10 −16 s −1 , see Fig. 9). Considering that our sightlines are located at much lower E(B-V) 10 (0.1 -0.7 mag), the values are reasonably consistent with the low E(B-V) portion by Indriolo & McCall (2012). Our results are consistent with the model H by Padovani et al. (2018), in which CRs (and consequently the CRIR) are attenuated as a function of A V (Padovani et al. 2018(Padovani et al. , 2022Gaches et al. 2022). However, model L , which corresponds to a "low" cosmic-ray spectrum based on Voyager-1 data, is found to underestimate the CRIR by almost an order of magnitude.

CONCLUSIONS
We analyze the abundances of CO, OH, and HCO + , and their abundance ratios in chemical equilibrium. We present a new approach to constrain the CRIR using these species. We calculated the column densities of HCO + and CO toward diffuse LOSs against quasars and compare them with 3d-pdr models. Our inferred values 2. Analyzing the chemical response of different molecules, we found that the abundance of CO increases with gas density and decreases with increasing FUV intensity, while the abundances of OH and HCO + mostly have an opposite trend to that of CO.
3. Our analytic expressions give an excellent abundance of OH when ζ 2 <10 −15 s −1 . This is potentially useful for hydrodynamical simulations as it reduces the computational expense of chemical networks. In the diffuse gas, the abundance of OH is proportional to the CRIR.
4. At a given N CO , the abundance ratio of OH/CO is monotonically increasing with the increase ζ 2 in the diffuse cloud, while the abundance ratio of HCO + /CO increases reaching a local maximum value at ζ 2 ≈ 10 −15 s −1 , before it decreases again. The downward trend of HCO + /CO at high ζ 2 is caused due to the increased destruction of the precursor of HCO + -CO + by electrons and atomic H at higher gas temperatures.
5. By comparing the observational values of HCO + /CO and chemical models, we find that the average ζ 2 in our sample is (8.0±6.4)×10 −15 s −1 . This value is ∼2 times higher than that of higher extinction regions, which is consistent with the hypothesis of decreasing ζ 2 as the increasing A V in theoretical studies.
With high-sensitivity measurements from HCO + and CO, we have demonstrated the possibility of using HCO + /CO ratios to constrain the CRIR in diffuse gas without knowing the exact molecular abundances relative to H 2 . We propose that the abundance ratios of OH/CO and HCO + /CO can be used to constrain the CRIR, especially with interferometry observations where high-resolution H 2 information is inaccessible. Due to the difficulties in obtaining accurate excitation and optical depth of OH with current facilities, we foresee that future instruments, such as SKA, will produce large samples of data sets for which our approach will be potentially useful. We thank Marco Padovani for the useful discussions on the CR attenuation models and for providing the latest data of these models. This work has been supported by the National Natural Science Foundation of China (grant No. 12041305), China Postdoctoral Science Foundation (grant No. 2021M691533), the Program for Innovative Talents, Entrepreneur in Jiangsu, the science research grants from the China Manned Space Project with NO.CMS-CSST-2021-A08. P. Z. acknowledges support from the National Natural Science Foundation of China (grant No. 12273010).  High confidence Low confidence High confidence in Luo et al. 2023 Low confidence in Luo et al. 2023 Figure 9. The red circles (high confidence fit) and black triangles (low confidence fit) show the inferred ζ2 in this work. Blue squares denote the measurements from H + 3 by Indriolo & McCall (2012). Brown crosses and gray stars are the high confidence fit and low confidence fit of ζ2 toward nearby star-forming cloud-IC 348. Black dashed, black solid, and black dash-dotted curves represent the polynomial fit to ζ2 at different E(B-V) for models L , H and with low-energy spectral slope α = −1.2 by Padovani et al. (2022).
The destruction of H 2 O is dominated by photodissociation and C + in the diffuse cloud (note that the destruction of H 2 O will be dominated by H + in dense cloud), thus, we have n(H 2 O) = n(H 2 O + )n(e)k R11 k pd (R12) + n(C + )k R13 .
In equilibrium, the formation and destruction of OH reach a balance (F[n(OH)] = D[n(OH)]).