Dissociative Recombination of Rotationally Cold OH⁺ and Its Implications for the Cosmic Ray Ionization Rate in Diffuse Clouds

Cosmic rays (CRs) are an important component of our Galaxy. In the local Galaxy, their total energy density (∼1 eV cm⁻³) is comparable to that of starlight, the magnetic field, and the cosmic microwave background (Grenier et al. 2015). CRs are the dominant ionization source for atomic and molecular hydrogen in the cold neutral medium (CNM). Through this ionization, they play a central role in initiating astrochemistry in the diffuse atomic and molecular and dense molecular phases of the interstellar medium (ISM; Hollenbach et al. 2012; Indriolo & McCall 2013; Gerin et al. 2016). In addition, the ionization provides coupling with the magnetic field, and in molecular cloud cores, slows down core collapse and inhibits both star and disk formation (McKee & Ostriker 2007; Padovani et al. 2018). Thus, the cosmic ray ionization rate (CRIR) is an important parameter for both chemical and dynamical models in these environments.

Observations have indicated that CR shielding decreases the ionization rate with increasing column density into a cloud. The rates estimated from low column density diffuse clouds (e.g., Neufeld & Wolfire 2017) are higher than those found in dense cloud interiors (Caselli et al. 1998; van der Tak & van Dishoeck 2000) by factors of ∼10. Here, we focus on diffuse clouds in order to obtain an ionization rate at cloud surfaces that are largely unshielded. The ionization rate at the cloud surface is one of the key parameters needed to model the astrochemistry in the ISM. In addition, comparing the surface rate to those estimated at larger columns acts as a guide to future researchers in understanding the details of CR shielding.

The CRIR in diffuse clouds is commonly inferred through observations of molecular cations such as ${{\rm{H}}}_{3}^{+}$, OH⁺, and ArH⁺ (Indriolo & McCall 2012; Gerin et al. 2016; Neufeld & Wolfire 2017). These are associated with specific local fractions of molecular hydrogen ${f}_{{{\rm{H}}}_{2}}\equiv 2n({{\rm{H}}}_{2})/{n}_{{\rm{H}}}$, where n(H₂) is the number density of molecular hydrogen, and n_H is the total hydrogen nuclei number density. ${{\rm{H}}}_{3}^{+}$ is found both in dense and diffuse molecular clouds (e.g., McCall et al. 2003; Brittain et al. 2004; Indriolo et al. 2007). OH⁺ probes the mostly atomic layer of diffuse clouds with ${f}_{{{\rm{H}}}_{2}}\sim 0.1$ (e.g., Indriolo et al. 2012). ArH⁺ is found in almost purely atomic diffuse clouds with ${f}_{{{\rm{H}}}_{2}}\lt 0.01$ (e.g., Schilke et al. 2014). Here, we focus on OH⁺ and the astrochemistry relevant to inferring the CRIR from the observed abundances.

The relative simplicity of OH⁺ chemistry in diffuse clouds makes it a powerful probe of the CRIR (Hollenbach et al. 2012). CR ionization of H forms H⁺, which can then undergo nearly resonant charge transfer with O to form O⁺. This can then be followed by an exoergic hydrogen abstraction reaction to form OH⁺ via

$$\begin{eqnarray}&&{{\rm{O}}}^{+}+{{\rm{H}}}_{2}\to {\mathrm{OH}}^{+}+{\rm{H}}.\end{eqnarray} \tag{ 1 }$$

The direct link of OH⁺ to the CR ionization of H makes its abundance roughly proportional to the CRIR over a wide range of parameter space. ⁷ OH⁺ can be destroyed by a hydrogen abstraction reaction, forming H₂O⁺ via

$$\begin{eqnarray}&&{\mathrm{OH}}^{+}+{{\rm{H}}}_{2}\to {{\rm{H}}}_{2}{{\rm{O}}}^{+}+{\rm{H}}\end{eqnarray} \tag{ 2 }$$

and also by dissociative recombination (DR) with free electrons via

$$\begin{eqnarray}&&{\mathrm{OH}}^{+}+{{\rm{e}}}^{-}\to {\rm{O}}+{\rm{H}}.\end{eqnarray} \tag{ 3 }$$

Recent laboratory studies have measured the hydrogen abstraction chain for diffuse cloud temperatures with an accuracy of 20% (Kovalenko et al. 2018; Kumar et al. 2018; Tran et al. 2018). Neufeld & Wolfire (2017) have identified the OH⁺ DR rate coefficient (i.e., kinetics) as one of the key uncertainties impacting the inferred CRIR using OH⁺ chemistry in diffuse clouds.

The reliability of the OH⁺ DR rate coefficient currently listed in astrochemical databases is a significant issue of concern. In diffuse clouds, OH⁺ is predicted to be in its lowest electronic, vibrational, and rotational level, neglecting its hyperfine structure. The KIDA database (Wakelam et al. 2012) has adopted the theoretical work of Guberman (1995) for the ground rotational level of OH⁺. A later theoretical study by Stroe & Fifirig (2018) expanded the calculations of Guberman (1995) by including core-excited Rydberg states into the existing theoretical framework and reproduced the previous results to within 20%. Quantum mechanical calculations are extremely challenging due to the many-body nature of the problem and the infinite number of intermediate states involved in the DR process. For example, the recent unified theoretical treatment of DR by Forer et al. (2023) for CH⁺ reproduces the experimental cross section results of Paul et al. (2022) only to within a factor of 2–5. Further theoretical improvements are required to reach experimental accuracy. As a result, laboratory measurements are still expected to be the most reliable means to generate DR kinetics data. The UMIST database (McElroy et al. 2013) and most astrochemical modelers follow this recommendation and have adopted the experimental results of Mitchell (1990), which are a factor of 6 larger than the calculations of Guberman (1995). However, the rate coefficient of Mitchell (1990) was derived from single-pass merged-beams measurements for which the OH⁺ ions were electronically, vibrationally, and rotationally excited. Amitay et al. (1996) studied the DR of electronically and vibrationally relaxed OH⁺ at the room-temperature Test Storage Ring (TSR). Their measurement, though, was only on a relative scale. Moreover, the rotational level population in the experiment was most likely close to room temperature. Recent experimental work for DR of HeH⁺ and CH⁺ has shown that for CNM temperatures the DR rate coefficient for internally excited ions can be over an order of magnitude larger or smaller than that for fully relaxed ions (Novotný et al. 2019; Paul et al. 2022). One aim of our work here is to generate DR kinetics data for internally cold OH⁺ to an accuracy of ∼20%, so that any remaining discrepancies between diffuse cloud chemical models and observations cannot be attributed to uncertainties in the underlying chemistry but rather begin to tell us about the astrophysics of diffuse clouds.

