Is the Gas-phase OH+H₂CO Reaction a Source of HCO in Interstellar Cold Dark Clouds? A Kinetic, Dynamic, and Modeling Study

The CRESU technique and the experimental setup employed have been previously described elsewhere (Antiñolo et al. 2016; Jiménez et al. 2015, 2016), and only a brief overview will be given below. A set of five Laval nozzles were used to get the jet temperatures between 22 and 107 K, and characterization of the obtained flows was discussed elsewhere (Canosa et al. 2016). The adiabatic gas expansion through the Laval nozzle from the reservoir to the reaction chamber provokes the cooling of the gas and forms a uniform jet in temperature and total gas density over several tens of cm, i.e., over several hundreds of μs hydrodynamic time, t_hydro (Table 3). OH radicals were generated in the gas jet by pulsed laser photolysis of gaseous H₂O₂ at 248 nm. The time evolution of OH radicals in the presence and absence of H₂CO was monitored by laser-induced fluorescence (LIF) ca. 310 nm, which was collected by a filtered photomultiplier tube. The excitation wavelength (282 nm) was delivered from a frequency-doubled dye laser pumped by the second harmonic of an Nd-YAG laser. The analysis of the exponential decays yields the pseudo–first-order rate coefficient, k', which is linearly related to [H₂CO] concentration under pseudo–first-order conditions (Appendix Figure 8). Varying the mass flow rate of a diluted mixture containing a well-known quantity of H₂CO and maintaining the OH-precursor concentration constant, k₁(T) was obtained from the slope of the plot of k' (or k'–k₀) versus [H₂CO] (Appendix).

A full-dimensional analytical PES has been developed to treat the reaction dynamics at temperatures below 100 K, as described in a recent publication (Zanchet et al. 2017). The intrinsic reaction coordinate (IRC) diagram in Figure 1 shows the quality of the analytical fit and describes the energy diagram of the reaction.

Zoom In Zoom Out Reset image size

QCT calculations have been performed using an extension of the code miQCT applied to N atoms (Dorta-Urra et al. 2015). The Hamilton equations are integrated using a step-adaptive Adams–Bashforth–Moulton predictor-corrector method (Shampine & Gordon 1975), using the Cartesian coordinates of the position vectors of all atoms. The total integral cross section, σ_vj(E), is calculated according to the following equation (Karplus et al. 1965),

$$\begin{eqnarray}&&{\sigma }_{\mathrm{vj}}(E)=\pi {b}_{\max }^{2}{P}_{{\rm{r}}}(E),\end{eqnarray} \tag{ 6 }$$

where b_max is the maximum impact parameter corresponding to the maximum distance at which trajectories are efficiently deviated to form products and P_r(E) is the reaction probability given by

$$\begin{eqnarray}&&{P}_{{\rm{r}}}(E)={N}_{{\rm{r}}}(E)/{N}_{\mathrm{tot}}(E),\end{eqnarray} \tag{ 7 }$$

where N_tot(E) is the maximum number of trajectories at collision energy E with b < b_max, while N_r(E) is the number of reactive trajectories leading to products.

A Monte Carlo sampling of initial conditions is done for each reagent independently (OH and H₂CO). To include the zero-point energy (ZPE), an adiabatic switching method is performed (Grozdanov & Sovolev 1982; Johnson 1987; Qu & Bowman 2016). In this method, a random selection of initial conditions is done according to analytical normal modes, and then the harmonic potential is adiabatically switched to the real one. Since the energy spreading is rather high, we select a particular trajectory with the closest energy to the anharmonic ground state eigenvalue of each fragment, which is propagated for 100,000 steps. In this case, OH and H₂CO are initially set at zero angular momentum by using the method of Nyman and coworkers (Nyman et al. 1990). Once the internal vibrational distribution is set, the relative orientation of the two fragments is randomly selected. The impact parameter is randomly selected up to a maximum distance ${b}_{\max }^{0}$, at which trajectories are efficiently deviated. This maximum value is given by the capture model (Levine & Bernstein 2001):

$$\begin{eqnarray}&&{b}_{\max }^{0}={3}^{1/2}{(A/2E)}^{1/3}.\end{eqnarray} \tag{ 8 }$$

The initial distance between the two reagents is set to 120 au. These long distances are required for the low energies considered here.

In the present system, the long-range interaction between the permanent dipole moments (d_X, X = OH and H₂CO) varies as:

$$\begin{eqnarray}&&V(R\to \infty )=\tfrac{-A}{{R}^{3}}\end{eqnarray} \tag{ 9 }$$

with A = 4 d_OH ${d}_{{{\rm{H}}}_{2}\mathrm{CO}}$.

Two different strategies are followed to analyze the collision dynamics.

• (1)  
  First, the individual cross section for H₂CO and OH in the ground vibrational and rotational states is calculated in the microcanonical ensemble according to Equation (6). This is done using a grid of collision energies from 0.1 to 1000 meV. For each translational energy, 10⁵ trajectories were run. Finally, the state-dependent rate, k_{v = 0,j = 0}(T), is obtained by numerical integration over a Boltzmann distribution. The final rate coefficient ${k}_{v=0,j=0}^{{\rm{e}}}$(T) is obtained by multiplying this k_{v = 0,j = 0}(T) by the electronic partition function, q_e, described in Section 3.
• (2)  
  Second, the thermal rate coefficients, k(T), are obtained from Equation (10) by running 10⁶ trajectories for each temperature in the macrocanonical ensemble, defining initial conditions for thermal rotational distributions of the two reagents,

  $$\begin{eqnarray}&&k(T)={q}_{e}{\left(\tfrac{8{kT}}{\pi \mu }\right)}^{\displaystyle \frac{1}{2}}\pi {b}_{\max }^{2}{P}_{r}(T),\end{eqnarray} \tag{ 10 }$$

where μ is the reduced mass of the OH + H₂CO system.

To evaluate the hypothesis that H₂CO could be a precursor of HCO radicals in cold interstellar clouds, we have run a chemical model of a cold dark cloud in which we adopt typical parameters, i.e., a kinetic temperature of 10 K, a volume density of H nuclei of 2 × 10⁴ cm⁻³, a visual extinction of 30 mag, a cosmic-ray ionization rate of H₂ of 1.3 × 10⁻¹⁷ s⁻¹, and the so-called "low-metal" elemental abundances (Agúndez & Wakelam 2013). We have run models using two chemical networks to have an idea of the robustness of the results. Concretely, we have used the network kida.uva.2014 from the KIDA (KInetic Database for Astrochemistry) database (Wakelam et al. 2015) and the rate 2012 network from the UMIST (University of Manchester Institute of Science and Technology) database (McElroy et al. 2013). In both chemical networks, we have implemented the reaction ${{\rm{H}}}_{2}\mathrm{CO}+\mathrm{OH}\to \mathrm{HCO}+{{\rm{H}}}_{2}{\rm{O}}$ with the extrapolated rate coefficient at 10 K obtained in this work (see Section 3.1).

The measured rate coefficients k₁(T) have been compiled in Table 2. They have been determined at ca. 22 K using three different Laval nozzles with total He densities ranging from 3.37 × 10¹⁶ to 1.67 × 10¹⁷ cm⁻³. At ca. 45 K, two Laval nozzles were employed for different gas densities, 6.90 × 10¹⁵ and 1.02 × 10¹⁷ cm⁻³(N₂ and N₂/He, respectively). All physical conditions for all Laval nozzles employed are listed in Table 3. In Table 2, k₁(T = 22 K) and k₁(T = 45 K) are presented at each gas density. As can be seen, no pressure dependence of k₁(T = 22 K) was observed in the investigated gas-density range, implying that reaction (1) at these temperatures is either a bimolecular reaction or a three-body process already in the high-pressure limit at the lowest gas density. At 45 K, a decrease of 27% is observed in k₁ when the gas density is varied by more than one order of magnitude. Within the experimental uncertainties, this decrease is not expected to be originated by a pressure effect. Such an effect on k₁(T) should be more noticeable at even higher temperatures for which the rate coefficient is lower and then separates even more from the high-pressure limit of a falloff curve. In other words, the high-pressure limit would be reached at higher total pressures. At the highest temperatures investigated in this work (107 K), it can be seen that the rate coefficient is similar to that at 89 K (see Table 2), whereas the total density (1.8 × 10¹⁷ cm⁻³) is 3.6 times higher. If the process was thermolecular, one should expect that k₁ would increase with gas density and when the temperature is lowered. Hence, we should have found a significantly higher value at 89 K than at 107 K, while it was found, within the experimental uncertainties, k₁(89 K) ≅k₁(107 K). For all these reasons, it can be argued that reaction (1) is essentially bimolecular in the temperature range of the present work. As shown in Table 2, an increase in k₁(T) is observed when lowering T.

When combining our data with those reported by Sivakumaran et al. (2003) at T < 300 K, the best fit is

$$\begin{eqnarray}&&{k}_{1}(T=22\mbox{--}300{\rm{K}})=(7.73\pm 1.82)\times {10}^{-12}\\ &&\quad \times {(T/300{\rm{K}})}^{-(1.03\pm 0.09)}\,{\mathrm{cm}}^{3}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 11 }$$

