© ESO, 2014

1. Introduction

Isotopic fractionation reactions have already been invoked by Watson (1976) and Dalgarno & Black (1976) to explain the enrichment of heavy isotopes of molecules in dark cold interstellar cloud environments. The exothermicity involved in the isotopic exchange reaction directly depends on the difference of the zero-point energies between the two isotopes, if one assumes that the reaction proceeds in the ground-rovibrational states of both the reactant and product molecule. This assumption has been questioned for the reaction H + HD⇌H₂D⁺ + H₂, where some rotational excitation in H₂ may reduce the efficiency of the reverse reaction (Pagani et al. 1992; Hugo et al. 2009).

In this paper we revisit some fractionation reactions involved in the ¹²C/¹³C, ¹⁶O/¹⁸O, and ¹⁴N/¹⁵N balance by reinvestigating the potential energy surfaces (PESs) involved in the isotopic exchange reactions. Within the Born-Oppenheimer approximation, a single nuclear-mass-independent PES is considered for all isotopic variants of molecules under consideration. The nuclear motions are introduced subsequently and isotopologues, molecules of different isotopic compositions and thus different masses, possess different rotational constants, different vibrational frequencies, and different ground-state (zero-point) vibrational energies, in other words, different thermodynamic properties (Urey 1947). Differences in zero-point energies can become important under cool interstellar cloud conditions where molecules rather undergo isotopic exchange (fractionation) than react chemically. This thermodynamic effect may result in isotopologue abundance ratios (significantly) deviating from the elemental isotopic ratios. Knowledge of the abundance ratios may in return provide valuable information on molecular processes at low collision energies.

As far as astrophysical models are concerned, ¹³C and ¹⁸O isotopic fractionation studies involving CO and HCO⁺ (Le Bourlot et al. 1993; Liszt 2007; Röllig & Ossenkopf 2013; Maret et al. 2013) are based on the pionneering paper by Langer et al. (1984), who referred to the experimental studies by Smith & Adams (1980) and used theoretical spectroscopic parameters for the isotopic variants of HCO⁺ reported by Henning et al. (1977). Lohr (1998) derived the harmonic frequencies and equilibrium rotational constants for CO, HCO⁺, and HOC⁺ at the configuration interaction (including single and double excitations) level of theory (CISD/6-31G**) and tabulated reduced partition function ratios and isotope exchange equilibrium constants for various isotope exchange reactions between CO and HCO⁺. Surprisingly, this paper has not received much attention in the astrophysical literature, and its conclusions have never been applied.

The studies of Langer et al. (1984) and Lohr (1998) led to qualitatively different conclusions regarding the following fractionation reaction: (1)which was found to be endothermic with ΔE/k_B = − 5 K by Langer et al. (1984) and exothermic with ΔE/k_B = 12.5 K by Lohr (1998), where k_B is the Boltzmann constant. To clear up this discrepancy, we carried out numerically exact calculations for the vibrational ground state of HCO⁺ using a PES previously developed by Mladenović & Schmatz (1998). Our calculations gave ΔE/k_B = 11.3 K for reaction (1), in good agreement with the harmonic value of Lohr (1998). In addition, we noticed that the ΔE/k_B values of Henning et al. (1977) for the reactions (2)and (3)were quoted as 17 ± 1 K and 7 ± 1 K by Smith & Adams (1980) and as 9 and 14 K by Langer et al. (1984). Reconsidering the original values of Henning et al. (1977), we found that Langer et al. (1984) permuted the zero-point energies for H¹³C¹⁶O⁺ and H¹²C¹⁸O⁺ in Table 2 of their paper. From the original spectroscopic parameters of Henning et al. (1977), we derive ΔE/k_B = 10.2 K for reaction (1), in good agreement with our result and the result of Lohr (1998).

The permutation of the zero-point vibrational energies of H¹³C¹⁶O⁺ and H¹²C¹⁸O⁺ affects the exothermicities and rate coefficients summarized in Tables 1 and 3 of the paper by Langer et al. (1984). These data are actually incorrect for all isotope fractionation reactions CO+HCO⁺, except for (4)The rate coefficients reported by Langer et al. (1984) are still widely used when including isotopes such as ¹³C and ¹⁸O into chemical (molecular) networks (Maret et al. 2013; Röllig & Ossenkopf 2013). With these points in mind, our goal is to provide reliable theoretical estimates for the zero-point vibrational energies first of H/DCO⁺ and to derive proper rate coefficients for the related fractionation reactions. Our improved results for the exothermicities and rate coefficients are summarized in Tables 2 and 5.

Henning et al. (1977) also reported spectroscopic parameters for various isotopic variants of N₂H⁺. This was our initial motivation to expand the present study to ion-molecule reactions between N₂H⁺ and N₂. ¹⁵N fractionation in dense interstellar clouds has been first considered by Terzieva & Herbst (2000), who referred to the experimental information of the selected ion flow-tube (SIFT) studies at low temperatures of Adams & Smith (1981).

The reactions discussed in this paper, CO+HCO⁺ and N₂ + HN, are the most obvious candidates for isotopic fractionation. In addition, they have been studied in the laborataory, which allows a detailed discussion. A similar reaction has been invoked for CN (Milam et al. 2009), but no experimental and/or theoretical information is available there.

