A&A 376, 663-666 (2001)
DOI: 10.1051/0004-6361:20010998

On the deuteration of C $\mathsf{_{3}}$H $\mathsf{_{2}}$ in dense interstellar clouds via
the reaction C $\mathsf{_{3}}$H $\mathsf{_{3}^{+}}$ + HD $\longrightarrow$ C $\mathsf{_{3}}$H $\mathsf{_{2}}$D⁺ + H $\mathsf{_{2}}$

D. Talbi¹ - E. Herbst²

1 - LETMEX, Université de Nice Sophia Antipolis, Parc Valrose, UFR Sciences, 06108 Nice Cedex 2, France
2 - Departments of Physics and Astronomy, The Ohio State University, Columbus, OH 43210, USA

Received 16 March 2001 / Accepted 6 July 2001

Abstract
Although deuterium isotope fractionation in cold, dark interstellar clouds is reasonably well understood via gas-phase chemistry, there are some discrepancies between observation and theory. For example, the observed abundance ratio between the deuterated species C₃HD and the cyclic molecule C₃H₂ is significantly higher than most theoretical values unless the exchange reaction C₃H₃⁺+ HD $\longrightarrow$ C₃H₂D⁺+H₂ is efficient. In this paper, we report quantum chemical and dynamical calculations on this reaction, which show it to possess a large activation energy barrier and to be very slow at all normal interstellar temperatures.

Key words: ISM: abundances - ISM: molecules - molecular processes

1 Introduction

Deuterium isotopic fractionation is a continuing topic of interest in the study of dense interstellar clouds. Although the deuterium to hydrogen elemental abundance ratio is only ${\approx} 10^{-5}$ (e.g. Wright et al. 1999), the abundance ratios of singly deuterated trace isotopomers such as H₂D⁺, DCO⁺, NH₂D, DCN, etc. to their respective normal species are measured to be 0.01-0.10 or even higher (depending on excitation model) in such clouds. Indeed, the doubly deuterated isotopomer NHD₂ has recently been detected in the dark cloud L134N with an abundance ratio NHD₂/NH₂D of ${\approx} 0.05$ (Roueff et al. 2000), which is similar to the value measured for NH₂D/NH₃ of ${\approx} 0.1$ (Tiné et al. 2000).

It is well known that for cold dark clouds, exothermic exchange reactions between the abundant ions H₃⁺, CH₃⁺ and C₂H₂⁺ and HD are very important in producing enhanced abundances of deuterated ions, which then undergo further reactions to produce other deuterium-containing species. Models for dark clouds including these and other reactions have been reported by Millar et al. (1989) and more recently by Roberts & Millar (2000a,2000b) and by Turner (2001). In these models, gas-phase reactions generally can account for the abundances of singly and doubly deuterated isotopomers, although the agreement is not uniform. One of the problems appears to be the large deuterium fractionation involving the cyclic molecule C₃H₂. Early measurements by Bell et al. (1988) obtained that C₃HD/C₃H $_{2}\approx 0.08{-}0.16$ depending upon position in TMC-1. More recent measurements by Turner (2001) in TMC-1 show that the ratio appears to be closer to 0.048. Calculated values tend to be on the order of ${\approx} 0.01$ (Millar et al. 1989; Roberts & Millar 2000a; Turner 2001) unless there is significant freeze out of molecules onto dust (Roberts & Millar 2000a,2000b), in which case all abundance ratios of deuterated to normal isotopomers increase. To complicate matters though, a somewhat higher value of 0.032 has been calculated by Turner (2001) for TMC-1 at early time.

In 1993, Howe & Millar (1993) suggested that if the exothermic left-to-right reaction of the reaction system

$$% {\rm C_{3}H_{3}^{+} + HD \Longleftrightarrow C_{3}H_{2}D^{+} + H_{2}},$$ (1)

occurs rapidly at low temperatures, the predicted amount of C₃HD can be enhanced to its observed value because dissociative recombination of C₃H₂D⁺ produces C₃HD directly. To the best of our knowledge, this suggestion has not yet been tested theoretically or in the laboratory. Unlike most ion-molecule reactions, reactions between ions and HD or H₂ can possess activation energy barriers (Henchman et al. 1988). In this paper, we report quantum chemical calculations on the reaction system in Eq. (1) and show that it indeed possesses a large activation energy barrier and thus cannot be invoked to reduce the likely disparity between theory and observation for the abundance ratio C₃HD/C₃H₂.

2 Quantum chemical calculations

