EDP Sciences
Free Access
Issue
A&A
Volume 495, Number 2, February IV 2009
Page(s) 513 - 521
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361:200810967
Published online 14 January 2009

A sensitivity study of the neutral-neutral reactions C + C$_{\sf 3}$ and C + C$_{\sf 5}$ in cold dense interstellar clouds

V. Wakelam1,2 - J.-C. Loison3 - E. Herbst4,5 - D. Talbi6 - D. Quan7 - F. Caralp3


1 - Université de Bordeaux, Laboratoire d'Astrophysique de Bordeaux, BP 89 33271 Floirac Cedex, France
2 - CNRS/INSU, UMR 5804, BP 89 33271 Floirac Cedex, France
3 - Institut des Sciences Moléculaires CNRS UMR 5255, Université Bordeaux 1, 33405 Talence, France
4 - Department of Physics, The Ohio State University, Columbus, OH 43210, USA
5 - Departments of Astronomy and Chemistry, The Ohio State University, Columbus, OH 43210, USA
6 - Université Montpellier II, Groupe de Recherches en Astronomie et Astrophysique du Languedoc, CNRS, UMR 5024, place Eugène Bataillon, 34095 Montpellier, France
7 - The Chemical Physics Program, The Ohio State University, Columbus, OH 43210, USA

Received 15 September 2008 / Accepted 19 December 2008

Abstract
Aims. Chemical networks used for models of interstellar clouds contain many reactions, some of them with poorly determined rate coefficients and/or products. In this work, we report a method for improving the predictions of molecular abundances using sensitivity methods and ab initio calculations.
Methods. Based on the chemical network osu.2003, we used two different sensitivity methods to determine the most important reactions as a function of time for models of dense cold clouds. Of these reactions, we concentrated on those between C and C3 and between C and C5, both for their effect on specific important species such as CO and for their general effect on large numbers of species. We then used ab initio and kinetic methods to determine an improved rate coefficient for the former reaction and a new set of products, plus a slightly changed rate coefficient for the latter.
Results. Putting our new results in a pseudo-time-dependent model of cold dense clouds, we found that the abundances of many species are altered at early times, based on large changes in the abundances of CO and atomic C. We compared the effect of these new rate coefficients/products on the comparison with observed abundances and found that they shift the best agreement from $3 \times 10^4$ yr to $(1{-}3)\times 10^5$ yr.

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

1 Introduction

The chemistry in astrochemical models is driven by large networks including thousands of reactions with poorly determined rate coefficients. For several years now, astrochemists have been working on the quantification of model uncertainties caused by uncertainties in rate coefficients. Using similar methods but different models and networks, Vasyunin et al. (2004,2008) and Wakelam et al. (2006a,2005) have studied the error propagation of such uncertainties in various environments such as dense clouds, diffuse gas, hot cores, and protoplanetary disks. This approach to quantify the model error is useful quantitatively comparing calculated and observed molecular abundances. In addition to this important aspect, these authors have developed sensitivity methods for identifying those reactions in complete networks that might be worth studying theoretically or experimentally because model results are very sensitive to their rate coefficients and uncertainties. Identifying these reactions is critical because detailed experimental or theoretical studies of a single reaction, especially at low interstellar temperatures, can take years. Guiding chemists and physicist as to which reactions are the most important for them to study is an efficient method for reducing error in chemical models.

``Important'' reactions can be defined in several ways. First, one can define important reactions for specific species; such reactions strongly affect the abundances of individual species through changes in their rate coefficients for formation or destruction processes. One must be precise with this definition, because even though some reactions dominate the formation or destruction of certain species, changes in their rate coefficients within certain ranges can have little effect. Quan et al. (2008), for instance, showed that although molecular oxygen is formed by the reaction  $\rm O + OH \rightarrow O_2 + H$, the rate coefficient of this reaction can be changed by several orders of magnitude without affecting the O2 abundance significantly. The reason is that the reaction is also the main route of OH destruction. When the rate coefficient is decreased, the OH abundance increases so that the rate of O2 formation stays constant. As a counterexample, the reaction  $\rm C + C_{3} \rightarrow C_{4} + h\nu$ is very important to the abundance of CO, as will be discussed below. Secondly, ``important'' reactions can be defined through the number of molecular abundances they influence strongly. As an example, the cosmic-ray ionization of molecular hydrogen is important for almost all of the species in the chemical models of dense interstellar clouds (Wakelam et al. 2006b; Vasyunin et al. 2008).

In this paper, we present a complete approach to improving chemical models through identifying sensitive reactions followed by an increase in our knowledge of their rate coefficients via theoretical means, and using the newly revised rate coefficients in model predictions. We direct our focus on two reactions: the radiative association reaction between $\rm C$ and $\rm C_3$ and the neutral-neutral reaction  $\rm C + C_5$. We start, in Sect. 2, with a discussion of the methods used for our sensitivity analyses and the determination of the important reactions in cold dense clouds. The theoretical study of the  $\rm C + C_3$ and  $\rm C + C_5$ rate coefficients is presented in Sect. 3 and the impact of new rate coefficients on the model predictions are studied in Sect. 4. Following our conclusion, an Appendix on other reactions of the $\rm C + C_{n}$ class is presented.

2 Selection of important reactions

In order to identify important reactions of both classes, we use two different sensitivity methods and compare the results of the two. Each sensitivity method is based on linear correlations but in the first one, only one reaction is modified at a time whereas in the second method all reactions are modified at the same time. The first one, which we will call ``individual perturbations'', is simply based on the study of the linear effect of individual rate coefficient perturbations on molecular abundances. We have already used this method previously (Wakelam et al. 2006a,2005). In the second method, we compute Pearson correlation coefficients.

Table 1:   Most important reactions for selected species at 105 yr.

2.1 Chemical model