In the Langevin model, the long-range contribution to the intermolecular potential is described by the isotropic interaction between the charge of the ion and the induced dipole of the neutral. Theoretical approaches based on this standard assumption may qualitatively explain the behaviour of the association rates. However, they generally provide rate coefficients that are higher than experimental results (Langer et al. 1984). The rate coefficients for ion-molecule reactions are quite constant at higher temperatures but increase rapidly at lower temperatures. The latter feature is an indication of barrierless PESs. The electrostatic forces are always attractive and can be experienced over large distances even at extremely low temperatures relevant for dark cloud enviroments. Short-range forces appear in closer encounters of interacting particles and may (prominently) influence the overall reaction rate. To explore the short-range effects we also undertake a study of linear proton-bound ionic complexes arising in the reactions involving HCO⁺, HOC⁺, and N₂H⁺ with CO and N₂, which are common interstellar species.

Our theoretical approach is described in Sect. 2. The specific aspects of the fractionation reactions of HCO⁺ and HOC⁺ with CO are reanalysed in Sect. 3.1 and the fractionation reactions N₂H⁺+N₂ in Sect. 3.2. We discuss the equilibrium constants and rate coefficients of CO+HCO⁺/HOC⁺ in Sect. 4.1, providing the astrochemical implications of the new exothermicities in Sect. 4.2. The isotope fractionation reactions N₂H⁺+N₂ are considered including the nuclear spin angular momentum selection rules in Sect. 4.3. The linear proton-bound cluster ions are analysed in Sect. 4.4. Our concluding remarks are given in Sect. 5.

2. Calculations

The global three-dimensional PES developed by Mladenović & Schmatz (1998) for the isomerizing system HCO⁺/HOC⁺ and by Schmatz & Mladenović (1997) for the isoelectronic species N₂H⁺ were used in the rovibrational calculations. These two PESs still provide the most comprehensive theoretical descriptions of the spectroscopic properties for HCO⁺, HOC⁺, and N₂H⁺ and are valid up to the first dissociation limit. Potential energy representations recently developed by Špirko et al. (2008) and by Huang et al. (2010) reproduce the experimental fundamental transitions within 11[6] and 4[3] cm^-1 for N₂H⁺[N₂D⁺], respectively, whereas the PES of Schmatz & Mladenović (1997) predicts the fundamental transitions for both N₂H⁺ and N₂D⁺ within 2 cm^-1.

The rovibrational energy levels of HCO⁺/HOC⁺ and N₂H⁺ are calculated by a numerically exact quantum mechanical method, involving no dynamical approximation and applicable to any potential energy representation. The computational strategy is based on the discrete variable representation of the angular coordinate in combination with a sequential diagonalization/truncation procedure (Mladenović & Bačić 1990; Mladenović & Schmatz 1998). For both molecular systems, the rovibrational states are calculated for the total angular momentum J = 0−15. These rovibrational energies are used to evaluate theoretical partition functions and to model rate coefficients for proton transfer reactions involving HCO⁺ and N₂H⁺.

To gain a first insight into dynamical features of ion-molecule reactions, additional electronic structure calculations were carried out for linear proton-bound cluster ions of HCO⁺, HOC⁺, and N₂H⁺ with CO and N₂. The PESs were scanned by means of the coupled cluster method with single and double excitations including perturbative corrections for triple excitations [CCSD(T)] in combination with the augmented correlation consistent triple ζ basis set (aug-cc-pVTZ). Only valence electrons were correlated. The ab initio calculations were carried out with the MOLPRO (Werner et al. 2012) and CFOUR (Stanton et al. 2012) quantum chemistry program packages.

3. Results

The PES of Mladenović & Schmatz (1998) provides a common potential energy representation for the formyl cation, HCO⁺, and the isoformyl cation, HOC⁺, where the local HOC⁺ minimum is 13 878 cm^-1 (166 kJ mol^-1) above the global HCO⁺ minimum. Inclusion of the zero-point energy reduces this separation by 640–650 cm^-1 for the hydrogen-containing isotopologues and by 570–580 cm^-1 for the deuterium variants. The angular motion is described by a double-minimum anharmonic potential with a non-linear saddle point at 26 838 cm^-1 (321 kJ mol^-1) above the HCO⁺ minimum, such that low-lying states of HCO⁺ and HOC⁺ are well separated.

The PES of Schmatz & Mladenović (1997) for N₂H⁺ (dyazenilium) has two equivalent colinear minima as a consequence of the S₂ permutation symmetry, separated by an isomerization barrier 17 137 cm^-1 (205 kJ mol^-1) above the energy of the linear geometries. Low-lying states of N₂H⁺ are, thus, localized in one of the two wells. The double-well symmetry and nuclear spin symmetries are lifted for mixed nitrogen isotope forms.

3.1. Reaction of CO with HCO⁺ and HOC⁺

The ground-state vibrational energies calculated in this work for isotopic variants of HCO⁺ and HOC⁺ are collected in Table 1. There we additionally show the harmonic zero-point energy estimates of Lohr (1998) and the anharmonic values of Martin et al. (1993) available only for three isotopologues, as well as the values obtained by Langer et al. (1984) and in the present work from the spectroscopic [CI(corr)] parameters of Henning et al. (1977). Our values for CO are computed at the theoretical level used to construct the PES for HCO⁺/HOC⁺ [CCSD(T)/cc-pVQZ].