where the uncertainties are ±2σ. The relative difference of the calculated k₁(T) from Equation (11) with respect to the experimental values in the range 22–107 K varies from 4% to 21%. In Figure 9, the confidence-interval bands of the fit are depicted at a 95% level of confidence in addition to the experimental data and the abovementioned fit. The k₁(T) values in the 200–300 K range are predicted by Equation (11) within an uncertainty of 1%–5% for most of the available measurements.

As can be seen in Figure 2, when comparing our results with all of the previous kinetic studies on the OH + H₂CO reaction between 200 and 1500 K, a U-shaped curvature is observed in the trend of k₁(T), with a minimum value around 300 K. The three-parameter expression recently proposed by Wang et al. (2015), which combines their high-temperature data with those at lower temperatures from Sivakumaran et al. (2003), is shown by a blue solid line in Figure 2. This expression represents perfectly experimental data for T > 200 K (within 6%–15%) but fails to predict k₁(T) at T < 120 K (difference >50%). For example, the predicted k₁(T = 107 K) is almost twice the experimental value. Therefore, we recommend using our expression for k₁(T) below 250 K and Wang's formula above. Our expression has been derived in the range 22 K < T < 250 K. Even though extrapolation can be risky, Equation (11) can be used to obtain an estimation of the rate coefficient at 10 K needed for the modeling of HCO and H₂CO abundances in cold interstellar clouds. This leads to a value of (2.60 ±0.60) × 10⁻¹⁰ cm³ s⁻¹.

