Theoretical Study of the Reaction of Carbon Monoxide with Oxygen

ARTICLE pubs.acs.org/JPCA

Theoretical Study of the Reaction of Carbon Monoxide with Oxygen Molecules in the Ground Triplet and Singlet Delta States Alexander Sharipov and Alexander Starik* Central Institute of Aviation Motors, Scientific Educational Center “Physical-Chemical Kinetics and Combustion”, Moscow, Russia ABSTRACT: Quantum chemical calculations were carried out to study the reaction of carbon monoxide with molecular oxygen in the ground triplet and singlet delta states. Transition states and intermediates that connect the reactants with products of the reaction on the triplet and singlet potential energy surfaces were identified on the base of coupled-cluster method. The values of energy barriers were refined by using compound techniques such as CBS-Q, CBS-QB3, and G3. The calculations showed that there exists an intersection of triplet and singlet potential energy surfaces. This fact leads to the appearance of two channels for the triplet COþO2(X3Σg-) reaction, which produces atomic oxygen in the ground O(3P) and excited O(1D) states. The appropriate rate constants of all reaction paths were estimated on the base of nonvariational transitionstate theory. It was found that the singlet reaction rate constant is much greater than the triplet one and that the reaction channel COþO2(a1Δg) should be taken into consideration to interpret the experimental data on the oxidation of CO by molecular oxygen.

1. INTRODUCTION Kinetics of CO oxidation is in the basis of the reaction mechanism of ignition and combustion of various hydrocarbons and is involved as a principal part in the mechanism of the reaction of syngas with air.1-5 The process COþO2fproducts is believed to be one of the most important reactions of chain initiation in this case. Kinetic data on the reaction of CO with triplet ground state oxygen O2(X3 Σg-) has been studied by different research groups and are well documented nowadays.6-17 However, novel challenges that address the reactions of different species with excited singlet oxygen molecules O2(a1Δg) and O2(b1Σgþ) come into play due to extensive studies of the processes in the atmospheres of different planets,18 in electric discharges,19 and in combustible systems activated by laser radiation20,21 or by nonequilibrium plasma.22,23 This dictates a necessity to analyze kinetics of the reaction COþO2fproducts more accurately taking into account the reaction channels both with normal ground state triplet oxygen and with excited singlet oxygen molecules. It is worth noting that the reaction of CO with singlet delta oxygen O2(a1Δg) is likely to be an initiation step in CO and syngas combustion assisted by electric discharge or resonance laser radiation and, possibly, is also important under slow CO oxidation in the middle Earth’s atmosphere. Unfortunately, until now there is no any analysis that concerns such a reactive system. Therefore, the aim of the present paper is to provide the detailed kinetic data on the reaction of CO oxidation when both ground triplet state oxygen and singlet excited delta oxygen participate in the reaction. 2. METHODOLOGY The strategy consisted of series of ab initio quantum chemical calculations of CO-O2(X3Σg-) and CO-O2(a1Δg) systems. r 2011 American Chemical Society

These calculations were used to examine the potential energy surface (PES) and kinetics of the reaction of CO with molecular oxygen. The geometry of reactants, transition states, and possible products of reactions was optimized by using an unrestricted coupled-cluster method with single and double excitations (UCCSD)24 and Dunning’s correlation consistent cc-pVDZ basis set.25 To define the energy values of critical points at PESs more accurately, four modern computational models were used: (i) UCCSD(T)/cc-pVTZ, (ii) CBS-Q,26 (iii) CBS-QB3,27 and (iv) G3.28 A frequency analysis using the CCSD/cc-pVDZ level of theory was performed on each stationary point to characterize it either as a minimum (all vibrational frequencies are real) or as a first-order saddle point (there exists a single imaginary frequency). The fact that reactants and products are connected by the transition state was confirmed by following minimum energy paths of the reactions using the Gonzalez-Schlegel intrinsic reaction coordinate (IRC) algorithm29 in both directions at CCSD/cc-pVDZ level. Electronic structure calculations were conducted by using the PC GAMESS version30 of the GAMESS (US) QC package31 and GAUSSIAN 09 package.32 For the each reaction channel the appropriate temperaturedependent rate constant was estimated with the usage of conventional nonvariational transition-state theory (CTST) formalism including Wigner’s tunneling correction Γ(T)33
 kb T Q TS Ea kðTÞ ¼ ΓðTÞg exp h QA QB kb T

ð1Þ

Received: October 29, 2010 Revised: January 31, 2011 Published: February 22, 2011 1795

dx.doi.org/10.1021/jp110345s | J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A ΓðTÞ ¼ 1 þ


 1 hν 2 24 kb T

ARTICLE

ð2Þ

Here QTS and QA, QB are the total partition functions for the transition state and the reagents denoted as A and B, respectively (calculated by using the harmonic oscillator and the rigid rotator approximation), g the reaction path degeneracy, Ea the activation energy which includes the zero-point energy correction to the electronic energy, kb is the Boltzmann constant, h is the Planck constant, and ν* is the imaginary frequency of the activated complex at the top of the barrier.