To represent cold dark clouds, we use a pseudo- time-dependent gas-phase model, which computes abundances as a function of time for a temperature of 10 K, a total hydrogen density of $2\times 10^4$ cm-3 and a visual extinction of 10. The model is known as nahoon (http://www.obs.u-bordeaux1.fr/amor/VWakelam/). For comparison with previous studies, the elemental abundances are the O-rich low-metal elemental abundances (see Table 1 in Wakelam & Herbst 2008). For this work, we used an older version of the osu database (osu.2003) since one goal of this work is to show how we corrected a mistake in this database. We focus on the specific time of 105 yr, since this so-called ``early time'' normally results in the best agreement between calculated abundances and observations (Smith et al. 2004).

2.2 Sensitivity methods

The idea of the individual perturbation approach is to look at the individual effect of each reaction. With a network of N reactions, we run the model N times and modify only one rate coefficient each run. One then obtains the chemical evolution of the fractional abundance: Xij(t) for each species j and run i . Note that the index i can stand for the run or the reaction with varied rate coefficient. We then compare the abundance of each species Xij(t) with the standard abundance $X^{\rm ref}_j(t)$ obtained without changing the rate coefficients. In our previous uses of this approach to select important reactions for hot-core (Wakelam et al. 2005) and dense-cloud (Wakelam et al. 2006a) chemistry, we multiplied each rate coefficient by its uncertainty factor. This approach assumes that the real uncertainty in the rate coefficient is not much larger than the considered uncertainty factor, which may not be true. Indeed, as we will show later, this approach can miss important reactions if the real rate coefficient error is much larger than the uncertainty factor. For this reason, here we multiply each rate coefficient by the same factor of two. This factor is currently the default uncertainty used for most classes of reactions in current reaction networks (http://www.udfa.net/; http://www.physics.ohio-state.edu/~eric/research.html). To quantify the effect of these variations in the reaction rate coefficients on the abundances, we compute a parameter labeled Rij(t) for each species and each reaction, where

\begin{displaymath}
R^i_j(t)=\frac{(X^i_j(t)-X^{\rm ref}_j(t))}{X^{\rm ref}_j(t)}\cdot
\end{displaymath} (1)

As an example, a value of Rij(t)=0.2 tells us that the multiplication of the reaction i by a factor of two increases the abundance of species j by 20%. We also did these calculations by multiplying the rate coefficients by a factor of 10 or dividing them by 2 but the list of important reactions is rather the same as for the case presented here.

The second method is based on linear correlations but is less time-consuming. Here we run about 2000 different models l, randomly varying all rate coefficients within an uncertainty range of a factor of two, with both increases and decreases considered. We then compute Pearson correlation coefficients for each species j and each reaction i using the following relation:

\begin{displaymath}P^i_j(t)=\frac{\sum^l (X^l_j(t)-\bar{X_j(t)})(k^l_i-\bar{k_i}...
...sum^l (X^l_j(t)-\bar{X_j(t)})^2 \sum^l (k^l_i-\bar{k_i})^2}},
\end{displaymath} (2)

where $\bar{X_j(t)}$ and $\bar{k_i}$ are the mean values of the abundance for species j and the rate coefficient i. Pearson correlation coefficients indicate the strength and direction of a linear relationship between two variables by dividing the covariance of the two variables by the products of their standard deviations (see for instance Press & Flannery 1992).

2.3 Results of sensitivity calculations

Table 2:   Reactions influencing the most species at 105 yr using R.

Table 3:   Reactions influencing the most species at 105 yr using P.

In this section, we show the results of both sensitivity methods to determine the most important reactions for the abundances of selected individual species and to determine the reactions that influence the abundances of the largest numbers of species.

In Table 1, we show the most important reactions for the species CO, O2, C4H and HC3N at an integrated time of 105 yr using the individual perturbation (R) and the correlation coefficient (P) methods. In Tables 2 and 3, we list the reactions that influence the most species using R and P, respectively, as a criterion. For each reaction, we count the number of species for which R or P is either larger than 0.1 or smaller than -0.1.

A glance at Table 1 shows that the important reactions governing the abundances of individual species are generally similar for both methods. There are two main surprises, however. First, for the case of CO, the values of R are much smaller than those for P. Presumably, the explanation for this dichotomy is that most of the carbon goes into CO at that time. As a result, a change in the rate coefficients does not change the CO abundance, which is then controlled by the elemental abundance of C. This effect also explains the small uncertainty in the CO abundance found by Wakelam et al. (2006a). Secondly, the most unexpected reactions, because they seem so unrelated to the listed species, are those between C and $\rm C_3$ and between C and $\rm C_5$. Note that the latter was assumed, erroneously, in osu.2003 to be a radiative association (see below). The importance of the C + $\rm C_3$ association seems to be associated with the fact that C3 acts as a catalyst to convert C and O into CO via the association reaction between C and C3 followed by

\begin{displaymath}{\rm O + C_{4} \longrightarrow CO + C_{3},}
\end{displaymath} (3)

which regenerates the C3. A similar cycle happens with C5 if it undergoes radiative association with C.

Considering the number of species affected, the most important reactions determined by both methods are quite similar although the order of importance can differ. As already noticed by several authors (Vasyunin et al. 2004; Wakelam et al. 2006b,a; Vasyunin et al. 2008), the H2 and He ionization rates by cosmic-rays appear to be important for many species since they are at the origin of formation and destruction pathways. Most of the other important reactions are ion-neutral and associative reactions and dissociative recombination reactions, such as those for H3O+ with electrons. Among the important reactions, $\rm C^+ + H_2 \longrightarrow CH_2^+ $, $\rm CH_3^+ + H_2 \longrightarrow CH_5^+$, $\rm H_3^+ + C \longrightarrow CH^+ + H_2$, $\rm H_3^+ + CO \longrightarrow HCO^+ + H_2 $, $\rm H_3^+ + O \longrightarrow OH^+ + H_2 $ have also been identified as important reactions for disk chemistry by Vasyunin et al. (2008). Most of these reactions are at the start of a chain of reactions that builds up hydrocarbons. It should be emphasized that our sensitivity methods need not tell us the details of the pathway of reaction cycles. If, for example, a destruction reaction leads to a loop of subsequent reactions that reforms the molecule being destroyed with some non-zero efficiency, the sensitivity method might just decide that the initial destruction reaction is less important, at least for the species being destroyed. The problem of enhanced molecular abundances due to artificial loops has been minimized whenever possible in the construction of the network by including a number of different product channels. For example, if the major destruction route of a neutral species leads to a protonated ion, and the major destruction route of the ion is dissociative recombination to form neutral molecules, the channel that loops back to the original neutral species plus a hydrogen atom will not be assumed to occur on more than 50% of the recombination reactions (unless, as is highly unusual, experimental results tell us the contrary). Let us look at the important reactions causing the destruction of HC3N in Table 1. One of them is the rather direct destruction via C+ ions. The products of this reaction do not lead back to HC3N in any direct fashion. If, on the other hand, there were a charge exchange reaction to form C3HN+, and this ion were hydrogenated by reaction with H2 to form, C3H2N+ + H, the hydrogenated ion could react with electrons to form HC3N + H. In the network, however, the probability of this neutral product channel is set to approximately 50% so that the loop is not particularly efficient.

The only neutral-neutral reactions that seem important for many species in Tables 2 and 3 are the (presumed) radiative associative reactions $\rm C + C_3$ and $\rm C + C_5$. The rate coefficients for these reactions used in osu.2003 are 10-10 and 10-11 cm3 s-1 based on Smith et al. (2004), without temperature dependence. These reactions have been looked at again by Smith et al. (2006) and it was realized that the reaction between C and $\rm C_5$ has the normal exothermic products $\rm C_3 + C_3$. Given the importance of these reactions, we decided to study them theoretically in some detail.

3 New calculations for C + C$_{\sf 3}$ and C + C$_{\sf 5}$ reactions

3.1 C + C 3

At low collision energies, the reaction of C(3P) atoms with linear C3 in its ground electronic state $(\rm X~^{1}\Sigma^{+}_{g})$ leads only to adduct formation, C4, because the exit channel C2 + C2 is endothermic. The lowest energy pathway to produce C4 is the linear one, and we focus on this pathway here. Under interstellar low-density conditions, the only way to stabilize the C4 adduct is by emission of a photon. The formation of the linear C4 adduct is strongly exothermic and spin-allowed to its triplet manifold of states, designated $^{3}\rm C_{4}$:

\begin{displaymath}\rm C\left(^3P\right) + ^1C_3 \longrightarrow ^{3}C_4 + h\nu.
\end{displaymath} (4)

There are a variety of values for this exothermicity. A high-temperature experiment yields a value of ( $5.21 \pm 0.21$) eV from the heats of formation of C4 and C3 (Gingerich et al. 1994), a PST calculation based on a photodissociation experiment yields a lower value of ( $4.84 \pm 0.15)$ eV (Choi et al. 2000), while high-level ab initio theoretical values (including our own work) computed at the RCCSD and RCCSD(T) levels range from 4.94-4.98 eV without zero-point correction (Martin & Taylor 1995). The high-temperature experimental value has been used in the RRKM calculations discussed below. The rhombic isomer of C4 in its ground singlet state has been calculated to to be nearly isoenergetic with the linear isomer (Masso et al. 2006), but the pathway from C and linear C3 to the first excited triplet state shows a high barrier at the CASSCF level of ab initio theory (see below).

In the $\rm C_{\infty,v}$ symmetry, the three electronic states of C4 correlating to the three spin-orbit states of C(3P0,1,2) lead to two different potential energy surfaces corresponding to the fundamental state ${\rm X~^3 \Sigma}^-_ g$ and the first excited $^3 \Pi _g$ state. (Note that the ground electronic state of C is of g symmetry). The energetic and structural parameters of these states of C4 and the second excited $^3\Pi_u$ state are listed in Table 4. These parameters and others such as the pathway to the rhombic isomer of C4 have been calculated at the Restricted Hartree-Fock (RHF) and Complete Active Space Self Consistent Field (CASSCF) level with the double-zeta cc-pVDZ basis set of Dunning (1989), using MOLPRO 2000, a package of ab initio programs designed by Werner, and Knowles (see http://www.molpro.net for more details). This level of theory was also used to determine if the three electronic curves are attractive along the linear reaction pathway. In these calculations, which underestimate the exothermicity of reaction (4), the C $_{\rm\infty v}$ (linear) approach of C toward the C3 axis, shown in Fig. 1, yields no barrier for the ground C $_4(^3 \Sigma ^-_g$) electronic state but shows a high barrier of 1.5 eV for the first excited C4($^3 \Pi _g$) state. Using $\Omega$ and MJ conservation, we obtain that only the C($\rm ^3P_0$) and C($\rm ^3P_1$, $M_{\rm J}=-1,+1$) fine-structure substates are reactive; these three substates form vibrationally excited linear C4 in its fundamental $^3 \Sigma^-_g$ state, while the other six substates lead to the formation of the $^3 \Pi _g$ excited state (through a high barrier). Any others approach on the triplet surface, as the perpendicular one leading to the rhombic C4, shows high barriers at the CASSCF/cc-pVDZ level. As the C4 is formed with a great excess of energy, the potential minima of several electronic states, not just the ground state, lie below the dissociation limit (the potential energy of ground-state reactants.). Stabilisation of C4 adduct can then occur via purely vibrational transitions from the ground-state adduct, but also via electronic transitions from excited state adducts if they can be formed without barriers (Herbst & Bates 1988). Here, we must therefore consider three mechanisms for adduct stabilisation:

  • vibrational radiative transitions from quasi-continuum vibrational levels in the ground state $\rm { (X~^3\Sigma}^-_g)$ adduct to lower bound vibrational levels;
  • ``inverse internal conversion'' from the ${\rm X~^3 \Sigma}^-_ g$ adduct to an excited electronic state followed by electronic photon emission to bound vibrational levels of the ground state (case 1);
  • direct radiative emission from quasi-continuum vibrational levels in the ground ${\rm X~^3 \Sigma}^-_ g$ state to lower vibrationally bound levels in excited states (case 2).
Other possibilities arise from the isomerisation of the triplet linear C4 state toward the higher-energy singlet rombric C4 state. However this isomerisation has to overcome a high barrier as well as intersystem crossing. We can reasonably neglect those possibilities.

Table 4:   Vdz/CASSCF optimized structural parameters.

 \begin{figure}
\par\includegraphics[width=8cm,clip]{0967fig1.eps}
\end{figure} Figure 1:

Co-linear VDZ-CASSCF one-dimensional cuts of the lowest triplet electronic states of linear C4 molecules along the C3-C stretch. The C3 bonds are optimized. The labeling of the C4 electronic states is based on D $_{\rm\infty h}$ symmetry.

Open with DEXTER

Of the three itemized mechanisms, the case 2 situation can be ignored. Both experiments and detailed calculations on our part show that the potential surfaces of all relevant low-lying excited states of linear C4 lie almost always above the ground state at all relevant geometries so that emission transition from the ground state are impossible or very inefficient. Therefore, the only possibility for stabilization via electronic emission is case 1, with the dominant channel being potential curve crossing from ${\rm X~^3 \Sigma}^-_ g$ to $^3\Pi_u$ at long C2-C2 distances (see Masso et al. 2006, for example) followed by an allowed $ u \rightarrow g$ electronic transition. The u-g potential curve crossing is not allowed when C4 lies in the symmetric D $_{\rm\infty h}$ linear configuration, but is allowed by vibrational excitation, which reduces the symmetry. Since the $^3\Pi_u$ state has an electronically allowed transition not only with the fundamental state ${\rm X~^3 \Sigma}^-_ g$ but also with the first electronically excited $^3 \Pi _g$ state, we have to include these two transitions.

If we designate unstable adducts (``complexes'') with an asterisk, the overall mechanism for association can now be written as:

$\displaystyle \rm C + C_3 \rightarrow{\textit{k}_{\rm ass}} \left(C_4\left(^3\Sigma^-_{\it g}\right)\right)^* \rightarrow{\textit{k}_{\rm diss}} C+ C_3$     (5)


$\displaystyle \rm\left(C_4\left(^3\Sigma^-_{\it g}\right)\right)^* \rightarrow{\textit{k}_{rad,v}} C_4\left(^3\Sigma^-_{\it g}\right) + h\nu$     (6)


$\displaystyle \rm\left(C_4\left(^3\Sigma^-_{\it g}\right)\right)^*\rightlefthar...
...extit{k}_{-cc}]{\textit{k}_{cc}} C_4\left(^3\Pi_{\it u}, ^3\Pi_{\it g}\right)^*$     (7)


$\displaystyle \rm C_4\left(^3\Pi^-_{\it u}\right)^* \rightarrow{\textit{k}_{rad,el}} C_4 \left(^3\Sigma^-_{\it g}, ^3\Pi_{\it g} \right) + h\nu,$     (8)

where ass and diss designate association and dissociation, rad, v designates vibrational emission, cc designates curve crossing, rad,el designates electronic emission, and k designates a rate coefficient. Assuming that curve crossing is more rapid than the assorted emission processes, we obtain the following expression for the rate coefficient of the overall radiative association ( $k_{\rm RA}$):

\begin{displaymath}k_{\rm RA} = k_{\rm ass} [k_{\rm rad,v} + k_{\rm rad,el} K_{\rm cc}] / (k_{\rm rad,v} + k_{\rm diss}),
\end{displaymath} (9)

where the equilibrium coefficient $K_{\rm cc}$ is equal to $k_{\rm cc}/ k_{\rm -cc}$ and pertains to the $^3\Pi_u$ only .

We first determine the microcanonical rate coefficient for dissociation, $k_{\rm diss}$, from the ground electronic state of linear C4 as a function of internal energy E above the C + C3 dissociation energy. We employ the so-called RRKM (Rice-Ramsperger et al.) unimolecular approach (Holbrook et al. 1996) using the FALLOFF system (Forst 1993). A so-called loose transition state is used with the canonical flexible transition state theory (CFTS) approach (Holbrook et al. 1996). The needed structural and thermochemical parameters of the studied system derive from quantum calculations. Thermal averaging of the microcanonical rate coefficient for dissociation of the C4 adduct obtained from the RRKM calculation allows us to calculate the thermal dissociation rate coefficient. The RRKM thermally averaged lifetime of the C4 adduct, $\rm\tau = \frac{1}{\textit{k}_{diss}(\textit{T})}$, is presented in Fig. 2 as a function of T. From the equilibrium constant governing formation and dissociation of the adduct, the rate coefficient for formation, which differs from  $k_{\rm ass}$ in that fine-structure correlations are not yet taken into account, can be obtained.

 \begin{figure}
\par\includegraphics[width=9cm]{0967fig2.eps}
\end{figure} Figure 2:

RRKM lifetime (in s) of the C4 adduct as a function of the internal energy (in K) above the C + C3 dissociation limit.

Open with DEXTER

The calculated capture rate constant obtained by this procedure at 298 K is equal to $2.0 \times 10^{-10}$ cm3 s-1 with very little temperature dependence in the 10-300 K range, in good agreement with theoretical treatments for reactions between atomic carbon and assorted alkenes and alkynes (Clary et al. 1994) and with the experimental results of Chastaing et al. (2001) for such reactions. This thermal rate coefficient, however, does not take into account the fraction of fine-structure states leading to the attractive ground state C4 surface. With this factor and the assumption that the fine-structure states are in thermal equilibrium, we obtain a refined association rate coefficient of $f(T) \times 2 \times 10^{-10}$ cm3 s-1 in the range 10-1000 K, where

\begin{displaymath}
f(T) = \left(1+2\rm {\rm e}^{-16.4/T}\right)/\left(1 + 3\times {\rm e}^{-16.4/T} + 5\times {\rm e}^{-43.4/T}\right).
\end{displaymath} (10)

The method of estimating $k_{\rm rad,v}$ equates this rate coefficient with the overall vibrational emission rate (Smith 1989; Herbst 1982):

\begin{displaymath}A(E_{\rm vib}) = \sum_i \sum_{n_i} P_{n_i}^{(i)} (E_{\rm vib})A_{n_i \rightarrow n_i -1}^{(i)}
\end{displaymath} (11)

where $A_{n_i \rightarrow n_i -1}^{(i)}$, the Einstein A coefficient for mode i spontaneously emitting from level ni to ni - 1, is equal to $n_i \times A_{1 \rightarrow 0}^{(i)}$), $P_{n_i}^{(i)} (E_{\rm vib})$ is the probability of finding mode i in level ni, and the sums are over all modes i and all energetically accessible quanta ni. The vibrational frequencies and Einstein coefficients, calculated at the B3LYP/vtz level, are listed in Table 4, leading to $k_{\rm rad,v} = 349$ s-1.

Table 5:   Vibrational properties of the ${\rm X~^3 \Sigma}^-_ g$ state of the C4 molecule.

To estimate the rate coefficient $k_{\rm rad,el}$ for stabilization by emission of electronic radiation, we use the results of ab initio calculations that the electronic potential curves of the C4($^3\Pi_u$), C4( $^3 \Sigma^-_g$) and C4($^3 \Pi _g$) states have very similar shapes. Therefore, we can approximately ignore nondiagonal Franck-Condon factors. With the additional assumption that the diagonal factors are unity, we can write that $\langle \nu'\vert\nu''\rangle=\delta_{\nu',\nu''}$, so that the spontaneous electronic emission rate coefficients are given by the expression (esu-cgs units):

