Theoretical Studies of the O (3P)+ C2 Reaction at Hyperthermal

Nov 18, 2012 - The O + C2 reaction has been investigated with the quasiclassical trajectory (QCT) method in conjunction with direct dynamics electroni...
2 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCC

Theoretical Studies of the O(3P) + C2 Reaction at Hyperthermal Energies Mausumi Ray, Biswajit Saha, and George C. Schatz* Department of Chemistry, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208, United States S Supporting Information *

ABSTRACT: The O + C2 reaction has been investigated with the quasiclassical trajectory (QCT) method in conjunction with direct dynamics electronic structure calculations using density functional theory (DFT) forces. Trajectory surfacehopping calculations have also been performed to study spin-forbidden reactions. Calculations were performed at collision energies of 1−5 eV so as to simulate conditions relevant to erosion of carbon-based materials on spacecraft in low Earth orbit (LEO). Since the energy difference between the electronic ground state (X1Σ+g ) and the first excited triplet state (a3Πu) of the C2 molecule is only 2.1 kcal/ mol, two reactions, O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu), have been studied. We present here the detailed mechanism, electronic branching, product energy disposal, and angular distribution for these reactions. The calculations show that the O(3P) + C2(a3Πu) reaction can occur on singlet, triplet, and quintet surfaces to give the spin-allowed electronically excited CO(1Σ) + C(1D), CO(3Π) + C(3P), and CO(3Π) + C(1D) products as well as the ground state product CO(1Σ) + C(3P), with CO(3Π) + C(3P) being the most important, while O(3P) + C2(X1Σ+g ) reacts on triplet surfaces to give primarily the CO(1Σ) + C(3P) product with only minor branching to spin-forbidden excited states. Reactions at 1 eV energy proceed on all surfaces through formation of the collision complex CCO, while the collision complex only forms briefly at 5 eV. The CO + C cross section for O(3P) reacting with C2(a3Πu) is three times smaller than with C2(X1Σ+g ). Angular distributions show that the product CO + C is more and more backward scattered as collision energy is increased as can be explained in terms of collision lifetime shortening at higher energies. Product energy disposal shows that for O(3P) + C2(X1Σ+g ) about 50% of the total available energy is deposited in relative translation, 10% is in CO rotation, and 40% is in CO vibration. For O(3P) + C2(a3Πu), about 50% of the available energy ends up as electronic excitation. The partitioning for each electronic state in the C2(a3Πu) reaction is strongly state dependent, but for CO(3Π) + C(3P) product vibrational excitation accounts for 60% of the energy in excess of electronic energy. and other combustion flames13 and during arc vaporization of carbon. The reactions of the dicarbon molecule in its (X1Σ+g ) electronic state are also relevant to formation of polycyclic aromatic hydrocarbons (PAH) and carbon-rich nanostructures such as fullerenes.14 As a result of the high heat of formation of C2, O + C2 and many other reactions with C2 are highly exothermic, leading to electronically excited products.15 Detailed kinetics information concerning reaction of C2 with combustion-relevant species such as atomic and molecular oxygen are therefore of general interest. Since the energy difference between the electronic ground state C2(X1Σ+g ) and the first excited electronic triplet state (a3Πu) of the C2 radical is only about 2.1 kcal/mol, both electronic states play an important role in C2 chemistry.16,17 Despite several investigations of the O + C2 reaction kinetics, the state of understanding of this system is still quite limited.15,18 Becker and co-workers investigated this reaction at room temperature and found that the reaction was

1. INTRODUCTION Atomic oxygen produced by photodissociation of oxygen molecules by solar radiation is the dominant species at low Earth orbital (LEO) altitudes. These oxygen atoms impact on orbiting spacecraft with ∼4.5 eV energy and an arrival flux approximately 1015 O atoms cm−2 s−1which are sufficient to cause carbon fiber composite materials, commonly used on spacecraft, to be oxidized at rates which can limit the durability of many spacecraft components.1 Thus, hyperthermal collisions between spacecraft and the O(3P) atoms in addition to exposure to vacuum ultraviolet (VUV) light, X-rays, and charged particles produce important chemical and physical changes in the exposed surfaces of spacecraft resulting in erosion and degradation of carbon based materials.1−4 A model for the reaction of hyperthermal oxygen with carbon-based materials is to consider O + carbon clusters. This process is important in other contexts as well, such as in combustion processes5−7 and soot formation.8−10 and the chemical evolution of extraterrestrial environments such as molecular clouds and circumstellar envelopes of dying carbon stars.11 The smallest carbon cluster, dicarbon (C2), in its electronic ground state 1Σ+g has been detected in hydrocarbon12 © 2012 American Chemical Society

Received: July 5, 2012 Revised: November 3, 2012 Published: November 18, 2012 26577

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