Zoom In Zoom Out Reset image size

The experimental determination of the reaction rate coefficient is done through the disappearance of the OH. One way to check if the measured disappearance of OH corresponds to the formation of the HCO product is to perform complete dynamical simulations of the ${{\rm{H}}}_{2}\mathrm{CO}+\mathrm{OH}\to \mathrm{HCO}+{{\rm{H}}}_{2}{\rm{O}}$ reactive collisions. This reaction presents a low barrier of 27 meV, with a ZPE of the order of 1 eV (see Figure 1), in agreement with previous studies. If pure TST is applied and the reaction coordinate is separated from the orthogonal vibrational modes, the reaction would not happen for collision energies below 27 meV (approximately 300 K), and the rates would tend to zero when decreasing temperature, in opposition to the behavior of the experimental results. One way to overcome this difficulty is to introduce tunneling in the TST model, giving rise to different approaches based essentially in the anharmonicity of the potential at the saddle point. These approaches are approximate, and the purpose of this work is to carry out full-dimension reaction dynamics studies that better describe the reaction mechanism at low temperatures.

QCT calculations have been done for the ${{\rm{H}}}_{2}\mathrm{CO}+\mathrm{OH}\to \mathrm{HCO}+{{\rm{H}}}_{2}{\rm{O}}$ reaction on the reactive force field (RFF) and the full (RFF+MB (many-body)) PESs described in Zanchet et al. (2017). These PESs were developed on purpose and properly describe the long-range interactions and all the energetics, especially the stationary points of the PES. The calculations performed in the microcanonical ensemble at fixed collision energy allow a better understanding of the reaction mechanism. The total reaction cross section, σ_{v = 0,j = 0}(E), calculated for the ground rotational and vibrational state of the two reactants shows a continuous increase with decreasing collision energy for both the total RFF+MB and RFF PESs (see bottom panel of Figure 3). Clearly, the increase of σ_vj is associated with the increase of b_max (middle panel of Figure 3), while P_r(E) remains nearly constant (top panel of Figure 3). The increase in the impact parameter is well described by the capture model (Levine & Bernstein 2001) using Equation (8). The agreement with the one found in the QCT simulations is nearly perfect until an energy of 100 meV, where the corresponding b_max is of the order of the size of the molecules. This confirms that the potential is able to capture trajectories at low energy with an increasing efficiency as collision energy decreases. This leads to trajectories that are trapped in the OH–H₂CO well for a long time, where the dynamics becomes ergodic and chaotic, as recently found in ultracold collisions (Croft & Bohn 2014).