\begin{displaymath}k_{\rm rad,el} = \left(64\pi^4{\rm e}^2a_0^2/3h\right)\nu^3\left\vert\mu^{\rm el}\right\vert^2,
\end{displaymath} (12)

where e is the electronic charge, $\nu$ is the transition frequency, a0 is the Bohr radius, and $\vert\mu^{\rm el}\vert$ is the electronic transition moment. From the $^3\Pi_u$ state, there are allowed transitions toward the $^3 \Sigma^-_g$ and the $^3 \Pi _g$ state (Larsson 1983). We obtain that $k_{\rm rad,el} \approx 6\times 10^4$ s-1. Finally, we obtain the equilibrium coefficient $K_{\rm cc}$ by assuming that the relative populations of the $^3\Pi_u$, $^3 \Pi _g$, and $^3 \Sigma^-_g$ states are equal to their relative densities of vibrational states at an internal C4 energy of 5.21 eV, the experimental exothermicity of reaction (4). This assumption leads to an equilibrium coefficient $K_{\rm cc}$ for the $^3\Pi_u$ state only of 0.15, and an overall radiative stabilization rate $k_{\rm rad,v} + k_{\rm rad,el}K_{\rm cc}$ of 9350 s-1.

Once the partial rate coefficients are all calculated, we can obtain the total rate coefficient for radiative association, $k_{\rm RA}$ as a function of temperature. The results are shown in Fig. 3; the value at 15 K is $3~\times~10^{-12}$ cm3 s-1. The value at 15 K if only vibrational stabilization is allowed is a much lower $1 \times 10^{-13}$ cm3 s-1, which serves as a reasonable lower limit. The value obtained with the $\langle \nu'\vert\nu''\rangle=\delta_{\nu',\nu''}$ assumption is certainly an extreme upper limit, because the diagonal Franck-Condon overlaps are certainly lower than unity, the electronic transition moment, assumed to be independent of energy, may well decrease with increasing energy for certain states (Boyé et al. 2002), and the internal conversion from the ground ${\rm X~^3 \Sigma}^-_ g$ of C4 to excited electronic states may not reach equilibrium. Estimation of the size these factors have on the association rate coefficient is not facile, but a reasonable guess is that they contribute to produce a factor of $\approx$0.25, leading to a reduced value at 15 K of $8 \times 10^{-13}$ cm3 s-1. We adopt as an error limit an uncertainty of $\pm$75%, allthough this error does not quite cover the extreme maximum and minimum values discussed above. From the RRKM calculations, we obtain a temperature dependence for $k_{\rm RA}$ of approximately T-1, in excellent agreement with the approach of Bates & Herbst (1988) for a system with two rotational degrees of freedom. We obtain finally a temperature-dependent association rate coefficient $k_{\rm RA}(T)$ of $(4.0 \pm 3.0) \times 10^{-14} \times (T/300)^{-1}$ cm3 s-1 in the range 10-300 K. Compared with the radiative association rate coefficient of $1 \times 10^{-10}$ cm3 s-1 in Smith et al. (2004) and osu.2003, our new value is lower by a factor of $\approx$100 at 10 K. This discrepancy between old and new values is larger than the formerly assumed uncertainty of an order of magnitude.

 \begin{figure}