calculation. Calculations at 10 different crossing-point geometries give SOC values that fall in the range 0.1−120 cm−1. We used a high-end SOC value of 100 cm−1for the present TSH calculation. However, we also performed calculations for sample trajectories with a SOC value of 70 cm−1 to check the sensitivity of the results to the coupling parameter. 2.B. Direct Dynamics Calculations. The O(3P) + C2 reaction has been investigated with the QCT version of Born− Oppenheimer molecular dynamics in which energies and forces were computed at each step of the trajectory with the DFT(B3LYP) method and a 6-31G(d,p) basis set. For calculating such B3LYP trajectories, we interfaced the QChem program package with a molecular dynamics code utilizing a standard fifth-order predictor, sixth-order corrector integrator to propagate Hamilton’s equations of motion. The reactants C2(X1Σ+g ) and C2(a3Πu) were prepared by choosing the appropriate spin symmetry and determining the equilibrium geometry. Then the rovibrational initial conditions were determined by running an intramolecular trajectory with zero-point energy starting from the corresponding equilibrium geometry. We utilized the dynamic reaction coordinate (DRC) subroutine of GAMESS30 for running this intramolecular trajectory, and the trajectory was integrated for many vibrational periods. The phase of the vibration motion was then sampled from this trajectory, and the attacking oxygen atom was randomly placed around the vibrating C2 molecule with an initial center-of-mass distance of 10a0 and a maximum sampled impact parameter (bmax) of 7a0. Reagent relative translational energies of 1−5 eV were chosen to simulate the LEO conditions. The trajectory was terminated when the product pair was sufficiently separated (12a0) or a maximum number of integration steps of 3000 was reached. The time step chosen for the dynamics was 0.24 fs. Upon completion of the trajectory propagation, an analysis was performed to derive various dynamical properties of the products. The reactive cross section for channel i was computed as

sufficiently exothermic to allow formation of electronically excited products.15 However, the details of the reaction mechanism and branching between possible products are unknown. Here we present a detailed direct dynamics study of the reaction of atomic O with both singlet and triplet C2 at hyperthermal energies between 1 and 5 eV using the quasiclassical trajectory (QCT) method and also including trajectory surface-hopping (TSH) calculations to study intersystem crossing (ISC) effects.19−22 Complete construction of a multidimensional potential energy surface (PES) for use in dynamical simulations is a demanding task, even for a three-atom system, as it requires high-quality electronic structure calculations over a broad range of molecular configurations. In the present study, we used the direct dynamics method wherein energy and forces were computed as the trajectory evolves with energy gradients determined using the density functional theory (DFT) method where the B3LYP functional was used in combination with a 631G(d,p) basis set. Our implementation of the TSH method follows a previous study of O + C2H4 wherein only gradients for the currently occupied surface are computed as the trajectory evolves and hops are only allowed when crossings between the singlet and the triplet surfaces are encountered.23 A decision to hop is then made on the basis of the Landau− Zener transition probability formula.

2. COMPUTATIONAL DETAILS 2.A. Electronic Structure Calculations. The reaction paths for both the O(3P) + C2(X1Σ+g ) and the O(3P) + C2(a3Πu) reactions were evaluated at stationary points along the PES with DFT using the 6-31G(d,p) basis set. Geometry optimizations and energy calculations were carried out with both the B3LYP and the BMK functional, where the latter was used because BMK-calculated barrier heights are often in better agreement with experiment than those from B3LYP.24,25 We ascertained that none of the equilibrium geometries exhibited imaginary frequencies and each transition state possessed only one imaginary frequency. Intrinsic reaction coordinate (IRC) calculations were performed to make sure the transition states connect the expected minima on the potential energy surfaces. Energies were also evaluated using the CCSD(T) method and a cc-pVTZ basis set to calibrate the DFT-calculated energetics at the DFT-optimized geometries. The Q-Chem program package was used for all these electronic structure calculations.26 The magnitude of the spin−orbit coupling (SOC) is small compared to the energy difference between the surfaces for most geometries. Thus, the only significant triplet−singlet (T− S) transitions occur when the surfaces are very close to a crossing. To simplify our trajectory calculations, surface hops were allowed only at points where the triplet and singlet surfaces cross. We did not actually calculate the spin−orbit matrix elements at each crossing point of each trajectory. Instead, the large majority of T−S crossings occurred in the vicinity of the low-energy CCO linear geometry where the SOC matrix elements do not vary much with geometry, so we used a constant value for the SOC valid in this region of the PES throughout our calculations. To verify this, we calculated the spin−orbit matrix elements at selected crossing-point geometries detected in sample trajectories using the Breit− Pauli method as implemented by Fedorov and Gordon27 in GAMESS.28 Orbitals were generated from a state-averaged CASSCF calculation using an aug-cc-pVDZ basis set.29 Eight electrons in eight active orbitals were included in the CASSCF

σi =

2 2πbmax N

N

∑ xkPik

(1)

k=1

where xk = bk/bmax with bk being the impact parameter of trajectory k and bmax being the maximum sampled impact parameter, N is the total number of trajectories, and Pik is 1 if the trajectory yields channel i and 0 otherwise. To study ISC effects in the O(3P) + C2 reaction, we used the direct dynamics quasiclassical variant of the TSH method with B3LYP forces. The trajectories were initiated on the ground triplet electronic state at reagent collision energies of 1−5 eV to simulate LEO collisions. Initially the trajectories used energy and gradients from the triplet surface, but the singlet surface energy was also calculated, and the singlet−triplet gap was evaluated to identify crossings. At intersections of singlet and triplet surfaces the transition probability was computed using the Landau−Zener (LZ) formula PLZ

⎡ 2 2πV13 ⎢ = 1 − exp⎢ − dE 3 dE1 ⎢⎣ ℏŻ dR − dR

⎤ ⎥ ⎥ ⎥⎦

(2)

where V13 is the SOC between the two surfaces, R is the coordinate perpendicular to the crossing seam, Ż is the nuclear coordinate velocity perpendicular to the crossing seam, and |(dE1)/(dR) − (dE3)/(dR)| is the magnitude of the difference 26578

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