Zoom In Zoom Out Reset image size

In the top panel of Figure 3, it is shown that P_r(E) changes smoothly with the collision energy and does not drop to zero even at 0.1 meV, energy well below the reaction barrier of the two PESs (of 27 meV). The difference between the reaction cross section found for the two PESs, RFF and RFF+MB, is explained by the change in the reaction probability, since the impact parameter remains nearly the same for both PESs. This change of the reaction probability is attributed to the many-body term introduced in the potential to better describe the calculated ab initio points. The RFF by itself is too isotropic, facilitating access from the complex to the saddle point, giving rise to an overestimation of the reaction cross section. This fact is corrected by the MB term.

The remaining point to understand is why the reaction probability is nonzero at energies below the top of the barrier. Considering the ZPE of reactants, the system has enough total energy to overcome this barrier. If the reactants' ZPE remains in the coordinates orthogonal to the reaction coordinate, the reaction could not take place, since there would not be enough energy in the reaction coordinate. Clearly, this is not the case, and there must be couplings between the reaction coordinate and the remaining degrees of freedom producing an energy transfer that explains the nonzero reaction probability even at 0.1 meV of collision energy. Such a situation is typical in non-intrinsic reaction coordinate (non-IRC) trajectories (Lourderaj et al. 2008), since they explore large regions of the configuration space far from the minimum energy path determined by the IRC, where the anharmonicity is larger. The concept of ZPE at the transition state is an approximation based on the separation between the reaction coordinate (in the dissociation continuum) and the orthogonal vibrations. This separation is adequate for trajectories going along the IRC, since the potential can be locally described by a parabolic separable potential. Out of this IRC the potential is no longer separable, hence the energy can be transferred among the different degrees of freedom. Moreover, the harmonic approximation made to estimate the ZPE at the saddle point is questionable, since the PES at the saddle point has a large anharmonicity that tends to reduce the ZPE.

The rate coefficient is calculated by integrating the cross section over a Boltzmann distribution. To account for the spin–orbit splitting of the OH(²Π) states into the F₁(²Π_1/2) and F₂(²Π_3/2) with an energy separation of 200.3 K, the rate coefficient is multiplied by the electronic partition function, q_e = [1+exp(−200.279/T)]⁻¹. This term arises from the assumption that only the ground spin–orbit state reacts. The rate coefficient thus obtained only corresponds to the ground rotational states of OH and H₂CO, and it is labeled as k^e_{v = 0,j = 0}(T) in Figure 4. In order to compare with the experimental data of this work, also in Figure 4, we have to consider the rotational populations of the reactants, which vary as a function of temperature. This is done by directly computing the thermal rate coefficient, k(T), in the macrocanonical ensemble according to Equation (10), which is depicted in Figure 4. Clearly, the computed rate coefficient k(T) is slightly lower than the experimental one, but it shows the same behavior. At T > 50 K, k(T) is a factor of ca. 2 lower; however, the difference between theory and experiment increases somewhat when lowering the temperature. We therefore conclude that the behavior of the simulation is semiquantitatively correct, explaining the experimental results obtained in this work.

Zoom In Zoom Out Reset image size