An extensive investigation of the singlet potential energy surface corresponding to an HD molecule approaching the c-C₃H₃⁺ion was initially undertaken at the RHF and MP2 levels in combination with a triple split-valence basis set to which diffuse orbitals have been added together with polarization functions (namely 6-311+G(2d,p)). The stationary points located on this surface were reoptimized using the density functional formalism (B3LYP functional; Becke 1993) and their character (either minimum energy point, for which all vibrational frequencies are real, or saddle point, characterized by one imaginary frequency) was confirmed by vibrational analysis done at the same B3LYP level. For accurate electronic energies, single point calculations were performed at the CCSD(T)/6-311+G(3df,2pd) theoretical level using the B3LYP/6-311+G(2d,p) optimized geometries (hereafter CCSD(T)/6-311+G(3df,2pd)//B3LYP/6-311+G(2d,p)). This is a coupled cluster singles and doubles method with a perturbative treatment of the triple excitations (Raghavachari et al. 1989), employing the previous triple-zeta basis set to which extra sets of polarization functions have been added. With such a basis set we reach a flexibility in the single determinant representation that approaches the Hartree-Fock limit. All energies were corrected for zero-point energy (ZPE) contributions calculated at the B3LYP/6-311+G(2d,p) level and scaled by a factor of 0.9613 (Frisch et al. 1996). When relative energies were small (for the case of the long-distance complexes in the entrance and exit channels) they were adjusted for the basis set superposition error (BSSE) following the counterpoise procedure (Boys & Bernardi 1970; Liu & McLean 1973, 1989). All the quantum chemical calculations of this study were performed by means of the Gaussian 98 package (Frisch et al. 1998). Uncertainties in calculated energies at the highest level of calculation are estimated to be 1-2 kcal/mol.

	Figure 1: CCSD(T)/6-311+G(3df,2pd//B3LYP/6-311+G(2d,p) lowest energy reaction pathway for the reaction system C3H3+ + HD $\leftrightarrow$ C3H2D+ + H2. ZPE corrections are included as well as BSSE corrections for the long distance complexes.
Open with DEXTER

	Figure 2: Structures of the (C3H4D+) stationary points optimized at the B3LYP/6-311+G(2d,p) level. The distances are in Å.
Open with DEXTER

The CCSD(T)/6-311+G(3df,2pd)//B3LYP /6-311+ G(2d,p) energy profile of the c-C₃H₃⁺ + HD reaction is shown in Fig. 1 and the structures of the corresponding (C₃H₄D⁺) stationary points, optimized at the B3LYP/6-311+G(2d,p) level, are given in Fig. 2. The energy profile reveals the existence on the reaction path of a high transition state, labelled (C₃H₂... H₂D⁺)$_{\rm TS}$ in Fig. 2, which can be described as a complex between c-C₃H₂ and H₂D⁺. This transition state corresponds to the minimum energy needed to go from reactants to products and is best thought of as a saddle point in the potential barrier. Its vibrational analysis shows one imaginary frequency of 1658 cm^-1 corresponding to a stretching of the two top hydrogens toward dissociation. The shallow minimum in the entrance channel corresponds to a long-distance complex between c-C₃H₃⁺ and HD (see Fig. 2). Its binding energy after ZPE and BSSE corrections are made is -0.4 kcal/mol, which is on the order of the uncertainty in our calculations. The presence of a minimum energy structure on this part of the potential energy surface is therefore not proved and needs more investigation. However, whether this loose complex exists or not is of no importance for the present study considering the height ( $80\,\pm\,2$ kcal/mol above the reactants) of the following transition state barrier. This argument is equally valid for the loose (C₃H₂D⁺-H₂) complex in the exit channel, which corresponds to another shallow minimum in the potential energy profile. Table 1 contains vibrational frequencies and rotational constants for the three intermediate stationary points.

 

Table 1: Computed harmonic frequencies and rotational constants of intermediate structures.

Species	Harmonic frequencies a	Rotational constants b
	(cm-1)	(GHz)
C3H $_{3}^{+}\cdot\cdot\cdot$HD	$a^{\prime} = 64$, 147, 312, 903, 920	30.26996
	1016, 1259, 1261, 1590, 3082	8.00741
	3113, 3151, 3653	6.33231
	$a^{\prime \prime} = 7$, 82, 749, 965, 986
(C3H $_{2}\cdot\cdot\cdot$H2D+)$_{\rm TS}$	$a^{\prime} = 1658i$, 178, 566, 908, 934	30.88821
	1110, 1289, 1406, 1588, 2379	11.11820
	2930, 3115, 3147	8.17545
	$a^{\prime \prime} = 189$, 236, 790, 897, 956
(C3H2D $^{+}\cdot\cdot\cdot$H2)	$a^{\prime} = 68$, 181, 355, 717, 900	30.43985
	984, 1238, 1252, 1546, 2329	9.90485
	3113, 3144, 4217	7.47315
	$a^{\prime \prime} = 9$, 81, 652, 885, 965

         ^a Harmonic frequencies have been calculated at the B3LYP/6-311+G(2d,p) level and scaled by 0.9613.
         Group theoretical designations refer to the symmetry of the vibrational modes. $i = \surd-1$.
         ^b Rotational constants have been calculated at the B3LYP/6-311+G(2d,p) level.

3 Rate coefficients

The rate coefficient for the reaction between C₃H₃⁺ and HD (or the reverse reaction) can be calculated with "activated complex'' theory (Weston & Schwarz 1972). In this method, an equilibrium is assumed to occur between the reactants and the transition state (see Fig. 1), while the weakly-bound entrance and exit channels complexes are ignored. Activated complex theory does not give us the probability that the transition state decays to products rather than reactants, but at low energies, the exothermic channel is preferred. With the assumption of decay into exothermic products, the rate coefficient for the left-to-right reaction in system (1) is given by the equation

$$% k_{1;\rm L \rightarrow R} = (kT/h) R_{\rm trans}R_{\rm rot... ... vib}R_{\rm elec} \Gamma_{\rm t} {\rm {\exp}}(-E_{\rm ts}/T)$$ (2)

where k is the Boltzmann constant, h is Planck's constant, T is the temperature, $E_{\rm ts}$ is the energy in K of the transition state relative to the energy of reactants (with zero-point energies included), and the R's are ratios between partition functions of the transition state and of the reactants for translational, rotational, vibrational, and electronic modes (Weston & Schwarz 1972). A tunneling correction $\Gamma_{\rm t}$ also appears as a factor (Herbst 1996).

The dominant term in this expression is the exponential one. From Fig. 1, it can be seen that $E_{\rm ts} = 80$kcal/mol ${\approx} 40\,000$ K. Thus even at temperatures approaching room temperature (300 K), the very large negative value of the exponent ( $-E_{\rm ts}/T = -134$) renders the rate coefficient effectively zero whatever the values of the other terms. Nevertheless, for completeness we have evaluated the other factors in Eq. (2) (Herbst & Talbi 1998). At 10 K we obtain

$$% k_{1;\rm L \rightarrow R}\;({\rm cm^{3}\; s^{-1}}) = 2.4\times 10^{-7} \exp(-4030)$$ (3)

while at 300 K

$$% k_{1;\rm L \rightarrow R}\;({\rm cm^{3}\; s^{-1}}) = 6.3\times 10^{-11} \exp(-134).$$ (4)

One can see that the process is not important in interstellar clouds, since a rate coefficient near the Langevin value (10^-9 cm³ s^-1) is required. Reaction system (1) should thus not be used in gas-phase fractionation models.

Acknowledgements

The Astrochemistry Program at The Ohio State University is supported by The National Science Foundation. Parts of the calculations reported here were supported by the IDRIS (France); the support is gratefully acknowledged.

References

 
• Becke, A. D. 1993, J. Chem. Phys., 98, 5648 In the text NASA ADS  
• Bell, M. B., Avery, L. W., Matthews, H. E., et al. 1988, ApJ, 326, 924 In the text NASA ADS  
• Boys, S. F., & Bernardi, F. 1970, Mol. Phys., 19, 553 In the text  
• Frisch, M. J., Trucks, G. W., & Cheeseman, J. R. 1996, in Recent Developments and Applications of Modern Density Functional Theory, Theoretical and Computational Chemistry, ed. J. Seminario (Amsterdam: Elsevier), 4, 679 In the text  
• Frisch, M. J., Trucks, G. W., Schlegel, H. B., et al. 1998, Gaussian 1998 (Revision A.5), (Pittsburgh: Gaussian, Inc.) In the text  
• Henchman, M. J., Paulson, J. F., Smith, D., Adams, N. G., & Lindinger, W. 1988, in Rate Coefficients in Astrochemistry, ed. T. J. Millar, & D. A. Williams (Dordrecht: Kluwer), 201 In the text  
• Herbst, E. 1996, in Atomic, Molecular, & Optical Physics Handbook, ed. G. W. F. Drake (Woodbury, NY: AIP Press), 429 In the text  
• Herbst, E., & Talbi, D. 1998, A&A, 338, 1080 In the text NASA ADS  
• Howe, D. A., & Millar, T. J. 1993, MNRAS, 262, 868 In the text NASA ADS  
• Liu, B., & McLean, A. D. 1973, J. Chem. Phys., 59, 4557 In the text  
• Liu, B., & McLean, A. D. 1989, J. Chem. Phys., 91, 2348 In the text NASA ADS  
• Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906 In the text NASA ADS  
• Raghavachari, K., Trucks, G. W., Pople, J. A., & Head-Gordon, M. 1989, Chem. Phys. Lett., 157, 479 In the text  
• Roberts, H., & Millar, T. J. 2000a, A&A, 361, 388 In the text NASA ADS  
• Roberts, H., & Millar, T. J. 2000b, A&A, 364, 780 In the text NASA ADS  
• Roueff, E., Tiné, S., Coudert, L. H., et al. 2000, A&A, 354, L63 In the text NASA ADS  
• Tiné, S., Roueff, E., Falgarone, E., Gerin, M., & Pineau des Forêts, G. 2000, A&A, 356, 1039 In the text NASA ADS  
• Turner, B. E. 2001, ApJS, in press In the text  
• Weston, Jr., R. E., & Schwarz, H. A. 1972, Chemical Kinetics (Englewood Cliffs, NJ: Prentice-Hall) In the text  
• Wright, C. M., van Dishoeck, E. F., Cox, P., Sidher, S. D., & Kessler, M. F. 1999, ApJ, 515, L29 In the text NASA ADS