changes are presented here at the B3LYP/6-31G(d) level. As shown in Figure 1a, first the O atom attacks one carbon atom of C2 to give a linear intermediate INT1 wherein one C−O bond of length 1.173 Å was formed while the C−C distance was elongated to 1.365 Å. Then, the O atom moved to the other carbon atom of C2, and again a linear C2O intermediate INT2 was formed via a transition state TS with a very high activation barrier of 72.7 kcal/mol at the B3LYP/6-31G(d,p) level. In TS, C−O and C−C distances were elongated to form the intermediate INT2 in the next step with C−O and C−C bond lengths of 1.427 and 1.415 Å, respectively. Finally, products CO and C were formed directly from INT2. The activation barrier calculated using BMK/6-31G(d,p) is 73.7 kcal/mol, which is similar to the B3LYP result, and both the B3LYP and the BMK results show similar trends. CCSD(T) calculations with the cc-pVTZ basis set were done for each DFT(B3LYP)- and DFT(BMK)-optimized stationary point to check the accuracy of the DFT energetics. The results (Figure 1a) show that the CCSD(T) energies calculated for all stationary points are much higher than that obtained with the DFT. However, the CCSD(T) calculation presents the same trend in the potential energy surface; moreover, the CCSD(T)calculated activation barrier of 72.5 kcal/mol is similar to the B3LYP result, indicating that the B3LYP method presents a reliable electronic structure method for the present system. Figure 1a and Table 1 show that the reaction is highly

between singlet and triplet energy gradients at the crossing point. These gradients were determined by interpolating the gradients to the crossing location using values obtained from time steps just before and just after the crossing. Thus, it is only at crossing points that gradient information is required on both surfaces. Once the LZ probability was calculated, it was then compared with a random number ζ between 0 and 1 to decide if a surface hop should occur.23 If PLZ > ζ, the electronic state was switched from triplet to singlet and the trajectory was continued on the singlet surface. No momentum adjustment had to be made on the new surface because hopping only occurred at the point of intersection where E1 = E2. If PLZ < ζ, the trajectory continued on the triplet surface. The treatment of a switch from the singlet to the triplet state followed exactly the steps described above. No attempt was made to restrict the frequency of successive hops.

3. RESULTS AND DISCUSSION 3.A. Potential Energy Surfaces. The reaction between O(3P) and C2(X1Σ+g ) is found to always involve formation of the intermediate C2O(3Π) on its way to give the ground state product CO(X1Σ) + C(3P). All triplet stationary points along the potential energy surface between reactants and products are shown in Figure 1a, with energies given relative to the energy of the reactants. Unless specified, the geometry and energy

Table 1. Reaction Energy (ΔE) for O(3P) + C2(1Σ) and O(3P) + C2(3Π) Reactionsa Calculations reaction O(3P) O(3P) O(3P) O(3P) O(3P)

+ + + + +

C2(1Σ) → CO(1Σ) + C(3P) C2(3Π) → CO(3Π) + C(3P) C2(3Π) → CO(3Σ) + C(3P) C2(3Π) → CO(1Σ) + C(1D) C2(3Π) → CO(3Π) + C(1D)

ΔE, kcal/mol −133.5 +23.6 −110.9 −66.9 +67.6

(−118.6) (+19.2) (−120.9) (−87.3) (+52.9)

ΔE is presented for both DFT(B3LYP) and CCSD(T) (in parentheses). a

exothermic for formation of both the spin-allowed CO(X1Σ) + C(3P) and the spin-forbidden CO(X1Σ) + C(1D) products, although the latter is 44 and 33.6 kcal/mol more endoergic in both B3LYP and CCSD(T) methods, respectively. In the case of the O(3P) + C2(X1Σ+g ) → CO(X1Σ) + C(3P) reaction, the 9-fold electronic degeneracy of the reactants (i.e., 3 P ignoring spin−orbit splitting) leads to three potential surfaces of triplet multiplicity. CASSCF(10,10) calculations with an aug-cc-pVDZ basis set at points along the reaction coordinate of Figure 1a show that all three surfaces have a barrier well below 5 eV, as shown in Figure S1 (Supporting Information), and hence, all of them are reactive at the highest energy we consider. As a result, the electronic degeneracy factor is taken to be 1 for this reaction, as shown in Table 2. We use this electronic degeneracy factor to evaluate the cross-section of the reaction, as described below. The reaction O(3P) + C2(a3Πu) can occur on potential surfaces with singlet, triplet, and quintet spin multiplicities. B3LYP/6-31G(d,p) calculations indicate that reaction for the lowest energy surfaces of each spin multiplicity passes through formation of stable intermediates. The relevant stationary points for the lowest singlet surface associated with O(3P) + C2(a3Πu) are shown in Figure 1b, and those for the lowest triplet and quintet surfaces are shown in Figures S2(a) and

Figure 1. Geometry and potential energy changes for (a) O(3P) + C2(1Σ+g ) and (b) O(3P) + C2(a3Πu) reactions. Energies (kcal/mol) are presented in black, blue, red, and green for calculations using B3LYP, BMK, CCSD(T)/B3LYP, and CCSD(T)/BMK, respectively. Note that the green and red results are nearly identical. Energies are relative to (a) the lowest triplet and (b) the lowest singlet. 26579

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

Table 2. Electronic Degeneracy Factors in the O + C2 Reaction reaction O(3P) O(3P) O(3P) O(3P) O(3P) O(3P) O(3P)

+ + + + + + +

C2(1Σ) → CO(1Σ) + C(3P) C2(3Π) → CO(1Σ) + C(1D) C2(3Π) → CO(3Π) + C(3P) C2(3Π) → CO(1Σ) + C(3P) C2(3Π) → CO(3Π) + C(3P) C2(3Π) → CO(3Π) + C(1D) C2(3Π) → CO(3Π) + C(3P)

potential surface