OH⁺ in diffuse clouds was first detected via absorption by rotational lines in the far-infrared (far-IR; Gerin et al. 2010; Neufeld et al. 2010; Wyrowski et al. 2010) and by electronic transitions in the near-ultraviolet (near-UV; Krełowski et al. 2010). Further near-UV observations found CRIR estimates similar to those derived from far-IR lines (Porras et al. 2013; Zhao et al. 2015). A subsequent IR survey by Indriolo et al. (2015) detected OH⁺, H₂O⁺, and H₃O⁺, providing one of the largest samples to date and demonstrating the usefulness of oxygen-bearing ions for probing the CRIR of hydrogen. The ratio of the OH⁺ to H₂O⁺ column densities was used to infer ${f}_{{{\rm{H}}}_{2}}$, an important parameter in determining the CRIR. Bacalla et al. (2019) extended and reevaluated the previous near-UV OH⁺ observations, implementing an update of the electronic-band line oscillator strengths and an updated chemical model. Their results suggest a somewhat higher CRIR compared to that from the far-IR studies. However, the near-UV studies lack corresponding H₂O⁺ observations, requiring those studies to assume a value for ${f}_{{{\rm{H}}}_{2}}$ in order to model the OH⁺ chemistry. In addition, both the far-IR and near-UV studies used a single-zone chemical model to interpret the OH⁺ observations, leading to the approximation that OH⁺ and H₂O⁺ exist in the cloud at the same location and set of physical conditions.

Some researchers have developed one-dimensional (1D) models as a function of visual extinction in diffuse clouds, so as to be able to more reliably interpret the simultaneous observation of multiple molecules within a diffuse cloud (Hollenbach et al. 2012; Neufeld & Wolfire 2016, 2017). These models predict that OH⁺, H₂O⁺, and H₃O⁺ peak in abundance at different depths into the cloud, contrary to the assumptions of single-zone models. However, a full comparison of the observed and modeled column densities requires reliable DR data for all involved species, most importantly for OH⁺ and H₂O⁺.

Here, we report our measurements for the DR rate coefficient of OH⁺ in its ground electronic, vibrational, and rotational level. Our results are applicable to the chemistry of OH⁺ in the ISM over a broad range of kinetic temperatures. In addition, we assess the impact of our data on models of diffuse cloud chemistry and the resulting implications for deriving the CRIR from astronomical observations. As for DR of H₂O⁺, the room-temperature storage ring experiments have measured DR of electronically and vibrationally relaxed H₂O⁺ ions, but the rotational level population was approximately room temperature (Jensen et al. 1999; Rosén et al. 2000). We are unaware of any DR measurements for rotationally relaxed H₂O⁺ and consider it a candidate for a future DR study.

DR measurements were performed on internally cold OH⁺ using the Cryogenic Storage Ring (CSR; von Hahn et al. 2016) facility at the Max Planck Institute for Nuclear Physics in Heidelberg, Germany. The methodology for performing DR measurements at CSR has been described in detail elsewhere (Novotný et al. 2019; Paul et al. 2022). Here, we provide only a brief overview, with an emphasis on those aspects that are most specific to the present results. Additional details can be found in the Appendices.