The isotopologues in Table 1 are arranged in order of increasing total molecular mass. For CO and H/DOC⁺, the zero-point energies decrease as the total molecular mass increases, which is not the case for H/DCO⁺. Inspection of the table shows that the substitution of the central atom by its heavier isotope (¹²C→¹³C in H/DCO⁺ and ¹⁶O→¹⁸O in H/DOC⁺) results in a more pronounced decrease of the zero-point energy than the isotopic substitution of the terminal atom (¹⁶O→¹⁸O in H/DCO⁺ and ¹²C→¹³C in H/DOC⁺). This feature shared by H/DCO⁺ and H/DOC⁺ in Table 1 is easy to rationalize since a central atom substitution affects all three vibrational frequencies.

The zero-point energy differences for the proton transfer reactions CO+HCO⁺/HOC⁺ are listed in Table 2. The reactions involving the formyl cation are labelled with F and the reactions involving the isoformyl cation with I. The deuterium variant of reaction F1 is denoted by F1(D) and similar for all other reactions. The reactions F1, F2, F3, F4, F5, and F6 are numbered as 1004, 3408, 3407, 3457, 3406, and 3458 by Langer et al. (1984).

In Table 2, our results, the values rederived from the spectroscopic parameters of Henning et al. (column HKD77), and the harmonic values of Lohr (column L98) all agree within less than 5 K. These three data sets predict the same direction for all listed reactions, whereas Langer et al. (column LGFA84) reported reaction F6 as endothermic. The replacement of our theoretical values for CO by the experimental values taken from Huber & Herzberg (1979) affects the zero-point energy differences by at most 0.4 K.

The general trend seen in Table 2 is that ¹³C is preferentially placed in H/DCO⁺ and ¹⁸O in H/DOC⁺. This is in accordance with Table 1, showing a stronger decrease of the zero-point energy upon isotopic substitution of the central atom. The substitution of the two ¹⁶O by ¹⁸O or the two ¹²C by ¹³C has nearly no influence on the exothermicities, as seen by comparing ΔE for reactions F1, F1(D), I1, I1(D) with ΔE for reactions F2, F2(D), I2, I2(D) and silimar for reactions F3, F3(D), I3, I3(D) versus F4, F4(D), I4, I4(D). Slightly higher exothermicities appear for reactions involving deuterium. The exothermicities for the reactions with the isoformyl isomers are lower than for the reactions with the formyl forms.

From the measured forward reaction k_f and backward reaction k_r rate coefficients, Smith & Adams (1980) calculated the experimental zero-point energy differences using (5)where K_e is the equilibrium constant. The total estimated error on k_f and k_r is reported to be ±25% at 80 K. Table 2 indicates that the new/improved theoretical values, and the experimental finding for reaction F1 agree within the experimental uncertainty. For reactions F3 and F6, we see that the theoretical results consistently predict a higher ΔE value for H¹²C¹⁸O⁺ reacting with ¹³C¹⁶O (reaction F6) than for H¹²C¹⁶O⁺ reacting with ¹²C¹⁸O (reaction F3), whereas the opposite was derived experimentally. Note that Smith & Adams (1980) reported for ¹³C¹⁶O reacting with H¹²C¹⁸O⁺ in addition to reaction F6 also a yield of 10% for the rearrangement channel (6)The latter transformation is not of a simple proton-transfer type (but bond-rearrangement type) and must involve a more complicated chemical mechanism probably including an activation energy barrier.

3.2. Reaction of N₂ with N₂H⁺

The zero-point vibrational energies calculated for N₂H⁺ are summarized in Table 3. In addition to the results obtained for the PES of Schmatz & Mladenović (1997), Table 3 also provides the values we derived from the spectroscopic parameters of Huang et al. (2010, column HVL10) and of Henning et al. (1977, column HKD77). The values for N₂ are taken from Huber & Herzberg (1979). The zero-point energy differences are given in Table 4. As seen there, our results agree with the values obtained from the spectroscopic parameters of Huang et al. (2010) within 0.4 K. The reactions involving diazenylium (or dinitrogen monohydride cation) are labelled with D in Table 4.

The ¹⁴N/¹⁵N substitution at the central-atom position lowers the zero-point energy more than the terminal-atom substitution (Table 3), such that ¹⁵N preferentially assumes the central position in N-N-H⁺ in all reactions of N₂H⁺ with N₂ in Table 4. The exothermicities are found to be slightly higher for the reactions involving deuterium.

In Table 4, the experimental (SIFT) results of Adams & Smith (1981) are listed as given in their paper. Note, however, that the elementary isotope fractionation reactions D2 and D3 involve common reactants, whereas reactions D4 and D5 have common products. The two reaction pairs are related by the ¹⁴N→ ¹⁵N substitution. Using thermodynamic reasoning, it is easy to verify that the following relationship (11)is strictly fulfilled for , where corresponds to reaction D6, (12)Note that the factor also provides the thermal population of ¹⁵N¹⁴NH⁺ relative to ¹⁴N¹⁵NH⁺.

4. Discussion