electronic degeneracy factors

triplet singlet singlet triplet triplet triplet quintet

1 5/54 1/9 1/6 1/3 5/9 5/54

S2(b) (Supporting Information). It is noted that the stationary points for the lowest triplet surface are the same as those for O(3P) + C2(X1Σg+) (Figure 1a) (as further discussed in the next section). As shown in Figure 1b, the linear C2O intermediate INT2 can form via formation of INT1 and TS. In INT1, two C−O bonds of 1.353 Å were formed while the C−C bond was elongated to 1.506 Å. In TS, one C−O distance was elongated to 2.001 Å while the other C−O distance was shortened to 1.322 Å to form only one C−O bond in INT2. The activation barrier for TS is 57.6, 65.8, and 37.7 kcal/mol in B3LYP/6-31G(d,p), BMK/6-31G(d,p), and CCSD(T)/ccpVTZ calculations, respectively. Finally, CO and C are directly formed from the linear intermediate INT2. As shown in Figure 1b and Table 1, the O(3P) + C2(a3Πu) reaction is highly exothermic for formation of the spin-allowed electronically excited CO(X1Σ) + C(1D) products. The other electronically excited spin-allowed products CO(3Π) + C(3P) are 23.6 and 19.2 kcal/mol endoergic according to the B3LYP and CCSD(T) calculations, respectively. Note that the lowest triplet excited state of CO is 3Π, and its vertical excitation energy from the ground singlet state X1Σ+ is 138 kcal/mol.31,32 In the triplet surface O(3P) + C2(a3Πu) reaction, as shown in Table 1 and Figure S2(a), Supporting Information, the ground state product CO(X1Σ) + C(3P) is highly exothermic while the other spin-allowed electronically excited products CO(3Π) + C(3P) and CO(3Π) + C(1D) are endoergic in both DFT and CCSD(T) calculations. In the quintet surface O(3P) + C2(a3Πu) reaction, the spin-allowed electronically excited product CO(3Π) + C(3P) is endoergic in both DFT and CCSD(T) calculations, as shown in Table 1 and Figure S2(b), Supporting Information. For the O(3P) + C2(a3Πu) → CO(1Σ) + CO(1D) reaction, the reactants have 54-fold electronic degeneracy and CASSCF(10,10)/aug-cc-pVDZ calculations at points along the reaction coordinate of Figure 1b show that 6 singlets, 6 triplets, and 1 quintet are reactive as they have barriers below 5 eV, as shown in Figure S3 (Supporting Information). To take this into account, for the singlet surface O(3P) + C2(a3Πu) → CO(1Σ) + CO(1D) reaction, the electronic degeneracy factor is 5/54 as the reactants have 54-fold electronic degeneracy and only 5 singlets can give products having 5-fold electronic degeneracy. In the same way, electronic degeneracy factors for the products of the O(3P) + C2(a3Πu) reaction on the triplet and quintet surface were determined. The degeneracy factors which are explained in detail in Figure S3, Supporting Information, are given in Table 2. Plots comparing the lowest triplet surface for O(3P) + C2(X1Σ+g ) and the lowest singlet, triplet, and quintet surfaces for O(3P) + C2(a3Πu) are shown in Figure 2. Note that singlet C2 (X1Σg) is lower in energy than triplet C2 (a3Πu) by 2.4 kcal/

Figure 2. Triplet vs singlet and quintet potential energy surfaces for the O + C2 reaction. Solid lines represent B3LYP/6-31G(d,p) calculation, while dotted lines represent CCSD(T)/cc-pVTZ calculation. In contrast to Figure 1, a single energy scale is used for multiple potential surface reactions.

mol at the CCSD(T)/cc-pVTZ level, but singlet C2 (X1Σg) is higher in energy than triplet C2 (a3Πu) by 22.6 kcal/mol at the DFT(B3LYP)/6-31G(d,p) level. Figure 2 shows that the triplet PES of O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu) is more stable than the singlet and quintet PES of O(3P) + C2(a3Πu) in both B3LYP and CCSD(T) calculations. B3LYP/6-31G(d,p) calculations show that the activation barrier (73.7 kcal/mol) on the triplet surfaces of the O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu) reactions is considerably higher than that on the singlet (57.6 kcal/mol) and quintet (16.2 kcal/mol) surfaces of the O(3P) + C2(a3Πu) reaction. However, all barriers are well below the available energy from the reactants, so these barriers will not control reactivity. 3.B. Dynamics. In this section, we present the results from direct dynamics QCT calculations for both O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu) reactions with Ecoll = 1−5 eV. For each Ecoll, 1000 trajectories were generated. Trajectories were simulated on the lowest triplet surface in the case of the O(3P) + C2(X1Σ+g ) reaction, while for the O(3P) + C2(a3Πu) reaction trajectories were simulated on the lowest singlet, triplet, and quintet surfaces. Note that the lowest triplet surfaces for O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu) refer to the same electronic state (in density functional theory), but because the geometries sampled by the reactants are distinct, the cross sections, product branching, and energy distributions are different. For the O(3P) + C2(X1Σ+g ) reaction, all reactive trajectories gave only the CO(X1Σ) + C(3P) product, while for the O(3P) + C2(a3Πu) reaction, about 65% of the reactive trajectories on the singlet surface gave the spin-allowed electronically excited product CO(X1Σ) + C(1D) and about 35% gave the other spin-allowed electronically excited product CO(3Π) + C(3P); about 65% of the reactive trajectories on the triplet surface gave the CO(X1Σ) + C(3P)product, 25% gave the CO(3Π) + C(3P) product, and about 10% gave the CO(3Π) + C(1D) product; all reactive trajectories on the quintet surface gave only the CO(3Π) + C(3P) product. All these results predict that substantial spin-allowed electronically excited products can be formed in the O + C2 reaction if a significant amount of triplet C2 is present in the source. 26580

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