\par\includegraphics[width=9cm]{0967fig3.eps}
\end{figure} Figure 3:

Radiative association rate coefficient for C4 production as a function of temperature.

Open with DEXTER

3.2 C + C 5

For the reaction between C in its ground $^3\rm P$ state and linear C5 in its ground $^1\Sigma^+_g$ state, the formation of the linear C6 adduct in its $^3 \Sigma^-_g$ ground state is strongly exothermic and spin-allowed. This case is then very similar to the C + C3 one, and to determine if the electronic curves are attractive along the entrance channel reaction pathway, we performed similar ab initio calculations. Our results are also similar, with no barrier for the C6( $^3 \Sigma^-_g$) electronic state but high barriers for the C6($^3\Pi$) electronic states. Unlike the case of the C4 adduct, however, the C6 adduct can fragment by an exothermic route; in particular

                   $\displaystyle {\rm C\left(^3P\right)}$ + $\displaystyle {\rm C_5\left(^1\Sigma +_{\it g}\right) \longrightarrow C_6\left(^3\Sigma -_{\it g}\right)}$  
  $\textstyle \longrightarrow$ $\displaystyle {\rm C_6\left(^1\Sigma\right) \longrightarrow C_3\left(^1\Sigma +_{\it g}\right) + C_3\left(^1\Sigma +_{\it g}\right)},$ (13)

