ARTICLE pubs.acs.org/JPCA
Fulvenallene Decomposition Kinetics Daniela Polino and Carlo Cavallotti* Department Chimica, Materiali e Ingegneria Chimica “G. Natta”, Politecnico di Milano, via Mancinelli 7, 20131 Milano, Italy
bS Supporting Information ABSTRACT: While the decomposition kinetics of the benzyl radical has been studied in depth both from the experimental and the theoretical standpoint, much less is known about the reactivity of what is likely to be its main decomposition product, fulvenallene. In this work the high temperature reactivity of fulvenallene was investigated on a Potential Energy Surface (PES) consisting of 10 wells interconnected through 11 transition states using a 1 D Master Equation (ME). Rate constants were calculated using RRKM theory and the ME was integrated using a stochastic kinetic Monte Carlo code. It was found that two main decomposition channels are possible, the first is active on the singlet PES and leads to the formation of the fulvenallenyl radical and atomic hydrogen. The second requires intersystem crossing to the triplet PES and leads to acetylene and cyclopentadienylidene. ME simulations were performed calculating the microcanonical intersystem crossing frequency using LandauZener theory convoluting the crossing probability with RRKM rates evaluated at the conical intersection. It was found that the reaction channel leading to the cyclopentadienylidene diradical is only slightly faster than that leading to the fulvenallenyl radical, so that it can be concluded that both reactions are likely to be active in the investigated temperature (15002000 K) and pressure (0.0550 bar) ranges. However, the simulations show that intersystem crossing is rate limiting for the first reaction channel, as the removal of this barrier leads to an increase of the rate constant by a factor of 23. Channel specific rate constants are reported as a function of temperature and pressure.
1. INTRODUCTION The gas phase kinetics of toluene has been the subject of much experimental and theoretical research in the last years.120 Toluene is in fact not only the simplest alkylated aromatic species, but it is also found in considerable concentration in crude oils as well as in jet fuels and in gasoline, and it is formed easily in the gas phase during the pyrolysis or combustion of hydrocarbons. A detailed understanding of its decomposition kinetics and of its secondary chemistry would thus be extremely useful to model the initial stages of combustion or pyrolysis of the mentioned fuels as well as to understand if and how it can contribute to the nucleation and growth of soot. Currently, there is a general consensus on which are the main products formed through the unimolecular decomposition of toluene: the benzyl and the phenyl radicals. Kinetic constants for these two competitive reaction channels have been determined experimentally by several authors,6,9,2023 with the first measures dating back to the 1950s.24 More recently, the experimental data have been interpreted through refined ab initio and kinetic calculations, accounting explicitly for intermolecular energy transfer and thus extending the experimental estimates to an ample pressure and temperature range.19 The secondary chemistry started by toluene decomposition, in particular for what concerns the benzyl radical, is understood at a lower level of detail. It is known that benzyl undergoes unimolecular decomposition above 1450 K and that a good fitting of shock tube and H-ARAS data can be obtained assuming that the r 2011 American Chemical Society
main decomposition products are atomic hydrogen and a C7H6 chemical species, whose structure could not be experimentally resolved.6,2528 However, it was also found that a good modeling of the main products of toluene pyrolysis could be obtained assuming that benzyl decomposes directly to the cyclopentadienyl radical and acetylene.13,20,28 This assumption was necessary to be able to reproduce the high concentration of acetylene experimentally measured. These evidences suggest that, if C7H6 is formed from benzyl decomposition, then it is likely that a fast route for its further decomposition to acetylene exists as well. Recently, we have proposed on the basis of theoretical calculations that it is most likely that the C7H6 chemical species formed from benzyl decomposition is fulvenallene.29 Successive experimental and theoretical studies seemed to confirm the reasonability of this hypothesis.15,3033 These studies clearly indicate that a detailed investigation of the reactivity on the C7H6 potential energy surface (PES) might help greatly to understand how it can contribute to the overall toluene secondary chemistry. This was the subject of a paper companion to the present study, in which we studied using first principle calculations the interconversion kinetics active on the C7H6 PES and focused, in
Received: March 24, 2011 Revised: August 5, 2011 Published: August 05, 2011 10281
dx.doi.org/10.1021/jp202756s | J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Summary of the reactions that can take place on the C7H6 potential energy surface.
particular, on unimolecular decomposition reactions.34 The overall PES so determined is summarized in Figure 1. The simulations confirmed previous experimental3538 and theoretical3941 studies, indicating that several ring enlargement and ring contraction reactions are active on this PES and involve phenylcarbene (23), spiroheptatriene (7), cycloheptatetrane (5), ethynylcyclopentadiene (24, 25, 26), and fulvenallene (1), among the most stable chemical species that can be accessed from transition states having activation energies smaller than the lowest dissociation threshold energy. Several possible decomposition channels were also identified. In particular, two emerged as competitive as their kinetic constants were comparable: the reaction of decomposition into cyclopentadienylidene and acetylene, which takes place on the triplet PES, and that of dissociation into the fulvenallenyl radical and atomic hydrogen. The latter reaction has recently been proposed by da Silva and Bozelli as an important decomposition channel for fulvenallene.42 However, it seems unlikely
that it is the only decomposition channel for fulvenallene, since at least two different H-ARAS studies25,27 found that only one H atom is produced for the decomposition of each benzyl molecule, while the vertical decomposition mechanism benzyl f fulvenallene + H f fulvenallenyl + H predicts the formation of two H atoms. In the present work we will thus investigate whether alternative decomposition channels with respect to that leading to fulvenallenyl are possible.
2. METHOD AND THEORETICAL BACKGROUND Pressure-dependent rate constants for the multichannel, multiwell reacting system here investigated were determined using the combined RRKM/master equation theories implemented in the stochastic kinetic Monte Carlo (KMC) code we have recently developed.43 The data used as input for the simulations are vibrational frequencies, rotational constants, and energy of stationary points 10282
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A taken from a previous study in which we investigated the C7H6 PES.34 Briefly, in that work we determined geometries, vibrational frequencies, and rotational constants of all the intermediates and transition states using density functional theory (DFT) with B3LYP functionals44,45 and the 6-31+G(d,p) basis set. Energies were calculated at the CCSD(T) level46 using the cc-pVDZ and ccpVTZ basis set47 and then extrapolated to the infinite basis set adopting the scheme suggested by Martin.48 All energies were corrected with zero point energies (ZPE) calculated at the B3LYP/ 6-31+G(d,p) level. An exception to this approach is given by the decomposition reaction of fulvenallene to fulvenallenyl radical and hydrogen, which does not have a saddle point, as its potential energy surface has an asymptotic behavior. For this reason, the kinetic constant of this reaction has been determined using the microvariational form of transition state theory on a PES calculated adopting a multireference ab initio approach, because it can take into account more properly the contribution offered by unoccupied orbitals to the molecular wave function. Hence, intermediate structures were calculated as function of the reaction coordinate at the CASSCF(10,10)49,50 level adopting the 6-311G(d,p) basis set, while energies were calculated at the CASPT251,52/cc-pvTZ level and scaled by the bond dissociation energy calculated at the CCSD(T)/CBS level, as discussed in more detail in our previous paper.34 Also, transition state TS30, the last reaction step that leads to C5H4 and C2H2, has a multireference character, as its T1 diagnostic calculated at the UQCISD(T)/6-311G(d,p) level is 0.402. Thus, its structural parameters, vibrational frequencies, and energies were calculated with the same procedure used for TS1. An advantage of this approach is that the activation energies of the two decomposition channels are calculated at the same level of theory. All DFT and QCISD(T) calculations were carried out adopting the Gaussian 03 computational suite,53 while CASSCF, CASPT2, and CCSD(T) calculations were performed adopting the Molpro 2008.1 computational suite.54 Additional details about the CASPT2 calculations are reported as Supporting Information. Microcanonical rate constants k(E,J) were computed from the convoluted rotational and vibrational density of states as a function of energy E and angular momentum J. All the vibrational and external rotational degrees of freedom were considered in the calculation of k(E,J) for reactions with a chemical barrier. The approach adopted to determine k(E,J) for the reactions with a loose transition state is the E,J model proposed by Miller et al.55 The degeneration of low vibrational frequencies into torsional rotors was treated in the hindered rotor approximation. To take into account this contribution within the microcanonical ensemble it has been necessary to consider these degrees of freedom as internal rotations in computing the rovibrational density of states. It is known that the overall density of states can be expressed, within the approximation of separable degrees of freedom, as the convolution of the density of states of individual modes. Thus, we computed the density of states for the hindered rotors by direct count of the quantized energy levels that were within the energy range E to E + ΔE, with ΔE equal to 1 cm1. The spectrum of eigenvalues for each hindered rotor was obtained by a quantum eigenstate calculation, that is, solving the 1D rotational Schr€odinger equation, as described in previous works.56,57 Moreover, for those species whose internal hindered rotations were more than one, we first summed the eigenvalues related to the degrees of freedom considered and successively we convoluted the density of states as mentioned above. Successively, this density was convoluted with the vibrational density of states using the BeyerSwinehart algorithm with an initialization
ARTICLE
vector that contains the density of states previously evaluated for internal rotations, as suggested by Astholz et al.58 The PESs of the hindered rotors are reported as Supporting Information together with the respective inertia moments. The reaction that takes place through transition state TS* and that involves the W1 and W7 wells is of key importance as it occurs in correspondence of the intersection between the singlet and triplet potential energy surface (PES) of C7H6. The procedure used to determine the conical intersection is described in our previous paper in which we investigated the C7H6 PES.34 TS* lies on the intrinsic reaction coordinate connecting W1 and W7, so that, if the crossing between the singlet and triplet PES takes place, the reacting system evolves toward the minimum triplet energy well found in proximity of TS*, which is W18. The possibility to compute the probability of intersystem crossing was recently included in our code59 by evaluating the k(E,J) microcanonical rate as suggested by Harvey:60 Z E E0 F6¼ðEI , JÞphop ðE E0 EI ÞdEI 0 ð1Þ kðE, JÞ ¼ hFðE, JÞ In this case, the microcanonical rate constant depends on phop, the probability of intersystem hopping, which can be calculated from LandauZener theory61 as the probability of hopping in a double pass through the Minimum Energy Crossing Point (MECP) as phop ðEÞ ¼ ð1 PLZ Þð1 þ PLZ Þ
ð2Þ
where the LandauZener transition probability was calculated as rffiffiffiffiffi! 2 2πHSO μ ð3Þ PLZ ¼ exp pΔF 2E In this expression, HSO is the off-diagonal spinorbit coupling element of the 4 4 Hamiltonian matrix between the singlet and triplet state, ΔF is the relative slope of the two surfaces at the MECP, μ is the reduced mass of the reacting moieties, and E is the system kinetic energy. The spin orbit matrix element HSO was determined, adopting singlet and triplet wave functions calculated using multireference configuration interaction62,63 (MRCI) theory, including six orbitals and six electrons in the active space and adopting the 6-31G basis set. The forces at the MECP were calculated at the B3LYP/6-31+G(d,p) level. Master equation (ME) calculations were performed adopting the 1D master equation formulation proposed by Miller at al.55 implemented in our computational code, as described in detail in Barbato et al.43 The code is based on the stochastic integration of the master equation, which is performed discretizing the investigated rovibrational internal molecular energy into finite bins whose reacting properties were determined with RRKM theory. In the present work, the single well integration protocol of the master equation implemented in our previous works was extended to describe the unimolecular reaction dynamics on a multiple well PES. This is accomplished straightforwardly by allowing the excited molecule to react to a different energetic state, differing from the reacting well for the reaction energy change, with updated collisional transitional probabilities between different energy levels and RRKM reaction rates. The computational protocol is designed so that reaction and transitional collisional properties for all the considered wells and transition states are calculated before starting the ME integration 10283
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
cycle using an integration energy step of 1 cm1 and are then averaged over discrete energy bins whose size is determined through a series of test calculations, which are needed to find a reasonable compromise between computational efficiency and numerical accuracy. The computational efficiency is, however, independent from the size (number of wells and transition states) of the PES, which makes this approach particularly apt to study complex reacting systems, in which many wells can be visited from an excited complex before a reaction takes place. Several calculations were performed to test the correct implementation of the interwell transition dynamics. In particular, we checked the capability to reach the high pressure limit and the internal energy conservation for each jump between different wells. Additionally, we introduced a new procedure of evaluation of the kinetic constants in order to obtain more accurate and robust solutions. Three different approaches were thus used to determine the kinetic constant. The first is the approach proposed in our previous study, which is based on a linear regression over the simulation time that each molecule takes to react. The second is the approach adopted by Barker and co-workers64 in the Multiwell program suite, for which the kinetic constant is calculated as the inverse of the sum of the time that each molecule has taken to react divided by the total simulation time. The third is the direct evaluation of the kinetic constant from RRKM rates weighted over the energy distribution function calculated stochastically from the time spent in each energy bin. The convergence criteria can thus be defined as the limit at which the kinetic constants calculated using the three different approaches reach a constant value, which was defined as the iteration at which the student’s t test has a 1% confidence limit. The simulations revealed that this condition is reached after about 109 transitions over the simulated PES, including collisional energy transfers. Interestingly, it was found that for this value the results obtained with the three approaches differ by less than 5%, which is probably the intrinsic numerical error of the adopted computational procedure. It is here important to highlight that all the simulations were performed assuming that the reactions leading to the formation of bimolecular products (i.e., TS1 and TS30) are irreversible sinks. The rate coefficients are then calculated assuming that the rate law can be expressed simply as dni ¼ kD ni dt
ð4Þ
where ni is the concentration of the reactant (fulvenallene in this case). The kinetic constant kD for each reaction channel is then calculated through linear regression over the stochastic reaction times. There are two subtle but important assumptions in calculating the rate coefficients using eq 4, which are the elimination of the backward reaction channel and the hypothesis that the contribution to the rate coefficient from stabilized intermediate wells is not significant. The rate coefficient determined through this approach represents thus a flux coefficient, that is, the coefficient that multiplied by the reactant concentration gives the number of molecules that per unit time and volume are converted to the reaction product. As clearly explained by Widom in a seminal paper published several decades ago,65 this coefficient is not equal to the phenomenological rate coefficient, which is the rate coefficient accessible to experimental measurements. The difference between phenomenological rate coefficients and reaction fluxes is subtle and counterintuitive, as
reviewed recently by Miller and Klippenstein.66 However, it is likely that this difference may not be too large for the case under investigation, because the products of the two decomposition reactions considered here are usually present in very low concentrations in reaction environments in which these reactions are active due to fast secondary reaction channels so that the irreversible sink assumption most probably holds and because the stabilization of intermediate wells is negligible, as discussed in the next section.
3. RESULTS AND DISCUSSION The aim of this study is to determine which are the main decomposition products of fulvenallene and to evaluate the respective reaction rates as a function of temperature and pressure using a combined RRKM/ME approach. The C7H6 PES on which the present calculations are based is illustrated in Figure 1 and comprehends 27 wells connected through 39 paths. A preliminary analysis showed that only a small portion of this PES is likely to be visited for a reaction pathway started from the fulvenallene well so that it can significantly be simplified without loss of accuracy. For this purpose, the rate of each reaction was compared with that of the first accessible dissociation channel, which is given by TS1 and leads to the formation of the fulvenallenyl radical and atomic hydrogen, assuming equilibrium between the reactant (fulvenallene) and the corresponding transition state. Only those paths with a high pressure rate constant faster or comparable with that of TS1 were included in the kinetic scheme, as the wells with lower reactivity would have hardly been visited because at least one dissociation reaction would proceed faster. It was thus found that, of the six possible dissociation channels, only two have comparable rates: TS1 and TS30. In fact, though other dissociation reactions to acetylene and a C5H4 isomer may have rates comparable to that of TS1 (notably, TS21 and TS26), the accessibility to the wells to which they are connected is hindered by significant energy barriers or low pre-exponential factors of reactions positioned previously on the PES. The PES considered in the simulations is illustrated in Figure 2. It is composed of 10 wells connected by 11 transition states. Two decomposition channels, proceeding through TS1 and TS30, are possible. Of all the wells reported in Figure 2, only W18 lies on the triplet PES. Two main decomposition channels, thus, appear as competitive: the dissociation to fulvenallenyl and atomic hydrogen through TS1 and the decomposition to cyclopentadienylidene and acetylene through TS30. Despite the higher activation energy, the second reaction channel is faster if the decomposition rate is computed assuming equilibrium between the involved wells because of its high pre-exponential factor. TS30 transition state has in fact a loose character and is, thus, characterized by five low vibrational frequencies, two of which degenerate into hindered rotors, as evidenced by an analysis of the rotational PESs. Details about the hindered rotor treatment are reported as Supporting Information. To check whether the approach used to determine the TS30 kinetic constant is consistent with experimental data, we calculated the kinetic constant of the inverse process, the addition of C2H2 to C5H4, and compared it to the reaction of C2H2 and C2H. We believe the two reactions are similar in several ways. First, because one of the reactants, C2H2, is the same, then because in both reactions C2H2 adds to highly unsatured radicals, 10284
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
Figure 2. Portion of the C7H6 PES considered in the simulations. Energies were calculated at the CCSD(T) level using the cc-pVDZ and cc-pVTZ basis sets and extrapolated to the complete basis set limit. ZPE corrections are included and were calculated using vibrational frequencies determined at the B3LYP/6-31+G(d,p).
C2H and C5H4, and finally for the fact that both reactions are almost barrierless. The kinetic constant measured for the addition of C2H to C2H2 at 1600 K by Shin and Michael is 1.6 1014 cm3/mol/s,67 while the one we calculated for the addition of C2H2 to C5H4 is 6.3 1013 cm3/mol/s. The similarity suggests that the transition state parameters determined for TS30, which are the same used to calculate the addition of C5H4 to C2H2, are reasonable. However, though TS30 appears to be faster than TS1, it should be remembered that this reaction takes place on the triplet PES, so that intersystem crossing may be limiting. Moreover, eventual fall off effects would favor the first decomposition channel. Given the high rate of the TS30 reaction channel, it is possible that also the decomposition to acetylene and the singlet state of C5H4, which lies 6.4 kcal/mol above the triplet and does not require intersystem crossing, may contribute to the overall C7H6 decomposition as the reaction is barrierless. In the proposed reaction scheme, the reaction channel leading to the decomposition of fulvenallene into cyclopentadienylidene and acetylene involves intersystem crossing from the singlet to the triplet state. The intersystem crossing frequency was determined using LandauZener theory, as described in the Method and Theoretical Background. This approach requires as input structure, energy, spinorbit coupling coefficient, force constants, and vibrational frequencies at the conical intersection between the singlet and the triplet PES. Conical intersections that may occur on the C7H6 PES were searched in proximity of the wells that may be reached from W1 scanning the PES at the B3LYP/6-31+G(d,p) level. The well-known singlet triplet crossing of phenylcarbene41 was not considered in the calculations as it is unlikely that phenylcarbene may be accessed from W1. Among the several conical intersections that were located, the one with the highest spinorbit coupling, TS*, was found on the intrinsic reaction coordinate connecting W1 to W7, after TS10. TS* was considered to be accessible both from W1 and from W7
in the simulations. The spinorbit matrix element HSO calculated for the corresponding MECP is 9.3 cm1. The difference between the forces acting on the two states is more difficult to determine, because an intrinsic reaction coordinate scan performed on both the triplet and the singlet PESs showed that the conical intersection takes place contextually to a rototraslation of the terminal CH group, so that it is difficult to define a reaction coordinate. The force difference was thus computed as the difference between the forces acting on the angle between the CH group and the C_C5H4 moiety on the singlet and triplet surfaces, as the IRC analyses showed that it is the coordinate that matches most closely the reaction path. The value so calculated is 9.4 kcal mol1 Å1. The reduced mass used in the calculations is that of CH. These parameters were used in a first set of simulations. However, to take into account the possibility of the existence of an alternative fast intersystem crossing channel and considering the uncertainty of many parameters, a second set of simulations was performed assuming that intersystem crossing is fast, which was obtained imposing a probability Fhop of 1 to eq 1. The master equation was solved assuming that the bath gas is pure Argon, adopting an exponential down model for the collisional energy transfer and using LennardJones collision parameters, which were assumed to be the same as those of toluene for all the C7H6 isomers.9,55,68 Vibrational frequencies, rotational constants, and activation energies of the relevant intermediates and saddle points located on the investigated potential energy surface were calculated from ab initio simulations, as described in the Method and Theoretical Background. Because energy transfer parameters are not well-known for these molecules, we assumed an empirical value for the mean downward transfer energy ÆΔEdownæ equal to 2000 cm1, which was used by da Silva et al.33 to study the benzyl decomposition reaction. The sensitivity of the calculated rate constant to this 10285
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
Figure 3. Overall kinetic constants for C7H6 decomposition. Simulations were performed (a) accounting for intersystem crossing using LandauZener theory and (b) assuming that intersystem crossing is fast.
Table 1. Kinetic Constant for the Decomposition of C7H6 Calculated Accounting for Intersystem Crossing Using LandauZener Theory pressure (bar) T (K)
0.05
1500
5.72 102
6.67 102
1.06 103
1.15 103
1.18 103
1.21 103
1.43 103
1600
2.52 10
3
3.33 10
4.81 10
5.71 10
3
6.77 10
3
7.45 10
7.92 103
1700
8.01 10
3
1.10 10
1.96 10
2.30 10
4
3.05 10
4
3.33 10
3.64 104
1800
2.09 10
4
3.06 10
5.93 10
7.28 10
4
1.05 10
5
1.15 10
1.34 105
1900
4.50 10
4
6.67 10
1.50 10
1.93 10
5
3.04 10
5
3.43 10
3.97 105
2000
8.37 10
4
1.29 10
3.16 10
4.31 10
5
7.38 10
5
8.61 10
1.07 106
0.1
0.5 3 4 4 4 5
1.0 3 4 4 5 5
5.0
10.0
50.0 3 4 5 5 5
Table 2. Kinetic Constant for the Decomposition of C7H6 Calculated Assuming that Intersystem Crossing is Fast pressure (bar) T (K)
0.05
1500
0.1
0.5
1.0
6.02 10
2
7.67 10
1.12 10
1600
2.72 10
3
3.41 10
1700
8.35 10
3
2 3
5.0
1.36 10
3
6.01 10
7.52 10
3 4
3 3
10.0
50.0
1.70 10
3
1.55 10
1.79 103
9.86 10
3
1.05 10
1.16 104
4.49 10
4
4
3 4
1.20 10
2.29 10
2.99 10
5.24 10
6.26 104
1800 1900
4
2.21 10 4.64 104
4
3.20 10 7.11 104
4
6.89 10 1.76 105
4
9.29 10 2.42 105
5
1.57 10 4.69 105
5
1.88 10 5.72 105
2.43 105 8.07 105
2000
8.58 104
1.37 105
3.65 105
5.31 105
1.14 106
1.46 106
2.23 106
4
4
parameter has been examined. Because we were interested in the initial product branching ratio of the C7H6 decomposition reaction, we used an “infinite sink” approximation for the two main products; that is, the reverse C2H2 + C5H4 and H + C7H5 reactions have been neglected. To solve the master equation, the reported values were first averaged over J, as proposed by Miller et al.55 and then in discrete reaction bins. To determine the effect of the bin size over the accuracy of the calculations, a sensitivity analysis was performed on this parameter. In particular, the bin size was systematically reduced from 100 cm1 to 50 and 25 cm1. The simulations revealed that the rate constants determined by solving the ME are not particularly sensitive to the size of the energy bin (i.e., the deviation from the 100 cm1 reference value is smaller than 5%), as long as the maximum energy allowed for a collisional energy jump is at least 9000 cm1. Also, we performed a sensitivity analysis on the ÆΔEdownæ parameter. It was found that a variation of (1000 cm1 at 1 atm leads to a
maximum change of the rate constant of approximately (50% at high temperatures (i.e., greater than 1900 K), while no significant changes are detected at low temperatures (i.e., lower than 1700 K). This behavior was expected as the variation of the ÆΔEdownæ parameter has an effect only when the rate constant is deep in the fall off regime. Also, the branching rate constants change both accordingly and this yields to a product branching ratio only slightly affected by the variation of ÆΔEdownæ. Unfortunately, as no experimental kinetic data are available for this reaction, it remains uncertain which is the most appropriate value. The global rate coefficients for C7H6 decomposition determined both for the case in which intersystem crossing was considered and when it was not are illustrated in Figure 3 and reported in Tables 1 and 2. Calculations were performed in the 0.0550 bar pressure range and in the 15002000 K temperature range. The reaction is in its fall off regime for pressures up to 12 bar when intersystem crossing was considered and up to 10 bar when it was not. 10286
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
Figure 4. Branching reaction constants for channels 1 and 2 calculated (a) accounting for intersystem crossing using LandauZener theory and (b) assuming that intersystem crossing is fast.
Figure 6. Product branching ratio as a function of temperature and pressure calculated accounting for intersystem crossing using Landau Zener theory (lowest points) and assuming that intersystem crossing is fast (upper points).
Figure 5. Energy distribution functions of W1 and W18 calculated at 1 bar and 1700 K and Boltzmann population of W1.
In addition to the total rate constant estimation, one of the attempts of the present study is to provide an insight into the mechanism and the product branching under different conditions. The channel specific rate constants for decomposition into cyclopentadienylidene and C2H2 (channel 1) and into the fulvenallenyl radical (channel 2) are reported in Figure 4. It can be observed that the two pathways have similar rate coefficients for different pressures and temperatures if intersystem crossing is explicitly considered in the calculations, while if intersystem crossing is assumed to be fast channel 1 is up to a factor three faster than channel 2. Furthermore, it is interesting to observe that the two channels exhibit a different pressure dependence. In fact, the homolytic dissociation reaction rate has a minor dependence with pressure compared to the other reaction. This can also be seen in Figure 5 where are represented the population densities of wells W1 and W18 compared with the Boltzmann population of W1. It is thus interesting to observe that the population density of W1 deviates from the Boltzmann population only at energies higher than the corresponding dissociation threshold energy, confirming that the direct decomposition to C7H5 and H is hardly affected by fall off behavior even at high temperature. The second reaction channel, instead, involves the transition from well W1 to well W18. The population of the most excited states of W18 is more depleted with respect to the threshold energy reached in TS3, which leads to a significant pressure dependence. The slight deviation of the W1 energy distribution function from the Boltzmann population that can be observed at low energies can be ascribed to the fact that simulations are started from energy 0, which biases slightly the energy distribution function. It is also interesting to observe that no significant population of energies below the TS* critical energy is observed for W18, which means that collisional stabilization is not significant.
Table 3. Kinetic Constants for the Two Fulvenallene Reaction Channels Calculated Considering Explicitly Intersystem Crossing (Ghop = f(PLZ) and under the Assumption that It is Fast (Fhop = 1)a C7H6 f C5H4 + C2H2 P (bar) Log10A Fhop = f(PLZ)
Fhop = 1
0.1
R
Ea
C7H6 f C7H5 + H Log10A
R
Ea
130.10 32.19 173.97 116.01 28.46 157.86
1
100.99 23.86 154.21 66.69
10
61.12
12.81 120.50 27.07
14.64 119.55 3.54 89.15
0.1
135.37 33.58 179.46 102.86 24.94 144.42
1
124.48 30.18 176.46 108.80 26.20 156.95
10
129.63 31.11 191.02 98.54
23.11 153.51
Kinetic constants expressed as k = A T(K)R exp(Ea(kcal/mol)/RT) in s1. a
In Figure 6 it is displayed the product branching ratio, defined as the ratio of the number of molecules produced by channel 1 over the number of molecules produced by channel 2, as a function of pressure and temperature. When intersystem crossing is explicitly considered, the branching ratio is almost constant with both temperature and pressure with a mean value of 1.06. In the assumption of fast intersystem crossing, instead, the branching ratio varies proportionally with pressure and temperature, ranging from 0.8 up to 3.8. This behavior is motivated by two factors: channel 2 has a reaction energy barrier smaller than that of channel 1, hence, the higher dependence of the branching ratio with temperature, while channel 1, because of fall off effects, has a greater dependence from pressure than channel 2. Channel specific rate constants were interpolated between 1500 and 2000 K and are reported at three different pressures as a function of temperature in Table 3.
4. CONCLUSIONS The decomposition kinetics of fulvenallene was investigated using RRKM/ME theory studying the reaction dynamics on a comprehensive PES in which have been included the most relevant 10287
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A intermediates and reactive paths. The structural and energetic parameters needed to carry out the simulations were determined through ab initio calculations. It was found that C7H6 can decompose mainly through two reaction channels: direct H-elimination to form the fulvenallenyl radical and decomposition to cyclopentadienylidene and C2H2. Our calculations demonstrate that both reactions are feasible and are likely to happen with almost equal probability if the intersystem crossing probability is explicitly considered in the calculations. However, if intersystem crossing is fast, the decomposition to cyclopentadienylidene and acetylene becomes significantly favored, in particular, at high temperatures and pressures. Further experimental studies focusing directly on the C7H6 high temperature reactivity would help greatly to achieve a better comprehension of the mechanism of this important reaction, in particular, considering that several C7H6 isomers, notably phenylcarbene, can be produced directly from the thermolysis or photolysis of phenyldiazomethane.69
’ ASSOCIATED CONTENT
bS
Supporting Information. Details about CASPT2 and CASSCF calculations. Details about hindered rotors calculations. Structures, vibrational frequencies, and inertia moments of the 10 wells and 10 saddle points of the C7H6 potential energy surface considered in the RRKM simulations. This material is available free of charge via the Internet at http://pubs.acs.org.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ REFERENCES (1) Ikeda, N.; Nakashima, N.; Yoshihara, K. J. Chem. Phys. 1985, 82, 5285. (2) Pamidimukkala, K. M.; Kern, R. D.; Patel, M. R.; Wei, H. C.; Kiefer, J. H. J. Phys. Chem. 1987, 91, 2148. (3) Tsukiyama, K.; Bersohn, R. J. Chem. Phys. 1987, 86, 745. (4) Brouwer, L. D.; Mullermarkgraf, W.; Troe, J. J. Phys. Chem. 1988, 92, 4905. (5) Ackermann, L.; Hippler, H.; Pagsberg, P.; Reihs, C.; Troe, J. J. Phys. Chem. 1990, 94, 5247. (6) Hippler, H.; Reihs, C.; Troe, J. Z. Phys. Chem. (Neue Folge) 1990, 167, 1. (7) Hippler, H.; Troe, J. J. Phys. Chem. 1990, 94, 3803. (8) D’Alessio, J.; Lazzaro, M.; Massoli, P.; Moccia, V. Opt. Laser Eng. 2002, 37, 495. (9) Eng, R. A.; Gebert, A.; Goos, E.; Hippler, H.; Kachiani, C. Phys. Chem. Chem. Phys. 2002, 4, 3989. (10) Sivararnakrishnan, R.; Tranter, R. S.; Brezinsky, K. Combust. Flame 2004, 139, 340. (11) Sivaramakrishnan, R.; Tranter, R. S.; Brezinsky, K. Proc. Combust. Inst. 2005, 30, 1165. (12) Oehlschlaeger, M. A.; Davidson, D. F.; Hanson, R. K. J. Phys. Chem. A 2006, 110, 9867. (13) Sivaramakrishnan, R.; Tranter, R. S.; Brezinsky, K. J. Phys. Chem. A 2006, 110, 9400. (14) Sivaramakrishnan, R.; Tranter, R. S.; Brezinsky, K. J. Phys. Chem. A 2006, 110, 9388. (15) Detilleux, V.; Vandooren, J. J. Phys. Chem. A 2009, 113, 10913. (16) Zhang, T. C.; Zhang, L. D.; Hong, X.; Zhang, K. W.; Qi, F.; Law, C. K.; Ye, T. H.; Zhao, P. H.; Chen, Y. L. Combust. Flame 2009, 156, 2071. (17) Granata, S.; Faravelli, T.; Ranzi, E. Combust. Flame 2003, 132, 533.
ARTICLE
(18) Colket, M. B.; Seery, D. J. Proc. Combust. Inst. 1994, 25, 883. (19) Klippenstein, S. J.; Harding, L. B.; Georgievskii, Y. Proc. Combust. Inst. 2007, 31, 221. (20) Smith, R. D. J. Phys. Chem. 1979, 83, 1553. (21) Cavallotti, C.; Mancarella, S.; Rota, R.; Carra, S. J. Phys. Chem. A 2007, 111, 3959. (22) Braun-Unkhoff, M.; Frank, P.; Just, T. Proc. Combust. Inst. 1989, 22, 1053. (23) Sivaramakrishnan, R.; Michael, J. V. Proc. Combust. Inst. 2011, 33, 225. (24) Szwarc, M. J. Chem. Phys. 1948, 16, 128. (25) Braun-Unkhoff, M.; Frank, P.; Just, T. Ber. Bunsen-Ges. Phys. Chem. 1990, 94, 1417. (26) Oehlschlaeger, M. A.; Davidson, D. F.; Hanson, R. K. J. Phys. Chem. A 2006, 110, 6649. (27) Sivaramakrishnan, R.; Su, M.-C.; Michael, J. V. Proc. Combust. Inst. 2011, 33, 243. (28) Rao, V. S.; Skinner, G. B. Proc. Combust. Inst. 1986, 21, 809. (29) Cavallotti, C.; Derudi, M.; Rota, R. Proc. Combust. Inst. 2009, 32, 115. (30) Detilleux, V.; Vandooren, J. Proc. Combust. Inst. 2011, 33, 217. (31) Li, Y. Y.; Cai, J. H.; Zhang, L. D.; Yuan, T.; Zhang, K. W.; Qi, F. Proc. Combust. Inst. 2011, 33, 593. (32) Zhang, L. D.; Cai, J. G.; Zhang, T. C.; Qi, F. Combust. Flame 2010, 157, 1686. (33) da Silva, G.; Cole, J. A.; Bozzelli, J. W. J. Phys. Chem. A 2009, 113, 6111. (34) Polino, D.; Famulari, A.; Cavallotti, C. J. Phys. Chem. A 2011, 115, 7928. (35) Wentrup, C. Adv. Heterocycl. Chem. 1981, 28, 231. (36) Joines, R. C.; Turner, A. B.; Jones, W. M. J. Am. Chem. Soc. 1969, 91, 7754. (37) Johnson, R. P. Chem. Rev. 1989, 89, 1111. (38) Baron, W. J.; DeCamp, M. R.; Hendrick, M. E.; Jones, M. J.; Levin, R. H.; Sohn, M. B. Carbenes from diazo compounds. In Carbenes; Jones, M., Jr., Moss, R. A., Eds.; Wiley: New York, 1973; Vol. I, pp 1. (39) Matzinger, S.; Bally, T.; Patterson, E. V.; McMahon, R. J. J. Am. Chem. Soc. 1996, 118, 1535. (40) Schreiner, P. R.; Karney, W. L.; Schleyer, P. V.; Borden, W. T.; Hamilton, T. P.; Schaefer, H. F. J. Org. Chem. 1996, 61, 7030. (41) Wong, M. W.; Wentrup, C. J. Org. Chem. 1996, 61, 7022. (42) da Silva, G.; Bozzelli, J. W. J. Phys. Chem. A 2009, 113, 12045. (43) Barbato, A.; Seghi, C.; Cavallotti, C. J. Chem. Phys. 2009, 130, 074108. (44) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (45) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (46) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479. (47) Woon, D. E.; Dunning, T. H. J. Chem. Phys. 1995, 103, 4572. (48) Martin, J. M. L. Chem. Phys. Lett. 1996, 259, 669. (49) Knowles, P. J.; Werner, H. J. Chem. Phys. Lett. 1985, 115, 259. (50) Werner, H. J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. (51) Finley, J.; Malmqvist, P. A.; Roos, B. O.; Serrano-Andres, L. Chem. Phys. Lett. 1998, 288, 299. (52) Roos, B. O.; Andersson, K.; Fulscher, M. P.; Malmqvist, P. A.; SerranoAndres, L.; Pierloot, K.; Merchan, M. Multiconfigurational perturbation theory: Applications in electronic spectroscopy. In Advances in Chemical Physics; John Wiley and Sons: New York, 1996; Vol. 93, Vol Xciii, p 219. (53) Frisch, M. J.; ; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A.Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; 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.; Ayala, P. Y.; Morokuma, K.; 10288
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
The Journal of Physical Chemistry A
ARTICLE
Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision B.05; Gaussian, Inc.: Pittsburgh, PA, 2003. (54) Werner, H.-J.; Knowles, P. J.; Lindh, R.; Manby, F. R.; Sch€utz, M.; Celani, P.; Korona, T.; Mitrushenkov, A.; Rauhut, G.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hetzer, G.; Hrenar, T.; Knizia, G.; K€oppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; Palmieri, P.; Pfl€uger, K.; Pitzer, R.; Reiher, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. MOLPRO, version 2008.1, a package of ab initio programs; http://www.molpro.net. (55) Miller, J. A.; Klippenstein, S. J.; Raffy, C. J. Phys. Chem. A 2002, 106, 4904. (56) Fascella, S.; Cavallotti, C.; Rota, R.; Carra, S. J. Phys. Chem. A 2004, 108, 3829. (57) Fascella, S.; Cavallotti, C.; Rota, R.; Carra, S. J. Phys. Chem. A 2005, 109, 7546. (58) Astholz, D. C.; Troe, J.; Wieters, W. J. Chem. Phys. 1979, 70, 5107. (59) Polino, D.; Barbato, A.; Cavallotti, C. Phys. Chem. Chem. Phys. 2010, 12, 10622. (60) Harvey, J. N. Phys. Chem. Chem. Phys. 2007, 9, 331. (61) Zener, C. Proc. R. Soc. A 1932, 137, 696. (62) Werner, H. J.; Knowles, P. J. J. Chem. Phys. 1988, 89, 5803. (63) Knowles, P. J.; Werner, H. J. Chem. Phys. Lett. 1988, 145, 514. (64) Barker, J. R. Int. J. Chem. Kinet. 2001, 33, 232. (65) Widom, B. Science 1965, 148, 1555. (66) Miller, J. A.; Klippenstein, S. J. J. Phys. Chem. A 2006, 110, 10528. (67) Shin, K. S.; Michael, J. V. J. Phys. Chem. 1991, 95, 5864. (68) Hippler, H.; Troe, J.; Wendelken, H. J. J. Chem. Phys. 1983, 78, 6709. (69) Wentrup, C. Top. Curr. Chem. 1976, 62, 173.
10289
dx.doi.org/10.1021/jp202756s |J. Phys. Chem. A 2011, 115, 10281–10289