Figure 3. Snapshots from representative trajectories showing two different mechanisms in the O(3P) + C2(1Σ+g ) reaction at Ecoll = 1 eV.

C2(a3Πu) reaction than in the O(3P) + C2(X1Σ+g ) reaction at all collision energies. 3.B.I. Excitation Functions. The excitation functions (cross section vs reagent translational energy) for the ground and excited state products are presented in Figure 4. Note that the electronic degeneracy of reactants and products has been included in the cross-section calculation, and the total reaction cross section σi calculated using eq 1 is multiplied by the

Upon inspecting the trajectories in detail we observed that mainly two types of reaction mechanisms are involved in the O(3P) + C2(X1Σ+g ) reaction. Representative snapshots from trajectories with Ecoll = 1 eV for these two mechanisms are shown in Figure 3. In both mechanisms (we assigned them as M1 and M2), the product CO + C is formed through formation of a collision complex. In M1, the O atom attacks one of the C atoms of C2 with a concomitant stretch of the C− C bond (see the 138 fs snapshot). Then the O atom migrates from the initially attacked C atom to the other C atom of the C2 molecule at around 143 fs. Finally, the CO + C product forms through rupture of the C−C bond of a linear CCO collision complex at around 173 fs. The CCO collision complex was found in the electronic structure calculation, as shown in Figure 1. In mechanism M1, the O atom attacks the C atoms sequentially and CO is formed from the second attacked carbon atom. In M2, the O atom attacks one of the C atoms of the C2 molecule and then a linear CCO collision complex is formed that after a few wide oscillations decomposes into C and CO. The carbon atom in the CO product is the one the O attacked initially. Mechanisms similar to M1 and M2 were observed in the trajectories at 5 eV; the only difference observed is that the collision complex CCO was formed very briefly and trajectories were short-lived compared to those at 1 eV collision energies. In the O(3P) + C2(a3Πu) reaction, the CO + C product was formed also via the collision complex CCO at 1 eV, while the collision complex was formed very briefly at 5 eV; therefore, mostly short-lived trajectories were observed. Trajectories are shorter lived in the O(3P) +

Figure 4. Excitation functions for the O(3P) + C2(1Σ) and O(3P) + C2(3Π) reactions. Standard deviation error bars are provided for all results. Excitation function for O(3P) + C2(3Π) → CO + C is the sum of four associated elementary processes O(3P) + C2(3Π) → CO(1Σ) + C(1D), O(3P) + C2(3Π) → CO(3Π) + C(3P), O(3P) + C2(3Π) → CO(3Σ) + C(3P), and O(3P) + C2(3Π) → CO(3Π) + C(1D). 26581

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

electronic degeneracy factor, as shown in Table 2. In the case of the O(3P) + C2(X1Σ+g ) → CO(X1Σ) + C(3P) reaction on the triplet surface, the total reaction cross-section σi calculated using eq 1 is multiplied by 1 as the electronic degeneracy factor is 1 in this case, as described in the potential energy surface section above. In the same way, electronic degeneracy factors were included in cross-section calculations for all products of the O(3P) + C2(a3Πu) reaction. For the O(3P) + C2(a3Πu) reaction, the CO(1Σ) + C(1D), CO(3Π) + C(3P), CO(3Π) + C(1D), and CO(1Σ) + C(3P) products are found from the trajectory calculations on the singlet, triplet, and quintet surfaces. As shown in Figure 4, the cross section for producing CO + C in the O(3P) + C2(a3Πu) reaction is about three times smaller than that for the O(3P) + C2(X1Σ+g ) reaction. The lower reactivity of the quintet surface plays an important role in this difference as this contributes significantly to the spin degeneracy. As collision energy increases, the CO + C cross section decreases rapidly for both reactions. For the O(3P) + C2(a3Πu) reaction, the largest cross section is found for the electronically excited CO(3Π) + C(3P) product and this cross section is about 1.5 times the ground state CO(1Σ) + C(3P) cross section but about twice that for the other electronically excited product CO(1Σ) + C(1D). The cross section for the CO(3Π) + C(1D) product is much smaller than the other three products from the O(3P) + C2(a3Πu) reaction. In the O(3P) + C2(a3Πu) reaction, the cross sections to give CO(1Σ) + C(3P) and CO(3Π) + C(3P) decrease rapidly as collision energy increases from 1 to 3 eV and then decreases slowly from 3 to 5 eV, while the cross section to give the CO(1Σ) + C(1D) and CO(3Π) + C(1D) products does not change much from 1 to 5 eV. From the above discussions, we predict that if the source contains a mixture of C2(X1Σ+g ) and C2(a3Πu) substantial electronically excited products, especially CO(3Π) and C(1D), will be found in the O + C2 reaction. 3.B.II. Product Energy Partitioning. The reaction mechanisms on the surfaces that we studied predominantly involve collision complex formation at a collision energy of 1 eV and direct reaction as energy was increased. This transition in reaction mechanisms affects product energy partitioning considered in this section. The average fraction of total available energy going into relative translation, CO vibration, and CO rotation as a function of the reagent translational energy is shown in Figure 5 for the various products that we considered. We focus on the high collision energy regime as the main interest of this paper lies in characterization of the reaction dynamics at hyperthermal energies due to interest in LEO structural damage. For the O(3P) + C2(X1Σ+g ) reaction, as shown in Figure 5, about 50% of the total available energy is deposited in relative translation, about 10% is in CO rotation, and about 40% is in CO vibration. The fraction in relative translational energy increases slowly with Ecoll, as shown in Figure 5a, while CO rotational energy increases (Figure 5b) and CO vibration decreases (Figure 5c). For the O(3P) + C2(a3Πu) reaction, the dominant CO(3Π) + C(3P) product is produced on singlet, triplet, and quintet surfaces, and we see in Figure 5 that CO vibration receives 60% of the available energy (in excess of electronic energy) at a collision energy of 1 eV. This fraction decreases with increasing collision energy, while relative translation grows, receiving 50% of the available energy at 5 eV. CO rotation receives only 4% of the available energy at all collision energies. The CO(1Σ) +