The equilibrium constant K_e for the proton transfer reaction (13)under thermal equilibrium conditions is given by (14)where Q(X) is the full partition function for the species X. Making the translation contribution explicit, we obtain (15)where the mass factor f_m is given by (16)for m(X) denoting the mass of the species X, whereas ΔE stands for the zero-point energy difference between the reactants and the products, (17)The zero-point energies E₀ are measured on an absolute energy scale. For isotope fractionation reactions, the internal partition function, Q_int, includes only the rovibrational degrees of freedom (no electronic contribution) and is given by the standard expression (18)where for a total angular momentum J is the rovibrational energy measured relative to the corresponding zero-point energy (J = 0). The factor (2J + 1) accounts for the degeneracy relative to the space-fixed reference frame and g for the nuclear spin (hyperfine) degeneracy, (19)in which α labels the constituent nuclei having the nuclear spin I_N,α. For the nuclei considered in the present work, we have I_N(H)= 1/2, I_N(D)= 1, I_N(¹²C)= 0, I_N(¹³C)= 1/2, I_N(¹⁶O)= 0, I_N(¹⁸O)= 0, I_N(¹⁴N)= 1, and I_N(¹⁵N)= 1/2.

Introducing the ratio (20)the equilibrium constant is compactly written as (21)where the partition function factor F_q is (22)For reactions proceeding in the ground-rovibrational states of the reactants and the products, the partition function ratios and are both equal to 1. Even then the corresponding partition function factor F_q of Eq. (22) is, strictly speaking, different from 1 because of the mass term f_m defined by Eq. (16). For the reactions F1–F6 in Table 2, for instance, the values are 0.998, 0.998, 0.997, 0.997, 0.995, and 1.002, respectively, which are different from 1 at most by 0.5%.The values are computed from Eq. (16) using the following atomic masses m(H) = 1.007825035, m(D) = 2.014101779, m(¹²C) = 12, m(¹³C) = 13.003354826, m(¹⁶O) = 15.99491463, and m(¹⁸O) = 17.9991603 u, as given by Mills et al. (1993).

The terms of Q_int in Eq. (18) decrease rapidly with energy and J. In the low-temperature limit relevant for dark cloud conditions, the discrete rotational structure of the ground-vibrational state provides the main contribution to Q_int. That said, the rotational energy cannot be treated as continuous and one must explicitly sum the terms to obtain Q_int. With increasing temperature, the rotational population in the ground-vibrational state increases and other vibrational states may also become accessible, leading to partition function factors F_q, which may show (weak) temperature dependences.

For a given PES, numerically exact full-dimensional strategies insure the determination of accurate level energies and therefore accurate partition functions and equilibrium constants. To predict/estimate rate coefficients, we may use kinetic models, such as e.g. the Langevin collision rate model for ion-molecule reactions. Uncertainties in the rate coefficients are thus defined by uncertainties in the model parameters. In the case of the system CO+HCO⁺, we employ the total rate coefficients from Table 3 of Langer et al. (1984) and the uncertainties of these quantities also provide the uncertainties of the rate coefficients derived in the present work.

4.1. Reaction of CO with HCO⁺

The equilibrium constants for HCO⁺ reacting with CO are given in Table 5. Our K_e values are obtained in accordance with Eq. (21) by direct evaluation of the internal partition functions Q_int from the computed rovibrational energies. The forward reaction k_f and backward reaction k_r rate coefficients are calculated using our ΔE values and the total temperature-dependent rate coefficients k_T given by Langer et al. (1984), where (23)such that The results for the deuterium variants are also listed in Table 5. Their rate coefficients k_f and k_r are calculated assuming the same total rate coefficients k_T as for the H-containing forms (due to nearly equal reduced masses). For the purpose of comparison, note that the Langevin rate for CO+HCO⁺ is k_L = 8.67 × 10^-10 cm³ s^-1.

The partition function factors F_q deviate from 1 by approximately 2% in Table 5. They also exhibit marginal temperature dependences. This reflects the influence of rotational and vibrational excitations in the reactants and the products. Only the rotationally excited ground-vibrational states contribute to Q_int at temperatures T< 200 K. The contribution of the bending ν₂ level is 0.5% at 200 K and 3.6–3.8% at 300 K, whereas the contributions from 2ν₂ are 0.1% at 300 K. To appreciate the effect of F_q, we employed the rate coefficients measured at 80 K by Smith & Adams (1980) to determine the ΔE value for reactions F1, F3, and F6 by means of Eq. (21). Using the F_q values from Table 5, we obtain ΔE/k_B of 13.8, 15.1, and 4.8 K, respectively. For F_q = 1, we find 12.3, 14.6, and 3.8 K, which are lower by 1.5 K (12%), 0.5 K (3%), and 1 K (26%) than the former F_q ≠ 1 results.

The equilibrium constants K_e reported by Langer et al. (1984) deviate from the present results and those of Lohr (1998) very prominently at low temperatures in Table 5. At 10 K, we see deviations of 43% and 217% with respect to our values for reactions F1 and F3, respectively, and the related rate coefficients k_f,k_r are accordingly different. An even larger discrepancy is seen for reaction F6 of Eq. (1), which was previously predicted to be endothermic. In accordance with this, the values of k_f and k_r derived by Langer et al. (1984) are given in the reverse positions as (k_r,k_f) for reaction F6 in Table 5.