with an overall reaction enthalpy of $\Delta H = -1.36$ eV. Although this reaction is formally spin-forbidden to the products in their singlet ground states, the adduct is likely to last long enough to allow intersystem crossing, as in the C + C2H2 reaction (Mebel et al. 2007), and the barrier on the singlet exit channel is certainly much lower than the exothermicity (Nelson et al. 1982). The first excited triplet state of C3 lies about 2 eV above the ground singlet state so is probably not accessible. Based on experiments, the assumption that the rate coefficient is independent of the number of hydrogen atoms on the C5 reactant, and using the fine structure factor in Eq. (10) , we estimate the rate coefficient to form the adduct and subsequent products to be $1.2 \pm(0.8) \times 10^{-10}$ cm3 s-1 for 10-300 K. This value is somewhat less than that obtained from the capture methodology of Clary et al. (1994). Compared with the value in osu.2003 ( $k_{\rm RA} = 1.4 \times 10^{-10}$ cm3 s-1), the rate coefficient for the reaction $\rm C + C_5$ is approximately the same but the products are C3 + C3 instead of C6.

The rate coefficients for other C + Cn reactions are discussed in the Appendix. These reactions are not among those the model results are very sensitive to.

4 Impact on predicted abundances

4.1 General impact

 \begin{figure}