Figure 5. Energy release in terms of average fraction of the available energy for the C2(1Σ) + O(3P) and C2(3Π) + O(3P) reactions.

C(3P) product is associated with the triplet surface, and about 55% of the total available energy goes to relative translation, about 10% to CO rotation, and about 35% to CO vibration at collision energy of 1 eV. These results vary weakly with Ecoll with trends that are similar for the same product that arises from O(3P) + C2(X1Σ+g ). For the other electronically excited product CO(1Σ) + C(1D), based on the singlet surface, about 60% of the total available energy goes to relative translation, about 10% to CO rotation, and about 30% to CO vibration for 1 eV collision energy, with vibration decreasing and rotation increasing with increasing collision energy, and translation is largely unchanged. 3.B.III. Angular Distributions. An important property for understanding the microscopic reaction mechanism is the product angular distribution. Figure 6 shows angular distributions for both reactions at different reagent translational energies, while angular distributions for various products from different potential surfaces are shown in detail in Figure S4 (Supporting Information). Note that the angle plotted is the angle of the outgoing diatomic with respect to the incoming oxygen atom direction. The figure shows a general tendency for angular distributions that are very broad with close to forward− backward symmetry for all products at some energies; however, there are some deviations worthy of comment. The angular distribution in Figure 6a for the product CO(X1Σ) + C(3P) from the O(3P) + C2(X1Σ+g ) reaction at Ecoll = 1 and 2 eV shows a preference for forward scattering, while it has a backward peak at Ecoll = 3 eV and almost forward and backward 26582

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

b for different reagent translational energies. The plots show that the largest impact parameters that contribute to reaction gradually decrease with increasing energy, thus reducing forward scattering. 3.B.IV. TSH Results. We carried out TSH calculations to study spin-forbidden effects in the O + C2 reaction dynamics for 1−5 eV collision energy. Inspection of the O(3P) + C2(X1Σ+g ) trajectories shows that a small number of trajectories starting on the triplet surface hop once to a singlet surface, leading to CO(X1Σ) + C(1D). Figure 7 shows the change in

Figure 7. Triplet and singlet potential energies as a function of time for geometries along a reactive trajectory for the C2(X1Σ) + O(3P) reaction from TSH calculations. Singlet−triplet crossing point where the reactive trajectory hops from triplet to singlet surface is shown in inset.

triplet and singlet potential energy with time for one such reactive trajectory at Ecoll = 1 eV. The singlet surface was initially well above the triplet, then it intersected with the triplet surface at around 97.8 fs, at which point hopping to the singlet surface occurred. The inset of Figure 7 shows an expanded view of the triplet−singlet crossing event. Then the reactive trajectory spent the rest of its time on the singlet surface until reaction was completed. As shown in Figure 8, the number of trajectory hops from triplet to singlet gradually drops from 0.8% to 0.3% as collision energy increases from 1 to 5 eV, as might be expected from shortening of the complex lifetime with increasing collision energy. We did another set of calculations for the O(3P) + C2(a3Πu) reaction where

Figure 6. Angular distribution expressed as normalized differential cross section (DCS, (2π/σ)(dσ/dΩ′)) for the C2(1Σ) + O(3P) and C2(3Π) + O(3P) reactions. For the CO(3Π) + C(3P) product in d, the spin-weighted differential cross section is presented.

symmetric scattering at Ecoll = 4 and 5 eV. With the O(3P) + C2(a3Πu) reaction, angular distributions for the CO(3Π) + C(3P), CO(X1Σ) + C(3P), and CO(1Σ) + C(1D) products (Figure 6b−d) have a forward peak at 1 eV, while it is more sideward and backward with increasing collision energies. All of these results are typical of reactions that proceed through intermediate complexes with lifetimes that are comparable to the rotational period.33 At Ecoll = 1 eV, one sees a preference for forward scattering in many of the results, which means that there is an excess population of complexes with lifetimes that allow for one-half the rotational period. At higher energies, the rotational period would be shorter and the excess forward scattering diminishes. At 5 eV there are some direct trajectories, but even these lead to broad angular distributions. Further insights are provided by the opacity functions that are presented in Figure S5 (Supporting Information). These show the variation of total reactivity with the impact parameter

Figure 8. Number (%) of trajectories (ntrjhop) that hops from triplet to singlet surface in the O(3P) + C2(X1Σ) reaction and singlet to triplet surface in the O(3P) + C2(a3Π) reaction as a function of collision energy (Ecoll). 26583

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