The quantitative differences may be attributed to two major factors. First, the reaction probabilities displayed in Figure 3 depend on the accuracy of the PES. Work is now in progress to improve the current PES. Second, quantum effects are expected to be important at low collision energies, while they are neglected in the present treatment. To analyze the qualitative effect of quantum dynamics, we can refer to a similar reaction, $\mathrm{OH}+{\rm{F}}\to {\rm{O}}+\mathrm{HF}$, with a similar energy profile, for which the cross section also increases significantly at low energies (Gómez-Carrasco et al. 2004a, 2004b, 2005). In that triatomic system, exact quantum wave-packet calculations were performed, finding that the reaction cross section increases with decreasing collision energy, being always larger than the QCT cross sections. This behavior was interpreted by an indirect mechanism at low energies due to the presence of resonances in both the quantum and classical simulations. It may be concluded that the QCT approach yielded qualitatively good results but tended to lead to lower rate coefficients at low temperatures. Similar effects are expected in the H₂CO + OH reaction, and the inclusion of quantum effects may yield to a better agreement with experiments, especially at low temperatures. Quantum effects may be the reason for the increasing difference, particularly the ZPE effects. The reaction mechanism consists of two steps: the capture to form the complex and the dissociation to either reactants or products. The lifetime of this complex is determined by the number of open channels. This number in the products is already rather large, since the reaction is exothermic. Then, the crucial quantity determining the reactivity is the number of open channels in the reactants. For low temperatures, this number reduces to 1, i.e., the ground rovibrational state of the two reactants. This is so in a quantum description. However, using a QCT description, there is no quantization. Therefore, the dissociation back to reactants is somehow overestimated, and this probably explains the disagreement between simulations and experiments, especially at low temperatures.

The total rate coefficient corresponds to the rearrangement channel HCO + H₂O. It should be noted that the products end in highly excited rovibrational states. In the particular case of HCO, with a rather weak bond, it may fragment, giving H + CO. Thus, the formation of CO + H + H₂O is included in the present model. However, its proportion is rather low, varying between 2% and 4% in the entire temperature range.

In order to evaluate the impact of the titled reaction on the abundances of HCO and H₂CO in cold dark clouds, we performed a chemical model using two different chemical networks, UMIST (McElroy et al. 2013) and KIDA (Wakelam et al. 2015), and set k₁(10 K) to 2.6 × 10⁻¹⁰ cm³ s⁻¹ (value from Equation (11)).