\par\includegraphics[width=8cm]{0967fig4.eps}
\end{figure} Figure 4:

Average fractional difference (AFD, see text) as a function of time. The dashed lines are for case (1), the dashed-dotted lines for case (2), and the solid lines for case (3). Black lines are for case (a) whereas grey curves are for case (b).

Open with DEXTER

To quantify the overall impact of change in the rate coefficients/products computed in the previous section, we sum Rij(t), the fractional difference for species j and reaction i as defined in Eq. (1), over the species j and divide by the number of species $n_{\rm j}$ to obtain the average fractional difference AFD(t) for reaction i:

\begin{displaymath}AFD(t) = \frac{1}{n_j} \sum_j \left\vert R^i_j(t)\right\vert.
\end{displaymath} (14)

Here, as in Eq. (1), we define the reference fractional abundance $X^{\rm ref}_j(t)$ for each species as that obtained with the old rate coefficient/products, but define the new values of Xij(t) as those obtained with the modified values for reactions i obtained here. An unimportant reaction at a specific time would give an AFD close to zero at that time. It should be noted that although we are using a method similar to Sect. 2.2, the goal here is completely different. In the sensitivity analysis Sect. 2.2, we multiply all reactions, one after another by a factor of 2 in order to compare the abundance sensitivity to individual reactions. Here, we change the rate coefficients for which new calculations have been done, for instance the C + C3 reaction is changed by two orders of magnitude, and we look how much the new rate coefficients have changed the model predictions. If we redo the sensitivity analysis as in Sect. 2.2 with the updated network, we may find different ``important'' reactions since the network has been changed.

For the dense cloud model discussed in Sect. 2, we consider three different cases for reactions i: (1) we change only the rate coefficient of the C + C3 reaction; (2) we change only the products and the rate coefficient of the C + C5 reaction; and (3) we perform the changes for both reactions. For each of these cases, we compute the AFD parameter by (a) summing over all species (452), and by (b) summing only over those species with a fractional abundance greater than 10-10 at 105 yr. For (b), the sum is made over 79, 80 and 77 species for cases (1), (2) and (3) respectively. Figure 4 shows the results for the six different cases. It is clear from the results that the most important effect is obtained by the combined effect of changing both C + C3 and C + C5 reactions. For case (3), the sum over all species yields an AFD of 15 at slightly more than 105 yr, while the sum over the most abundant species only yields an AFD of 13 at 105 yr. If only one of the two reactions is changed, C + C3 has the greater effect whereas C + C5 affects the species at early times.

4.2 Impact on specific species

 \begin{figure}
\par\includegraphics[width=7.5cm,clip]{0967fig5.eps}
\end{figure} Figure 5:

Abundances with respect to total hydrogen of several important species. Black dash-dotted lines (0) are calculated abundances when the standard network osu.2003 is used; grey dashed lines (1) when only C + C3 is modified; grey dash-dotted lines (2) when only C + C5 is modified, and black solid lines (3) when both C + C3 and C + C5 are modified.

Open with DEXTER

In Fig. 5, we show how the abundances of selected species are affected by the modifications of the C + C3 and C + C5 reactions at intermediate times; at other times there is little to no effect. For carbon monoxide, the reduction of the C + C3 reaction rate coefficient results in the lowering of its abundance in the period 103-105 yr by at most one order of magnitude. Changing the products of the C + C5 reaction on the contrary slightly increases the CO abundance. For neutral atomic carbon, the gradual decrease in abundance with increasing time occurs much later when the rate coefficient of C + C3 is reduced; at its maximum, there is a difference in abundance of about 3 orders of magnitude. The change in C + C5 speeds up the destruction of C. In fact, the only species shown in Fig. 5 for which the modification of this reaction has a strong effect are C5 and C6H.