for the O(3P) + C2(a3Πu) reaction varies with electronic state, and we note that for reactions as a whole, 57−34% of the available energy ends up in electronic excitation as collision energy increases from 1 to 5 eV. For the CO(1Σ) + C(3P) and CO(1Σ) + C(1D) products, 60% of available energy (in excess of electronic energy) goes to product translation and 40% to product internal energy. For the CO(3Π) + C(3P) product, internal energy receives most (65%) of the available energy (in excess of electronic excitation) for Ecoll = 1 eV while 50% of available energy goes into product translation and 50% into product internal energy at 5 eV. For this CO(3Π) + C(3P) product, product vibration receives 90% of available internal energy while 10% goes to the product rotation. For the other products in the O + C2 reaction, product vibration receives 60− 80% of available internal energy at Ecoll = 1 eV, while energy deposition in product vibration decreases and that in product rotation increases with increasing collision energy from 1 to 5 eV. The angular distributions for O + C2 are in all cases dominated by backward−forward symmetry but with specific nonstatistical trends that arise from the short lifetime of the intermediate complexes. Trajectory surface-hopping calculations show that less than 1% of the trajectories hop from the triplet to singlet surface at 1 eV and the number gradually drops as collision energy increases in the O(3P) + C2(X1Σ+g ) reaction. For the singlet to triplet surface hopping in the O(3P) + C2(a3Πu) reaction, less than 2% trajectories hop and the number drops at hyperthermal energies. These results indicate that intersystem crossing is of minor importance for the O + C2 reaction, as makes sense given the short lifetime of the complexes. We identified a number of important features of the O + C2 reaction under LEO conditions related to reactivity, electronic branching, and product angular distributions and energy partitioning. While the C2 radical is not representative of species found on the surface of spacecraft, it is likely to be found in secondary reactions associated with LEO materials erosion and also in hydrocarbon combustion processes such as might arises in rocket exhaust. However, further experimental work is necessary to confirm our predictions about spin-allowed electronically excited products, angular distributions, and product energy partitioning. We hope our present study will stimulate further experimental and theoretical work to study the interactions of hyperthermal O(3P) atoms with carbon materials used for coating spacecraft and in hydrocarbon combustion.

trajectories were initiated from the singlet surface and then allowed to hop to the triplet surface. In this case also only a small number of trajectories hopped from the singlet to the triplet surface and spent the rest of the time on the triplet surface to give CO(1Σ) + C(3P). The number of trajectory hops from the singlet to the triplet increased from 0.3% to 2.0% with increasing collision energy from 1 to 3 eV but then dropped as collision energy increased from 3 to 5 eV. All these results indicate that ISC plays an unimportant role in the O + C2 reaction dynamics.

4. DISCUSSION AND CONCLUSION We studied the O + C2 reaction dynamics at hyperthermal collision energies that would be important for low Earth orbit (LEO) material erosion, and we also provided insights into the thermal kinetics and dynamics of this system. Quasiclassical trajectory (QCT) direct dynamics electronic structure calculations were performed for both O(3P) + C2(X1Σ+g ) and O(3P) + C2(a3Πu). Energies and forces were computed using the DFT(B3LYP) method with a 6-31G(d,p) basis set, and other electronic structure methods were used to validate the results. Intersystem crossing effects in the O + C2 reaction were investigated using a simplified version of the quasiclassical trajectory surface-hopping (TSH) method in which triplet− singlet transitions were only allowed at points where the two surfaces cross. The singlet C2 reaction O(3P) + C2(X1Σ+g ) occurs only on the triplet potential surface, while the triplet C2 reaction O(3P) + C2(a3Πu) occurs on singlet, triplet, and quintet potential surfaces. Electronic structure calculations with the DFT(B3LYP) method and 6-31G(d,p) basis set show that both the O(3P) + C2(X1Σ+g ) and the O(3P) + C2(a3Πu) reaction proceed through formation of two intermediates and one transition state in their electronic ground states but with no net barrier to reaction for either reaction. In addition, these reactions are highly exothermic, and this combined with the reagent hyperthermal energies allows for the spin-allowed formation of not only the ground state CO(X1Σ) + C(3P) but also, for the O(3P) + C2(a3Πu) reaction, formation of the electronically excited products CO(X1Σ) + C(1D), CO(3Π) + C(3P), and CO(3Π) + C(1D) . The QCT calculations show that in the absence of intersystem crossing only the ground state CO(X1Σ) + C(3P) products are formed in the O(3P) + C2(X1Σ+g ) reaction on triplet surfaces, while both the ground state CO(X1Σ) + C(3P) and the spin-allowed electronically excited products CO(3Π) + C(3P), CO(1Σ) + C(1D), and CO(3Π) + C(1D) are formed in the O(3P) + C2(a3Πu) reaction on singlet, triplet, and quintet surfaces as determined by degeneracy factors in Table 2. The CO + C cross section for O(3P) reacting with C2(a3Πu) is three times smaller than with C2(X1Σ+g ). The cross section for these reactions reflects formation of short-lived complexes with no activation energy and with the lower reactivity of the quintet surface playing a role in the reduced reactivity of C2(a3Πu). A large cross section at all collision energies strongly suggests that CO(3Π) + C(3P) is one of the major products of the O(3P) + C2(a3Πu) reaction. From these results, we predict that ground state product CO(1Σ) + C(3P) as well as the electronically excited products especially CO(3Π) + C(3P) can be found in the O + C2 reaction. We find that 50% of the available energy appears in product translation and 50% in product internal energy for the O(3P) + C2(X1Σg+) → CO(1Σ) + C(3P) reaction. Energy partitioning

■

ASSOCIATED CONTENT

S Supporting Information *