3. PESS AND CO þ O2 REACTION PATHWAYS The triplet and singlet potential energy surfaces for the COO2 system were investigated. In PES calculations, 6 internal coordinates were introduced: interatomic distances r1, r2, r3, bond-angles R and β, as well as the torsion angle γ (see Figure 1). The results of the calculations of reaction pathways that connect the reactants and products of the reaction COþO2(X3Σg-) are schematically depicted in Figure 2. Herein and hereafter the following designations for critical points at PESs are used: R for reactants, SP for saddle points, P for products, and IM for local minima; the left superscripts denote the multiplicity. The appropriate energy values of critical points estimated at the different levels of theory are summarized in Table 1. The geometrical parameters of critical points optimized at the CCSD/cc-pVDZ, CCSD(T)/cc-pVDZ, and CCSD(T)/ccpVTZ levels of theory are presented in Table 2. Rotational temperatures and vibrational frequencies of saddle points and intermediates are summarized in Table 3. It was revealed that the reaction of CO with O2(X3Σg-) can occur through two different channels with planar transition states. Shown in the diagrams of parts a and b of Figure 2 are

Figure 1. Internal coordinates used in PES calculations.

the saddle points and intermediates of planar structures (Cs symmetry) for which the torsion angle γ is equal to 0 and 180
, respectively. In fact, the saddle points, which connect reactants and products, have two isomeric forms: cis (γ = 0
) and trans (γ = 180
). In each case the pathway connecting the reactants and products passes through the cascade of two successive saddle points (3SP1cis and 3SP2cis; 3SP1trans and 3SP2trans) separated by the shallow local minima (3IMcis and 3IMtrans). It should be emphasized that these minima exist on the surfaces of pure electronic energy only, and they vanish if we take into account the energy of zero-point motion (see Figure 2). Thereby, only 3SP2cis and 3SP2trans saddle points are really significant for the reactivity of the system under study. To visualize this topic, the relaxed potential scans in the saddle points vicinity were performed. In these calculations, the r2 bond length and β angle were varied, and for each fixed pair of r2 bond length and β angle the partial optimization of remaining coordinates were conducted. The resulting relaxed PESs calculated at the CCSD/cc-pVDZ level of theory are shown in Figure 3. As was mentioned above, 3SP2cis and 3SP2trans structures are characterized mainly by different values of torsion angle γ (see Table 2). The torsion potential for 3SP structure is presented in Figure 4. This potential was obtained for 3SP2cis structure with all parameters except the γ angle fixed. The energy of 3SP2cis saddle point is less than the energy of 3SP2trans by the value of ∼6000 K. There are no saddle points of the first order except the stated above. So, one can conclude, that the pathway running across the 3 SP2cis structure is the only preferable reaction pathway, and saddle point 3SP2cis is, indeed, a transition state for the reaction CO þ O2 ðX3 Σg - Þ ¼ CO2 þ Oð3 PÞ

ðR1Þ

The results of the calculations of pathways that connect the reactants and products of the reaction CO þ O2(a1Δg) are schematically depicted in Figure 5. In this case, the reaction pathway connecting reactants and products has 2-fold degeneracy, because there exist two different minimum energy paths with identical energy barriers. This fact is visualized in Figure 6, where the torsion potential for 1SP structure is presented. One can see that there are two optical isomers of 1SP structure separated by the low inversion barrier of ∼200 K height. To gain more insight in the geometry of 1SP structure the relaxed potential scans in the singlet saddle point vicinity were performed in the same manner as for the triplet saddle points. The resulting relaxed PES performed at the CCSD/cc-pVDZ level of theory assuming γ = 80
is shown in Figure 7.

Figure 2. Relative energies of reactants, saddle points, intermediates and products on the triplet PES at the CCSD/cc-pVDZ level of theory: closed symbols are the pure electronic energy; open symbols denote the energy values with ZPE correction. All values are given with respect to the electronic energy of CO þ O2(X3Σg-) system. 1796

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A

ARTICLE

Table 1. Energy Values of Critical Points Obtained at Different Levels of Theory CCSD/cc-pVDZ

CCSD(T)/cc-pVDZ

-263.0201632

-263.0410061

-263.2848388

3

-262.9287724

-262.9523165

-263.2023282

-262.929795 -262.9252937

-262.9525115 -262.9489697

-263.2028451 -263.19 76865

SPtrans1 3 trans

-262.9304047

-262.957244

-263.2035362

-262.9305276

-262.9518753

-263.2036173

3

-262.903097

-262.9274488

-263.1750774

-263.0388292

-263.0580963

-263.3011 844

SPcis1

3

cis

IM SPcis2

3 3

IM

SPtrans2

3

P

-262.9669383

-262.9923476

-263.2371865

1

-262.9028435

-262.9490193

-263.1986631

1

-262.9489259

-262.9749124

-263.2198591

1

R SP P

R

β

γ

1.315

126.59

110.79

0.0

1.353

129.85

111.83

0.0

1.311

1.510

138.83

113.01

0.0

1.166

1.497

1.312

121.72

112.42

180.0

3

1.172

1.443

1.328

122.92

111.38

180.0

3

1.188

1.285

1.597

137.32

105.95

180.0

1

1.155

1.514

1.296

133.78

113.86

76.86

1

1.155

1.514

1.296

133.74

113.91

283.14

CCSD(T)/cc-pVTZ

3

R

