Experimental and theoretical study of the low-temperature kinetics of the reaction of CN with CH₂O and implications for interstellar environments

Abstract

Rate coefficients for the reaction of CN with CH₂O were measured for the first time below room temperature in the range 32–103 K using a pulsed Laval nozzle apparatus together with the Pulsed Laser Photolysis–Laser-Induced Fluorescence technique. The rate coefficients exhibited a strong negative temperature dependence, reaching (4.62 ± 0.84) × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ at 32 K, and no pressure dependence was observed at 70 K. The potential energy surface (PES) of the CN + CH₂O reaction was calculated at the CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ level of theory, with the lowest energy channel to reaction characterized by the formation of a weakly-bound van der Waals complex, bound by 13.3 kJ mol⁻¹, prior to two transition states with energies of −0.62 and 3.97 kJ mol⁻¹, leading to the products HCN + HCO or HNC + HCO, respectively. For the formation of formyl cyanide, HCOCN, a large activation barrier of 32.9 kJ mol⁻¹ was calculated. Reaction rate theory calculations were performed with the MESMER (Master Equation Solver for Multi Energy well Reactions) package on this PES to calculate rate coefficients. While this ab initio description provided good agreement with the low-temperature rate coefficients, it was not capable of describing the high-temperature experimental rate coefficients from the literature. However, increasing the energies and imaginary frequencies of both transition states allowed MESMER simulations of the rate coefficients to be in good agreement with data spanning 32–769 K. The mechanism for the reaction is the formation of a weakly-bound complex followed by quantum mechanical tunnelling through the small barrier to form HCN + HCO products. MESMER calculations showed that channel generating HNC is not important. MESMER simulated the rate coefficients from 4–1000 K which were used to recommend best-fit modified Arrhenius expressions for use in astrochemical modelling. The UMIST Rate12 (UDfa) model yielded no significant changes in the abundances of HCN, HNC, and HCO for a variety of environments upon inclusion of rate coefficients reported here. The main implication from this study is that the title reaction is not a primary formation route to the interstellar molecule formyl cyanide, HCOCN, as currently implemented in the KIDA astrochemical model.

---

1 Introduction

The current list of rate coefficients and product branching ratios determined experimentally for reactions at very low temperatures, typical of astrochemical environments, is still relatively sparse. Theoretical calculations of these kinetics parameters require knowledge of the potential energy surface (PES), in particular the energy of transition states and intermediates along the reaction coordinate. At very low temperatures, only reactions with small or no activation barriers are likely to proceed quickly. Calculated rate coefficients are highly sensitive to small changes to the PES, which can be on the order of the uncertainties of the ab initio theory being used (estimated to be 3.0–4.5 kJ mol⁻¹ at the CCSD(T)/aug-cc-pVTZ level of theory¹), making kinetics predictions less reliable.² Extrapolation of parameterisations to very low temperatures determined from rate coefficients at higher temperatures is fraught with difficulty, and has been shown to be incorrect by orders of magnitude.³ In some cases, this disagreement can be due to a change in the chemical mechanism which governs the temperature dependence of rate coefficients, k(T), at low temperatures.^3,4

There is a desire to understand what aspects of the reaction potential energy landscape control the behavior of the reaction rate coefficient and the branching to products as a function of temperature. As discussed in a recent perspective,³ both PESs with small barriers to reaction and/or submerged barriers to reaction can cause sharp increases in rate coefficients at low temperatures until in some cases the collision limit is reached. For the case of small barriers to reaction under low temperature reaction conditions, weakly bound complexes in the entrance channel for reaction are formed with less energy compared with higher temperatures. This causes their lifetime for dissociation back to reactants to become sufficiently long that quantum mechanical tunneling through reaction barriers to products can become competitive, leading to a dramatic increase in the rate coefficient with decrease in temperature.³ Two examples are the reaction of OH with Complex Organic Molecules (COMs) such as CH₃OH,^3,5–8 and the reaction of C(³P) with H₂O.⁹ For the case of submerged barriers, an increase in the rate coefficient with a decrease in temperature is observed consistent with a negative activation energy to reaction. Two examples are the reaction of C₂H with O₂,¹⁰ and the reaction of O(¹D) with CH₄.¹¹ If a calculated PES is used in conjunction with kinetic models to predict k(T), any uncertainty in the height or width of relatively small potential energy barriers can lead to large uncertainties in the prediction of reaction rate coefficients. Therefore, these reaction-specific PES attributes are important to closely examine alongside robust experimental data spanning a broad temperature range, especially including low temperatures.

The current study focuses on the low temperature reaction of the cyano radical with formaldehyde, CN + CH₂O (Reaction 1), which was suggested by astrochemists to be the primary formation route to the interstellar molecule formyl cyanide, HCOCN,¹² a suggestion reinforced recently by the ab initio quantum calculations of the PES and master equation calculations of the rate coefficients by Tonolo et al.¹³ Chemical intuition and results from past studies of similar reactions^14–16 suggest, however, that CN will likely attack CH₂O via a hydrogen abstraction mechanism, forming hydrogen cyanide (HCN) and the formyl radical (HCO). Both reactants and predicted hydrogen abstraction products have already been observed in several astrochemical environments, such as Asymptotic Giant Branch (AGB) stellar winds,^17,18 the InterStellar Medium (ISM),^19–26 and dark clouds,^27–30 providing further motivation for the current study at low temperatures. Additionally, the CN + CH₂O reaction is relevant to combustion processes,³¹ planetary atmospheres,³² and Titan's atmosphere.^33,34 The inclusion of new kinetic data for neutral–neutral reactions at very low temperatures can make a significant difference for the abundance of key species in interstellar environments.^3,35