In Figure 5, we show the calculated fractional abundances of HCO and H₂CO as a function of time using the UMIST and KIDA networks. It is seen that the calculated abundance of HCO at times between 10⁵ and 10⁶ yr is in good agreement with the range of values derived from observations. Even though the reaction between H₂CO and OH is very fast at low temperatures, the chemical model indicates that it accounts for a few percent in the production of HCO radicals in cold dark clouds. This result is consistently obtained no matter which chemical network is adopted. According to the chemical model, the main HCO-forming routes are reactions involving atomic oxygen (reactions (2) and (3)) and the dissociative recombination of the molecular ion H₂COH⁺ (reaction (4)). The relative importance of these reactions varies depending on the adopted chemical network. This is illustrated in Figure 6, where we show the contributions of the five most important HCO-forming reactions as a function of time and compare them with the HCO production rate given by reaction (1). At times between 10⁵ and 10⁶ yr, reaction (2) dominates when the UMIST network is adopted, while if the KIDA network is used, the main route to HCO becomes reaction (4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The rate coefficient of the reaction between CH₂ radicals and O atoms has been measured at room temperature (Böhland et al. 1984), although it has not been studied at low temperature and it is not clear which is the branching ratio of the channel leading to HCO. According to Böhland et al. (1984), the HCO radicals formed in the O + CH₂ reaction should contain enough excess energy to dissociate into CO + H. It is therefore uncertain whether reaction (2) is actually an important source of HCO radicals. The UMIST network adopts a rate coefficient of 5 × 10⁻¹¹ cm³ s⁻¹ estimated in high-temperature experiments, while in the KIDA network a smaller value, 2 × 10⁻¹² cm³ s⁻¹, is adopted based on a discussion by Loison et al. (KIDA datasheet 281; Wakelam et al. 2015). The importance of reaction (2) strongly depends on its rate coefficient at low T. This explains the big difference between the UMIST and KIDA results shown in Figure 6. For ${\rm{O}}+{{\rm{C}}}_{2}{{\rm{H}}}_{4}^{+}$ (reaction (3)) and H₂COH⁺ + e⁻ (reaction (4)), the kinetics is better constrained from experiments. For the first reaction, its rate coefficient has been measured to be 2.4 × 10⁻¹⁰ cm³ s⁻¹ at room temperature, and the branching ratio for the channel yielding HCO has been measured as 45% (Scott et al. 2000). Since the neutral reactant is nonpolar, the rate coefficient should be independent of temperature, and thus the value measured at room temperature should pertain at low temperatures as well. Therefore, the value adopted at 10 K in both networks was k₃ = 1.08 × 10⁻¹⁰ cm³ s⁻¹.

The dissociative recombination of H₂COH⁺, as with most reactions of this type, is fast at low temperatures and has been found to preserve the C–O bond in most of the collisions (Hamberg et al. 2007). The branching ratios of the channels leading to H₂CO, HCO, and CO are not known. However, they are assumed to be equal by the KIDA and UMIST networks. The values for k₄ used were 1.10 × 10⁻⁶ cm³ s⁻¹ (KIDA) and k₄ = 3.08 × 10⁻⁶ cm³ s⁻¹ (UMIST).

The chemistry of HCO radicals in prestellar cores was discussed by Bacmann & Faure (2016) based on simple steady-state arguments and different assumptions on the main source and sink reactions for HCO. These authors evaluated which HCO abundance would result if this radical was exclusively formed (i) by reaction (1) or (ii) by the dissociative recombination of H₂COH⁺ (reaction (4)), in both cases assuming that the main sinks of HCO are reactions of proton transfer from abundant ions that can easily give a proton (e.g., HCO⁺, ${{\rm{H}}}_{3}^{+}$). They concluded that the ionic pathway scenario (reaction (4)) results in HCO abundances that agree with the values derived from observations, while the neutral–neutral reaction OH + H₂CO needs to have a rate coefficient of 4 × 10⁻¹⁰ cm³ s⁻¹ to produce HCO in the observed quantities. This is rather close to the value used in the present model extrapolated from our experimental results, k₁(10 K) = 2.6 × 10⁻¹⁰ cm³ s⁻¹. According to the chemical model (see Figure 6), reaction (1) would need to be more than one order of magnitude faster than indicated by our measurements to compete with other reactions in the synthesis of HCO radicals in cold dense clouds. We stress, however, that this conclusion is based on a chemical model that is subject to significant uncertainties because the low-temperature kinetics of most of the reactions included is not precisely known, as exemplified earlier with, e.g., the O + CH₂ reaction.

The formyl radical is particularly abundant in PDRs, where its abundance is enhanced by two orders of magnitude with respect to cold dense clouds (Gerin et al. 2009). It is, however not clear how this radical is formed in these environments. Gerin et al. (2009) carried out a pure gas-phase model and found that HCO abundances of the order of magnitude of the observed ones can be obtained if HCO is efficiently produced in the photodissociation of H₂CO and in the reaction between CH₂ radicals and O atoms (reaction (2)). These two routes, however, lack sufficient empirical support, and it remains to be seen if alternative HCO-forming routes—as, e.g., reaction (1)—would need to be invoked to explain the relatively large abundances of HCO observed in PDRs such as the Horsehead Nebula.

In summary, despite the uncertainties in the rate coefficients and branching ratios of various reactions important for the chemistry of HCO, it seems that, in the present state of knowledge and chemical uncertainties, in cold dense clouds the reaction OH + H₂CO is less efficient than alternative gas-phase routes involving O atoms and the precursor ion H₂COH⁺ in the formation of HCO radicals. However, it remains to be explored whether the dramatic enhancement in the reactivity of OH and H₂CO could be responsible for the formation of HCO radicals in other relatively cold interstellar media, such as diffuse clouds and PDRs.