The deuterium variants in Table 5 are associated with slightly lower F_q values and somewhat higher low-temperature K_e, resulting in somewhat faster foward reactions and slower backward reactions.

4.2. Astrochemical implications

We investigated the role of these new derived exothermicities under different density conditions relevant to cold dark interstellar clouds. We display in Table 6 steady-state results for isotopic ratios of CO, HCO⁺ and DCO⁺ for two chemical models performed at a temperature of 10 K with a cosmic ionization rate ζ of 1.3 × 10^-17 s^-1 per H₂ molecule with the old ΔE values by Langer et al. (Model A: LGFA84) and the present ΔE values listed in Table 2 (Model B). The ratios of the principal isotope to the minor isotope obtained for Model A, R_A, and for Model B, R_B, are compared using the relative difference δ, (26)The chemical network contains 288 chemical species including ¹³C and ¹⁸O containing molecules as well as deuterated species and more than 5000 reactions. We assumed that the elemental ¹²C/¹³C and ¹⁶O/¹⁸O isotopic ratios are 60 and 500, so that any deviation relative to these values measures the amount of enrichment/depletion with respect to the elemental ratios. For the ¹³C¹⁸O-containing molecules the value of 30 000 is the reference. The zero-point energies of other isotopic substitutes do not pose any problem because the reactions involved in the interstellar chemical networks are significantly exothermic and the solutions of the chemical equations are independent of these quantities.

The isotopic fractionation reactions are introduced explicitly in the chemical network, whereas the other reactions involving isotopologues are built automatically from the reactions involving the main isotope in the chemical code. The adopted method has first been presented in Le Bourlot et al. (1993), where statistical arguments were used to derive the various branching ratios in the chemical reactions. The procedure is limited to three carbon-containing molecules (oxygen-containing molecules have a maximum of two oxygen atoms in our chemical network) and does not disitinguish between C¹³CC- or ¹³CCC-containing species. A similar approach has recently been applied by Röllig & Ossenkopf (2013) for photon-dominated region models. However, Röllig & Ossenkopf (2013) used the old (LGFA84) exothermicity values. We also explicitly introduce the relation given by Langer et al. (1984) that k_f + k_r = k_T. The forward reaction k_f and reverse reaction k_r rate coefficients involved in the isotopic exchange reaction are then evaluated from the total rate coefficient k_T as follows (27)and (28)These expressions have also been included in the study of fractionation in diffuse clouds presented by Liszt (2007).

The results summarized in Table 6 show that CO/¹³CO has the elemental value, whereas rarer isotopologues are very slightly depleted. The results for Models A and B are also very similar because no differences were used for the reaction rate coefficients between ¹³C⁺ and CO. However, more significant are the differences for the results for the isotopic ratio of HCO⁺, which directly arise from the variations of the exothermicities found in the present work. We also introduced a fractionation reaction for the deuterated isotope, whose rotational frequencies have been measured in the laboratory (Caselli & Dore 2005) and are detected in the interstellar medium (Guelin et al. 1982; Caselli et al. 2002). As the exothermicity of the deuterated isotopologues is somewhat higher, the isotopic ¹³C ratio is somewhat lower than in the hydrogenic counterpart.

The general trend seen in Table 6 is that the new Model B predicts lower fractional abundances x for H¹²C¹⁶O⁺ (up to 2%) and D¹²C¹⁶O⁺ (up to 5%), lower relative abundances R_B of H¹²C¹⁸O⁺ (7–21%), and higher relative abundances of the ¹³C-containing isotopologues (up to 40% for the hydrogenic forms and up to 75% for the deuterated forms) than Model A.

4.3. Reaction of N₂ with N₂H⁺

Molecular nitrogen is a homonuclear diatomic molecule with a X ground-electronic state with the three naturally occurring isotopologues: ¹⁴N₂, ¹⁴N¹⁵N, and ¹⁵N₂. Whereas ¹⁴N is a spin-1 boson, ¹⁵N is a spin-1/2 fermion, such that the two symmetric forms ¹⁴N₂ and ¹⁵N₂ follow different nuclear spin statistics. In the states with a higher nuclear spin degeneracy (ortho states), we have g = (I_N + 1)(2I_N + 1), whereas g = I_N(2I_N + 1) holds for the states with lower nuclear spin degeneracy (para states). To properly account for this effect, we evaluated the internal partition functions separately for even and odd J values, where Σ′ denotes summation in steps of 2. Multiplying each term by the appropriate nuclear spin (hyperfine) degeneracy factor, we obtain the partition function for N₂ as ¹⁴N¹⁵N is not a homonuclear diatomic molecule, such that (33)The equilibrium constants K_e and rate coefficients for the isotopic variants of N₂H⁺ reacting with N₂ are shown in Table 7. There we assumed the total rate coefficient k_T given by the Langevin collision rate (in SI units) (34)where e is the elementary charge, μ_R the reduced mass for the collision, and α(N₂) the polarizability of N₂ (α(N₂) = 1.710 ų, Olney et al. 1997), giving thus k_T = k_L = 8.11 × 10^-10 cm³ s^-1. The rate coefficients k_f and k_r are determined from k_T and K_e with the help of Eqs. (24) and (25), respectively. Spectroscopic parameters of Trickl et al. (1995) and Bendtsen (2001) were used for the states of ¹⁴N₂, ¹⁴N¹⁵N, and ¹⁵N₂.