Previous experimental work on the CN + CH₂O reaction measured overall rate coefficients between 297–673 K by Yu et al.³⁶ and 294–769 K by Chang and Wang³⁷ (for ease of reading these two works are referred to from now on as YCW^36,37) using the Pulsed Laser Photolysis–Laser-Induced Fluorescence (PLP–LIF) technique in resistively-heated flow cells. In these studies, k₁(T) was found to decrease slightly with decreasing temperature, from approximately 4 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ near 770 K to approximately 1.5 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ at room temperature. The uncertainties quoted for these studies were on average 2.8% for Yu et al.,³⁶ and 1.2% for Chang and Wang.³⁷ It was predicted that a hydrogen abstraction process, forming HCN, was likely to dominate the chemical mechanism. Indeed, a prior theoretical study (QCISD/6-31G**//UHF/6-31G**) found only a very small (∼2.7 kJ mol⁻¹) barrier to this H atom abstraction reaction.³⁸ In comparison, a very recent calculation by Tonolo et al.¹³ (CCSD(T)/CBS+CV) reports a small submerged barrier for the abstraction pathway forming HCN. Further aspects of the reaction potential energy surface were explored in this calculation, with particular attention paid to addition reaction pathways forming CN–CH₂O adducts, with bonds formed between the C of formaldehyde and either end of CN. However, the high-pressure limit rate coefficient (2 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹) calculated using transition state theory and based on this new PES¹³ is close to an order of magnitude larger than the previously measured room temperature studies of YCW. Exothermic product channels for the reaction of CN with CH₂O include:^13,39,40

Enthalpies for reactions (1a)–(1c) are taken from Ruscic and Bross 2019;³⁹ for reaction (1d) from Born et al.;⁴⁰ but for reaction (1e) the zero-point energy corrected energy change calculated by Tonolo et al.¹³ is given. In the work presented here, low-temperature reaction rate coefficients were measured for the first time below 294 K for the CN + CH₂O reaction using the PLP–LIF technique in a pulsed Laval nozzle apparatus. Using CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ, an ab initio PES was calculated for CN + CH₂O and utilized in the Open Source MESMER software package⁴¹ to calculate temperature-dependent rate coefficients and product branching ratios. Using a fitting procedure, the parameters within MESMER were optimized to give the best agreement with experimental data over the range 32–769 K. These MESMER rate coefficients were then used to develop a parameterisation via a modified Arrhenius (MA) equation and extrapolated to 4 K. The fits of the rate coefficients were then incorporated into astrochemical simulations of cold dense clouds, AGB stars, and hot molecular cores.

2 Methods

2.1 Experimental measurement of rate coefficients at low temperatures

Rate coefficient measurements for the reaction of CN with CH₂O were performed in a pulsed Laval nozzle apparatus using the PLP–LIF technique, as shown schematically in Fig. 1.

Fig. 1 Schematic of the pulsed Laval nozzle apparatus and the PLP–LIF technique. Reprinted and adapted with permission from D. E. Heard, AcChR, 2018, 51, 2620³. Copyright 2018 American Chemical Society.

The generation of uniform, thermalized, cold, gas flows for kinetics measurements with this apparatus was described previously in detail and will only be described briefly here.^14,42–47 Gaseous formaldehyde (CH₂O) is a difficult reagent to work with and particular care was given to its preparation from paraformaldehyde. Preparation of mixtures with the bath gas, its handling, and quantification of its concentration are necessary for accurate kinetics studies. Generation of the CH₂O reagent gas was performed by controllably heating paraformaldehyde powder (Sigma-Aldrich, 95%) in an evacuated 500 mL glass bottle (Duran) with a heat gun (Steinel, model HL1810S) to ∼70 °C. During heating of paraformaldehyde, the nascent CH₂O gas was then passed through a −10 °C cold trap. The cold trap was submerged in ethanol (VWR, 99.96%) which was chilled with a refrigerated immersion probe (LaPlant, model 100CD). The cold trap was utilized to trap any water or other condensable byproducts generated during the heating of paraformaldehyde. The purified CH₂O was then allowed to flow into evacuated cylinders until reaching a pressure of ∼200 torr (∼26.7 kPa). The cylinders were then filled with argon (BOC, 99.998%) or nitrogen (BOC, 99.998%) to ∼6 atm (∼608 kPa) total (generating ∼4.4% mixtures of CH₂O/bath gas). The gas in the cylinders was then allowed to mix for >12 hours. This method of generating CH₂O is similar to methods utilized previously by our group and other groups.^14,48,49 To generate the cyano radical, CN, the vapor from the solid precursor cyanogen iodide (ICN, Acros Organics, 98%, vapor pressure ∼1 torr (∼133 Pa) at 298 K.^50,51) was entrained in ∼2 atm (∼203 kPa) of Ar or N₂ bath gas.

The controlled mixing of reagent gases was accomplished through the use of Mass Flow Controllers (MFC) (MKS, type 1179A) such that final mixtures of gases were ∼0.1–1% CH₂O, ∼0.005% ICN, and ∼99% Ar or N₂ bath gas with a total flow rate of ∼2–5 slm. The gases were then allowed to mix in a mixing manifold prior to being pulsed through the Laval nozzle. Since CH₂O was found to slowly repolymerize, depositing solid paraformaldehyde in the gas lines and causing the MFC calibration to slowly drift, the concentration of CH₂O was determined via UV absorption measurements by sampling gas from the tubing between the mixing manifold and the pulsed valves for each concentration of CH₂O utilized in a kinetics experiment.¹⁴ Measurements of the concentration of CH₂O were made using a custom-made 1 m path length UV/Vis absorption cell filled to ∼1.2 atm (∼122 kPa) with the gas-mixture as measured by a capacitance manometer (MKS, 0–100 PSIA (∼689 kPa)). The light source for absorption measurements was a UVB lamp (EXOTERRA, UVB200) with continuous output around ∼290–350 nm. Absorption measurements were performed with a UV/Vis spectrometer (Ocean Optics, HR4000CG-UV-NIR) with 0.75 nm resolution. Absorption spectra were integrated for ∼2 seconds and 4 spectral traces were averaged. Representative averaged UV absorption spectra for CH₂O are shown in Fig. S1 in the ESI.† Using concentrations determined from this spectrometer, a typical half-life of the generated ∼4.4% CH₂O in an N₂ bath in the cylinders was determined to be approximately 4 days.

After the mixing manifold, the gas mixtures were sent to two solenoid valves (Parker, Series 9) to be pulsed at 5 Hz into the pre-expansion region upstream of the Laval nozzle. Gas in the pre-expansion reservoir underwent a controlled expansion through a custom-made, axisymmetric, converging-diverging Laval nozzle such that the flow into the vacuum chamber at 0.3–2 torr (40–267 Pa) was a wall-less reactor. The vacuum chamber was continuously evacuated by two Roots blowers operating in parallel: a Roots blower (Leybold RUVAC 251) backed by a rotary pump (Leybold D65B) and a Roots blower (Edwards EH250) backed by a rotary pump (Edwards ED660) and the pressure in the vacuum chamber was monitored by a capacitance manometer (Leybold, type CTR90, 0–10 torr (0–1.3 kPa)). Flow temperatures between 32–103 K were achieved by switching between Laval nozzles of Mach numbers between 2.49 and 5.00 and/or switching between Ar and N₂ bath gases. Impact pressure measurements using a Pitot tube were utilized to determine the Mach number along the flow, and via the Rayleigh equations, the density and rotational–translational temperature of the flow,⁴ as illustrated for temperature in Fig. S2 (ESI†).

In order to initiate reactions of CN + CH₂O, the ICN precursor was photolyzed at 266 nm (∼30 mJ per pulse) using the fourth harmonic of a Nd:YAG laser (Quantel Q-Smart 850), which was directed co-linearly with the supersonic flow along the axis of the Laval nozzle (“Pump laser” in Fig. 1), generating a uniform density of CN radicals. Although the CN radical is initially generated with rotational excitation,^52–54 collisions with the bath gas at the densities used ensure that rotational relaxation occurs on a timescale that is short compared with the timescale for removal of CN due to its reaction with CH₂O. The time evolution of the CN radicals produced are monitored using laser-induced fluorescence following excitation of the R₁(2) rotational line of the B²Σ − X²Σ (1,0) vibronic transition at 357.893 nm, which is generated by a Nd:YAG (Quantel, Q-smart 850) pumped dye laser (Sirah, Cobra Stretch, with Pyridine 2 dye). This “probe” laser beam (Fig. 1) crossed the gas flow perpendicularly at a fixed position ∼10–25 cm from the nozzle exit depending on the stable flow length of the particular nozzle used. Fluorescence from CN (B²Σ − X²Σ) was focused through a series of lenses and through a bandpass filter at 400 nm with a FWHM of 40 nm (Thorlabs, FB400-40) onto a Channel PhotoMultiplier (CPM) (PerkinElmer, C1952P). The gain of the CPM was controlled with a custom-built time-dependent high voltage gating module in order to block scattered light from the pump laser at 266 nm. Digitization and integration of the CPM signal were performed on an oscilloscope (LeCroy, Waverunner LT264), which is then transferred and saved to a computer for further analysis via a LabVIEW program. The timing of the experiment was also controlled via LabVIEW communication with a digital delay generator (BNC, Model 555), which randomly varied the order of the time delays between the pump and probe laser after each gas pulse.

2.2 Computational methods

Theoretical approaches were used in this work to calculate the PES of the CN + CH₂O reaction in order to further explore the reaction mechanisms responsible for the behavior of the reaction rate coefficients as a function of temperature. Geometric structures of stationary points (reactants, products, intermediates, and transition states (TSs)) were first optimized at the BHandHLYP/aug-cc-pVDZ level of theory^55–58 and further refined using M06-2X/aug-cc-pVTZ.⁵⁹ Higher-level single-point energy calculations were performed at the CCSD(T)/aug-cc-pVTZ level^60,61 to obtain more accurate energies. Vibrational frequency calculations were performed to evaluate zero-point vibrational energies (ZPVE), where TSs were found to have only one imaginary vibrational frequency. The vibrational frequency scaling factors for BHandHLYP/aug-cc-pVDZ and M06-2X/aug-cc-pVTZ are taken to be 0.9589 and 0.956, respectively.^62,63 Intrinsic reaction coordinate (IRC) calculations were carried out for all TSs located during the PES search, unless otherwise specified, to verify that they are indeed saddle points on the minimum energy pathways connecting the respective local minima. In order to further explore the long-range reaction PES as the two reactants approach each other, relaxed scans were performed along the reaction entrance channels at the BHandHLYP/aug-cc-pVDZ and M06-2X/aug-cc-pVTZ levels of theory. For a relaxed scan, one of the geometric parameters is selected as the scan coordinate while all the other geometric parameters, unless being specifically frozen (in which it is referred to as a partially constrained scan), are allowed to be optimised to give the minimum energy geometry. All electronic structure calculations were carried out using the Gaussian 09 program.⁶⁴

From the generated PES, statistical rate theory calculations were performed using the MESMER software program⁴¹ in order to obtain the rate coefficients of the system. The stationary points of the ab initio calculations provide the energies, rotational constants, and vibrational constants required by the MESMER input file. The energy wells along the PES are divided into energy grains, where each grain couples the reactant, intermediate, and product species to one another via the microcanonical rate coefficients, k(E). The individual grains can be populated/depopulated by exchange with other grains via collisional energy transfer with the buffer gas. The microcanonical rate coefficients were calculated with either Rice, Ramsperger, Kassel, and Marcus (RRKM) theory⁶⁵ for reactions involving a defined transition state or the inverse Laplace transformation (ILT) method⁶⁶ for barrierless reactions. Collisional energy transfer probabilities were described using the exponential-down model.⁶⁷ Corrections for quantum mechanical tunneling were also included using the Eckart expression.⁶⁸ The set of coupled differential equations that describe each of the energy grains is known as the energy grained master equation (EGME) and can be described by:

where p is the population density vector containing populations of each grain from each well, and M is the transition matrix that describes the population evolution due to collision energy transfer and reaction. The solution to eqn (2) is

p = Ue^ΛtU⁻¹p(0) (3)

where p(0) contains the initial conditions for each grain, U is the matrix of eigenvectors obtained from the diagonalization of M, and Λ is the diagonal matrix of corresponding eigenvalues, where the smallest are the chemically significant eigenvalues (CSE). MESMER solves the EGME and obtains the phenomenological rate coefficients from the CSE using the procedure described by Bartis and Widom.⁶⁹ Using the time dependence of the concentration of all species calculated by MESMER, the branching yield of different products can be determined. In addition, MESMER has a built-in fitting feature using available experimentally measured rate coefficients, whereby input parameters, e.g. the energies of the stationary points, can be adjusted to best fit to the experimental data. The input file for MESMER simulations used in this work are included in the online ESI.†

3 Results and discussion

3.1 Experimental results

The temperature-dependent rate coefficients for the reaction of CN + CH₂O were obtained by first measuring the time-evolution of the relative transient LIF signal of CN as the radical reacted with an excess abundance of CH₂O under pseudo-first-order conditions. The integrated LIF signal at each time delay between the photolysis laser (pump) and the dye laser (probe) was collected at least 5 times in order to obtain averaged traces of the temporal evolution of CN. Any background signal due to the probe-laser only was subtracted using the average value of data recorded prior to the pump-laser pulse.

The initial photolytic production of CN and fast collisional relaxation of rotationally-excited^52–54 CN into the CN (X²Σ, N = 2) laser-probed state led to a rapid rise in the CN LIF signal (≲1 μs) which was not resolved in our experiments. Therefore, traces were analyzed after ∼5 μs following the initial rise in the signal, after which there was an exponential decay of CN due to diffusion, reaction with CH₂O, and reaction with other possible species. Diffusion of CN out of the volume of supersonic flow through which the pump laser passed is given by:

The reaction of CN with CH₂O (reaction 1, bimolecular rate coefficient k₁), and potential reactions with other species (for example, with the precursor, other photolysis products, or reagent impurities) is given by:

where X_i represents other non-reagent species i. The observed pseudo-first-order rate coefficient for the loss of CN, k_obs, is thus given by:where N is the total number of species other than CH₂O in the cold gas flow and k_int represents the intercept of a linear fit of a plot of k_obsversus [CH₂O], with k_int expected to be dominated by k_diff. Each CN integrated LIF decay trace, I(t), which is proportional to [CN]_t, was fitted with a single exponential decay function:where t is the pump–probe time delay, and I₀ is the fitted value of LIF signal at t = 0 for exponential decay fits of the signal. Fits of eqn (7) to the data were performed starting at a pump–probe time delay of ∼5 μs. This ensured the completion of all significant relaxation of CN into the rovibronic state that was being probed that could otherwise interfere with the analysis. Examples of the background corrected, averaged, temporal evolution of the CN integrated LIF signal together with such fits are shown in Fig. 2.

Fig. 2 Averages of ≥5 temporal decays of the CN integrated LIF signal utilized to determine the pseudo-first-order rate coefficients, k_obs, for loss of CN at 32 K, both in the absence of CH₂O and for [CH₂O] = 6.06 × 10¹³ molecule cm⁻³, at a total density of (3.24 ± 0.24) × 10¹⁶ molecule cm⁻³ in Ar bath gas, together with exponential fits (eqn (7)) to the data. No data could be recorded for longer pump-probe time delays due to the limited reaction time within the low-temperature uniform supersonic flow prior to the breakup of the flow.

For each CH₂O concentration, the experiment was repeated at least five times and the k_obs values obtained were averaged to give k̄_obs. Second-order (bimolecular) plots of k̄_obsversus [CH₂O] were then generated at each temperature, examples of which are given in Fig. 3 (with the intercept k_int subtracted for clarity) with a linear least-squares fit of eqn (6) used to determine k₁ from the gradient.

Fig. 3 Variation of the intercept-subtracted, average loss rate of CN, k̄_obs − k_int, with [CH₂O], together with linear-least squares fits of eqn (6) to the data for T = 32 K, total density (Ar bath gas) = (3.24 ± 0.24) × 10¹⁶ molecule cm⁻³, and T = 84 K, total density (N₂ bath gas) = (7.56 ± 0.63) × 10¹⁶ molecule cm⁻³. Error bars represent one standard deviation of the value of k̄_obs obtained from at least 5 temporal decays of CN (Fig. 2).

For the linear fits of the second-order plots such as shown in Fig. 3, care was taken to include only the range of [CH₂O] for which k̄_obs is linear with [CH₂O] to avoid any possible influence from CH₂O dimers. Evidence for dimerization of CH₂O was found in previous low-temperature studies at higher [CH₂O] performed using a Laval nozzle, namely for rate coefficient determinations of the reactions OH + CH₂O,¹⁵ CH + CH₂O,¹⁴ and NH₂ + CH₂O.⁷⁰ The range of values of the intercept (at [CH₂O] = 0 molecule cm⁻³) was k_int = 2200–7600 s⁻¹, and the variation of k̄_obs with [CH₂O] without the subtraction of k_int can be found in Fig. S3 of the ESI.† We have been very cautious in the maximum [CH₂O] used to obtain k_obs for a given temperature, always using less than the [CH₂O] at which second-order plots have been seen to become non-linear in similar studies of radical + CH₂O reactions.^14,15 The values of k₁(T) determined in this work for the CN + CH₂O reaction, together with corresponding experimental conditions, are given in Table 1. At 70 K, k₁ was determined at two overall densities providing evidence that k₁ is independent of pressure.

Table 1 Rate coefficients and experimental conditions for kinetics studies of the CN + CH₂O reaction

T (K)	Bath gas	N total (10^16 molecule cm^−3)	k 1(T)^b (10^−11 molecule^−1 cm^3 s^−1)
a Pitot tube measurements of impact pressures were utilized to determine T and Ntotal. Here 1σ fluctuations of these values in the cold flow along the axis of the nozzle are reported. b The error of each k1(T) value represents the standard error in the fitted value of the gradient of k̄obsversus [CH2O].
32 ± 2	Ar	3.24 ± 0.24	3.57 ± 0.53
32 ± 2	Ar	3.24 ± 0.24	4.62 ± 0.84
40 ± 4	Ar	8.36 ± 1.19	2.26 ± 2.11
53 ± 4	Ar	7.04 ± 0.74	3.18 ± 1.19
56 ± 6	Ar	7.58 ± 1.11	1.56 ± 1.94
70 ± 11	Ar	11.18 ± 2.54	1.51 ± 0.49
70 ± 2	N2	2.91 ± 0.20	1.30 ± 0.52
84 ± 3	N2	7.56 ± 0.63	0.99 ± 0.27
92 ± 6	N2	4.99 ± 0.74	1.80 ± 0.29
103 ± 10	N2	6.80 ± 1.57	1.45 ± 0.19

---

The measurement uncertainty for k₁ for CN + CH₂O given in Table 1 is larger than reported in our determination of the rate coefficient for the CH + CH₂O reaction.¹⁴ The uncertainty quoted is a statistical error only from linear fits of the type shown in Fig. 3, and adding systematic errors did not increase the overall error significantly. The primary reason is that k₁(T) is about 10–100 times smaller than k_CH+CH2O(T) across the range of temperatures studied (0.99 × 10⁻¹¹–4.62 × 10⁻¹¹versus 4.86 × 10⁻¹⁰ – 11.15 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹, respectively). Hence, for the same [CH₂O] and available reaction time (controlled by the length of the stable flow for a given Laval nozzle), the LIF signal from CN does not decay as much as for CH, resulting in a relatively larger uncertainty in the fit of eqn (7) to the data to determine k_obs. The values of k₁(T) determined in this work for T = 32–103 K for the CN + CH₂O reaction are shown in Fig. 4, and also in Fig. S8(a) and S10 of the ESI,† together with experimental measurements from two other groups^36,37 reported over the temperature range 294–769 K.

Fig. 4 Measured values of k₁(T) versus temperature from this work (black squares), Chang and Wang³⁷ (blue triangles), and Yu et al.³⁶ (green triangles) (collectively YCW). The error bar of each k₁(T) from this work are those given in Table 1, and the YCW values are assigned an error equal to 0.064 × k₁(T), see text for details. The Laval fit (black line) are the MESMER simulations which are a best-fit to the experimental Laval measurements only, with ab initio parameters (Table 2) including a submerged transition state. For details of the Laval + Lit 1 and Laval + Lit 2 fitting scenarios and parameters see Table 2. The inset shows the extrapolation of Laval (black line) and Laval + Lit 2 (red line) models to lower temperatures, where significant deviation in these models is only observed below 30 K. See text and Table 2 for further details of these scenarios.

Fig. 4 shows that k₁(T) exhibits a weak positive temperature dependence above room temperature but a strong negative temperature dependence below ∼100 K, with a likely minimum somewhere between 100–200 K, suggestive of a change in reaction mechanism between the low and high temperature regimes. The values of k₁(T) at low temperatures are less precise than those previously reported at higher temperatures,^36,37 due to the challenges for this experiment as discussed above. It is noted though that for the previous data reported at higher temperatures, the absolute concentration of CH₂O was not determined by UV absorption spectroscopy (in contrast to the low temperature work reported here), rather manometric/flow methods were used to determine [CH₂O].

3.2 Ab initio calculations of the CN + CH₂O potential energy surface

The overall potential energy surface for the CN + CH₂O reaction is shown in Fig. 5, with energies obtained using CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ and shown relative to the CN + CH₂O entrance channel.

Fig. 5 PES of the CN + CH₂O reaction obtained at the CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ level of theory. All energy values are in kJ mol⁻¹ and are corrected with scaled ZPVE. The red line indicates the van der Waals complex VDW found in this work. The blue and green paths show the subsequent pathways from VDW leading to HCN + HCO (P1) and HNC + HCO (P2) products, respectively.

The geometries of the stationary points, obtained at the M06-2X/aug-cc-pVTZ level of theory, are shown in Fig. S4 (ESI†). The optimized Cartesian coordinates, vibrational frequencies, energy values, and additional scans along the reaction potential energy surface can be found in the ESI.†

A weakly bound van der Waals complex VDW (−13.3 kJ mol⁻¹) is identified following the approach of CN to CH₂O. As seen in Fig. 5, this complex leads to four different product pathways, namely H atom abstraction to form HCN (P1), H atom abstraction to form HNC (P2), and addition of CN onto the O atom (P3–4) eventually forming either NCO + ³CH₂ (P3) or HC(O)CN + H (P4). The formation of HCN + HCO (P1, −163.4 kJ mol⁻¹) is accessible through the submerged barrier TS_VDW/P1 (−0.62 kJ mol⁻¹) while the formation of HNC + HCO (P2, −103.6 kJ mol⁻¹) involves a small positive energy barrier relative to the CN + CH₂O entrance channel (TS_VDW/P2, 3.97 kJ mol⁻¹). The error of these two calculated barrier heights, even with the high level of theory used here, are such that they could both be either positive or submerged barriers if calculated at a different level of theory. At the CCSD(T)/aug-cc-pVTZ level of theory, the error of the energy is estimated to be 3.0–4.5 kJ mol⁻¹ (250–450 cm⁻¹).¹

The addition of CN onto the O atom of CH₂O involves surmounting a large barrier TS_VDW/P3-4 (32.9 kJ mol⁻¹), leading to the formation of the intermediate H₂C–O–CN (Int1, −119.2 kJ mol⁻¹). H₂C–O–CN can dissociate to form NCO + ³CH₂ (P3, 174.4 kJ mol⁻¹) or undergo cyclization through TS_1/2 (−44.2 kJ mol⁻¹) to form the cyclic intermediate Int2 (−81.2 kJ mol⁻¹). By going through TS_2/3 (−62.5 kJ mol⁻¹), the ring opens to form the intermediate H₂C(O)CN (Int3, −153.4 kJ mol⁻¹). Breaking one of the CH bonds gives the products HC(O)CN + H (P4, −66.7 kJ mol⁻¹). Although the IRC calculation did not converge successfully for TS_3/P4 (−41.0 kJ mol⁻¹), judging from the vibrational mode of the imaginary frequency, it is likely that it is the TS connecting Int3 and P4. A dashed line connecting TS_3/P4 reflects the incomplete mapping of the IRC along this coordinate.

Relaxed scans at the BHandHLYP/aug-cc-pVDZ and M06-2X/aug-cc-pVTZ level of theory attempted to map out the approach of the two reactants, CN and CH₂O. Two different reactant approaches were investigated: the CN radical approaching from the oxygen side of CH₂O, and CN approaching from the hydrogen side. For each scan, the distance between the two reactants was fixed while all other coordinates were allowed to optimize. When CN approaches from the oxygen side of CH₂O, as shown in Fig. S5 (ESI†), it is favorable for CN to orient the carbon towards CH₂O at the beginning of the approach due to the dipole–dipole attraction. The potential energy decreases smoothly as the distance between CN carbon and CH₂O oxygen decreases, eventually reaching the potential energy well where the van der Waals structure VDW is located.

In the case where CN approaches from the hydrogen side of CH₂O, as shown in Fig. S6 (ESI†), it is favorable at first for CN to orient itself such that the nitrogen side points towards CH₂O. The potential energy decreases with decreasing separation between the two moieties until encountering a fairly “flat” region of the potential energy surface. Here, with a distance of approximately 3.3 Å between nitrogen of CN and carbon of CH₂O, or a center of mass distance between the moieties of approximately 4.5 Å apart, the overall energy of the system (BHandHLYP/aug-cc-pVTZ, uncorrected for ZPVE) reaches −5 kJ mol⁻¹ relative to the entrance channel. In order to further explore this “flat” region of the PES, a relaxed scan from this point has been done using the O=C⋯N angle as the scanning parameter, as shown in Fig. 6.

Fig. 6 Potential energy curve (ZPVE corrected) calculated at the M06-2X/aug-cc-pVTZ level for the rotation of CN around CH₂O for O=C–N angle from 55° to 180°. Energy is relative to the sum of the energies (ZPVE corrected) of the two separate reacting species, CN and CH₂O, at the M06-2X/aug-cc-pVTZ level.

While the CN radical rotates around CH₂O, the potential energy first experiences a fairly flat region of the potential and then falls smoothly into the potential energy well corresponding to the van der Waals structure VDW. Thus, it is suggested that both ways of approach can eventually lead to the van der Waals structure VDW.

A direct H abstraction mechanism of CN from CH₂O to form HCN + HCO was studied in a previous calculation using the QCISD/6-31G**//UHF/6-31G** level of theory.³⁸ This prior work suggested that the reaction mechanism involved overcoming a small positive ∼2 kJ mol⁻¹ (emerged) barrier relative to the CN + CH₂O entrance channel. The current study and a very recent study¹³ have not identified this pathway on the reaction coordinate, perhaps because of the lower level of theory used for geometry optimization (UHF/6-31G**) in the previous study.¹³ In the current study, instead of direct abstraction, an indirect channel to form HCN + HCO was found which involves the van der Waals complex VDW and a small submerged barrier TS_VDW/P1. A very recent theoretical study¹³ by Tonolo et al. also identified the VDW structure on the reaction PES (−13.0 kJ mol⁻¹ relative to reactants, including ZPE, CCSD(T)/CBS + CV) leading to hydrogen abstraction to form HCN + HCO through a submerged barrier (−1.26 kJ mol⁻¹). This pathway is consistent with our current work, with less than 0.5 kJ mol⁻¹ difference in energy in the complex and barrier. Our work also details for the first time the presence of the HNC pathway from the VDW complex.

The work from Tonolo et al. identifies several stationary points along the reaction coordinate not used in the current study. For example, Tonolo et al. claim a barrierless route to forming a C–C bond directly from the reactants, resulting in a tetrahedral intermediate which they label 1C in a deep potential well (−153 kJ mol⁻¹). Although we agree that this structure can be formed (Int3 in Fig. 5), we were unable to connect this structure to the reactants (relaxed scans, BHandHLYP/aug-cc-pVDZ and M06-2X/aug-cc-pVTZ) without surmounting a large barrier at least 50 kJ mol⁻¹ (CCSD(T)/aug-cc-pVTZ//M06-2X/aug-cc-pVTZ) higher in energy than the reactants, as shown in Fig. S7 (ESI†). It is likely that the overestimated rate coefficients from Tonolo et al. compared to previous and current experimental work are directly related to not identifying a substantial barrier to C–C bond formation in forming their intermediate 1C. This structure 1C and the remainder of the addition pathways Tonolo et al. identified coming from this structure, including roaming mechanism pathways, were not considered in our search for reactions which would occur at temperatures relevant to the current study. This impacts on the relative importance stated for the pathway leading to the formation of formyl cyanide (HCOCN), although we identified a new pathway to HC(O)CN + H through VDW which does not involve roaming and has a barrier (TS_VDW/P3-4).

In comparison with the PES of OH + CH₂O, which has been studied extensively in previous work,^2,16,71–75 the results from the current work on CN + CH₂O show that both systems share a similar shape of the PES in terms of the energy profile or mechanism for the H abstraction reaction. Following the approach of the two reacting species, a pre-reaction complex is formed followed by a transition state with a small barrier to form products. The difference in energy which determines whether the barrier is positive or submerged relative to the reactants energies is in the kJ mol⁻¹ range, and hence within the uncertainty of most calculation methods, and so whether the transition state is submerged or not will depend on the level of theory used.

For example, the energy values of the transition state for H abstraction channel of OH + CH₂O computed with more robust methods tend to give lower values,² decreasing the barrier from being slightly emerged to slightly submerged. The latest value reported from Machado et al.² obtained at CCSD(T)/CBS level is approximately −5.7 kJ mol⁻¹.

3.3 Calculation of rate coefficients using the MESMER package

The energy values, see Fig. 5 (and given in Tables S1–S3, ESI†), Cartesian coordinates of the stationary points (Tables S4 and S5, ESI†), vibrational frequencies (Tables S6 and S7, ESI†), and rotational constants (Tables S8 and S9, ESI†) of the species obtained from the ab initio calculations were used as the inputs for the master equation solver, MESMER,⁴¹ in order to calculate the rate coefficients for CN + CH₂O. From simulations over a wide range of temperatures (4–1000 K) and pressures (10¹⁵–10¹⁹ molecule cm⁻³) it was observed that HCN and HNC accounted for greater than 99.99% of the products, under all conditions. The reaction occurs initially via van der Waals complex (VDW) formation followed by transition states TS_VDW/P1 and TS_VDW/P2 to form HCN + HCO and HNC + HCO, respectively. The HNC channel never accounts for more than 1% yield, as shown in Fig. S8(b) (ESI†). However, uncertainties on calculated energy barriers (estimated to be 3.0–4.5 kJ mol⁻¹ at the CCSD(T)/aug-cc-pVTZ level of theory¹) can result in a reverse situation for which TS_VDW/P2 (Fig. 5) can be slightly submerged and TS/VDW/P1 somewhat emerged, and so there is a possibility that the yield of HNC may be significantly larger than this, and the yield of HCN correspondingly smaller.

As CN + CH₂O did not show a pressure dependence for k₁(T), either experimentally near 10¹⁷ molecule cm⁻³ or from MESMER simulations between 10¹⁵–10¹⁸ molecule cm⁻³, the gas density was set at 10¹³ molecule cm⁻³ (so making sure in a pressure independent region) for the MESMER simulations. Above 10¹⁹ molecule cm⁻³ and T < 50 K, a pressure dependence was evident and the VDW species was populated. No pressure dependence was observed in the current calculations, implying that the van der Waals complex is not significantly stabilised, even at the lowest temperatures.

Therefore, to a very good approximation, Fig. 5 can be reduced to just the formation of HCN and HNC as shown in Fig. 7 without losing chemical information; reducing the system to just HCN formation would still describe the system to better than 99% of the product yield.

Fig. 7 A simplified PES used for MESMER fitting scenarios, see Table 2, together with molecular structures at key stationary points. Initial van der Waals complex, VDW red, formation is followed by transition states to products, HCN + HCO (P1, blue) and HNC + HCO (P2, green) respectively. See text for details.

The inverse Laplace transformation (ILT) method⁶⁶ was used to calculate the microcannonical rate coefficients for the barrierless formation of the VDW complex from CN + CH₂O. The ILT approach overcomes the problem of explicitly assigning a transition state for a barrierless process, whose position would be varying as a function of temperature. The ILT method is especially convenient when experimental data are available, as for this study. The ILT parameters for VDW formation were assigned by:

where A^∞_ILT,VDW and n^∞_ILT,VDW are the ILT parameters which describe the rate coefficient for the formation of the VDW at the high pressure limit (in cm³ molecule⁻¹ s⁻¹ and unitless respectively). A reference temperature in the expression of 30 K was chosen so that the A factor is representative of the lowest experimental data point at 32 K, making the ILT more suited with the low temperature regime relevant to conditions in the ISM. The MESMER input file for these calculations is given at the end of the ESI.†

As well as performing simulations, MESMER can adjust the important rate controlling parameters in eqn (9), as well as the transition state energies in order to best fit to available experimental data. During such data fitting, MESMER uses the Marquardt algorithm⁷⁶ to adjust the parameters in order to minimize χ²:

where k_calc(T) is the MESMER calculated rate coefficient, and k_expt(T) and σ are the experimental rate coefficients and their associated error. For the MESMER fitting exercise, the present low temperature results used the experimental errors given in Table 1, whereas the room temperature and above rate coefficients from YCW^36,37 (denoted Laval + Lit 1 in Table 2) were assigned an error of 0.064 × k₁(T). The MESMER fitting scenarios used in this work and the optimised kinetic parameters are summarised in Table 2.

Table 2 Fitting scenarios and optimized kinetics parameters from MESMER fitting

Name of fitting scenario	Laval	Laval + Lit 1	Laval + Lit 2
A ^∞ILT,VDW and n^∞ILT,VDW are the ILT parameters defined in eqn (9) and the other parameters are the energy of the transition state (TS_VDV/P1, see Fig. 5 and 7) and its imaginary frequency.a Units are cm^3 molecule^−1 s^−1.b Units are cm^−1.c N is the number of degrees of freedom, the number of data points minus the number of fitting parameters. The reported errors are 1 sigma.
Temperature range (K)	32–103	32–769	32–769
A ^∞ILT,VDW ×10^−11^a	5.29 ± 0.8	2.5 ± 0.8	6.58 ± 0.7
n ^∞ILT,VDW (unitless)	−0.77 ± 0.17	0.016 ± 0.005	−0.10 ± 0.02
TS_VDW/P1 (kJ mol^−1)	−0.62, fixed	−0.62, fixed	4.0 ± 0.9
Imaginary frequency^b (cm^−1)	215, fixed	215, fixed	806 ± 24
χ ^2/N^c	1.09	4.24	0.80

---

Initially, MESMER data fitting was carried out with energies fixed at the ab initio calculated values, as shown in Fig. 7, with only the ILT parameters in eqn (9) being adjusted. When fitting to just the low temperature Laval nozzle data (Laval scenario in Table 2) an excellent visual fit with the data was obtained (see Fig. 4). χ²/N is a measure of the goodness of fit, where χ² is given by eqn (10) and N is the number of degrees of freedom, which is equal to the number of data points minus the number of fitting parameters. A χ²/N value close to 1.0 represents a good fit, and from Table 2 it can be seen that the Laval model is close to 1.0, and indicates the experimental errors in Table 1 are realistic.

When MESMER simultaneously fits to the low and high-temperature data of YCW^36,37 (denoted Laval + Lit 1 scenario in Table 2), the fit to the data is much worse, with χ²/N = 4.24, and as can clearly be seen in Fig. 4. Therefore it is concluded that MESMER calculations based on our ab initio calculated values is not capable of reproducing all the experimental data shown in Fig. 4. The reason is because the ab initio value for the energy of the transition state (TS_VDW/P1) is submerged, i.e. is negative with respect to the reagent energies, and a reaction mechanism proceeding via a submerged transition state is going to predict k₁ to decrease with increasing temperature, which is at odds with the high temperature literature data, which requires the transition state to have a positive energy with respect to the reagents. The simplest way to account for the variation of k₁ across the full range of temperatures is to increase the transition state energy so that it is positive, making it consistent with the high-temperature literature data trends and to increase the imaginary frequency of the tunnelling coordinate in the transition state so that quantum mechanical tunnelling overcomes the effect of the positive barrier, leading to an increase in k₁ at lower temperatures, as is observed experimentally. When the model scenario Laval + Lit 2 (Table 2), where both the energy and the imaginary frequency of both transition states for the channels producing HCN and HNC products (TS_VDW/P1 and TS_VDW/P2, respectively in Fig. 7) are adjusted in the MESMER simulations using a best-fit procedure, the k₁ experimental data across the full range of temperatures can be fitted well, evidenced by a value of χ²/N = 0.80. The optimum MESMER adjustments from the fitting were an increase in the energy from −0.62 to + 4.0 kJ mol⁻¹ and an increase in the imaginary frequency from 215 to 806 cm⁻¹ for TS_VDW/P1. We stress that the 4.59 kJ mol⁻¹ab initio energy difference between the two transition states leading to HCN and HNC products (TS_VDW/P1 and TS_VDW/P2 in Fig. 7, respectively) is maintained, with the energy of the transition state for HNC (TS_VDW/P2) now being +8.59 kJ mol⁻¹, and the imaginary frequency for the transition-state leading to HNC was also increased to ∼806 cm⁻¹. This is reasonable given the similar calculated geometries of the two transition states, but with just the C and N flipped, as shown in Fig. S4 (ESI†), and the tunnelling motion is therefore roughly the same.

The MESMER simulation model scenario Laval + Lit 2 of the fit to the data is shown in Fig. 4, and is much improved over Laval + Lit 1. In order for the Laval + Lit 2 scenario to give the best fit to the k₁ data, the transition state energy needed to be increased by ∼4.6 kJ mol⁻¹ (from −0.62 kJ mol⁻¹ to +4.0 kJ mol⁻¹ for TS_VDW/P1, which is still the TS via which the reaction occurs to form HCN products) from our ab initio calculated values at the CCSD(T)/aug-cc-pVTZ level of theory.¹ However, this adjustment seems reasonable given it is close to the estimated uncertainty at the CCSD(T)/aug-cc-pVTZ level of theory of ∼3.0–4.5 kJ mol⁻¹.¹ Similarly, even though the imaginary frequencies has been increased significantly, from 205 cm⁻¹ to 806 cm⁻¹, the adjusted value of ∼800 cm⁻¹ again is not unreasonable.⁵ As both transition states emanate from the same pre-reaction complex leading to H abstraction products, it seems reasonable to adjust the energies of each of these by the same amount whilst maintaining the same difference between them. The motivation for changing the initial ab initio results to the fitted results using MESMER is to show that the experimental data can be fitted well if the potential energy surface, specifically for these two transition states, is modified. Just changing the properties of the two initially formed transition states which are formed from the same pre-reaction complex is the most straightforward way to do this.

From Fig. 4 it can be seen that the minimum in the rate coefficient k₁ occurs around 150 K, which is the point at which the controlling influences of quantum mechanical tunnelling and the 4.0 kJ mol⁻¹ barrier for TS_VDW/P1 become balanced. Also included in Fig. 4 is the simulation of model Laval down to 7 K, where at 10 K the predicted value of k₁ is approximately twice that of the Laval + Lit 2 model. This significant difference emphasizes how sensitive k₁ is to the parameters used in the models when extrapolating down to very low temperatures, i.e. <10 K, as both the Laval and Laval + Lit 2 models give almost the same quality fits to the experimentally measured k₁ using the Laval nozzle.

3.4 Parameterisation of the temperature dependent rate coefficient k₁(T) for use in astrochemical modelling

In order to provide an analytical expression for the results of the master equation analysis which is suitable for use in astrochemical modelling, the best-fit parameters of the Laval and Laval + Lit 2 scenarios (see Table 2) were used as input for MESMER simulations to generate values of k₁(T) over a very wide temperature range from 4–1000 K, as tabulated in Table S10 (ESI†). Fig S8(a) (ESI†) shows all of the experimental data for k₁(T) together with the MESMER simulations using the Laval + Lit 2 model. As CN + CH₂O did not show a pressure dependence for k₁(T), neither experimentally near 10¹⁷ molecule cm⁻³ nor from MESMER simulations between 10¹⁵–10¹⁸ molecule cm⁻³, the gas density was set at 10¹³ molecule cm⁻³ for the MESMER simulations; above 10¹⁹ molecule cm⁻³ and T < 50 K a pressure dependence was evident and the VDW species begins to become stabilized. From 50–1000 K, a 30 cm⁻¹ grain size was sufficient to calculate a converged rate coefficient, but it was reduced down to 2 cm⁻¹ for the 4–50 K simulations. It was not possible to run MESMER simulations at lower temperatures. Over the range of 4–1000 K, the MESMER simulations indicate that HCN + HCO are the only significant products formed. The fractional yield of HNC is less than 0.33% for all temperatures between 4–1000 K, because the energy of TS_VDW/P2 (forming HNC) is 4.59 kJ mol⁻¹ above that of TS_VDW/P1 (forming HCN), and the imaginary frequencies of the both transition states is 806 cm⁻¹. The temperature dependence of the fractional yield of the HCN and HNC products is shown in Fig. S8(b) (ESI†). The results of k₁ for this simulation (Table S10 in the ESI†) are the basis for our recommended values for the overall CN + CH₂O rate coefficient, from which parameterisations are developed for use in astrochemical models.

The equation utilized commonly as a parameterisation in several astrochemical models (which use, for example, the UMIST Rate12 (UDfa, https://udfa.ajmarkwick.net) or KIDA (https://kida.obs.u-bordeaux1.fr⁷⁸) databases) to represent k(T) is a single modified Arrhenius (MA) equation:

where α, β and γ are best-fit parameters. Although Table S10 (ESI†) can be used to look up the MESMER simulated values of k₁(T) at a given temperature in the range 4–1000 K, a parameterisation in the form of eqn (11) is convenient to calculate k₁(T) using an expression traditionally used by astrochemical modellers. However, a single MA eqn (11) is not able to adequately fit k₁(T) for CN + CH₂O over the full range of temperatures 4–1000 K. It is especially inadequate for the k₁ values around 150 K, where the T dependence changes from negative to positive, see Fig. 4, and also at the lowest temperatures where k₁ is rapidly increasing. Running MESMER simulations down to the lowest temperatures means that no extrapolation of fits using eqn (11) is necessary and avoids the extra error associated with this. However, given the use of the single MA in astrochemical databases, a piecewise parameterisation using eqn (11) is presented here. Our temperature ranges were chosen for the piecewise least-squares fits of MA eqn (11) to the MESMER simulated values of k₁(T), and Fig. S10, ESI† shows these fits to the Laval + Lit 2 scenario for low (4–20 K), (20–100 K), (100–300 K) and (300–1000 K) temperature ranges. In addition, a fit to the Laval scenario for the low (3–80 K) temperature range was carried out and is shown in Fig. S10 (ESI†).

In the parameterization given by eqn (11) each MESMER data point was assigned a 5% error. The values for the best-fit parameters from eqn (11) fits to the Laval + Lit 2 simulation data are given in Table 3. From Table 3, it is noted that χ²/N is always less than 1, which implies that on average this representation is better or within 5% of the MESMER simulation data, see Fig. S10 (ESI†). Also included in Table 3 are the parameterization values for the Laval model for the range T = 3–80 K. But it is the parameterization of the Laval + Lit 2 model that represents our recommendation of k₁(T) for use in astrochemical modelling.

Table 3 Best-fit parameters from a single modified Arrhenius equation (eqn (11)) fitted to MESMER simulated data over various temperature ranges

T range (K)	α (cm^3 molecule^−1 s^−1)	β (Unitless)	γ (K)	χ ^2/N^b	MESMER fitting scenario
a We recommend using the Laval + Lit 2 fitting scenario for astrochemical modelling. b N is the number of degrees of freedom, which is equal to the number of MESMER simulated data points minus the number of fitting parameters.
3–80^a	(3.05 ± 0.14) × 10^−12	−1.11 ± 0.02	−2.18 ± 0.23	0.70	Laval
4–20	(1.18 ± 0.08) × 10^−11	−0.58 ± 0.03	0.53 ± 0.24	0.55	Laval + Lit 2
20–100	(3.72 ± 0.15) × 10^−12	−1.09 ± 0.04	5.2 ± 1.8	0.12	Laval + Lit 2
100–300	(5.99 ± 0.14) × 10^−11	−2.19 ± 0.04	−313.4 ± 7.7	0.11	Laval + Lit 2
300–1000	(6.26 ± 0.30) × 10^−11	−0.02 ± 0.03	398 ± 16	0.03	Laval + Lit 2

---

As a check of whether the parameterisation for k₁(T) gives values that are realistic at the very lowest temperatures, classical capture theory (CCT) calculations have been carried out for the CN + CH₂O reaction to calculate the rate coefficient at the collision limit k_coll(T), using the molecular parameters given in Table S11 (ESI†). As shown in Fig. S9 (ESI†) for the CN and CH₂O collision pair, k_coll(T) is dominated by large dipole–dipole interactions below 600 K. As shown in Fig. 8, k_coll at ∼3 K is equal to 1.3 × 10⁻⁹ cm³ molecule⁻¹ s⁻¹, and from the MESMER simulations k₁(T) is still calculated to be a factor of 10 or so less than k_coll(T), and does not get close nor indeed exceed k_coll(T) for all temperatures relevant to conditions used in astrochemical models.

Fig. 8 Plot of the classical capture theory (CCT) rate coefficient (blue curve) for the CN + CH₂O reaction versus the k₁(T) from the Laval + Lit 2 model (black filled circles), see Table S10 (ESI†), where errors have been propagated via the covariance matrix obtained from fitting to the experimental data, and the MA extrapolation to the lowest temperature (red curve). Note that the predicted values of k₁(T) increase markedly at low T, but never approach the values from CCT.

4. Astrochemical modelling

Three different astrochemical environments were considered in order to explore the impact of the newly measured and evaluated rate coefficients (k₁(T)) for the reaction of CN with CH₂O over a wide range of temperatures. These were a dark cloud model for a total density of 2 × 10⁴ cm⁻³, at T = 5, 10 and 30 K, hot cores/corino models at higher densities and temperatures, and C-rich and O-rich AGB outflow models, further details of which are given in the following sections. The network for the model utilised the UMIST Rate12 kinetic database⁷⁷ with updates to rate coefficients to include recent measurements of SiH with O₂⁷⁹ and CH with CH₂O.¹⁴ Further details are given in Section 3.4, together with the parameterisation used for k₁(T).

4.1 Dark cloud models

We calculated time-dependent chemical kinetic models using low-metal abundances appropriate for cold, dark interstellar clouds with a number density of n(H₂) = 10⁴ cm⁻³ and temperatures of 5, 10 and 30 K, comparing results using the new Laval + Lit2 MA parameterisation to those obtained with the rate coefficient of the CN + CH₂O reaction in the UMIST Rate12 (UDFA, https://udfa.ajmarkwick.net⁷⁷) database. The latter adopts an energy barrier of 826 K (6.9 kJ mol⁻¹) over the range 297–2500 K and, in the lack of low-temperature information, the rate coefficient extrapolates to zero at low temperatures. We find that the inclusion of reaction (1), even at its faster low-temperature fit, makes no difference to the abundances of CN, CH₂O, HCO, HCN and HNC (see Fig. S11, ESI†). This arises because each of these common interstellar molecules are formed through a large variety of reactions, many of which involve highly abundant atoms and radicals such as O, N and CH₂, as well as ion-neutral routes.

The inclusion of the title reaction does, however, impact on the abundance of HCOCN which our ab initio calculations rule out as a product of the CN + CH₂O reaction. Previous modeling¹³ that included this product channel with a rate coefficient of up to 2 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹ had shown that the calculated HCOCN fractional abundance was similar to that observed, 3.5 × 10⁻¹¹, in the cold dark cloud TMC-1, (n(H₂) = 10⁴ cm⁻³, T = 10 K).⁸⁰ Our results indicate that, since CH₂O is an abundant reactant, the total production rate of HCOCN is likely to be an order of magnitude less than the observational requirement.⁸⁰

4.2 Hot core/corino models

We also investigated higher density, higher temperature models appropriate for hot core/corino sources. In particular, we modeled two of the hot sources in the Orion Molecular Cloud, the Orion Hot Core (n(H₂) = 5 × 10⁷ cm⁻³, T = 225 K) and the Plateau (n(H₂) = 10⁶ cm⁻³, T = 125 K). In these sources, initial abundances are molecular and reflect those in dust grain ice mantles which evaporate in the hot gas surrounding a massive young star. Our initial fractional abundances are taken from Doddipatla et al.⁸¹ Although these values are specific to these objects, our results are likely appropriate to other hot core/corino sources. Fig. S12 (ESI†) shows the abundances of HCO and HCN calculated for the Orion Hot Core, those for the Plateau source show a similar behaviour. Again, we find no discernible difference between our newly calculated abundances and those using the UDfA (UMIST Rate12⁷⁷) rate coefficient. The major reason for this is that neither CN nor CH₂O are abundant species in these evaporating interstellar ice mantles nor are they produced in significant fractions in the hot gas chemistry. We note in passing that the difficulty mentioned above, in reproducing the observed HCOCN abundance in hot sources, is likely to persist until kinetic data for any new significant HCOCN production reactions are measured or calculated.

4.3 Asymptotic giant star (AGB) outflows

The AGB outflow model is based on the publicly available UMIST CSE model (McElroy et al.⁷⁷). Following Van de Sande et al.,⁸² the temperature of the outflows is parameterised as where T* and R* are the stellar temperature and radius, respectively, and ε is the power law exponent. We assumed a value of R* = 5 × 10¹³ cm, T* = 2300 K and ε = 0.7. The outflow is assumed to be spherically symmetric with an expansion velocity of 15 km s⁻¹, while varying the mass-loss rate between 10⁻⁵, 10⁻⁶, and 10⁻⁷ M_sun per year. Similar to West et al.,¹⁴ we considered both an O-rich and C-rich outflow, with the abundance of parent species taken from Agúndez et al.⁸³

We find that changing the parameterisation of the reaction rate to the new measured and calculated rates does not affect the chemistry of both C-rich and O-rich outflows. The abundances and column densities of all species involved (CH₂O, CH, HCO, HCN and HNC) are not changed. Fig. S13 and S14 (ESI†) show the abundances of HCO and HCN calculated for the high mass-loss rate O-rich and C-rich outflow, respectively. This behaviour is also seen in the lower mass-loss rate outflows.

5. Conclusions

Low-temperature rate coefficients, k₁(T), for the reaction CN + CH₂O have been measured for the first time below room temperature in a pulsed Laval apparatus using the PLP–LIF technique. A negative temperature dependence of k₁(T) was observed below ∼100 K, in contrast to the positive temperature dependence of k₁(T) reported previously in experiments performed above 300 K. Both ab initio calculations of the potential energy surface (PES) for the reaction and MESMER rate-theory simulations making use of this PES were performed in order to explore the overall mechanism that is consistent with the observed k₁(T). Two low energy pathways were found on the PES through a VDW complex (binding energy of 13.3 kJ mol⁻¹) with either a small positive barrier (3.97 kJ mol⁻¹) to form HNC + HCO, or a submerged barrier (−0.62 kJ mol⁻¹) to form HCN + HCO. This calculated PES can explain the negative temperature dependence of the rate coefficients at low temperatures but not the previously measured literature rate coefficients at room temperature and above. For the formation of formyl cyanide, HCOCN, a large activation barrier of 32.9 kJ mol⁻¹ was calculated. Increasing both the energy of the transition state for the formation of HCN + HCO products from −0.62 to +4.0 kJ mol⁻¹ (which is within the accuracy of the theory used) and also its imaginary frequency to 806 cm⁻¹ in the MESMER simulations, enabled the temperature dependence of k₁(T) across the entire experimental data to be reproduced very well. The energy and imaginary frequency of the other transition state for formation of HNC + HCO products were changed by the same amounts so their relative values were the same, ensuring that only the HCN + HCO remains the important one for this reaction. The postulated mechanism is similar to previous reactions of OH with volatile organic compounds where the formation of a weakly bound van der Waals complex, whose lifetime is extended at low temperature, is then followed by quantum mechanical tunneling through a small activation barrier to products.³

Using the MESMER package the product branching ratio was calculated and the formation of HNC was not found to be important as a reaction channel. MESMER was then used to simulate the rate coefficients k₁(T) from 4–1000 K, and this dataset was used to recommend a parameterization from best-fit modified Arrhenius expressions for k₁(T) for use in astrochemical modelling. Even at 4 K the simulated values of k₁(T) were a factor of 10 less than the collision limit. This parameterization of k₁(T) was used as input to astrochemical models of dark clouds, hot core/corinos, and Asymptotic Giant Star (AGB) stellar outflow environments using the UMIST Rate12 (UDfa) network of reactions. For a range of temperatures, the models yielded no significant changes in the abundances of HCN, HNC, and HCO upon inclusion of their formation using the rate coefficients reported here. It was found that the new rate coefficients for the reaction CN + CH₂O do not have a significant effect on the modeled abundances of reagents and products of this reaction due to several other competing reactions with large rate coefficients. However, we note that the uncertainty at the CCSD(T)/aug-cc-pVTZ level of theory, (∼3.0–4.5 kJ mol⁻¹) is similar to the difference in energy (4.59 kJ mol⁻¹) calculated for the transition states forming HCN+HCO and HNC+HCO products. Hence, although HCN+HCO are the most likely products, it is possible that the HNC channel may also be the dominant one, as within the calculation errors the HNC+HCO transition state might be the lower one. In order to examine the impact of the scenario in which HNC is the dominant product of the reaction, the rate coefficient k₁ was maintained at the same value, but it was assumed that HNC+HCO formed 100% of the products. For a 10 K rate coefficient and 100% formation to HNC, the title reaction provides less than 1% of the overall HNC production rate for cold environments such as molecular clouds. For the Orion Hot Core case at higher temperatures, the HNC+HCO reaction provides around 1% of the total HNC production rate. So, assuming 100% HNC production rather than very little production from the CN+CH₂O reaction made a negligible difference to the modelled HNC abundance.

The main astrochemical implication of this paper is that the title reaction cannot contribute to the formation of formyl cyanide, HCOCN, in interstellar clouds, as has been suggested previously and currently implemented in astrochemical models.

Author contributions

DEH and LD conceived the experiment, NAW and ER performed the experiments, LHDL and MAB performed the theoretical calculations, TJM and MVdS performed the astrochemical calculations, NAW, DEH and JHL wrote the manuscript with contributions from all authors.

Conflicts of interest

There is no conflict of interest to declare.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreements no. 646758 and 948525). TJM acknowledges the receipt of a Leverhulme Emeritus Fellowship and is grateful to the STFC for support through grant ST/T000198/1. MVdS acknowledges support from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 882991. We would like to thank the mechanical and electronics workshop in the School of Chemistry at the University of Leeds for support. The ab initio and MESMER work was undertaken on ARC3, part of the High-Performance Computing facilities at the University of Leeds, UK. We would like to thank Professor John Plane for helpful discussions on classical capture theory. The authors would like to thank Dr Tom Speak for valuable input on ab initio and MESMER calculations. The authors would like to thank the two reviewers, whose extensive reports and very constructive suggestions have led to a considerable improvement in the paper.

References

1. G. Czakó, I. Szabó and H. Telekes, J. Phys. Chem. A, 2014, 118, 646–654.
2. G. d S. Machado, E. M. Martins, L. Baptista and G. F. Bauerfeldt, J. Phys. Chem. A, 2020, 124, 2309–2317.
3. D. E. Heard, Acc. Chem. Res., 2018, 51, 2620–2627.
4. A. Potapov, A. Canosa, E. Jiménez and B. Rowe, Angew. Chem., Int. Ed., 2017, 56, 8618–8640.
5. R. J. Shannon, M. A. Blitz, A. Goddard and D. E. Heard, Nat. Chem., 2013, 5, 745–749.
6. K. Acharyya, E. Herbst, R. L. Caravan, R. J. Shannon, M. A. Blitz and D. E. Heard, Mol. Phys., 2015, 113, 2243–2254.
7. M. Antiñolo, M. Agúndez, E. Jiménez, B. Ballesteros, A. Canosa, G. El Dib, J. Albaladejo and J. Cernicharo, Astrophys. J., 2016, 823, 25.
8. A. J. Ocaña, S. Blázquez, A. Potapov, B. Ballesteros, A. Canosa, M. Antiñolo, L. Vereecken, J. Albaladejo and E. Jiménez, Phys. Chem. Chem. Phys., 2019, 21, 6942.
9. K. M. Hickson, J.-C. Loison, D. Nuñez-Reyes and R. Méreau, J. Phys. Chem. Lett., 2016, 7, 3641–3646.
10. D. Chastaing, P. L. James, I. R. Sims and I. W. M. Smith, Faraday Discuss., 1998, 109, 165–181.
11. Q. Y. Meng, K. M. Hickson, K. J. Shao, J. C. Loison and D. H. Zhang, Phys. Chem. Chem. Phys., 2016, 18, 29286–29292.
12. A. J. Remijan, M. Hollis, F. J. Lovas, W. D. Stork, P. R. Jewell and D. S. Meier, Astrophys. J., 2008, 675, L85–L88.
13. F. Tonolo, J. Lupi, C. Puzzarini and V. Barone, Astrophys. J., 2020, 900, 85.
14. N. A. West, T. J. Millar, M. Van de Sande, E. Rutter, M. A. Blitz, L. Decin and D. E. Heard, Astrophys. J., 2019, 885, 134.
15. A. J. Ocaña, E. Jiménez, B. Ballesteros, A. Canosa, M. Antiñolo, J. Albaladejo, M. Agúndez, J. Cernicharo, A. Zanchet, P. del Mazo, O. Roncero and A. Aguado, Astrophys. J., 2017, 850(28), 12.
16. A. Zanchet, P. del Mazo, A. Aguado, O. Roncero, E. Jiménez, A. Canosa, M. Agúndez and J. Cernicharo, Phys. Chem. Chem. Phys., 2018, 20, 5415–5426.
17. H. Olofsson, Astrophys. Space Sci., 1997, 251, 31–39.
18. K. Willacy, Astrophys. J., Lett., 2003, 600, L87.
19. A. McKellar, Publ. Astron. Soc. Pac., 1940, 52, 187–192.
20. L. E. Snyder, D. Buhl, B. Zuckerman and P. Palmer, Phys. Rev. Lett., 1969, 22, 679–681.
21. B. Zuckerman, D. Buhl, P. Palmer and L. E. Snyder, Astrophys. J., 1970, 160, 485.
22. A. Mann and D. Williams, Nature, 1980, 283, 721–725.
23. R. I. Kaiser, Chem. Rev., 2002, 102, 1309–1358.
24. I. W. M. Smith, A. M. Sage, N. M. Donahue, E. Herbst and D. Quan, Faraday Discuss., 2006, 133, 137–156.
25. L. E. Snyder, J. Hollis and B. Ulich, Astrophys. J., 1976, 208, L91–L94.
26. L. E. Snyder and D. Buhl, Astrophys. J., 1971, 163, L47.
27. E. Churchwell and J. H. Bieging, Astrophys. J., 1983, 265, 216–222.
28. P. Palmer, B. Zuckerman, D. Buhl and L. E. Snyder, Astrophys. J., 1969, 156, L147.
29. M. Schenewerk, L. Snyder, J. Hollis, P. Jewell and L. M. Ziurys, Astrophys. J., 1988, 328, 785–796.
30. J. Harju, Astron. Astrophys., 1989, 219, 293–302.
31. A. M. Dean and J. W. Bozzelli, Gas-phase combustion chemistry, Springer, 2000, pp. 125–341.
32. G. Visconti, Fundamentals of Physics and Chemistry of the Atmosphere, Springer, 2001.
33. C. A. Nixon, R. K. Achterberg, N. A. Teanby, P. G. Irwin, J.-M. Flaud, I. Kleiner, A. Dehayem-Kamadjeu, L. R. Brown, R. L. Sams and B. Bézard, Faraday Discuss., 2010, 147, 65–81.
34. V. A. Krasnopolsky, Icarus, 2009, 201, 226–256.
35. I. R. Cooke and I. R. Sims, ACS Earth Space Chem., 2019, 3, 1109–1134.
36. T. Yu, D. L. Yang and M. Lin, Int. J. Chem. Kinet., 1993, 25, 1053–1064.
37. Y.-W. Chang and N. S. Wang, Chem. Phys., 1995, 200, 431–437.
38. W. Lin Feng, Y. Wang, S. Wen Zhang and X. Yong Pang, Chem. Phys. Lett., 1997, 266, 43–46.
39. B. Ruscic and D. H. Bross, Active Thermodynamical Tables (ATcT) values based on version 1.122g of the Thermochemical Network, https://atct.anl.gov/Thermochemical%20Data/version%201.122/, 2019.
40. M. Born, S. Ingemann and N. M. M. Nibbering, J. Chem. Phys., 1996, 100, 17662–17669.
41. D. R. Glowacki, C.-H. Liang, C. Morley, M. J. Pilling and S. H. Robertson, J. Phys. Chem. A, 2012, 116, 9545–9560.
42. S. E. Taylor, A. Goddard, M. A. Blitz, P. A. Cleary and D. E. Heard, Phys. Chem. Chem. Phys., 2008, 10, 422–437.
43. R. J. Shannon, S. Taylor, A. Goddard, M. A. Blitz and D. E. Heard, Phys. Chem. Chem. Phys., 2010, 12, 13511–13514.
44. R. L. Caravan, R. J. Shannon, T. Lewis, M. A. Blitz and D. E. Heard, J. Phys. Chem. A, 2015, 119, 7130–7137.
45. J. C. Gómez Martín, R. L. Caravan, M. A. Blitz, D. E. Heard and J. M. C. Plane, J. Phys. Chem. A, 2014, 118, 2693–2701.
46. R. J. Shannon, R. L. Caravan, M. A. Blitz and D. E. Heard, Phys. Chem. Chem. Phys., 2014, 16, 3466–3478.
47. B. R. Rowe, A. Canosa and D. E. Heard, Uniform Supersonic Flows in Chemical Physics, World Scientific, 2022.
48. V. Sivakumaran, D. Holscher, T. J. Dillon and J. N. Crowley, Phys. Chem. Chem. Phys., 2003, 5, 4821–4827.
49. A. M. Oliveira, J. H. Lehman, A. B. McCoy and W. C. Lineberger, J. Chem. Phys., 2016, 145, 124317.
50. D. R. Stull, Ind. Eng. Chem., 1947, 39, 517–540.
51. E. Lemmon, M. McLinden, D. Friend, P. Linstrom and W. Mallard, NIST Chemistry WebBook, NIST Standard Reference Database Number 69 National Institute of Standards and Technology, Gaithersburg, 2011.
52. I. Nadler, D. Mahgerefteh, H. Reisler and C. Wittig, J. Chem. Phys., 1985, 82, 3885–3893.
53. W. Marinelli, N. Sivakumar and P. Houston, J. Phys. Chem., 1984, 88, 6685–6692.
54. A. P. Baronavski and J. McDonald, Chem. Phys. Lett., 1977, 45, 172–176.
55. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
56. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007–1023.
57. R. A. Kendall, T. H. Dunning Jr. and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
58. D. E. Woon and T. H. Dunning Jr., J. Chem. Phys., 1995, 103, 4572–4585.
59. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
60. J. A. Pople, R. Krishnan, H. B. Schlegel and J. S. Binkley, Int. J. Quantum Chem., 1978, 14, 545–560.
61. R. J. Bartlett and G. D. Purvis, Int. J. Quantum Chem., 1978, 14, 561–581.
62. M. L. Laury, S. E. Boesch, I. Haken, P. Sinha, R. A. Wheeler and A. K. Wilson, J. Comput. Chem., 2011, 32, 2239–2347.
63. NIST Computational Chemistry Comparison and Benchmark Database, http://cccbdb.nist.gov/.
64. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Revision D.01 edn., Gaussian, Inc., Wallingford CT, 2013.
65. R. A. Marcus, J. Chem. Phys., 1952, 20, 359–364.
66. J. W. Davies, N. J. B. Green and M. J. Pilling, Chem. Phys. Lett., 1986, 126, 373–379.
67. J. Troe, J. Chem. Phys., 1977, 66, 4745–4757.
68. W. H. Miller, J. Am. Chem. Soc., 1979, 101, 6810–6814.
69. J. T. Bartis and B. Widom, J. Chem. Phys., 1974, 60, 3474–3482.
70. K. M. Douglas, D. I. Lucas, C. Walsh, N. A. West, M. A. Blitz and D. E. Heard, Astrophys. J., Lett., 2022, 937, L16.
71. B. D'Anna, V. Bakken, J. Are Beukes, C. J. Nielsen, K. Brudnik and J. T. Y. Jodkowski, Phys. Chem. Chem. Phys., 2003, 5, 1790–1805.
72. H.-Y. Li, M. Pu, Y.-Q. Ji, Z.-F. Xu and W.-L. Feng, Chem. Phys., 2004, 307, 35–43.
73. V. Vasudevan, D. F. Davidson and R. K. Hanson, Int. J. Chem. Kinet., 2005, 37, 98–109.
74. S. Xu, R. S. Zhu and M. C. Lin, Int. J. Chem. Kinet., 2006, 38, 322–326.
75. Y. Zhao, B. Wang, H. Li and L. Wang, THEOCHEM, 2007, 818, 155–161.
76. D. W. Marquardt, J. Soc. Ind. Appl. Math., 1963, 11, 431–441.
77. D. McElroy, C. Walsh, A. J. Markwick, M. A. Cordiner, K. Smith and T. J. Millar, Astron. Astrophys., 2013, 550, A36.
78. V. Wakelam, E. Herbst, J.-C. Loison, I. Smith, V. Chandrasekaran, B. Pavone, N. Adams, M.-C. Bacchus-Montabonel, A. Bergeat and K. Béroff, Astrophys. J., Suppl. Ser., 2012, 199, 21.
79. T. Yang, A. M. Thomas, B. B. Dangi, R. I. Kaiser, A. M. Mebel and T. J. Millar, Nat. Commun., 2018, 9, 774.
80. J. Cernicharo, C. Cabezas, Y. Endo, M. Agúndez, B. Tercero, J. R. Pardo, N. Marcelino and P. De Vicente, Astron. Astrophys., 2021, 650, L14.
81. S. Doddipatla, C. He, S. J. Goettl, R. I. Kaiser, B. R. L. Galvão and T. J. Millar, Sci. Adv., 2021, 7, eabg7003.
82. M. Van de Sande, J. O. Sundqvist, T. J. Millar, D. Keller, W. Homan, A. de Koter, L. Decin and F. De, Ceuster, Astron. Astrophys., 2018, 616, A106.
83. M. Agúndez, J. Cernicharo, M. Guélin, C. Kahane, E. Roueff, J. Kłos, F. J. Aoiz, F. Lique, N. Marcelino, J. R. Goicoechea, M. G. García, C. A. Gottlieb, M. C. McCarthy and P. Thaddeus, Astron. Astrophys., 2010, 517, L2.

---

Footnotes

† Electronic supplementary information (ESI) available: The experimental method (Sections S1–S3), computational methods using Gaussian (Section S4), MESMER (Section S5) and capture theory (Section S6), parameterization of the rate coefficient (Section S7) and results of the astrochemical modelling (Section S8). See DOI: https://doi.org/10.1039/d2cp05043a

‡ Current address: Andor Technology Ltd., Belfast, BT12 7AL, UK.

§ Current address: School of Chemistry, University of Birmingham, Edgbaston, B15 2TT, UK.

---