As discussed in Sect. 2.3, with the use of the osu.2003 network, these two reactions are both involved in catalytic cycles to convert C and O into CO. By radiative association reactions with atomic carbon, C3 and C5 form C4 and C6, which then react with atomic oxygen to give CO and regenerate C3 and C5. In this work, we have shown that the C + C3 reaction is less efficient than the value used in osu.2003 by two orders of magnitude at 10 K and the more likely products of the C + C5 reaction are two C3 molecules. Even after changing these reactions, this catalytic process remains the most efficient one to form CO at 104 yr. Changing only the C + C5 reaction products makes this process more efficient by introducing more C3 into the cycle. For this reason, CO is slightly increased whereas C is decreased. Obviously C3 is more abundant because C + C3 consumes less C3 and C + C5 produces it. Contrary to C3, the C5 abundance decreases if the $\rm C + C_5 \rightarrow C_6$ - $\rm O + C_6 \rightarrow C_5 + CO$ cycle is broken because C5 is not locked in this loop anymore and is destroyed to give C3. The effect of these changes on C-rich molecules such as the ones shown in Fig. 5 is less evident to explain. From a general point of view, larger abundances of atomic oxygen and carbon between 104 and 105 yr destroy these species more efficiently than producing them. The radical C6H for instance is destroyed by C to give C7 and C7 in turn reacts with C to produce C3 and C5 (see Table A.1). Similarly C4H2 is destroyed by the reaction $\rm C_4H_2 + C \rightarrow C_5H + H$ and C5H is destroyed by O, producing CO and C4H. Large cyanopolyynes, such as HC7N, are reduced in abundance. On the other hand, large O-bearing molecules, such as CH3OH, are increased in abundance or not changed.

4.3 Comparison with observations

 \begin{figure}
\par\includegraphics[width=7cm,clip]{0967fig6.eps}
\end{figure} Figure 6:

Comparison between model and observations in two different clouds (black: L134N and grey: TMC-1). The solid line represents the comparison with the model using the original values of the C + Cn rate coefficients whereas the dashed lines represents the comparisons with the updated model.

Open with DEXTER

To study the impact of these new rate coefficients on the agreement between models and observations, we used the method described in Wakelam et al. (2006a). We modeled the chemical evolution of two different clouds: L134N and TMC-1 (Cyanopolyyne Peak), using the Nahoon chemical model. The list of reactions is osu.2003 with the original C + Cn rate coefficients and osu.2003 with the new values of the C + Cn rate coefficients listed in Table A.1. For the uncertainties in the rate coefficients of these reactions, we assumed a flat distribution (because the mean values of the rate coefficients are not necessarily the favorite values) with the uncertainty ranges described in Table A.1. For the other reactions, we used log-normal distributions (see Wakelam et al. 2004). The physical conditions are taken to be homogeneous in the observed regions: a gas temperature of 10 K, an H2 density of 104 cm-3 and a visual extinction of 10. The elemental abundances are the low-metal elemental abundances (see Table 1 in Wakelam & Herbst 2008). For TMC-1, we decreased the oxygen abundance to have a C/O elemental ratio of 1.2 in order to better reproduce the observations (see Wakelam et al. 2006a, and references therein). The observed abundances are given in Tables 3 and 4 of Garrod et al. (2007). In addition, we included the observational upper limits on the SiO abundances in these two clouds by Ziurys et al. (1989); the upper limits for this species are $2.6\times 10^{-12}$ and $2.4\times 10^{-12}$ for L134N and TMC-1 respectively. In total, we have observed abundances of 42 and 53 species for these two clouds to compare with our calculations.

The results of the comparisons are shown in Fig 6 in three panels, for each of which a different criterion is utilized. These criteria are based on the Monte Carlo approach to determination of the uncertainties in the computed abundances based on uncertainties in the rate coefficients, and the overlap of Gaussian distributions representing both the computed and observed abundances (Wakelam et al. 2006a). In the top panel, we use the fraction of reproduced molecules. In the middle panel, we use the so-called ``distance of disagreement'' D in which, for species not reproduced by the first criterion, we sum over the remaining species the distances between logarithms of the observed and model fractional abundances. Finally, in the bottom panel, we use the minimum distance of disagreement (over all times), divide by the distance of disagreement, and multiply by the fraction of reproduced molecules. See Wakelam et al. (2006a) for details on these parameters. The time of best agreement is roughly the same for all three parameters, although the sharpest agreement occurs for the bottom panel. It can be seen in all panels that the general agreement with the observations is slightly better using the new values of the C + Cn rate coefficients. However, the more interesting effect is the shift in time of the best agreement. With the original model, the best agreement was obtained around $3 \times 10^4$ yr whereas now we calculate the best agreement to be between 105 and $3\times 10^5$ yr. Note that to compare our results with the results of Wakelam et al. (2006a), we have only updated the C + Cn rate coefficients whereas many other reactions have been updated in the latest version of the osu database (osu.3.2008).

5 Conclusion