The nuclear spin degeneracy affects the equilibrium constants for the reactions involving either ¹⁴N₂ or ¹⁵N₂. At higher temperatures, K_e in Table 7 approaches 1/2 for reactions D2 and D3 having as a product, and 2 for reactions D4 and D5 having as a reactant. For reaction D1, (35)the effects from nuclear spin statistics cancel out and as the temperature increases. From Eqs. (31) and (32), the ortho-to-para ratio is given by R₁₄ = 6Q_evenJ/ 3Q_oddJ for ¹⁴N₂ and by R₁₅ = 3Q_oddJ/Q_evenJ for ¹⁵N₂. We may note that R₁₄ assumes a value of 2.41 (2.01) and R₁₅ a value of 2.60 (2.99) at 5 K (10 K). At high temperature equilibrium, we have R₁₄ = 2 and R₁₅ = 3.

Adams & Smith (1981) employed normal nitrogen (ratio 2:1 of ortho vs. para ¹⁴N₂) in the SIFT experimental study of N₂H⁺+N₂. To measure the forward reaction and backward reaction rate coefficients at a given temperature, they interchanged the ion-source gas and reactant gas. Using mass-selected samples, these authors, however, were unable to distinguish between the isotopomers ¹⁴N¹⁵NH⁺ and ¹⁵N¹⁴NH⁺, such that their results provide the overall yield of these cations (no information on the relative yields). This applies to the competing reactions D2 and D3 on one side and the competing reactions D4 and D5 on the other side. ¹⁴N¹⁵NH⁺ and ¹⁵N¹⁴NH⁺ are expected to be differently fractionated (see Table 4).

To simulate the experimental conditions of Adams & Smith (1981), we introduced the overall forward k₂₃ and overall reverse k_-23 rate coefficients for reactions D2 and D3, and the overall forward k₄₅ and overall reverse k_-45 rate coefficients for reactions D4 and D5, Here we explicitly assumed an equilibrium distribution between ¹⁴N¹⁵NH⁺ and ¹⁵N¹⁴NH⁺. The term is the state-distribution normalization factor.

The variation of the rate coefficients with the temperature is displayed in Fig. 1. The common feature seen there is that the forward reaction becomes faster and the backward reaction slower with decreasing temperature. We also see that k_f and k_r exhibit a very weak temperature dependence for T> 50 K. For reactions D1 and D6, k_f and k_r approach the same value (k_L/ 2 in our model) at higher temperatures in Fig. 1a and Table 7. The high temperature limits of k_i and k_−i for i = 2−5 are, however, different because of the nuclear spin restrictions, as clearly seen in Figs. 1b,c. The reverse rate coefficients k_-2 and k_-3 in Fig. 1b become even higher than k₂ and k₃ for T> 14.7 K and T> 3.4 K, respectively, inverting thus the reaction direction. Herbst (2003) also found k_-3>k₃ at T = 10 K for reaction D3 assuming a different rate-coefficient model.

Fig. 1 Temperature dependence of the rate coefficients for N2H+ + N2 for reactions D1 and D6 in a), reactions D2 and D3 in b), and reactions D4 and D5 in c).

Fig. 2 Minimum energy paths for the linear approach of CO to HCO+/HOC+ in a), for the linear approach of N2 to HN in b) and for the formation of the mixed linear cluster ions N2··· HCO+ and N2H+··· OC in c). The coordinate displayed on the x-axis is shown with the dotted line in the chemical formulas and X = O, C.

For the overall state-averaged rate coefficients in Fig. 1 and Table 7, we have k₂₃>k_-23 and k₄₅>k_-45 for all temperatures shown. This is in accordance with the SIFT experiment of Adams & Smith (1981). The rate coefficients k_±23 and k_{± 45} appear 30% higher than the experimental finding, reported with an error of ±25% at 80 K. Note, however, that the ratios k₂₃/k_-23 and k₄₅/k_-45 agree within 6% with the corresponding experimental values. Due to the nuclear spin angular momentum selection rules, the high-temperature limits (for ) of k_{± 23} and k_±45 are different from the high-temperature limits of k_±1 and k_{± 6}.

From Eqs. (36)–(39) and the relationship of Eq. (11), we easily obtain (40)and (41)Following the procedure of Adams & Smith (1981), we may model the temperature dependence of the latter ratios as e ^ΔE_ij/k_BT (compare with Eq. (5)). Using our results from Table 7 for the overall forward and overall reverse rate coefficients calculated at the temperatures of the SIFT experimental study, T = 80 K and T = 292 K, we derive ΔE₂₃,ΔE₄₅ = 6.5 K for both reaction pairs.

Adams & Smith (1981) estimated the zero-point energy difference of 9 ± 3 K for reactions D2 and D4 (see Table 4). In accordance with the analysis presented here, we see, however, that the results of Adams & Smith (1981) should be attributed to the reaction pairs {D2, D3} and {D4, D5}. This also explains a large discrepancy seen in Table 4 between the theoretical estimates and experimental finding for reaction D4.