OH⁺ ions were generated in a gas-discharge source, producing internally excited ions that were injected into CSR (see Appendix A). The stored ions rapidly relaxed to their ground electronic and vibrational states. Rotationally, the ions relaxed via radiative interactions toward thermal equilibrium with the cryogenic blackbody temperatures of CSR (O'Connor et al. 2016; Meyer et al. 2017; Kálosi et al. 2022). The resulting rotational distribution was much colder than that for previous OH⁺ studies using room-temperature storage rings (Amitay et al. 1996; Strömholm et al. 1997; Hechtfischer et al. 2019). Additional rotational cooling was provided through electron–ion rotational-level-changing collisions using the electron beam described below. Using the collisional-radiative model of Kálosi et al. (2022) for the populations of the rotational levels, labeled by quantum number N, we estimate that ≈80% of the ions were in the N = 0 ground level, and ≈20% were in the N = 1 level (see Appendix C).

DR measurements were performed using an electron–ion merged-beams configuration that is discussed in more detail in Appendix A. We merged a magnetically guided electron beam onto the stored ions in one of the straight sections of CSR. The relative collision energy between the merged beams was controlled by tuning the nominal laboratory-frame electron-beam energy in the interaction region, E_e. DR events resulted in neutral fragments (see reaction outlined in Equation (3)), at nearly the same laboratory-frame velocity as the initial ions. The resulting neutral reaction products were collected by a particle-counting detector downstream from the interaction region. The recorded count rate is proportional to the merged-beams rate coefficient

$$\begin{eqnarray}&&{\alpha }^{\mathrm{mb}}=\langle \sigma v\rangle ,\end{eqnarray} \tag{ 4 }$$

where σ is the energy-dependent DR cross section, v is the electron–ion collision velocity, and the angled brackets designate the average over the velocity distribution in the electron–ion overlap region. The velocity distribution is determined by several experimental factors: the perpendicular and parallel temperature components T_⊥ and T_∥, respectively, relative to the bulk electron-velocity vector; the merging geometry; and the variable laboratory-frame energy along the interaction region. The resulting velocity distribution is modeled following the approach of Novotný et al. (2013), as extended by Kálosi et al. (2022). This velocity distribution differs significantly from the Maxwell–Boltzmann (i.e., kinetic-temperature-dependent) distribution used in astrochemical models. Extraction of the underlying cross section and subsequent conversion of the data is required before it can be applied in kinetic models (e.g., Paul et al. 2022). The velocity distribution can also be translated into a collision-energy distribution as a function of the detuning energy E_d, defined as the nominal center-of-mass collision energy,

$$\begin{eqnarray}&&{E}_{{\rm{d}}}={\left(\sqrt{{E}_{{\rm{e}}}}-\sqrt{{E}_{0}}\right)}^{2},\end{eqnarray} \tag{ 5 }$$

where E₀ is the laboratory-frame electron-beam energy at matched electron–ion velocities, taking into account that the reduced mass of the collision system is essentially equal to that of the electron. We use the full width at half maximum (FWHM) of the collision-energy distribution as a measure of the energy resolution in the experiment. The evaluation of the experimental data and the absolute scaling of the merged-beams results are discussed in Appendix B.

Measurements were performed using a nearly pure beam of OH⁺. To achieve this, we operated CSR in the recently developed isochronous mode, i.e., as a time-of-flight mass spectrometer with a mass resolution of ∼10⁻⁵ u (Grieser et al. 2022). A particular concern was potential contamination with NH₃ ⁺, which differs in mass from OH⁺ by ≈2 × 10⁻² u. For the present results, we find that <0.1% of the ion beam was due to NH₃ ⁺, using the methods developed by Grieser et al. (2022).

The measured merged-beams rate coefficient α^mb(E_d) for DR of OH⁺ is plotted in Figure 1, showing a rich structure of features. We can compare our results with those of Amitay et al. (1996), which were measured at the room-temperature TSR. At the top of Figure 1, we have plotted the energy-dependent energy resolution ΔE for the CSR and TSR results. At the lowest energies, the energy resolution for the CSR results is almost an order of magnitude finer than that for the TSR data. In order to qualitatively compare our work with the TSR results, we have scaled their relative values with a common factor at all energies to best match our data in the 1–3 eV region. Between ∼0.1 and 1 eV, we find reasonable agreement in the shape of the CSR and TSR DR results, taking into account the differences in the energy resolution of the two data sets. Below ∼0.1 eV, the CSR results show structure that is not seen in the TSR data. We attribute these differences to the improved energy resolution of the CSR data and to the well-defined and lower internal excitation of the ions, as discussed below. Above ≈5 eV, approximately corresponding to the dissociation energy of OH⁺ (Hechtfischer et al. 2019), the CSR and TSR results diverge. We attribute this to the opening of the dissociative excitation (DE) channels

$$\begin{eqnarray}&&{\mathrm{OH}}^{+}+{{\rm{e}}}^{-}\to \left\{\begin{array}{l}{{\rm{O}}}^{+}+{\rm{H}}+{{\rm{e}}}^{-}\\ {\rm{O}}+{{\rm{H}}}^{+}+{{\rm{e}}}^{-}.\end{array}\right.\end{eqnarray} \tag{ 6 }$$

The CSR DR detection system is currently not capable of distinguishing between DR and DE, while the TSR experiments used a mass-sensitive detection method that enabled Amitay et al. (1996) to select for DR and discriminate against DE events.

Zoom In Zoom Out Reset image size

The astrochemical relevance of our DR measurements is due to three advances over the TSR results, which we discuss in no particular order. First, the order of magnitude finer energy spread in the present experiment enables us to access collision energies relevant to diffuse cloud conditions. This improved energy resolution is most clearly demonstrated by the resonant feature observed below 20 meV. Second, we have measured α^mb(E_d) on an absolute scale, to be contrasted with the relative values from TSR, and have converted it into a kinetic-temperature-dependent rate coefficient, α^k(T_k), for use in chemical models (see Appendix E). Here, T_k is the kinetic temperature of the gas, which characterizes the velocity distribution for all particles. Last is our ability to generate internally cold OH⁺ ions, similar to the conditions expected in diffuse clouds. Our model for the level populations predicts an N = 0 relative population during the measurement of (80% ± 5%) and N = 1 of (19% ± 4%), where here and throughout all uncertainties are given at a one-sigma confidence level (see Appendix C). This contrasts with the room-temperature TSR results, where these two levels are predicted to have relative populations of 8% and 20%, respectively. In diffuse clouds, the OH⁺ is predicted to be nearly 100% in the N = 0 level. The critical density for collisions of OH⁺ with H and e⁻ for the N = 0 → N = 1 transition are calculated to be on the order of ≈10⁷ and ≈10⁴ cm⁻³, respectively, which are orders of magnitude larger than the estimated H and e⁻ densities in diffuse clouds. Our storage-time dependent DR measurements, similar to the studies for CH⁺ of Paul et al. (2022), enable us to probe the changing contributions of the OH⁺ N = 0, and N = 1 levels (see Appendix D). We find that the differences are most noticeable for E_d < 30 meV but are negligible when converting α^mb to α^k.

We have generated α^k(T_k) from our α^mb(E_d) results, following the DR cross-section-extraction method of Novotný et al. (2013) with further improvements from Paul et al. (2022). We have determined α^k from T_k = 10 to 20,000 K, as plotted in Figure 2. Simple fitting formulae for the rate coefficient are given in Appendix E. Also shown in the figure are the literature values from the single-pass merged-beams experiment of Mitchell (1990) and the theoretical calculation of Guberman (1995). These DR rate coefficients have been adopted by the astrochemistry databases UMIST (McElroy et al. 2013) and KIDA (Wakelam et al. 2012), respectively. Figure 2 also shows the results of Stroe & Fifirig (2018), which agree with that of Guberman (1995) to within 20%. At typical diffuse cloud temperatures of 40–130 K (Shull et al. 2021), our kinetic temperature rate coefficient is a factor of ∼5 times larger than the experimentally derived value of Mitchell (1990) and a factor of ∼33 times larger than the theoretical results of Guberman (1995).

Zoom In Zoom Out Reset image size

We also note that the kinetic temperature rate coefficients of Mitchell (1990) and Guberman (1995) exhibit a nearly ${T}_{{\rm{k}}}^{-1/2}$ behavior. However, our experimentally derived value has a slope that varies significantly with temperature between ${T}_{{\rm{k}}}^{-1/2}$ and ${T}_{{\rm{k}}}^{-1}$. A variable slope with temperature has also been seen by Novotný et al. (2019) for HeH⁺ and by Paul et al. (2022) for CH⁺. Together with those works, our findings demonstrate that the DR kinetic temperature rate coefficient can exhibit a temperature dependence that differs significantly from the theoretically derived ${T}_{{\rm{k}}}^{-1/2}$ behavior for direct DR of diatomic ions (Guberman 1995).

We have investigated the astrophysical impact of our DR rate coefficient α^k using the single-zone model of Bacalla et al. (2019) and an updated version of the 1D model of Neufeld & Wolfire (2017). These translate the observed OH⁺ column densities N(OH⁺) into an inferred CRIR. Both studies have adopted similar reaction networks and identified the same dominant destruction pathways for OH⁺ in diffuse clouds, i.e., reactions outlined in Equations (2) and (3).

To demonstrate the impact of the present DR rate coefficient on the inferred CRIR values, we have calculated the destruction rate of OH⁺ due to DR and hydrogen abstraction as a function of ${f}_{{{\rm{H}}}_{2}}$, similar to the methods of Indriolo et al. (2015) and Bacalla et al. (2019). The resulting inferred CRIR of atomic hydrogen ζ_H is proportional to the sum of destruction rates

$$\begin{eqnarray}&&{\zeta }_{{\rm{H}}}\propto \left({x}_{{\rm{e}}}{k}_{3}+\displaystyle \frac{{f}_{{{\rm{H}}}_{2}}}{2}{k}_{2}\right),\end{eqnarray} \tag{ 7 }$$

where x_e = n(e⁻)/n_H is the electron fraction, where n(e⁻) is the electron number density, and k₂ and k₃ are the rate coefficients for reactions outlined in Equations (2) and (3), respectively. For further details on Equation (7), see Indriolo et al. (2015). The ratio of the destruction rates

$$\begin{eqnarray}&&R=\displaystyle \frac{\left({x}_{{\rm{e}}}{\alpha }^{{\rm{k}}}+\tfrac{{f}_{{{\rm{H}}}_{2}}}{2}{k}_{2}\right)}{\left({x}_{{\rm{e}}}{k}_{3}^{\mathrm{old}}+\tfrac{{f}_{{{\rm{H}}}_{2}}}{2}{k}_{2}\right)},\end{eqnarray} \tag{ 8 }$$

calculated with the present (k₃ = α^k, given in Appendix E) and the previous (${k}_{3}={k}_{3}^{\mathrm{old}}$) DR rate coefficients, multiplicatively increases the CRIR values inferred by Bacalla et al. (2019). ⁸ The resulting factor is shown as a function of ${f}_{{{\rm{H}}}_{2}}$ in Figure 3, using the model parameters of Bacalla et al. (2019). They adopted the DR rate coefficient of Mitchell (1990), which we use in Equation (8) for ${k}_{3}^{\mathrm{old}}$. For the value of k₂, we use 1 × 10⁻⁹ from Jones et al. (1981), which has been incorporated into the UMIST database and adopted by Indriolo et al. (2015) and Bacalla et al. (2019). This value does not significantly differ from the results of Tran et al. (2018) or Kumar et al. (2018), and we use it to stay consistent with the aforementioned astronomical studies. We also show the multiplicative increase using the data of Guberman (1995) for ${k}_{3}^{\mathrm{old}}$.

Zoom In Zoom Out Reset image size

We find a factor of 1.6 increase for the inferred CRIR values of Bacalla et al. (2019), who assumed a single value of ${f}_{{{\rm{H}}}_{2}}\approx 0.17$. However, their ${f}_{{{\rm{H}}}_{2}}$ value is likely to be an overestimate when contrasted with the IR observations of Indriolo et al. (2015), who used the ratio of the OH⁺ to H₂O⁺ column densities to infer ${f}_{{{\rm{H}}}_{2}}$ and found ${f}_{{{\rm{H}}}_{2}}\lt 0.1$ for the majority of observations. The above quoted values of ${f}_{{{\rm{H}}}_{2}}$ represent the average for those parts of a cloud where OH⁺ is most abundant. We do not use values for ${f}_{{{\rm{H}}}_{2}}$ determined from observed H and H₂ column densities (e.g., Winkel et al. 2017) because those measurements are dominated by regions where the OH⁺ abundance is insignificant. The impact of our present results on the interpretation of the Bacalla et al. (2019) observations is even greater when values of ${f}_{{{\rm{H}}}_{2}}\lt 0.17$ are considered. Complementary observations of H₂O⁺ for the sightlines in their work do not yet exist but would enable one to better determine the appropriate value of ${f}_{{{\rm{H}}}_{2}}$ and thereby make a more reliable comparison between the near-UV and far-IR observations. We note that the ${f}_{{{\rm{H}}}_{2}}$ values inferred by Indriolo et al. (2015) directly depend on the DR rate coefficient of H₂O⁺. The uncertainty in this rate coefficient represents a multiplicative scaling factor with larger values leading to an increase in the inferred ${f}_{{{\rm{H}}}_{2}}$ and smaller values leading to a decrease.

Moving on to 1D diffuse cloud models, we find that there is a significant impact from our DR rate coefficient on the CRIR inferred for diffuse atomic clouds. We have redetermined the CRIR implied by the OH⁺ and H₂O⁺ observations presented by Indriolo et al. (2015). Our modeling results are compared directly to the reported column densities for OH⁺, H₂O⁺, and atomic H (N(OH⁺), $N({{{\rm{H}}}_{2}{\rm{O}}}^{+})$, and N(H), respectively), and are independent of the single-zone model and from the ${f}_{{{\rm{H}}}_{2}}$ values of Indriolo et al. (2015). To facilitate a comparison with observations of ${{\rm{H}}}_{3}^{+}$, we report below the primary CRIR per H atom, ζ_p (H). The inferred total CRIR ζ_H additionally includes secondary ionizations by energetic electrons and depends on the ionization fraction and ${f}_{{{\rm{H}}}_{2}}$ (Dalgarno et al. 1999). Here, the new OH⁺ DR rate coefficient obtained in the present study was incorporated into the Neufeld & Wolfire (2017) model, along with several other updates: these include the use of (1) new estimates of the rates for several photoprocesses and their depth dependence (Heays et al. 2017); (2) a revised treatment of the heating rate in diffuse clouds ⁹ ; and (3) new measurements of the rate coefficients (Kovalenko et al. 2018; Tran et al. 2018) for the reactions of O⁺, OH⁺ and H₂O⁺ with H₂. Interpreted with this updated version of the Neufeld & Wolfire (2017) grid of 1D models for diffuse clouds, the OH⁺ and H₂O⁺ column densities measured by Indriolo et al. (2015) imply an average of the logarithmic values ${\mathrm{log}}_{10}\,{\zeta }_{p}({\rm{H}})/[1\,{{\rm{s}}}^{-1}]\,=(-15.18\pm 0.06)$. This average value corresponds to the geometric mean ζ_p (H) = (6.6 ± 1.0) × 10⁻¹⁶ s⁻¹. The stated uncertainties are purely statistical in nature and represent the standard error on the mean. Our new estimate of the CRIR is a factor of 3 times larger than that obtained by Neufeld & Wolfire (2017). A factor of 2 increase is directly attributable to the larger DR rate coefficient obtained in the present study, in agreement with the expectations indicated in Figure 3, with the remaining factor of 1.5 being the result of all the other updates to the model.

The changes are represented graphically in Figure 4, where example model results are shown in the plane of the two observable quantities: the column density ratios N(OH⁺)/N(H₂O⁺) and N(OH⁺)/N(H). The diamond symbols with error bars show the values measured by Indriolo et al. (2015), while the red and blue curves show the model predictions as a function of ζ_p (H)/n₅₀, where n₅₀ = n_H/[50 cm⁻³], and the total visual extinction through the cloud, A_V(tot). Here, the blue curves are loci of constant ζ_p (H)/n₅₀ in units of 10⁻¹⁶ s⁻¹; and the red curves are loci of constant A_V(tot) in magnitudes. Figure 4(a) shows the present model predictions, while Figure 4(b) shows the Neufeld & Wolfire (2017) predictions (as shown previously in Neufeld & Wolfire 2017; Figure 7). The downward displacement of the blue curves in Figure 4(a) with respect to Figure 4(b) is the result of our modifications to the Neufeld & Wolfire (2017) model and is the explanation for why our revised estimate of the CRIR is 3 times as large.

Zoom In Zoom Out Reset image size

We note that our results depend on the H⁺ + O charge exchange rate coefficient, and there is some disagreement in the literature on the predicted value. The data used here are from Stancil et al. (1999), which are a factor of ∼6 larger than those of Spirko et al. (2003). Here, using the new OH⁺ DR rate coefficient and the Spirko et al. (2003) calculations, our model yields unphysical astronomical results with cloud columns much greater than A_V(tot) of 3 and an even larger CRIR than presently inferred. If the Spirko et al. (2003) calculations are correct, we would need to reevaluate our code and update our conclusions. Further theoretical and experimental charge exchange studies are highly desirable to help resolve the issue. Whatever the resolution, our improved OH⁺ DR rate coefficient will be critical for reliably determining the CRIR from observations of OH⁺ and H₂O⁺.

Our new estimate of the CRIR derived from OH⁺ and H₂O⁺ observations in diffuse atomic clouds is now 2.5 times larger than that inferred from an analysis of ${{\rm{H}}}_{3}^{+}$ observations in diffuse molecular clouds using the updated model, whereas the previous study of Neufeld & Wolfire (2017) concluded that the CRIR derived from OH⁺ and H₂O⁺ was in good agreement with that derived from ${{\rm{H}}}_{3}^{+}$. This difference suggests that CRs are attenuated as they penetrate from the atomic into the molecular regions of diffuse clouds, where ${{\rm{H}}}_{3}^{+}$ is most abundant. This, in turn, could constrain the low energy spectrum of CRs, since the effects of shielding are most pronounced at low energies (e.g., Silsbee & Ivlev 2019). Such shielding effects have been hypothesized from the fact that the CRIR inferred for dense molecular clouds is lower by an order of magnitude or more than that determined in diffuse molecular clouds from observations of ${{\rm{H}}}_{3}^{+}$ (e.g., Indriolo & McCall 2012). To confirm our finding, future determinations of the DR rate coefficients for rotationally cold ${{\rm{H}}}_{3}^{+}$ and H₂O⁺ will be needed to refine our estimates of the CRIR determined in these different environments.

Accurate CRIR values inferred from OH⁺ observations for the outer layers of diffuse clouds require reliable DR data for OH⁺ in its ground electronic, vibrational, and rotational level. Here, we have reported the first DR measurement for OH⁺ in its lowest energy levels. Using these data, we have generated a kinetic temperature rate coefficient suitable for diffuse cloud chemical models (i.e., for OH⁺ in its ground rotational level). Our results are valid for gas kinetic temperatures of 10–20,000 K. Fit parameters for the analytical representations of the kinetic temperature rate coefficient are given in Tables 1 and 2. The rate coefficient reported here is suitable for modeling a range of molecular environments in the ISM, such as diffuse cloud and photodissociation regions, and is especially important for the on-going analysis of OH⁺ observations along diffuse sightlines, e.g., from the SOFIA HyGal survey (Jacob et al. 2022).

We have explored the astrophysical implications of our experimental results in two parts: first, by adopting the single-zone diffuse cloud model of Bacalla et al. (2019) and calculating the multiplicative scaling factor for their previously inferred values of the CRIR due to our present DR rate coefficient; second, by incorporating our experimentally generated DR rate coefficient into the model of Neufeld & Wolfire (2017), along with additional updates to their model. Our combined experimental and modeling results point to a significant increase in the CRIR estimated from observations of OH⁺, which has important implications for the shielding of CRs in the ISM. Within diffuse clouds, we find that the CRIR determined from OH⁺ and H₂O⁺ observations is larger by a factor of 2.5 than that inferred from observations of ${{\rm{H}}}_{3}^{+}$, suggesting that significant CR shielding is possible in the transition from diffuse atomic to diffuse molecular clouds.

Financial support by the Max Planck Society is acknowledged. A.K., D.A.N., D.P., D.W.S., and M.G.W. were supported in part by the NASA Astrophysics Research and Analysis program under 80NSSC19K0969. L.W.I. was supported by the Deutsche Forschungsgemeinschaft (DFG, grant No. 431145392).

OH⁺ was generated in a standard Penning ion source using an H₂/O₂ gas mixture. A current of up to a few hundred nA was extracted, accelerated to an energy of ≈280 keV, and mass-to-charge selected using a series of dipole magnets, prior to injection into CSR. The number of injected ions was varied between 1 × 10⁵ and 2 × 10⁶, in order to remain within the linear counting regime of our neutral particle detector. The injected ions were stored on a closed orbit in CSR and merged with a magnetically confined electron beam in one of the straight sections of the ring. The storage ring lattice was set up in a configuration where the momentum dispersion, i.e., the coupling between ion orbit and momentum, was near zero in the electron–ion overlap section, which improves the electron cooling capabilities of CSR.

The electron beam can be used to cool the stored ions through phase-space cooling (explained below) and through internal-state-changing collisions (see Appendix C), as well as a target for collision studies. Here, we generated a beam current of 15 or 3.75 μA, in both cases with a circular density profile that was nearly uniform. We will refer to these cases as high or low electron current DR measurements, respectively. The electrons were initially electrostatically accelerated to an energy of ≈30 eV and magnetically guided into CSR, following a step-wise decreasing magnetic field. This produced an adiabatic expansion of the beam, increasing its size and simultaneously reducing the energy spread perpendicular to the magnetic field lines. The effective interaction region of the merged-beams setup is defined by a set of biased drift tubes that control the laboratory-frame energy of the electron beam. Within this region, we applied the lowest guiding field (10 mT). This resulted in an expansion factor of 30 and an effective beam diameter of 12.4 ± 0.6 mm for the high current and 10.2 ± 0.9 mm for the low current measurements. For the perpendicular electron-beam temperature, we take the ${k}_{{\rm{B}}}{T}_{\perp }={2.0}_{-0.5}^{+1.0}$ meV estimate for an expansion factor of 20 of Paul et al. (2022) and extend its lower limit (i.e., 2.0 ± 1.0 meV) to account for the factor of 1.5 times larger magnetic expansion in the present experiment.

Elastic collisions between the electrons and ions enabled us to reduce the size and energy spread of the stored ion beam by a process known as phase-space or electron cooling (Poth 1990). We achieved electron cooling by matching the average laboratory-frame electron-beam velocity to that of the stored ions. Observing the revolution frequency and momentum spread of the stored ions by means of a Schottky pickup (von Hahn et al. 2016), we have verified the condition of matched velocities and derived the corresponding laboratory-frame electron-beam energy E₀ = 9.029 ± 0.013 eV, also called the cooling energy. During collision experiments, we applied electron cooling immediately after ion injection. For the high current measurements, within 10 s, this resulted in the reduction of the injected ion-beam horizontal and vertical FHWM to <4.7 ± 0.6 and <3.5 ± 0.4 mm, respectively, in the interaction region. For the low current measurements, these ion-beam FWHM limits were reached within 30 s.

Using the electron beam as a collision target, we probed the DR energy dependence by detuning the electron-beam energy E_e from matched velocities. This was done by varying the voltage on the drift tubes. For monoenergetic beams, the resulting collision energy is given by Equation (5). We measured for predefined sets of E_d values and recorded the corresponding DR product rates. The neutral DR products were collected using a position-sensitive particle-counting detector downstream from the interaction region. We operated the detector in fragment-imaging mode, which enabled us to determine the kinetic energy released in the observed DR process and to image the ion-beam profile. Data for a given E_d were collected for 25 ms, followed by electron cooling for 100 ms, a reference energy at E_d = 10 eV for 25 ms (used for consistency checks), and electrons off for 25 ms. We also implemented an ∼5 ms waiting time between each of these steps so that the system could stabilize. The experimental merged-beams rate coefficient α^mb is proportional to the measured DR rate, as is explained next.

The measured absolute merged-beams rate coefficient is given by

$$\begin{eqnarray}&&{\alpha }^{\mathrm{mb}}({E}_{{\rm{d}}})=\displaystyle \frac{{R}_{{\rm{e}}}({E}_{{\rm{d}}})}{{\eta }_{{\rm{d}}}({E}_{{\rm{d}}})\xi {N}_{{\rm{i}}}{n}_{{\rm{e}}}({E}_{{\rm{d}}}){\hat{l}}_{0}/{C}_{0}}.\end{eqnarray} \tag{ B1 }$$

R_e(E_d) is the electron-induced count rate due to DR and, above 5 eV, partly also due to DE. The quantity η_d(E_d) is our detection efficiency for DR events. The parameter ξ is the fraction of ions enclosed by the electron beam in the effective interaction region. Here, ξ > 0.997 was determined from the profiles of the ion and electron beam for all the measurements. N_i is the number of stored ions. The electron density is n_e(E_d). The effective electron–ion overlap length is ${\hat{l}}_{0}\,=0.77\pm 0.01{\rm{m}}$, and its definition is described in the supplemental material of Kálosi et al. (2022). C₀ = 35.12 ± 0.05 m is the circumference of the ion-beam orbit.

R_e(E_d) is given by the difference between the count rates for the measurement and electrons-off steps. The electrons-off step measures the background count rate due to residual gas induced collisions and the intrinsic dark rate of the neutral detector. Reaction studies have shown that turning the electron beam on and off does not cause fluctuations of the background rate (Á. Kálosi et al. 2023, in preparation). The dark rate of the neutral detector was measured before ion injection.

We determined η_d(E_d) using the approach of Paul et al. (2022) to compare single-to-double particle hits on the detector. This parameter is dependent on the kinetic energy released in the DR process and the branching ratios between open dissociation channels. Both of these can vary as a function of E_d. At matched velocities, we obtained η_d(0 eV) = 0.690. We probed the energy dependence η_d(E_d) at a few dedicated energies between 1 and 10 meV. The values ranged between 0.675 and 0.690. In the 10 to few hundred meV range, in order to have sufficient counting statistics, we subdivided and averaged our data over 4 ranges and found values within the above limits. Finally, we characterized η_d around 1.6 eV, below which the O(¹D) + H(n = 1) DR channel dominates and above which the O(³P) + H(n = 2) channel opens (Strömholm et al. 1997). Probing at selected energies from 1.5 to 1.9 eV, η_d(E_d) increased monotonically from 0.640 to 0.720. For our work here, we used a single effective η_d = 0.68 ± 0.01 for the entire energy range studied. This introduces a negligible systematic error for E_d ≳ 1 eV compared to the statistical counting uncertainties at those energies.

The merged-beams rate coefficient was initially analyzed on a relative scale using a proxy for the storage-time dependent ion number N_i(t). In specific, we used the residual gas induced count rate R_g(t), determined from the electrons-off rate minus the dark rate. R_g results from collisions of the stored ions with residual gas, generating neutral products through either fragmentation or charge exchange. These measurements were followed by a calibration of R_g(t) to N_i(t), expressed as a proportionality factor S_b = R_g(t)/N_i(t). N_i was measured using beam bunching combined with a capacitive current pickup (PU-C; von Hahn et al. 2016). The induced voltage conversion to current for the PU-C was calibrated by Paul et al. (2022) with a 10% systematic uncertainty. S_b is specific to the present OH⁺ campaign as it is proportional to residual gas density and collisional rate coefficient for OH⁺ on the residual gas. We applied ion-beam bunching for durations of 0.5 s at various storage times and for ion numbers in the range of 2 × 10⁵ to 1 × 10⁶ to verify the linear behavior of S_b. The combined uncertainty due to counting statistics and the PU-C voltage-signal quality contributes an additional 10% uncertainty to the determination of S_b. In addition, we also monitored the relative DR rate coefficient over the energy range E_d < 30 meV. These DR measurements were performed over the course of the OH⁺ DR campaign. From a comparison of the corresponding relative rate coefficients, we were able to verify that the residual gas density remained constant during the campaign, to within the statistical counting uncertainties of the measurements.

The energy-dependent n_e(E_d) was determined from the measured electron-beam current and radius, accounting for the laboratory-frame electron energy versus E_d. The high current measurement value at matched velocities was n_e(0 eV) = (4.4 ± 0.4) × 10⁵ cm⁻³. For the low current measurements, n_e(0 eV) = (1.6 ± 0.3) × 10⁵ cm⁻³.

The total systematic uncertainty of the absolute scaling for our merged-beams DR rate coefficient α^mb was determined by three dominant contributions. These are the uncertainties of the PU-C voltage-to-current calibration factor, the simultaneous measurement of residual gas induced count rate and bunched-beam induced voltage, and the electron-beam density. All contributing uncertainties were treated as random sign errors and added in quadrature, resulting in a total systematic uncertainty of 17% for the high current measurements and 25% for the low current measurements.

We have constructed a detailed collisional-radiative model in order to predict the rotational populations for the stored OH⁺ ions. A detailed description of the underlying rate equations is given in the supplemental material of Kálosi et al. (2022), who constructed a similar model for CH⁺.

The OH⁺ from our ion source are electronically, vibrationally, and rotationally excited. The low-lying excited electronic states of OH⁺ decay on timescales <1 ms, except for the first excited a ¹Δ state, which has a calculated radiative lifetime to the X ³Σ⁻ ground electronic state on the order of ∼30 ms (Strömholm et al. 1997). Werner et al. (1983) calculated the vibrational transition probabilities within the X ³Σ⁻ ground electronic state, finding radiative lifetimes <5 ms. Given the timescale of our experiments on the order of several tens to hundreds of seconds, we can safely assume that the ions for our measurements were in their ground electronic and vibrational states.

The rotational structure of X ³Σ⁻ OH⁺ is best represented by the quantum number N ≥ 0, analogous to the rotational ladder of a ¹Σ molecule. Each level with N > 0 is split into three J sublevels by spin–spin and spin–rotation interactions. Each J sublevel lies <3 cm⁻¹ (4.3 K) away from the N level rotational energy (Hechtfischer et al. 2019). The three J sublevels for a given N have practically equal radiative lifetimes. This substructure was omitted in the final model as radiative cooling in the CSR radiation field, when including the sublevels, leads to equilibrium J populations similar to an ≈17 K Boltzmann distribution. This temperature is sufficiently high compared to the level splittings so that the sublevels for each N are statistically populated with weights 2J + 1. This enabled us to neglect the sublevels in our final model for the time evolution of the N level populations.

The approximate rotational energies of OH⁺ X ³Σ⁻(v = 0) were calculated using the B₀ and D₀ spectroscopic constants of Hodges & Bernath (2017) as

$$\begin{eqnarray}&&{E}_{N}={B}_{0}[N(N+1)]-{D}_{0}{\left[N(N+1)\right]}^{2}.\end{eqnarray} \tag{ C1 }$$

Levels up to N = 19 were included in the model. We used the theoretical 2.26 D dipole moment of Werner et al. (1983) to calculate Einstein coefficients for spontaneous emission for rotational transitions using Equation (S1) from Kálosi et al. (2022), with J replaced by N.

Our radiative cooling model also accounts for the ambient radiation field in CSR, following the approach of Meyer et al. (2017). Their model consists of two components: one for the thermal radiation of the cryogenic chambers and the second for room-temperature leaks from various openings in CSR. Here, we estimated the cryogenic component to be T_low ≈ 6 K, based on the measured chamber temperatures. The room-temperature component is fixed to 300 K. Previous work has found its fraction to be ε = (1.0 ± 0.3) × 10⁻² (Kálosi et al. 2022). Of these two components, the room-temperature fraction dominates the accuracy of our rotational population predictions when approaching equilibrium with the CSR radiation field.

Rotational-level-changing collisions between electrons and ions have been experimentally demonstrated for CH⁺ to affect the rotational level populations in our experimental setup (Kálosi et al. 2022). To examine the role of these collisions for OH⁺, we adopted the theoretical electron-impact rotational excitation cross section of Hamilton et al. (2018) for ΔN = 1 and 2 transitions of OH⁺. The corresponding deexcitation cross sections were obtained by applying the principle of detailed balance. We then calculated rotational-level-changing merged-beams rate coefficients for all E_d corresponding to those used for our DR rate coefficient measurements. For a given transition, the importance of rotational-level-changing collisions can be shown by the ratio of the Einstein ${A}_{N^{\prime} \to N^{\prime\prime} }$ rate to the N-level-changing merged-beams rate coefficient ${\alpha }_{N^{\prime} \to N^{\prime\prime} }$, which give the critical electron density

$$\begin{eqnarray}&&{n}_{{\rm{c}}}={A}_{N^{\prime} \to N^{\prime\prime} }/{\alpha }_{N^{\prime} \to N^{\prime\prime} }.\end{eqnarray} \tag{ C2 }$$

Our DR measurements here were carried out for ions primarily in the N = 0 and 1 levels. For the N = 1 → N = 0 transition, we calculated A_1→0 = 1.89 × 10⁻² s⁻¹ and a critical electron density of (4.5 ± 0.8) × 10³ cm⁻³ at matched velocities. This value is comparable to the typical ring-averaged electron density of (9.6 ± 0.9) × 10³ cm⁻³ for the high current measurements, meaning that both collisions and radiative interactions are important in our experiment for the N = 0 and 1 levels.

C.1. Results for High Electron Current DR Measurements

The high electron current measurements were designed to collect statistically significant DR signal compared to background at all collision energies of interest, within the two weeks of the measurement campaign. The complete population model, including both radiative cooling in the CSR radiation field and rotational-level-changing collisions, was used to optimize the internal-state preparation scheme for the high current DR measurements, for which the measured data are shown in Figure 1. The scheme consisted of three ion storage phases. First, following injection, we applied electron cooling at matched velocities for 52 s. During this phase, rotational-level-changing collisions accelerated the rotational cooling of the stored ions compared to only radiative cooling. Next, we turned off the electron beam for 104 s to equilibrate the rotational level populations with the CSR radiation field. During this phase, the FWHM of the ion-beam profile slowly increased, due to ion intrabeam heating processes. Finally, we measured for predefined sets of E_d for 25 s. The α^mb data shown in Figure 1 have been evaluated after excluding the first 3 s of the measurement phase where the ion-beam FWHM was decreasing, due to electron cooling.

The modeled relative populations are shown in Figure 5 for the levels of interest (N ≤ 2) as a function of the storage time from injection. The initial rotational level populations are represented by a Boltzmann distribution within the ground vibrational state. The final population results for levels N ≤ 2 are insensitive to the initial populations for an assumed initial rotational temperature >300 K. In the initial cooling phase, the time needed for the N = 0 population to approach its end value is reduced by a factor of 3 compared to only radiative cooling. After turning off the electron current, the N = 0 population continues to grow as the system approaches equilibrium with the CSR radiation field. In the measurement phase, the populations are determined by the combined effects of collisions and radiation. The additional cycling between the selected E_d and accompanying electron cooling, reference, and electron-off steps results in a slight increase of the effective rotational temperature.

Zoom In Zoom Out Reset image size

The model uncertainties were determined primarily by varying the experimental parameters T_⊥, T_low, and ε. Two additional parameters can affect the model predictions, the dipole moment and the electron-impact rotational excitation cross sections for OH⁺. The magnitudes of these two parameters primarily affect the timescale of the rotational population evolution. The dipole moment is expected to be accurate within 7%, based on a comparison of the same level of theory and a measurement for OH⁻ (Meyer et al. 2017). The theoretical model used by Hamilton et al. (2018) to calculate the cross sections for rotational-state-changing collisions has been experimentally benchmarked for CH⁺ by Kálosi et al. (2022), and the agreement with theory was better than ∼40%. For the present high electron current results, when the populations are near equilibrium, the uncertainties in the dipole moment and the collisional cross sections do not contribute significantly to the uncertainty budget.

The model uncertainty determination uses a Monte Carlo approach. The dipole moment, electron-impact rotational excitation cross sections, and ε are drawn from normal distributions. For the cross sections, we assumed a systematic uncertainty of 40% at all energies and for all levels. The initial rotational temperature, T_low, and T_⊥ are drawn from uniform distributions to include the physical limits on these parameters. For the initial rotational temperature, an upper limit of 5000 K was chosen. For T_low, the lower and upper limits were set to 4 and 8 K, respectively. At each storage-time value, we evaluated the 16th and 84th percentiles of the population distributions, which is equivalent to ±one-sigma for a normal distribution. The shaded areas in Figure 5 are given for these percentiles.

Using our model, we estimate the average relative populations of rotational levels in the measurement time window to be 0.80 ± 0.05 for N = 0, 0.19 ± 0.04, for N = 1, and ≈0.01 for N = 2. The uncertainties for the two lowest N levels are anticorrelated, and the sum of their populations totals to ≈0.99.

C.2. Results for Low Electron Current DR Measurements

The low electron current measurements were designed to reduce the influence of rotational-level-changing collisions. The corresponding ring-averaged electron density was (3.5 ± 0.7) × 10³ cm⁻³ at matched velocities. These low current measurements were performed by cycling between the selected E_d and the accompanying electron cooling, reference, and electron-off steps, starting immediately after injection. The relative population results of the model are shown in Figure 6 for N ≤ 2. The initial rotational level populations are represented by a Boltzmann distribution within the ground vibrational state with an assumed initial rotational temperature of >300 K.

Zoom In Zoom Out Reset image size

The model uncertainties were determined using the Monte Carlo approach described above. Compared to the results shown in Figure 5, the uncertainty of the predictions is smaller at storage times where the N = 0 population is increasing, which enables us to more reliably extract the level-specific ${\alpha }_{N}^{\mathrm{mb}}$ for N = 0 than when using a higher electron current.

The determination of level-specific ${\alpha }_{N}^{\mathrm{mb}}$ values is enabled by the time evolution of the rotational level populations for stored ions. Here, we have determined ${\alpha }_{N}^{\mathrm{mb}}$ using the data from our low electron current DR measurements.

The measured α^mb(t) for a given storage-time window t is the sum

$$\begin{eqnarray}&&{\alpha }^{\mathrm{mb}}(t)=\displaystyle \sum _{N}{\bar{p}}_{N}(t){\alpha }_{N}^{\mathrm{mb}},\end{eqnarray} \tag{ D1 }$$

where the rotational-level-specific rate coefficients ${\alpha }_{N}^{\mathrm{mb}}$ are weighted by the average populations ${\bar{p}}_{N}(t)$ in the time window. The results shown in Figure 6 were used to calculate ${\bar{p}}_{N}(t)$. By including more storage-time windows than the number of contributing rotational levels, we can perform a least squares fit to Equation (D1) and extract ${\alpha }_{N}^{\mathrm{mb}}$.

The earliest storage-time window included in the analysis spans from ≈25 to ≈33 s, followed by several windows of the same length. Earlier storage times where the ion-beam FWHM was decreasing due to electron cooling were excluded. Based on our collisional-radiative model results shown in Figure 6, a rotational equilibrium was reached after ≈150 s of ion storage. The final storage-time window was extended to cover the times from ≈145 to ≈250 s. Only levels with >1% relative population (N ≤ 2) were included in the fit. In Figure 7(a), we compare the measured α^mb for two selected storage-time windows in the experiment: from ≈25 to ≈33 s and ≈145 to ≈250 s.

Zoom In Zoom Out Reset image size

In order to include the uncertainty of the collisional-radiative model in the level-specific results, we used our Monte Carlo model to generate a set of ${\bar{p}}_{N}(t)$, as described above. A fit by Equation (D1) was performed for each ${\bar{p}}_{N}(t)$ within the set. In the selected time windows, the N = 2 population was small but nonnegligible. This led to an uncertainty for ${\alpha }_{2}^{\mathrm{mb}}$ that was much larger than the measured rate coefficient at any storage time. Due to a strong correlation between the N = 2 and 1 populations, this uncertainty propagated into the ${\alpha }_{1}^{\mathrm{mb}}$ result as well. The combination of the large ${\alpha }_{2}^{\mathrm{mb}}$ uncertainty and the small differences between the measured α^mb for the sampled storage-time windows (see Figure 7(a)) led to an uncertainty for the ${\alpha }_{1}^{\mathrm{mb}}$ result that was also larger than the measured rate coefficient at all storage times. However, for ${\alpha }_{0}^{\mathrm{mb}}$, each fit yielded values ${\alpha }_{0,i}^{\mathrm{mb}}$ with uncertainties ${\delta }_{0,i}^{\mathrm{mb}}$ smaller than the fitted values. The final ${\alpha }_{0}^{\mathrm{mb}}$ shown in Figure 7 was calculated as the unweighted mean of the individual fits ${\alpha }_{0,i}^{\mathrm{mb}}$, and its statistical uncertainty was calculated as the unweighted mean of the individual fit uncertainties ${\delta }_{0,i}^{\mathrm{mb}}$. The variance of the individual fits ${s}_{0}^{\mathrm{mb}}$ was much smaller than any of the individual fit uncertainties ${\delta }_{0,i}^{\mathrm{mb}}$.

The final ${\alpha }_{0}^{\mathrm{mb}}$ of the fits by Equation (D1) for the Monte Carlo generated ${\bar{p}}_{N}(t)$ is plotted in Figures 7(b) and (c). Figure 7(b) compares the extracted ${\alpha }_{0}^{\mathrm{mb}}$ to the data from the final storage-time window of the low electron current measurements. Figure 7(c) compares the extracted ${\alpha }_{0}^{\mathrm{mb}}$ to the data from the high current measurements, i.e., from Figure 1. The present analysis shows that the measured rate coefficient for our high electron current DR measurements well approximates ${\alpha }_{0}^{\mathrm{mb}}$. When converting the measured data into α^k, the observed differences become negligible. Hence, α^k(T_k) shown in Figure 2 is well suited for chemical kinetics models for ground rotational level OH⁺.

Here, we provide two different analytic representations for the OH⁺ DR kinetic temperature rate coefficient α^k and its one-sigma error band shown in Figure 2. First, we follow Novotný et al. (2013) and fit α^k with an optimized function, with a small number of fit parameters, that accounts for typical DR features, i.e., broad peaks or dips from resonances. The fit function is given by

$$\begin{eqnarray}\begin{array}{rcl}{\alpha }^{{\rm{k}}}({T}_{{\rm{k}}})[{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}] & = & A{\left(\displaystyle \frac{300}{{T}_{{\rm{k}}}[{\rm{K}}]}\right)}^{n}\\ & & +{T}_{{\rm{k}}}{\left[{\rm{K}}\right]}^{-1.5}\displaystyle \sum _{r=1}^{4}{c}_{r}\exp \left(-\displaystyle \frac{{T}_{r}}{{T}_{{\rm{k}}}[{\rm{K}}]}\right),\end{array}\end{eqnarray} \tag{ E1 }$$

and the parameters are listed in Table 1. The maximum relative deviation of the fit is 0.5% below 6000 K and 1.4% from there to 20,000 K.

Second, we give the Arrhenius–Kooij (AK) representation, as is typically used in astrochemistry, combustion chemistry, and other chemical models and databases. Following the approach of Paul et al. (2022), we provide a set of piecewise-joined fit functions on several temperature intervals, since we cannot model the experimental results with a single AK fit function. This approach introduces discontinuities in the temperature dependence of the analytical rate coefficient between the temperature intervals, which can be avoided by using the representation by Equation (E1). The temperature-interval fit function is given by

$$\begin{eqnarray}&&{\alpha }^{{\rm{k}}}({T}_{{\rm{k}}})[{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1}]=A{\left(\displaystyle \frac{{T}_{{\rm{k}}}[{\rm{K}}]}{300}\right)}^{\beta }{e}^{-\tfrac{\gamma }{{T}_{{\rm{k}}}[{\rm{K}}]}},\end{eqnarray} \tag{ E2 }$$

and the parameters for each temperature interval are listed in Table 2. The maximum relative deviation of the fit is 1%. Following KIDA conventions, the relative uncertainty is described by the log-normal factor $F=\exp \left({\rm{\Delta }}{\alpha }^{{\rm{k}}}/{\alpha }^{{\rm{k}}}\right)$. This quantity is fitted on the same temperature intervals as for the AK fits by the fit function

$$\begin{eqnarray}&&F({T}_{{\rm{k}}})={F}_{0}\exp \left(g\left(\displaystyle \frac{1}{{T}_{{\rm{k}}}[{\rm{K}}]}-\displaystyle \frac{1}{300}\right)\right).\end{eqnarray} \tag{ E3 }$$

Due to the asymmetric error bands of α^k, the log-normal factor is calculated as the average of the upper and lower error bands. The continuity of both the α^k and error fits are guaranteed on the border of each temperature range.