CASSCF potential energies for O(3P) + C2(X1Σ+g ) reaction on the triplet surface; geometry and potential energies for the O(3P) + C2(a3Πu) reaction on the lowest triplet and quintet surfaces; CASSCF potential energies for O(3P) + C2(a3Πu) reaction; angular distribution in O + C2 reaction at different potential surfaces; change of opacity function P(b) with impact parameters for different collision energies; discussion of opacity function; movie files showing mechanisms M1 and M2 corresponding to Figure 3. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. 26584

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585

The Journal of Physical Chemistry C

Article

Notes

(30) Taketsugu, T.; Gordon, M. S. J. Phys. Chem. 1995, 99, 14597− 14604. (31) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure; Van Nostrand Reinhold: New York, 1979; Vol. 4. (32) Shi, D.-H.; Li, W.-T.; Sun, J.-F.; Zhu, Z.-L. Int. J. Quantum Chem. 2012, DOI: 10.1002/qua.24036. (33) Troya, D.; Lakin, M. J.; Schatz, G. C.; Harding, L. B.; González, M. J. Phys. Chem. B 2002, 106, 8148−8160.

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by AFOSR Grant FA9550-10-10205 and by Lockheed Martin Corporation. We acknowledge Dr. Wenfang Hu for her help with the TSH code, Dr. Jeffrey T. Paci for his help with the product angular distribution code, and Dr. György Lendvay and Dr. Marcin Ziółkowski for helpful discussions.

■

REFERENCES

(1) Nicholson, K. T.; Minton, T. K.; Sibener, S. J. Prog. Org. Coat. 2003, 47, 443−447. (2) Murr, L. E.; Kinard, W. H. Am. Sci. 1993, 81, 152−165. (3) Murad, E. J. Spacecraft Rockets 1996, 33, 131−136. (4) Troya, D.; Schatz, G. C. Int. Rev. Phys. Chem. 2004, 23, 341−373. (5) Kern, R. D.; Xie, K.; Chen, H. Combust. Sci. 1992, 85, 77−86. (6) Baukal, C. E. Oxygen-enhanced combustion; CRC Press: New York, 1998. (7) Pilling, M. J.; Greenhalgh, D. A.; Hayhurst, A. N.; Lindstedt, R. P.; Smith, D. B.; Smith, I. W. M.; Wolfrum, J. Combustion Chemistry: Elementary Reactions to Macroscopic Processes. London, 2001; The Faraday Division, Royal Society of Chemistry: London, 2002, 119, 1− 475. (8) Zhang, Q. L.; O’Brien, S. C.; Heath, J. R.; Liu, Y.; Curl, R. F.; Kroto, H. W.; Smalley, R. E. J. Phys. Chem. 1986, 90, 525−528. (9) Gerhardt, P.; Löffler, S.; Homann, K. H. Chem. Phys. Lett. 1987, 137, 306−310. (10) Kroto, H. W.; McKay, K. Nature 1988, 331, 328−331. (11) Kaiser, R. I.; Balucani, N. Int. J. Astrobiol. 2002, 1, 15−23. (12) Weltner, W.; van Zee, R. J. Chem. Rev. 1989, 89, 1713−1747. (13) Nelson, H. H.; Helvajian, H.; Pasternack, L.; McDonald, J. R. Chem. Phys. 1982, 73, 431−438. (14) Richter, H.; Mazyar, O. A.; Sumathi, R.; Green, W. H.; Howard, J. B.; Bozzelli, J. W. J. Phys. Chem. 2001, 105, 1561−1573. (15) Becker, K. H.; Donner, B.; Freitas Dinis, C. M.; Geiger, H.; Schmidt, F.; Wiesen, P. Z. Phys. Chem. 2000, 214, 503−517. (16) Martin, M. J. Photochem. Photobiol. A 1992, 66, 263−289. (17) Reisler, H.; Mangir, M. S.; Wittig, C. J. Chem. Phys. 1980, 73, 2280−2286. (18) Kruse, T.; Roth, P. In Twenty-Seventh International Symposium on Combustion; The Combustion Institute, The University of Colorado at Boulder, Boulder, Colorado, USA, 1998; pp 193−200. (19) Garton, D. J.; Minton, T. K.; Troya, D.; Pascual, R.; Schatz, G. C. J. Phys. Chem. A 2003, 107, 4583−4587. (20) Troya, D.; Pascual, R. Z.; Garton, D. J.; Minton, T. K.; Schatz, G. C. J. Phys. Chem. A 2003, 107, 7161−7169. (21) Troya, D.; Pascual, R. Z.; Schatz, G. C. J. Phys. Chem. A 2003, 107, 10497−10506. (22) Garton, D. J.; Minton, T. K.; Hu, W.; Schatz, G. C. J. Phys. Chem. A 2009, 113, 4722−4738. (23) Hu, W.; Lendvay, G.; Maiti, B.; Schatz, G. C. J. Phys. Chem. A 2008, 112, 2093−2103. (24) Boese, A. D.; Martin, J. M. L. J. Chem. Phys. 2004, 121, 3405− 3416. (25) Vandeputte, A.; Reyniers, M.-F.; Marin, G. Theor. Chem. Acc. 2009, 123, 391−412. (26) Shao, Y.; Fusti-Molnar, L.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T. B.; Slipchenko, L. V.; Levchenko, S. V.; O’Neill, D. P.et al.. Q-Chem, Version 3.1 ed.; Q-Chem, Inc.: Pittsburgh, PA, 2007. (27) Fedorov, D. G.; Gordon, M. S. J. Chem. Phys. 2000, 112, 5611− 5623. (28) 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. J. J. Comput. Chem. 1993, 14, 1347−1363. (29) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007−1023. 26585

dx.doi.org/10.1021/jp3066629 | J. Phys. Chem. C 2012, 116, 26577−26585