Table 2. Structural Parameters of Critical Points (Ångstroms, Degrees) Optimized at Different Levels of Theory r1, Å

r2, Å

r3, Å

3

1.163

1.526

3

1.174

1.409

3

1.178

3

CCSD/cc-pVDZ SPcis1 IMcis SPcis2 SPtrans1 IMtrans SPtrans2 SP SP

CCSD(T)/cc-pVDZ 3

SP

Now let us conduct an analysis of the mutual arrangement of triplet and singlet PESs corresponding to the CO þ O2(X3Σg-) and CO þ O2(a1Δg) systems. For this purpose IRC calculations were performed on the both triplet and singlet PESs. To study possible intersystem crossing (IC), we projected the triplet minimal energy path obtained in the course of IRC calculations on the singlet PES and vice versa. The resulting profiles are presented in Figure 8 and 9. From Figure 8a one can deduce that at γ = 0
the triplet and singlet PESs do not come close to one another. This leads to the conclusion that in this case the nonadiabatic transitions can not play a significant role. The calculation of PES sections at γ = 180
(see Figure 8b) revealed the triplet-singlet IC. Taking into account the triplet and singlet surfaces intersection we are led to infer that 3SP2trans structure is a transition state for two reaction channels ( 3 Oð PÞ 3 ðR2, R2'Þ CO þ O2 ðX Σg Þ f Oð1 DÞ These reactions produce oxygen atoms in the ground O(3P) (R2) and excited O(1D) (R20 ) electronic states. Because of significantly higher energy barrier than that of the R1 process, the reaction paths R2 and R20 can be neglected during evaluation of the rate constant of the overall reaction CO þ O2(X3Σg-)fproducts. Therefore, the determination of the ratio of the rate constants of R2 and R20 reactions k(R20 )/k(R2) does not provide any significant benefit from the practical point of view, and hereafter we will designate this reaction path as overall R2 reaction channel. Nevertheless, the fact that in the course of the reaction of CO with oxygen molecule in the ground electronic state O2(X3Σg-) the excited O(1D) atom forms, is quite interesting from the theoretical viewpoint, and it would be desirable to receive an experimental evidence of this fact. As was stated previously, the 1

CO þ O2 ða ΔgÞ ¼ CO2 þ Oð DÞ 1

ðR3Þ

reaction has 2-fold degenerated MEP running across the optical isomers of 1SP structure, and IC occurs in this case too (see Figure 9). This means that COþO2(a1Δg) system after passing through 1SP point has a probability to leave the path leading to CO2þO(1D) products. This fact should be taken into account when estimating the rate constant of the R3 process.

cis

1.173

1.497

1.334

127.38

111.06

0.0

3

1.777

1.430

1.357

129.32

111.63

0.0

3

1.184

1.321

1.506

138.13

112.90

0.0

3

1.168

1.668

1.336

123.52

111.36

180.0

3

1.181

1.405

1.330

125.12

108.93

180.0

3

1.194

1.291

1.602

136.80

105.56

180.0

1

SP 1 SP

1.178 1.178

1.413 1.413

1.377 1.377

131.20 131.11

106.44 106.41

93.51 266.49

3

1.162

1.513

1.317

126.80

110.81

0.0

3

1.171

1.416

1.348

129.59

111.65

0.0

3

1.176

1.302

1.507

139.96

113.19

0.0

3

1.165

1.492

1.312

121.52

112.89

180.0

3

1.170

1.443

1.327

122.63

111.86

180.0

3

1.186

1.273

1.583

138.88

106.12

180.0

1

1.677 1.677

1.434 1.434

1.388 1.388

128.35 128.32

100.39 100.37

91.10 268.90

1

IMcis SPcis2 SPtrans1 IMtrans SPtrans2

CCSD(T)/cc-pVTZ SPcis1 IMcis SPcis2 SPtrans1 IMtrans SPtrans2

SP 1 SP

Table 3. Rotational Temperatures (K) and Vibrational Frequencies (cm-1) of the Structures Optimized at the CCSD/ cc-pVDZ Level of Theory θrot1

θrot2

θrot3

f1

f2

f3

f4

f5

f6

1.345 1.342

0.325 0.340

0.262 0.271

i674 302

212 317

310 639

718 803

1070 970

1985 1954

3

1.290

0.317

0.254

i940

288

450

777

1005

1986

3

3.842

0.235

0.222

i439

164

376

652

1130

1984

3

3.887

0.241

0.227

183

321

467

764

1107

1966

3

2.969

0.232

0.216

i1780

326

332

741

943

1964

1

1.956

0.261

0.238

i668

82

425

635

1168

1998

1

1.958

0.260

0.238

i668

81

426

635

1170

1997

3

cis

SP 1 IMcis

3

SPcis2 SPtrans1 IMtrans SPtrans2 SP SP

4. REFINEMENT OF ACTIVATION BARRIERS As possible reaction pathways in the CO þ O2 system were established, it is necessary to determine their activation barriers more accurately than it was done while PES investigation. Table 4 summarizes the values of activation barriers obtained at different levels of theory for the processes R1-R3. The first line values were obtained by means of the coupled-cluster method with single and double excitations (CCSD/cc-pVDZ); in the second line the CCSD values were complemented with the noniterative 1797

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A

ARTICLE

Figure 3. The contour plots of PES sections obtained through the relaxed potential energy scan in the vicinity of triplet saddle points of cis (a) and trans (b) configurations. The values of energy are given in electronvolts with respect to the electronic energy of CO þ O2(X3Σg-) system.

triple corrections (CCSD(T)/cc-pVDZ) during reoptimization, and in the third line, the energy values obtained by using the CCSD(T) method in extended basis set (cc-pVTZ) are given. In addition, three compound methods involving a sequence of predefined ab initio molecular orbital calculations developed to achieve high target accuracy ((2 kcal/mol) were used: CBS-Q complete basis set method of Petersson et al.,26 its version CBSQB3,27 and G3.28 Furthermore, the adiabatic (with reoptimization of internuclear distance) energies (Te) of O2(a1Δg) and O(1D) species as well as the heat of reaction R1, estimated by means of different methods and borrowed from Gurvich et al.,34 are given here. One can conclude that all methods involved reproduce the heat of reaction R1 calculated by using the data from Gurvich et al.34 with an accuracy better than 2 kcal/mol. The electronic energies of O2(a1Δg) and O(1D) species (Te(O2) and Te(O)) are overestimated by the methods in use, but compound methods CBS-Q, CBS-QB3, and G3 reproduce the values of Te(O2) and Te(O) with reasonable accuracy. It is worth noting that during unrestricted coupled-cluster calculations for singlet critical points, pure singlet states were obtained. For triplet critical points, except 3SP2, the mild spin contamination was observed, i.e., spin-squared operator ÆS2æ was only slightly greater than those of the pure triplet states (ÆS2æ 2.0 e 0.03). For 3SP2 cis and trans species, the greater values of spin contamination were calculated (ÆS2æ -2.0 ≈ 0.2). As a matter of fact, this may result in the slight overprediction of the computed total energy of 3SP2 structures and, therefore, leads to the overestimation of Ea for triplet channels during CC calculations. However, as it was shown recently,35,36 highly correlated methods such as UCCSD(T) are rather insensitive to the spin contamination. In addition, CBS-Q and CBS-QB3 methods,27 in contrast to G3, include the empirical corrections for the spin contamination for open-shell species, but nevertheless, all these three compound techniques give similar values of Ea for triplet channels (see Table 4). Thereby, one can conclude that the spin contamination does not affect the computed values of activation barrier on the triplet PES seriously. Furthermore, the coupled-cluster (CCSD/cc-pVTZ) calculations were examined by the T1 diagnostic37,38 to find out whether the single determinant Hartree-Fock reference wave function describes the system satisfactorily. Hence, the larger the T1 value is the less reliable the results of the single-reference coupled cluster wave function are. In line with our calculations, the values of Euclidean norm of the T1 vector for reactants and

Figure 4. Torsion potential for the 3SP2 structure at the CCSD/ cc-pVDZ level of theory.

Figure 5. Relative energies of reactants, saddle points, intermediates and products on the singlet PES at the CCSD/cc-pVDZ level of theory: closed symbols are the pure electronic energy; open symbols denote the energy values with ZPE correction. All values are given with respect to the electronic energy of CO þ O2(X3Σg-) system.

products are clearly under the limit of reliability (T1 < 0.02) proposed by Lee and Taylor.38 For saddle points and intermediates, larger values of T1 were observed. Thus, for critical points of the triplet cis reaction path, T1 values are within the range 0.030.04. For critical points on the triplet trans and singlet reaction paths, even greater values of T1 were obtained (0.03-0.05). This indicates that the system has a non-negligible multireference character, and the application of multireference electron correlation procedures may be useful for the system under study. However, the multireference study is beyond the scope of the present work. Note that in the work of Rienstra-Kiracofe et al.,39 1798

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A when the reaction of C2H5 with O2 was studied by CCSD(T) method, the values of 0.04-0.05 for T1 were accepted as reliable.

5. RATE CONSTANTS FOR CO þ O2 REACTION PATHWAYS The pre-exponential factors A(T) in the rate constant expressions for the R1, R2, and R3 processes were derived in line with the relationships 1 and 2. The appropriate partition functions

Figure 6. Torsion potential for 1TS at the CCSD/cc-pVDZ level of theory.

Figure 7. The contour plot of PES section obtained through the relaxed potential energy scan in the vicinity of singlet saddle point. The values of energy are given in electron-volts with respect to the electronic energy of CO þ O2(a1Δg) system.

ARTICLE

were calculated by using the harmonic oscillator and rigid rotator approximations, and the moments of inertia and normal-mode frequencies were taken from the results of CCSD/ cc-pVDZ calculations. The values of path degeneracy g for R1 and R2 reactions were taken equal to unity, whereas for the R3 reaction it was taken equal to 2. It should be noted, that significant errors can appear if the harmonic oscillator approximation is used to calculate the partition function for the low frequency mode ( f2 = 82 cm-1) of 1SP structure relevant to the torsion motion. This degree of freedom of the 1SP structure should be treated as a hindered internal rotator. As can be clearly seen from Figure 6, where the appropriate torsion potential is shown, the torsion barrier V0 is approximately 1500 K. In the present work we used the approximation by Truhlar,40 which gives a smooth approximation from free rotator to harmonic oscillator pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3Þ Q hin ¼ Q ho tanhð πV0 =kT Þ where Qho is the partition function when the normal mode under consideration is treated as a harmonic oscillator. The comparison of different ways to treat the low frequency mode with f2 = 82 cm-1 is given in Figure 10. One can see that treating this mode as a harmonic oscillator overpredicts the partition function of transition state by 1.3 times at T = 5000 K. As was mentioned above, the reactive system after passing through the transition state has a certain probability to leave the singlet term via the singlet-triplet intersection point. In accordance with Landau-Zener model,41 the transition probability from the singlet state to the triplet state via the intersection of two PESs can be expressed in the following form rffiffiffiffiffi! 2 2πV12 μ p ¼ 1 - exp ð4Þ pjΔFj 2E where |ΔF| is the magnitude of the gradient difference of two diabatic PESs (V1 and V2) at the crossing point (|ΔF| = (∂V1/ ∂qr - ∂V2/∂qr)1/2), V12 the spin-orbit coupling matrix element between two interacted electronic states, μ the reduced mass of the bond (mode) whose direction is defined as the reaction coordinate, and E the energy of the nuclei motion along the reaction coordinate qr at the IC point. A comprehensive study of spin-orbit coupling effects is beyond the scope of the present work. As usual, the value of V12 for similar molecular system does

Figure 8. The electronic energy profile as a function of the reaction coordinate on the triplet PES (solid line) and the energy profile resulting from the projection of triplet reaction path on the singlet PES for γ = 0 and 180
(dashed line). All calculations are performed at the CCSD/cc-pVDZ level of theory. 1799

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A

ARTICLE

not exceed 0.01 eV.42,43 We can apply this value of V12 as an upper estimation of spin-orbit coupling matrix element. The value of the energy of nuclei motion E = 0.14 eV was applied assuming that the system passes through the saddle point with nearly zero velocity. As is seen from Figure 9, the transition from the saddle point to the crossing point is accompanied by the energy release into the nuclei motion along the reaction coordinate. Its value is depicted in Figure 9 as ΔE. These assumptions lead us to the value of transition probability p ≈ 0.07. The corresponding temperature dependencies of preexponential factors were calculated in line with the eqs 1-4. However, at low temperature they exhibit a nonmonotonic temperature dependence mainly due to tunneling effect. By assumprion of a negligible rate of the reaction under study at T < 1000 K, we approximated these dependencies for the reactions R1, R2, and R3 in the temperature range T = 8005000 K related to experimental conditions in the following form AðTÞ ¼ 7:625  106 T 1:670 ðfor the R1 reactionÞ AðTÞ ¼ 5:176  106 T 1:728 ðfor the R2 reactionÞ and AðTÞ ¼ ð1 - pÞ  7:279  107 T 1:600 ðfor the R3 reactionÞ

It is worth noting that these approximations are valid only at T g 800 K and that the usage of them at lower temperature may cause a significant underestimation of rate constants.

6. COMPARISON WITH EXPERIMENT Though the reaction CO þ O2 = CO2 þ O has been studied for a long time in a wide temperature range, the temperature dependence of its rate constant has not been well established until now.44 The measurements of rate constants of the CO þ O2 = CO2 þ O reaction conducted by different researchers are related to the temperature range T = 1000-3500 K.6,7,9-12,15-17 Furthermore, the data relevant to the high temperature range T = 2600-4600 K obtained from the measurements of the rate constant for the backward reaction CO2 þ O = CO þ O28,13,14 on the base of detailed balancing principle are presented too. When analyzing the experimental data concerning the reaction of CO with O2, we must keep in mind that the reacting mixture contains oxygen molecules in different electronic states, and reactions of CO with O2 molecule both in the ground triplet electronic state O2(X3Σg-) and in the singlet excited states O2(a1Δg) and O2(b1Σgþ) contribute to the total rate of reaction. Therefore, for the adequate description of the experiments, the rather complicated kinetic models involving electronically excited species should be applied. Assuming that the excitation processes of the O2(b1Σgþ) and O2(a1Δg) states due to electronic translational exchange O2 ða1 Δg Þ þ M ¼ O2 ðb1 Σg þ Þ þ M and pooling process O2 ðX3 Σg - Þ þ O2 ðb1 Σg þ Þ ¼ O2 ða1 Δg Þ þ O2 ða1 Δg Þ occur rather rapidly,45,46 we can consider the distribution of oxygen molecules over the electronic states as equilibrium. Therefore, in our calculations we need to take into consideration the contributions of both channels R1 and R3 assuming rapid quenching of excited atomic oxygen O(1D) that forms in the course of R3 reaction. In this case, the rate constant of overall reaction CO þ O2 = CO2 þ O can be expressed as follows 1 ktot ¼ kðR1Þ  γX ðX3 Σg Þ þ kðR3Þ  γa ða Δg Þ

Figure 9. The electronic energy profile as a function of the reaction coordinate on the singlet PES (solid line) and the energy profile resulting from the projection of singlet reaction path on the triplet PES (dashed line). All calculations performed at the CCSD/cc-pVDZ level of theory.


 Ei gi exp T 1 γi ði ¼ X3 Σg , a Δg Þ ¼ Qel

ð5Þ

ð6Þ

where gi is the degeneracy of the ith state, Ei the energy of the ith

Table 4. Values of Activation Barriers (Es) for the Processes R1-R3 as Well as the Energies of O2(a1Δg) (Te(O2)) and O(1D) (Te(O)) and Heat of Reaction R1, Obtained at Various Levels of Theory and Borrowed from Gurvich et al.34 Ea, K triplet (cis)

triplet (trans)

singlet

Te(O2), K

Te(O), K

ΔHf(1), K

CCSD/cc-pVDZ

30410

37280

20600

16760

28390

4970

CCSD(T)/cc-pVDZa

28650

36360

14040

15350

26265

4280

CCSD(T)/cc-pVTZa

27980

34970

12530

15000

25680

4240

CBS-Q

26910

34020

12570

14820

24790

4470

CBS-QB3

26950

33470

13660

14420

24400

4240

G3 Gurvich et al.34

26420

33850

11550

14860 11390

23870 22820

4960 4070

a

a

The ZPE correction obtained at the CCSD/cc-pVDZ level was incorporated into all coupled cluster energy values. 1800

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A

Figure 10. Temperature dependence of the partition function for the torsion mode of 1SP structure obtained on the base of harmonic oscillator and free rotator assumptions as well of as a smooth approximation from free rotator to harmonic oscillator calculated in line with Truhlar.40

state, and Qel the electronic partition function of oxygen. Figure 11 compares the effective rate constant for the CO þ O2 = CO2 þ O reaction predicted by using the eqs 5 and 6 with the state-resolved rate constants for R1 and R3 processes with Ea values obtained from CCSD(T)/cc-pVTZ calculations. As is seen, the rate constant estimated in line with eqs 5 and 6 equations is much greater than that when O2 molecule is in the X3Σg- state only (by a factor of 5-20). The importance of this effect is not surprising because the excited a1Δg state is a lowlying (Te = 0.98 eV) one, and already at T = 3000 K the mole fraction of singlet delta oxygen molecule in the total amount of oxygen achieves 1.5%. Therefore, this methodology should be applied to the analysis of available experimental data, when the characteristic times of the excitation of a1Δg and b1Σgþ states of O2 molecule are smaller than the characteristic time of the oxidation process. Let us analyze the experimental data on the rate constant for the reaction CO þ O2 = CO2 þ O. These data are presented in Figure 12. One can see that at T > 2500 K the measurements of different works (except Clark et al.8) are close to each other. However, at the temperatures below 2500 K more than 10-fold discrepancy is observed. It should be noted that all measurements were conducted behind the incident or reflected shock wave by using the shock-tube technique. The experimentalists are inclined to interpret this discrepancy in terms of impurity effects mainly due to the presence of H2O in the shock-tube in trace amounts.15,17 However, our results show that the other reason of this discrepancy may exist. It can be associated with the features of nonequilibrium processes that occur behind the shock front during experiment. In fact, the finite rate of equilibration between O2(X3Σg-) and O2(a1Δg) states in the postshock region can violate the validity of assumptions 5 and 6 owing to slow excitation of O2(a1Δg) state via electronic translational relaxation. This leads to the different contribution of the R3 reaction path to the overall rate of reaction CO þ O2 = CO2 þ O depending on the parameters behind the shock front and, consequently, different concentration of O2(a1Δg) molecule in the postshock region. Shown in Figure 12 are also the approximations of rate constant for the overall CO þ O2 = CO2 þ O process obtained applying different calculated values of Ea (see Table 4). One can conclude that when CCSD(T)/cc-pVDZ (curve 1), CBS-QB3 (curve 2), and CCSD(T)/cc-pVTZ (curve 3) values of Ea are

ARTICLE

Figure 11. Temperature dependence of the rate constant for the reaction of CO with O2(X3Σg-) (closed squares) and O2(a1Δg) (open squares) for the activation barriers calculated at the CCSD(T)/ cc-pVTZ level of theory. The solid curve is the overall rate constant for Boltzmann distribution of O2 molecules over the electronic states.

Figure 12. Temperature dependence of the rate constant for overall CO þ O2 = CO2 þ O reaction measured by different researchers (symbols)6-17 and estimated in line with expressions 5 and 6 assuming different calculated Ea values: (1) CCSD(T)/cc-pVDZ, (2) CBS-QB3, (3) CCSD(T)/cc-pVTZ, and (4) G3. The approximation of Tsang and Hampson48 is given by dotted curve.

taken, our predictions coincide rather well with measurements almost in the whole temperature range. Only at high temperature (T > 3500 K) the computed curves have tendency to deviate from most experiments. The curve 4 overestimates the experimental data at more broad temperature range (T > 2000 K) noticeably, but the discrepancy diminishes with the temperature decrease. Thus, we attempted to assess whether such disagreement was related to the fact that the nonvariational transition state theory formalism had been utilized instead of variational one. Therefore, in the present work for the R1 reaction path (see Figure 8a) we performed a normal-mode analysis along the reaction path at the CCSD/cc-pVDZ level of theory. The reaction rate constants estimated by canonical variational theory47 (CVT) kCVT(T) were compared with the rate constants predicted by CTST with the same activation barriers kCTST(T), and the ratio r(T) was determined rðTÞ ¼

kCVT ðTÞ kCTST ðTÞ

In line with our calculations, at T = 300-5000 K the value of r(T) is in the range 1.00-0.99. Therefore, one can conclude that 1801

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A application of variational theory can not appreciably improve the pre-exponential factors calculated above. For comparison, the approximation of Tsang and Hampson48 applied in the modern kinetic mechanisms for syngas-air mixture combustion2,3 is also presented in Figure 12. It is worth noting that this approximation was proposed for the temperature range T = 300-2500 K. One can see from plots in Figure 12 that the usage of this approximation at T > 3000 K leads to the noticeable underestimation of both of measurements13,14 and of theoretical predictions of the present work. Because the CBS-QB3 method reproduces the excitation energies of O2(a1Δg) and O(1D) and heat of reaction slightly better than the CCSD(T)/cc-pVTZ one (see Table 4), just CBSQB3 values of Ea may be recommended as the most reliable for the rate constant evaluation for the R1-R3 processes. On the basis of CBS-QB3 state-weighted calculations, the following approximation for the CO þ O2 = CO2 þ O reaction rate constant applicable at T = 800-5000 K can be proposed ktot ðTÞ ¼ 4:323  107 T 1:618 expð- 25018=TÞ It is necessary to emphasize that the assumptions 5 and 6 are not always valid when the postshock wave conditions considered, but the contribution of singlet oxygen must be taken into account, if we inquire about the rate constant for thermally equilibrium chemistry. According to our calculations, the expressions for R1, R2 R20 , and R3 reactions can be written in the form kðTÞ ¼ 7:625  106 T 1:670 expð- 26950=TÞ kðTÞ ¼ 5:176  106 T 1:728 expð- 33470=TÞ and kðTÞ ¼ ð1 - pÞ  7:279  107 T 1:600 expð- 13660=TÞ respectively, where the value of p was estimated as ∼0.07.

6. CONCLUSIONS Ab initio quantum chemistry calculations were conducted to examine the potential energy surfaces and to investigate the kinetics of CO þ O2(X3Σg-) and CO þ O2(a1Δg) systems. During the triplet and singlet PESs investigations, different reaction channels representing the planar and spatial types of the interaction of CO and O2 molecules were observed. The critical points of triplet and singlet PESs for these types of interaction were determined that allows us to calculate the energetic barriers for each reaction path configuration. The calculations of PES sections along the reaction pathways both for triplet CO þ O2(X3Σg-) and for singlet CO þ O2(a1Δg) systems revealed the intersections of triplet and singlet PESs. The appropriate rate constants for CO þ O2(X3Σg-) and CO þ O2(a1Δg) reactions were determined taking into account intersystem crossing. The existence of such an intersection leads to the appearance of the channel for triplet CO þ O2(X3Σg-) reaction which produces excited atomic oxygen O(1D). The rate constants for all reaction paths were calculated on the base of nonvariational transition-state theory. It was found that the singlet reaction rate constant is much greater than the triplet one. At T = 1000 K, the ratio of these rate constants is as large as 5500. This ratio decreases noticeably with increasing temperature. The calculated values of overall rate constant for the reaction CO þ O2 = CO2 þ O involving the triplet and singlet

ARTICLE

channels are in reasonable agreement with experimental data. It was shown that the contribution of the singlet reaction path CO þ O2(a1Δg) to the overall reaction of CO with oxygen molecule must be taken into account to interpret the experimental data on the rate of CO þ O2 = CO2 þ O2 reaction.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work was supported by Russian Foundation for Basic Research (Grants 10-08-90035 and 11-08-90443) and by the Federal Purpose Program “Scientific and science-educational manpower of innovational Russia” (State Contract No. 2.740.11.0074). ’ REFERENCES (1) Chaos, M.; Dryer, F. L. Combust. Sci. Technol. 2008, 180, 1053. (2) Li, J.; Zhao, Z.; Kazakov, A.; Chaos, M.; Dryer, F. L.; Scire, J. J. Int. J. Chem. Kinet. 2007, 39, 109. (3) Sun, H.; Yang, S. I.; Jomaas, G.; Law, C. K. Proc. Combust. Instit. 2007, 31, 439. (4) Simmie, J. M. Prog. Energy Combust. Sci. 2003, 29, 599. (5) Leung, K. M.; Lindstedt, R. P. Combust. Flame 1995, 102, 129. (6) Sulzmann, K. G. P.; Myers, B. F.; Bartle, E. R. J. Chem. Phys. 1965, 42, 3969. (7) Drummond, L. J. Aust. J. Chem. 1968, 21, 2631. (8) Clark, T. C.; Garnett, S. H.; Kistiakowsky, G. B. J. Chem. Phys. 1969, 51, 2885. (9) Brabbs, T. A.; Belles, F. E.; Brokaw, R. S. Symp. Int. Combust. Proc. 1971, 13, 129. (10) Dean, A. M.; Kistiakowsky, G. B. J. Chem. Phys. 1971, 54, 1718. (11) Gardiner, W. C.; McFarland, M.; Morinaga, K.; Takeyama, T.; Walker, B. F. J. Phys. Chem. 1971, 75, 1504. (12) Rawlins, W. T.; Gardiner, W. C. J. Phys. Chem. 1974, 7, 497. (13) Baber, S. C.; Dean, A. M. J. Chem. Phys. 1974, 60, 307. (14) Korovkina, T. D. High Energy Chem. 1976, 10, 75–76. (15) Thielen, K.; Roth, P. Ber.-Bunsen. Phys. Chem. 1983, 87, 920. (16) Koike, T. Bull. Chem. Soc. Jpn. 1991, 64, 1726. (17) Sutherland, J. W.; Patterson, P. M.; Klemm, R. B. Tech. report BNL-46985; Brookhaven National Laboratory: Upton, NY, 1992. (18) Capitelli, M.; Ferreira, C. M.; Gordietz, B. F.; Osipov, A. I. Plasma kinetics in atmospheric gases, Springer Series on Atomic, Optical, and Plasma Physics; Springer-Verlag: Berlin, 2000; p 31. (19) Ionin, A. A.; Kochetov, I. V.; Napartovich, A. P.; Yuryshev, N. N. J. Phys. D: Appl. Phys. 2007, 40, R25. (20) Starik, A. M.; Kuleshov, P. S.; Titova, N. S. J. Phys. D.: Appl. Phys. 2009, 42, 175503. (21) Starik, A. M.; Titova, N. S.; Bezgin, L. V.; Kopchenov, V. I. Combust. Flame 2009, 156, 1641. (22) Starik, A. M.; Lukhovitskii, B. I.; Titova, N. S. Combust. Explos. Shock Waves 2008, 44, 249. (23) Starik, A. M.; Loukhovitski, B. I.; Sharipov, A. S.; Titova, N. S. Combust. Theory Model. 2010, 14, 653. (24) Bartlett, R. J. J. Phys. Chem. 1989, 93, 1697. (25) Kendall, R. A.; Dunning, T. H., Jr.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (26) Ochterski, J. W.; Petersson, G. A.; Montgomery, J. A., Jr. J. Chem. Phys. 1996, 104, 2598. (27) Montgomery, J. A., Jr.; Frisch, M. J.; Ochterski, J. W.; Petersson, G. A. J. Chem. Phys. 1999, 110, 2822. (28) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V; Pople, J. A. J. Chem. Phys. 1998, 109, 7764. 1802

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803

The Journal of Physical Chemistry A

ARTICLE

(29) Gonzalez, C.; Schlegel, B. H. J. Chem. Phys. 1989, 90, 2154. (30) Granovsky, A. A. PC GAMESS version 7.1.5, http://classic. chem.msu.su/gran/gamess/index.html. (31) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (32) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, Jr., J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision A.1; Gaussian, Inc., Wallingford CT, 2009. (33) The theory of rate processes; Glasstone, S., Laidler, K. J., Eyring, H., Eds.; McGraw-Hill: New York, 1941. (34) Thermodynamics Properties of Individual Substances; Gurvich, L. V., Veyts, I. V., Alcock, C. B.; Hemisphere Pub. Co., 1989. (35) Szalay, P. G.; Vazquez, J.; Simmons, C.; Stanton, J. F. J. Chem. Phys. 2004, 121, 1. (36) Krylov, A. I. J. Chem. Phys. 2000, 113, 6052. (37) Lee, T. J.; Rice, J. E.; Scuseria, G. E.; Schaefer, H. F., III. Theor. Chim. Acta. 1989, 75, 81. (38) Lee, T. J.; Taylor, P. R. Int. J. Quantum Chem., Symp. 1989, S23, 199. (39) Rienstra-Kiracofe, J. C.; Allen, W. D.; Schaefer, H. F., III. J. Phys. Chem. A 2000, 104, 9823. (40) Truhlar, D. G. J. Comput. Chem. 1991, 12, 266. (41) Quantum Mechanics: Non-Relativistic Theory; Landau, L. D., Lifshitz, E. M., Eds.; Pergamon Press: Oxford, 1977. (42) Chernyi, G. G.; Losev, S. A.; Macheret, S. O.; Potapkin, B. V. Physical and chemical processes in gas dynamics: cross sections and rate constants. Progress in Astronautics and Aeronautics AIAA; Reston, V.A, 2002, p 196. (43) Fedorov, D. G.; Koseki, S.; Schmidt, M. W.; Gordon, M. S. Int. Rev. Phys. Chem. 2003, 22, 551. (44) NIST Chemical Kinetics Database, NIST Standard Reference Database 17, Version 7.0 (Web Version), Release 1.4.3, Data version 2008.12, National Institute of Standards and Technology, Gaithersburg, Maryland, 20899-8320. Web address:http://kinetics.nist.gov/ (45) Herron, J. T.; Green, D. S. Plasma Chem. Plasma Process. 2001, 21, 459. (46) Schurath, U. J. Photochem. 1975, 4, 215. (47) Garrett, B. C.; Truhlar, D. G. J. Phys. Chem. 1979, 83, 1052. (48) Tsang, W.; Hampson, R. F. J. Phys. Chem. Ref. Data. 1986, 15, 1087.

1803

dx.doi.org/10.1021/jp110345s |J. Phys. Chem. A 2011, 115, 1795–1803