This paper presents an approach to improve chemical models, which consists of first identifying important reactions, then obtaining improved rate coefficients for these reactions, and finally studying the impact of these new rate coefficients on the predictions of the model. In this example of the method, we found that under cold dense cloud conditions, the abundances of many species are sensitive to the rate coefficients of the C + C3 and C + C5 reactions. Using detailed ab initio and kinetic calculations, we have shown that the rate coefficient of the $\rm C + C_3 \rightarrow C_4$ previously used in the osu.2003 database is incorrect by two to three orders of magnitude at 10 K and has an inverse temperature dependence. We also found that the rate coefficient for the $\rm C + C_5$ reaction was roughly correct in the osu.2003 database, but the products are most likely $\rm C_3 + C_3$ instead of C6, as assumed earlier based on faulty thermodynamic information. Newly obtained rate coefficients for other reactions in the C + Cn series are also reported, although the model shows little sensitivity to these. With the new rate coefficient for the association reaction between C and C3 and the slightly changed rate coefficient and new products for the $\rm C + C_5$ reaction, we found that the abundances of many carbon chain species are strongly modified at the so-called early time of 105 yr. The main reason for this effect is that the $\rm C + C_3 \rightarrow C_4$ / $\rm O + C_4 \rightarrow CO + C_3$ and $\rm C + C_5 \rightarrow C_6$ / $\rm O + C_6 \rightarrow CO + C_5$ reactions comprise a very efficient catalytic cycle to transform C and O into CO. Slowing down the first cycle and removing the second one, both decreases the abundance of CO and enhances the abundances of atomic C and O at early times. These changes in turn cause an alteration in the early-time abundances of many other species. The agreement with observed abundances in the two dense clouds L134N and TMC-1 (CP) is not significantly modified by these new rate coefficients/products but the best agreement is obtained at later times using the new values. The newly calculated rate coefficients and/or products will be included in the latest network from the Ohio State group (http://www.physics.ohio-state.edu/~eric/research.html).

Acknowledgements
V. W. acknowledges the french national program PCMI for partial support to this work. E. H. thanks the National Science Foundation for its support of his research program in astrochemistry.

Appendix A: Other C + C n reactions

The calculated rate coefficients and products for C + Cn are listed in Table A.1 for n = 2-8 along with estimated uncertainties and the older values used in the osu.2003 network. The applicable temperature range is generally 10-300 K, although for radiative association reactions, the dependence can extend to higher temperatures. For the C($\rm ^3P$) + C2( ${\rm ^1\Sigma }^+_g$) reaction, the spin-allowed pathway leads to a C3( $\rm {a~^3\Pi} _u$) adduct, which is not the ground electronic state of C3 but is highly thermodynamically accessible, lying 5.16 eV below the reactants. Preliminary ab initio calculations (CASCCF, vdz basis, full valence active space) give no barrier for the production of the $\rm {a^3\Pi}_u$ state. To determine the rate coefficient for radiative association, we considered only vibrational relaxation, Computing vibrational frequencies for C3 at the B3LYP/vtz level, we obtained a relaxation rate of $k_{\rm rad,\nu} = 220$ s-1. We then estimated the association (capture) rate coefficient to be $2 \times 10^{-10}$ cm3 s-1 based on experiments of analogous systems (e.g. C + C2H2) and calculated the dissociation rate coefficient by the RRKM method. We obtained the radiative association rate coefficient, shown in Table A.1, to be $(3 \pm 2) \times 10^{-16} \times (T/300)^{-1.0}$ cm3 s-1, a rather low value due to the rapid dissociation of the small C3 adduct.

Table A.1:   Rate coefficients for the C + Cn reactions.

For the C($\rm ^3P$) + C4( $\rm ^3\Sigma^-_g$), C6( $\rm ^3\Sigma^-_g$), C8( $\rm ^3\Sigma^-_g$) reactions, the formation of the Cn+1 adducts is strongly exothermic and spin-allowed. The possible potential energy surfaces in the entrance channels for all of these reactions are $^{1,3,5} \Sigma^+$ and $^{1,3,5} \Pi$. Preliminary ab initio calculations (CASCCF, vdz basis) show that the three $^{1,3,5} \Sigma^+$ states have no barrier in the entrance channel, but that the three $^{1,3,5} \Pi$ states have high barriers along the reaction pathway. Exothermic fragment channels exist for all of these adducts. The rate coefficients for these reactions are once again estimated from experimental values of analogous systems, which in general are smaller than the capture values, derived from the formalism of Clary et al. (1994), and further reduced by a factor that allows for the multiplicity of the potential energy surfaces that correlate with the reagents and also by a factor for the multiplicity of the exit channels if relevant. Correlation involving the fine-structure states of atomic C is ignored. The assigned rate coefficients and estimated uncertainties are shown in Table A.1.

The C($\rm ^3P$) + C7( ${\rm ^1\Sigma }^+_g$) reaction leading to the formation of the C8 adduct in the $\rm ^3\Sigma^-_g$ ground electronic state is stronglyexothermic and spin-allowed. By comparison with C + C5 we assume that there is no barrier in the entrance channel for the formation of the ground C8 electronic state and, although this reaction is formally spin-forbidden to the products in their singlet ground states, the adduct is once again likely to last long enough to allow intersystem crossing. Once again, we use analogous experimental results rather than the formalism of Clary et al. (1994) to obtain the capture rate, and reduce this via Eq. (10).

References

All Tables

Table 1:   Most important reactions for selected species at 105 yr.

Table 2:   Reactions influencing the most species at 105 yr using R.

Table 3:   Reactions influencing the most species at 105 yr using P.

Table 4:   Vdz/CASSCF optimized structural parameters.

Table 5:   Vibrational properties of the ${\rm X~^3 \Sigma}^-_ g$ state of the C4 molecule.

Table A.1:   Rate coefficients for the C + Cn reactions.

All Figures

  \begin{figure}
\par\includegraphics[width=8cm,clip]{0967fig1.eps}
\end{figure} Figure 1:

Co-linear VDZ-CASSCF one-dimensional cuts of the lowest triplet electronic states of linear C4 molecules along the C3-C stretch. The C3 bonds are optimized. The labeling of the C4 electronic states is based on D $_{\rm\infty h}$ symmetry.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{0967fig2.eps}
\end{figure} Figure 2:

RRKM lifetime (in s) of the C4 adduct as a function of the internal energy (in K) above the C + C3 dissociation limit.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=9cm]{0967fig3.eps}
\end{figure} Figure 3:

Radiative association rate coefficient for C4 production as a function of temperature.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=8cm]{0967fig4.eps}
\end{figure} Figure 4:

Average fractional difference (AFD, see text) as a function of time. The dashed lines are for case (1), the dashed-dotted lines for case (2), and the solid lines for case (3). Black lines are for case (a) whereas grey curves are for case (b).

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7.5cm,clip]{0967fig5.eps}
\end{figure} Figure 5:

Abundances with respect to total hydrogen of several important species. Black dash-dotted lines (0) are calculated abundances when the standard network osu.2003 is used; grey dashed lines (1) when only C + C3 is modified; grey dash-dotted lines (2) when only C + C5 is modified, and black solid lines (3) when both C + C3 and C + C5 are modified.

Open with DEXTER
In the text

  \begin{figure}
\par\includegraphics[width=7cm,clip]{0967fig6.eps}
\end{figure} Figure 6:

Comparison between model and observations in two different clouds (black: L134N and grey: TMC-1). The solid line represents the comparison with the model using the original values of the C + Cn rate coefficients whereas the dashed lines represents the comparisons with the updated model.

Open with DEXTER
In the text


Copyright ESO 2009

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.