In recent studies of Bizzocchi et al. (2010, 2013), ¹⁴N¹⁵NH⁺ and ¹⁵N¹⁴NH⁺ were both detected in a prototypical starless core L1544 of low central temperature and an abundance ratio of 1.1 ± 0.3 was derived. Note that the ratio correlates with describing reaction D6 of Eq. (12). As seen in Table 7, we obtain of 1.22–1.02 for T = 40−292 K; the additional calculation at T = 30 K gave . Also note that the earlier model of Rodgers & Charnley (2004) has led to the ratio of 1.8–2.3, which correlates with our value of 2.27 (1.72) at T = 10 K (15 K).

4.4. Ionic complexes

Ion-molecule reactions were additionally examined using electronic structure calculations, carried out for the linear approach of the neutral CO and N₂ to the linear cations HCO⁺, HOC⁺, and N₂H⁺. The corresponding minimum-energy paths (MEPs) are displayed in Fig. 2. The MEPs are obtained optimizing three intramolecular distances for various monomer separations. Our calculations were performed at the CCSD(T)/aug-cc-pVTZ level of theory employing the standard MOLPRO and CFOUR optimization/threshold parameters.

The lower MEP in Fig. 2c is related to the reaction (42)which is considered to be the main destruction path for N₂H⁺ when CO is present in the gas phase at standard abundances [CO]/[H₂] ~ 10^-4 (Snyder et al. 1977; Jørgensen et al. 2004). For this reaction, Herbst et al. (1975) reported a rate coefficient of 8.79 × 10^-10 cm³ s^-1 at 297 ± 2 K. No reverse reaction was detected (Anicich 1993). For reactions involving HOC⁺, Freeman et al. (1987) measured a rate coefficient k of 6.70 × 10^-10 cm³ s^-1 for the following reaction (43)whereas Wagner-Redeker et al. (1985) reported k as 6.70 × 10^-10 cm³ s^-1 for (44)The Langevin collision rate is k_L = 8.67 × 10^-10 cm³ s^-1 for reactions (42) and (43) involving CO and k_L = 8.11 × 10^-10 cm³ s^-1 for reaction (44) involving N₂.

The common feature in Fig. 2 is the formation of a linear proton-bound ionic complex, which is 2000–7000 cm^-1 more stable than the separated monomers. The properties of the complexes are summarized in Table 8, where we give the geometric parameters r_i, the equilibrium rotational constants B_e, the harmonic wavenumbers ω_i for the main and deuterated isotopologues, and the harmonic zero-point energies E₀. The corresponding results for the constituent monomers are listed in Table 9. Note that the monomer values E₀ in Table 9 are harmonic and therefore different from the anharmonic results of Table 1. The coordinates r_i(i = 1−4) for A–B–H–C–D denote r₁ = r(A− B),r₂ = r(B− H),r₃ = r(H−C) and r₄ = r(C − D) in Table 8 and similar in Table 9. The dipole moments μ_z and the quadrupole moments Θ_zz in Table 8 and 9 are given with respect to the inertial reference frame with the origin in the complex centre of mass, where the position of the first atom A of A–B–H–C–D or A–B–C along the z axis is chosen to be the most positive.

The ionic complexes N₂HN and COHOC⁺ have linear centrosymmetric equilibrium structures. The complex OCH⁺···CO is asymmetric with a barrier height to the centrosymmetric saddle point OCHCO⁺(TS), seen at 358 cm^-1 in Fig. 2a. In the mixed-cluster ions, the proton is bound either to CO, when N₂··· HCO⁺ is formed, or to N₂, when N₂H⁺···OC is formed. Comparison of Tables 8 and 9 shows that the geometric parameters experience prominent changes (up to 0.01–0.02 Å) upon complexation. In this fashion, the ionic (molecular) complexes differ from van der Waals complexes, in which the monomers preserve their geometric parameters to a great extent.

The transformations in Fig. 2 are all of the proton transfer type. The neutral CO may approach H⁺ of the triatomic cation either with C or O since both C and O possess lone electron pairs. The proton attachment from the C side leads to a more stable complex. As seen in Fig. 2a, the complex OCH⁺···OC is 1785 cm^-1 above OCH⁺···CO and 9615 cm^-1 below COHOC⁺. We also see that N₂···HCO⁺ is 6996 cm^-1 more stable than N₂H⁺···OC. In all cases, the energy separation between the HCO⁺- and HOC⁺-containing complexes is smaller than the separation between free HCO⁺ and HOC⁺, seen to be 13 820 cm^-1 in Fig. 2. The results of Fig. 2 are consistent with the fact that the proton tends to localize on the species with higher proton affinity. The experimental proton affinity is 594 kJ mol^-1 (49 654 cm^-1) for CO on the C end and 427 kJ mol^-1 (35 694 cm^-1) for CO on the O end (Freeman et al. 1987). The experimental proton affinity of 498 kJ mol^-1 (41 629 cm^-1) was determined for N₂ (Ruscic & Berkowitz 1991).

The harmonic wavenumbers for the ionic complexes occurring in the course of reactions F1–F6 are provided in Table 10. In addition to the spectroscopic properties, we also give the harmonic zero-point energies of the complexes E₀, the reactants , and the products , as well as the dissociation energies including the harmonic zero-point energy correction in the direction of the reactants, , and in the direction of the products, , where D_e is the classical dissociation energy.

In Table 10, the vibrational mode ω₂, which is predominantly the diatom CO stretching vibration, is the most sensitive to isotopic substitutions. Compared with ω of free CO, ω₂ exhibits a blue-shift of 93 cm^-1 for the main isotopologue (Table 9 vs. Table 8). The modes ω₁ and ω₃, highly sensitive to the H→D substitution (Table 8), can be considered as the H-C-O stretching modes. The intermolecular stretching mode is ω₈. The zero-point-corrected dissociation energies in Table 10 are approximately 240 cm^-1 lower than the electronic dissociation energy of 4876 cm^-1 (Fig. 2). The harmonic ΔE^h/k_B values in Table 10 and anharmonic ΔE/k_B values in Table 2 agree within 0.5 K.

The proton-bound complexes OCH⁺··· CO and N₂··· HCO⁺ have large dipole moments μ_e of 1.16 ea₀ (2.94 D) and 1.39 ea₀ (3.53 D) (Table 8). For OCH⁺··· CO, the most intense infrared transitions are expected for ω₃ (with harmonic intensity of 2440 km mol^-1) and ω₁ ( km mol^-1), whereas the intermolecular stretch ω₈ has km mol^-1. The fundamental (anharmonic) transitions (ν₁,ν₂,ν₃,ν_4,5,ν_6,7,ν₈,ν_9,10) are calculated to be (2267, 2236, 1026, 1136, 346, 186, 208) for the main isotopologue (in cm^-1). The most intense infrared active transitions for N₂··· HCO⁺ are ω₁ ( km mol^-1), ω₃ ( km mol^-1), and ω₈ ( km mol^-1). For this complex, the fundamental vibrational (ν₁,ν₂,ν₃,ν_4,5,ν_6,7,ν₈,ν_9,10) transitions are determined to be (2357, 2321, 1876, 1045, 127, 186, 113) (in cm^-1). The anharmonic transitions are calculated from the cubic and semi-diagonal quartic force field in a normal coordinate representation by means of vibrational second-order perturbation theory, as implemented in CFOUR (Stanton et al. 2012).

Regarding the CCSD(T)/aug-cc-pVTZ method used here, we may note that our value of 358 cm^-1 in Fig. 2a for the barrier height of OC+HCO⁺ agrees reasonably well with previous theoretical results of 382 cm^-1 (the CCSD(T)/cc-pVQZ approach of Botschwina et al. 2001) and 398 cm^-1 (the CCSD(T)/aug-cc-pVXZ approach of Terrill & Nesbitt 2010 at the complete basis-set limit). A classical dissociation energy was previously determined to be 4634 cm^-1 for OCH⁺··· CO and 5828 cm^-1 for N₂HN at the complete basis-set limit (Terrill & Nesbitt 2010). The use of larger basis sets would ultimately be needed for converging theoretical results to stable values. Our primary goal here is the acquisition of first information relevant for the physical behaviour of the ionic complexes involving HCO⁺ and N₂H⁺. For these initial explorations of the PESs, the CCSD(T)/aug-cc-pVTZ approach is of satisfactory quality. A more detailed analysis of various basis-set effects, including the basis-set superposition error in systems with significantly deformed monomers, is being prepared and will be presented elsewhere.

5. Conclusion

Ion-molecule reactions are common in interstellar space, and investigating them helps to quantitatively understand the molecular universe (Watson 1976). We studied the isotope fractionation reactions of HCO⁺/HOC⁺ with CO and N₂H⁺ with N₂, as well as the linear proton-bound complexes formed in the course of these reactions. For OCH⁺+CO, we pointed out inaccuracies of previous exothermicity values that are commonly employed in chemical networks. The new exothermicities affect particularly prominently the rate coefficients derived at temperatures of dark interstellar cloud environments, which markedly changes the abundance ratios of the ¹³C- and ¹⁸O-containing formyl isotopologues.

The linear proton-bound cluster ions are found to be strongly bound (2000–7000 cm^-1). The ionic complexes OCH⁺··· CO and OCH⁺··· N₂ have sizeable dipole moments (2.9–3.5 D) and rotational constants of approximately 2000 MHz. If stabilized by means of collision and/or radiative processes, their high rotational population may facilitate the detection of these ions at low temperatures.

Acknowledgments

M.M. is grateful to Geerd H. F. Diercksen for sending her a copy of the MPI/PAE Astro 135 report. Marius Lewerenz is acknowledged for helpful discussions. Mila Lewerenz is thanked for helping with the literature search.

References

All Tables

All Figures

	Fig. 1 Temperature dependence of the rate coefficients for N2H+ + N2 for reactions D1 and D6 in a), reactions D2 and D3 in b), and reactions D4 and D5 in c).
In the text

	Fig. 2 Minimum energy paths for the linear approach of CO to HCO+/HOC+ in a), for the linear approach of N2 to HN in b) and for the formation of the mixed linear cluster ions N2··· HCO+ and N2H+··· OC in c). The coordinate displayed on the x-axis is shown with the dotted line in the chemical formulas and X = O, C.
In the text