Article pubs.acs.org/JPCA

Oxidation of Graphene-Edge Six- and Five-Member Rings by Molecular Oxygen Ravi I. Singh,‡ Alexander M. Mebel,*,† and Michael Frenklach*,‡,§ †

Department of Chemistry and Biochemistry, Florida International University, Miami, Florida 33199, United States Department of Mechanical Engineering, University of California at Berkeley, Berkeley, California 94720-1740, United States § Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States ‡

S Supporting Information *

ABSTRACT: To gain qualitative and quantitative understanding of oxidation processes of large polycyclic aromatics, soot particles, and graphene edges, a theoretical study is reported for the pyrenyl−O2 reaction system. First, possible reaction pathways and their energetics were investigated using high-level ab initio calculations. The results were utilized in RRKM−master equation calculations of rate coefficients and relative product yields at temperatures and pressures relevant to combustion. Finally, the deduced oxidation mechanisms of six- and fivemember rings and the computed rate coefficients were employed in kinetic Monte Carlo simulations of oxidation of a graphene “molecule” evolving in flame-like environments. Among the major findings from the latter simulations are the following: The oxidation system exhibits two basic pathways, thermal decomposition and regeneration of oxyradicals. Their competition is temperature-dependent, with the former dominating at higher and the latter at lower temperatures. The overall oxidation of the graphene substrate is computed to be time-dependent, with the initial rates consistent with the known experimental data. aromatic rings: “free migration” over the growing graphene edge16 and “embedded migration” within the grown layer of graphene.17 The migrating five-member rings turned out to transform (isomerize) into six-member rings.19 The evolution of graphene edges with the newly established five-member-ring chemistry demonstrated bending of graphene sheets.22 All of these studies were focused on PAH and graphene-edge growth chemistry. With the recent attention given also to the oxidation, new information began to appear on elementary reactions of aromatics oxidation.23−29 However, the very initial attempts to examine graphene-edge oxidation with the combined detailed growth and oxidation models encountered fundamental problems. The modeling, described in more detail later in the text, was carried out using detailed kinetic Monte Carlo (KMC) simulations coupled with molecular-dynamics (MD) geometry optimization steps.22 Performing KMC-MD simulations of a pregrown graphene “sheet” with an excess of molecular oxygen showed that degradation of the graphene moiety quickly comes

1. INTRODUCTION One of the persistent problems associated with utilization of fossil fuels is the formation of soot during combustion. Soot leaving exhaust signifies incomplete combustion of fuel carbon and hence a reduced thermal efficiency.1 The emitted soot is considered to be a significant pollutant with adverse effects on human health2,3 and global climate.4,5 In light of these facts, formation of soot particles has been studied extensively, with establishing the fundamental mechanism of its formation being one of the major goals.6,7 Structurally, it has long been suggested that soot particles are composed of turbostratically arranged polycyclic aromatic molecules (or, in other words, graphene flakes).8,9 Discovery of fullerenes10 led Kroto and co-workers to propose that a spherical soot particle could form through a continuously grown fullerene.11 While such a single-molecule mechanism was argued to be implausible to explain particle inception,12 it is conceivable that a curved graphene-edge front, once formed, can propagate surface growth retaining particle sphericity.13 Mechanistic modeling of soot-particle surface chemistry began with postulating chemical similarity between growth of particle constituent graphene edges and gas-phase polycyclic aromatic hydrocarbons (PAH).14,15 Being able to explain much of the soot formation phenomena,6 the initial models limited the aromatic extension to the planar growth. A much richer chemistry was discovered in subsequent theoretical studies,16−21 revealing migration of five-member © XXXX American Chemical Society

Special Issue: 100 Years of Combustion Kinetics at Argonne: A Festschrift for Lawrence B. Harding, Joe V. Michael, and Albert F. Wagner Received: January 27, 2015 Revised: April 13, 2015

A

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

detected because apparently it was unstable under the experimental conditions. The oxidation kinetics of the naphthyl radical was studied by Marinov et al.,43 who estimated the rate coefficient of the 1C10H7 + O2 reaction in their kinetic modeling of a laminar premixed n-butane flame, and by Lin and co-workers,44 who reported thermal rate coefficients for the 2-naphthyl + O2 system at 299−444 K and concluded that at these low temperatures the reaction proceeds via the associationstabilization mechanism to the naphthyl peroxy radical C10H7OO, which may then lose the oxygen atom yielding the naphthoxy radical, C10H7O. Mebel and co-workers 23,24 reported the most detailed ab initio studies of the potential energy surfaces (PESs) for the C6H5 + O2 and C10H7 + O2 reactions, and these results were recently utilized in RRKMmaster equation (RRKM-ME) computations of temperatureand pressure-dependent rate coefficients and product branching ratios under combustion conditions, at temperatures of 1500, 2000, and 2500 K and pressures of 0.01−10 atm.29 The latter study found that the dominant reaction channel in all cases is elimination of the oxygen atom from peroxy complexes formed at the initial O2 addition step, leading to the phenoxy and naphthoxy radical products, with its contribution increasing with temperature. Chemically activated phenoxy and naphthoxy radicals either decompose to the cyclopentadienyl + CO and indenyl + CO products, respectively, thus completing the conversion of the six-member ring to a five-member ring upon oxidation, or thermally equilibrate, with the relative yields of the two processes strongly depending on temperature and pressure. At the lowest temperature of 1500 K, the reactions also yield significant amounts of pyranyl/benzopyranyl + CO (up to 24%). Although the deduced trends in the oxidation kinetics of phenyl and naphthyl radicals appeared to be similar, the prior studies29 concluded that the size of the molecule and the position of the radical site in it may affect the relative product yields. To understand the oxidation mechanism of larger PAH molecules as well as soot particles and graphene sheets, a fuller set of reactions needs to be considered, including not only oxidation of aromatic radicals with O2 but also creation of the radical sites and competing reactions on those sites. The information in the literature concerning rate coefficients and the reaction mechanism of oxidation of larger PAH radicals is rather limited. Edwards et al.28 investigated the oxidation reaction of phenanthryl radical with OH using B3LYP and CBS-QB3 calculations of pertinent potential energy surfaces and RRKM-ME calculations of rate coefficients. They found that the reaction proceeds by OH addition to the radical site followed by H atom migration/elimination and by oxyradical decomposition via expulsion of the CO group. Although the oxyradical decomposition turning a six-member ring into a fivemember ring should be common for the oxidation processes both with O2 and OH, rate coefficients for the other steps involved in phenanthryl oxidation by molecular oxygen are not available. Raj et al.45 studied structural effects on the oxidation of soot particles in a combined experimental and theoretical study. Experimentally, they found that curved soot particles are oxidized easier than planar ones. To explain this observation and to understand the role of the PAH structure in the reactivity toward O2, these authors carried out theoretical calculations using density functional B3LYP and Hartree−Fock methods of the reaction pathways of 4-pyrenyl and 1-

to a stop. Inspection revealed that the graphene edge becomes “non-reactive” because of accumulation of embedded fivemember rings. Our initial model did not include direct oxidation of five-member rings based on the presumption of rapid thermal desorption of five-member rings migrating over the graphene edge and, consequently, infrequency of their embedding into the graphene layer. The numerical simulations disproved this assumption. One of the reasons for this is that the source of five-member rings in an oxidative environment is not only acetylene adsorption to a zigzag site or thermal isomerization of six-member-ring complexes but also oxidation of six-member rings. Another problem that appeared in the initial KMC-MD simulations was the formation of unrealistic compounds. For instance, starting with a pyrene molecule, the initial oxidation converts one of its six-member rings to an embedded fivemember ring

Following the same reaction mechanism leads to the formation of two adjacent five-member rings

which defies the so-called isolated pentagon rule.30−33 Thus, the primary purpose of the present study is to investigate theoretically the mechanism and products of oxidation of five-member rings embedded into a graphene layer, in comparison to the oxidation of six-member rings. Oxidation of six-member aromatic rings with molecular oxygen have been previously studied both experimentally and theoretically in small molecules (phenyl radical) and prototype PAH systems (naphthyl radicals). For instance, rate coefficients for the phenyl + O2 reaction have been measured in two relatively low-temperature ranges of 297−47334 and 418−81535 K and computed using the canonical variational transition-state theory in conjunction with a density functional potential.36 At low temperatures, the peroxy radical C6H5OO formed via the addition/stabilization mechanism was determined to be the major product,34 whereas at temperatures above 1000 K the reaction was inferred to produce phenoxy radical, C6H5O + O, with possible minor contributions from o- and p-benzoquinone, o/p-C6H4O2, plus atomic hydrogen.37 Phenoxy radical can undergo unimolecular decomposition to the cyclopentadienyl radical and CO;34,38 thus, the oxidation with O2 essentially converts the six-member ring of phenyl to a five-member ring. The only products of the C6H5 + O2 reaction observed in crossed molecular beam experiments under single-collision conditions at collision energies of 5−26 kcal/mol were C6H5O and O.39−41 However, the most recent experimental study of the reaction in a combustion-simulating chemical reactor at 873 and 1003 K with product identification by tunable vacuum ultraviolet photoionization gave the phenoxy radical, cyclopentadienyl, and ortho-benzoquinone as major primary products, whereas secondary products included para-benzoquinone, phenol, cyclopentadiene, 2,4-cyclopentadienone, vinylacetylene, and acetylene.42 2,4-Cyclopentadienone, C5H4O, was suggested to originate from fragmentation of another hypothetical primary product, pyranyl radical C5H5O, which was not B

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

starting from molecular orbitals obtained from restricted openshell Hartree−Fock (ROHF) calculations. As in the previous studies of the phenyl + O2 and naphthyl + O2 reactions,23,24 wave functions of all species exhibit a rather mild multireference character; hence, the CCSD(T) energies are expected to be generally reliable. The B3LYP and MP2 calculations were performed using the Gaussian 09 program package,50 whereas the MOLPRO 2010 program package51 was employed for RHF-RCCSD(T) calculations. 2.2. Entrance Channel Reaction Rate Coefficients. Calculations of thermal rate coefficients for the barrierless C16H9 + O2, C15H9O|i1 + O, and C15H9|iN + O2, N = 2−5, association channels were carried out using variable reaction coordinate-transition state theory (VRC-TST) in the VariFlex code developed by Klippenstein and co-workers.52,53 The notation of “CnHm|iN” used above and throughout the text indicates Nth isomer of species CnHm, with isomer structures identified in the corresponding PES figures. A similar notation, “CnHm|tsN”, is used to identify specific transition states (TS). Briefly, the VRC-TST approach separates the transition state’s degrees of freedom into “conserved” and “transitional” modes. The conserved modes marginally change over the reaction process and are treated as rigid rotors or harmonic oscillators by directly computing the corresponding sums of states. The transitional modes are treated classically through a phase-space integral. Then, the partition function of a transition state is calculated through a convolution of these two types of modes, by varying the distance along the reaction coordinate and recalculating the contribution of the transitional modes to the partition function until a minimum reaction rate is found.52 In a study by Kislov et al.,29 the potential along the minimal energy reaction path (MEP) for the addition of O2 to phenyl (C6H5) was obtained by performing single-point energy calculations at the multireference second-order perturbationtheory CASPT2(19,14)/aug-cc-pVDZ level. For the larger systems in this study, the CASPT2(19,14) potential for phenyl + O2 was scaled to reproduce the accurate C16H9−O2 and C15H9−O2 bond dissociation energies in the C16H9O2 and C15H9O2 radicals. The MEP for the C15H9O + O barrierless channel was computed at the B3LYP/6-311G** level of theory and scaled to match the oxygen atom elimination energy calculated at the G3(MP2,CC) level. The computed stretching potentials for VariFlex calculations were provided in a tabulated form and then interpolated by the VariFlex built-in procedure. All density functional and ab initio calculations were carried out using the Gaussian 0950 and MOLPRO 201051 program packages. Because of an artificial restriction in the VariFlex program, the reaction rate could not be variationally calculated for bond lengths shorter than the distance from the center of mass for each fragment to the bonding point. This limitation is a problem for all six of the barrierless reactions in the present study because the C16H9 and C15H9 species are relatively large. For example, the distance from the center of mass of C15H9|i2 radical to the bonding oxygen atom is 3.0 Å. This means that if the true reaction rate minimum occurs for a bond length lower than 3.0 Å, then the calculated rate constant will be overestimated. To correct this, a smaller molecule (phenyl) was chosen as a surrogate of the larger one for the barrierless reactions. The bond dissociation energies, potentials, number of electronic states, and symmetry numbers for the surrogate reaction were altered to match the corresponding values of the actual reaction. For the C15H9|i2 + O2 reaction, the true

corannulenyl with O2, using these radicals as models of planar and curved PAHs, respectively. The calculations suggested a mechanism in which the six-member ring of the pyrenyl radical is oxidized by an oxygen molecule and is eventually converted into a five-member ring, similar to the reactions of phenyl and naphthyl radicals with molecular oxygen. The oxidation of a five-member ring with a second O2 molecule follows and results in the full destruction of the five-member ring with the formation of the phenanthryl radical. The oxidation mechanism of 1-corannulenyl was found to be similar, but with significantly lower activation energies, which explained the experimentally observed higher reactivity of curved PAH structures. However, the energetic parameters computed along the reaction pathways at the B3LYP//HF level may not be reliable, which may significantly influence the evaluation of rate coefficients for the oxidation reactions. To achieve better qualitative and quantitative understanding of oxidation processes of large PAH radicals, soot particles, and graphene edges, here we report a theoretical study of the reactions of pyrenyl, which is used as a model of large planar PAH radicals, with two oxygen molecules. First, we compute possible reaction pathways using high-level ab initio calculations that are capable of providing chemical accuracy of 1−2 kcal/mol for the energetics of various species including transition states. Then, the results of quantum calculations are utilized in RRKM-ME calculations of rate coefficients and relative product yields at temperatures and pressures relevant to combustion. Finally, the deduced oxidation mechanisms of sixand five-member rings and the rate coefficients are employed in KMC simulations of oxidation of a graphene “molecule” evolving in a flame-like environment.

2. COMPUTATIONAL METHODOLOGY 2.1. Potential Energies. The geometries of various species (reactants, products, intermediates, and transition states) involved in the C16H9 + O2 and C15H9 + O2 reactions, H2/ H-assisted isomerization of the C15H9 radical, and rearrangements of C14H9 were optimized using the hybrid density functional B3LYP method with the 6-311G** basis set.46 The same B3LYP/6-311G** approach was used to compute vibrational frequencies and zero-point energy (ZPE) corrections. All the stationary points were characterized as local minima or transition states according to the number of imaginary frequencies. The intrinsic reaction coordinate (IRC) calculations were performed to confirm all connections between transition states and local minima. Next, we applied the G3(MP2,CC)//B3LYP version47,48 of the original Gaussian 3 (G3) scheme49 for high-level singlepoint calculations to obtain more accurate energies of all of species. Within the G3(MP2,CC)//B3LYP scheme, final total energies at 0 K were computed at the B3LYP optimized geometries as follows: E[G3(MP2,CC)] = E[CCSD(T)/6‐311G**] + ΔE(MP2) + ZPE[B3LYP/6‐311G**]

where ΔE(MP2) = E[MP2/G3large] − E[MP2/6-311G(d,p)] is the basis set correction. In the CCSD(T) and MP2 calculations, the energies were computed from restricted RHF-RCCSD(T) and unrestricted UMP2 energies, respectively, where RHF-RCCSD(T) here denotes partially spinadapted open-shell coupled-cluster singles and doubles theory augmented with a perturbation correction for triple excitations C

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. KMC Oxidation Reactions

D

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. continued

E

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. continued

minimum reaction rate occurred at a bonding distance of 2.1 Å. Figure S1 in the Supporting Information illustrates the surrogate reaction technique used for the VRC calculations.

Next, the microcanonical rates for the surrogate reactions were calculated for the full range and for the restricted range of bonding distances. The ratio needed to correct the restricted F

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

2.4. KMC Simulations. The KMC methodology followed that used previously.22,61 Briefly, the KMC simulations tracked a single graphene “molecule” evolving in a flame-like environment, but maintaining the gas phase in a constant state, i.e., at constant values of temperature, pressure, and species mole fractions, which allowed us to investigate unambiguously the influence of individual factors on the simulation outcome. At each time step, a reaction event was selected stochastically and then applied. The simulation was treated as a series of Markov processes, i.e., each subsequent simulation step was dependent only on the current simulation state and not on the previous states. The selection of the reaction event and specific graphene-edge site was done by application of the Gillespie algorithm57,58 adapted for surface processes.62,63 Ensemble averages of substrate size and five-member ring fraction, f R5, were calculated for each set of simulations; f R5 was defined as in Whitesides and Frenklach22

range values was then calculated. This correction factor was then applied to the microcanonical rates of the original system. These scaled microcanonical rates for the barrierless reactions were used as inputs for the MultiWell simulations. In contrast to the four other C15H9 + O2 reactions, the C15H9|i1 + O2 reaction has an energy barrier of 1.6 kcal/mol. The bimolecular rate coefficient of this reaction was calculated in VariFlex as well as using TST and employing rigid rotor harmonic oscillator−tight transition state treatment in evaluating the partition function of the TS. 2.3. RRKM-ME Calculations. The rate coefficients of the chemically activated reaction systems were computed using the MultiWell 2014.1b code.54−56 Potential energy surfaces of C16H9 + O2 and C15H9|iN + O2 reaction systems, N = 1−5, are depicted in Figures 1−5. MultiWell solves the one-dimensional time-dependent energy-transfer master equations for a multiwell, multichannel unimolecular reaction system using the Monte Carlo stochastic method.57,58 Microcanonical rate coefficients for all elementary reactions studied here were computed with MultiWell at the RRKM level. Key parameters required for the MultiWell simulations, such as reaction barriers, frequencies, and moments of inertia, were taken from the quantum-chemical calculations of the present study. Following Gilbert and Smith,59 vibrational frequencies were examined by graphically visualizing the associated normal mode vibrations to identify internal rotational modes and, accordingly, all internal rotors were treated as one-dimensional hindered rotors. The microcanonical rate constants for the barrierless entrance reactions (computed by VariFlex as described in the previous section) were read into MultiWell and directly incorporated into the ME simulation. Reaction rates were computed at seven different temperatures (1000, 1250, 1500, 1750, 2000, 2250, and 2500 K) and five pressures (0.01, 0.1, 1.0, 10, and 100 atm). Argon was chosen as the bath gas collider. The exponential-down model with ⟨ΔE⟩down = 260 cm−1 was used to describe the collisional energy transfer.26 Lennard-Jones parameters were estimated from an empirical correlation.60 The exact count, with an energy grain size of 10 cm−1 for the first segment of the double array and a maximum energy of 500 000 cm−1, was employed to determine the density of states. For each set of initial conditions, the number of trials was varied to keep statistical error below 5%. The high-pressure-limit rate coefficients for product formation were calculated by using the high-pressurelimit rate coefficients of all well-to-well channels involved in the reactant-to-product reaction path and numerically solving a differential equation system of the corresponding species rate equations. The high-pressure-limit rate constant at a given temperature was obtained from the computed reactant and product profiles. Details and an example of these calculations are given in section 2 of the Supporting Information. The chemical-activated rate coefficients were derived from the accumulated species fractions of the products, as was reported previously.18 The accumulated species fractions were found by running a MultiWell simulation started with a chemical-activated initial energy distribution. The species fractions were found after the average energy of the initial adduct converges, indicating the end of the chemically activated simulation. Simulation times ranged from 10−9 s, at T = 2500 K and P = 100 atm, to 8 × 10−5 s at T = 1500 K and P = 0.01 atm. Product formation rates were obtained by multiplying the product fractions by the high-pressure reaction rate coefficient of the adduct formation reaction.

fR5 =

NR5 32 12 NR5 + NR6

(1)

where NR5 and NR6 are the number of five- and six-member rings, respectively, and the coefficients normalize eq 1 to have f R5 = 1 for the Buckminster fullerene, C60. Simulations were run with an ensemble of 100 individual runs, each with a different starting random seed, thus reducing the stochastic fluctuations in the averaged data. For each set of conditions, the stochastic error (one standard deviation) did not exceed 5% for the number of carbon atoms and 10% for the five-member ring fraction, the level of accuracy sufficient for the objectives of an exploratory analysis of the present work. An example plot of the results with error bars is presented in section 3 of the Supporting Information. To quantify the reaction rates, the number of reaction events was counted and averaged over the ensemble of simulations. The averaged reaction counts are reported in Tables S15−S32 of the Supporting Information. The set of surface reactions employed in the KMC model contained 45 growth and 45 oxidation reactions. The surface growth reactions and their corresponding rate coefficients can be found in Whitesides and Frenklach.22,61 The surface oxidation reactions and their corresponding rate coefficients are specified in Table 1. Many of these reactions are simple O, H, and OH abstractions and additions, and their rate coefficients were assigned by analogy to similar benzene-ring reactions. The rate coefficients of the ring-oxidation steps were taken from recent studies, and that of reaction 90, oxidation of a five-member ring on a boat edge, is the one calculated in the present study for the C15H9|i1 + O2 → C15H9O|i10 + O reaction. The expressions listed in Table 1 are per-site rates, i.e., the rate coefficients referring to individual sites of a graphene edge.

3. RESULTS AND DISCUSSION 3.1. Potential Energy Surfaces. High-level ab initio calculations were carried out using the pyrene molecule as a model system for which we considered oxidation processes with two oxygen molecules. For all of the molecular structures discussed below, Table S1 in the Supporting Information contains the vibrational frequencies, rotational constants, unsymmetrical hindered rotor parameters, and the Cartesian coordinates for all optimized structures. Because the reactions of PAH molecules with O2 are normally preceded by formation of their radicals via direct H abstraction by other radicals G

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. Potential energy diagram for oxidation of pyrenyl radical by O2. Potential energies calculated at the G3(MP2,CC) level are shown in kilocalories per mole relative to the C16H9|i1 + O2.

Figure 2. Potential energy diagram for H2/H assisted isomerizations of C15H9 radicals. Potential energies calculated at the G3(MP2,CC) level are shown in kilocalories per mole relative to C15H9|i1.

only the reaction channels which were proved to be dominant, at least at combustion temperatures, for the reactions of phenyl and naphthyl radicals with O2,29 namely, addition of molecular oxygen to the radical site followed by elimination of an oxygen atom and consequent unimolecular decomposition of the remaining oxypyrenyl radical. The reaction channels involving

present in the system or by the pyrolitic cleavage of a C−H bond, we begin our study with the reaction of the pyrenyl radical with O2. 3.1.1. C16H9 (Pyrenyl) + O2. The potential energy diagram for the reaction of the pyrenyl with molecular oxygen is depicted in Figure 1. It should be noted that here we consider H

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 3. Potential energy diagram for oxidation of C15H9 radicals by O2: the C15H9 + O2 → C15H9OO → C15H9O + O pathways. Potential energies calculated at the G3(MP2,CC) level are shown in kilocalories per mole relative to C15H9|i1 + O2.

of only 5.8 kcal/mol relative to i2. The first step appears to be rate-determining. The decomposition reaction of the oxypyrenyl radical is 18.8 kcal/mol endothermic, and in the product C15H9|i1 one six-member ring of the original pyrenyl radical is replaced with a five-member ring. The further fate of the radical C15H9|i1 may be two-fold; it may either be oxidized by a second oxygen molecule or convert to a stable closed-shell molecule by recombining with a hydrogen atom or by abstracting H from H2 or any other molecular source available. The H addition or H abstraction reaction results in C15H10|i1 (Figure 2), a PAH molecule with three six-member rings and one five-member ring. The weakest C−H bond in C15H10|i1 is at the CH2 group in the five-member ring (80.9 kcal/mol), whereas the strengths of C−H bonds in the six-member rings are in the range of 111−112 kcal/mol. The i1 isomer is the most stable among the considered PAH C15H9 radicals, by ∼30 kcal/mol. The C15H9|i1 + H2 → C15H10| i1 + H reaction has a barrier of 25.2 kcal/mol and is 20.0 kcal/ mol endothermic, whereas for the other C15H9 isomers similar reactions exhibit higher barriers of ∼38 kcal/mol and are ∼10 kcal/mol exothermic. In the reverse direction, the direct H abstraction from C15H10 by a hydrogen atom has barriers of 5.0 kcal/mol to form C15H9|i1 and 17−18 kcal/mol to produce the other isomers. Even if the i1 isomer of C15H9 is stabilized by H atom addition, the activation of the C10H15 species will preferably occur at the five-member ring by H abstraction from the CH2 group. Therefore, the C15H9|i1 radical will most likely undergo further oxidation by reacting with molecular oxygen.

insertion of one of the oxygen atoms into the aromatic ring followed by decomposition of the resulting intermediate with the PAH•-O2 stoichiometry (here, C16H9O2•) are ignored, keeping in mind that they may be responsible for up to 15− 24% of the reaction products in the temperature range of 1500−2500 K. The addition of O2 to pyrenyl is calculated to occur without a barrier and to form the peroxy pyrenyl radical C16H9O2|i1 with the exothermicity of 46.2 kcal/mol (Figure 1), which is similar to those computed earlier for peroxy phenyl and naphthyl radicals. The next reaction step is emission of an oxygen atom by cleaving the O−O bond. Just like for C6H5OO and C10H7OO, this process goes via the O−O bond-cleavage transition state, C16H9O2|ts1, to a weakly bound van der Waals C16H9O···O complex, which decomposes to the oxypyrenyl radical, C16H9O|i1, without a barrier. The bond-cleavage transition state resides 1.7 kcal/mol below the final product but nevertheless represents a kinetic bottleneck on the path of the O atom emission, at least at combustion temperatures, and will therefore be used in rate-constant calculations. Overall, the C16H9 (pyrenyl) + O2 → C16H9O (oxypyrenyl) + O reaction is predicted to be exothermic by 17.0 kcal/mol. The most favorable mechanism for unimolecular decomposition of the oxypyrenyl radical is elimination of the CO molecule from the oxidized ring resulting in the contraction of this six-member ring to a five-member ring. The reaction occurs via two steps: the first one is the ring contraction via C16H9O|ts1 and a barrier of 60.4 kcal/mol leading to i2, which has an out-of-ring CO group, and the second is emission of the CO group via a barrier I

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 4. Potential energy diagram for oxidation of C15H9|i1 by O2: the C15H9 + O2 → C15H9OO → C14H9 + CO2 pathways. Potential energies calculated at the G3(MP2,CC) level are shown in kilocalories per mole relative to C15H9|i1 + O2.

highest barriers on the pathways of CO elimination are ∼70 kcal/mol for the i4 and i5 pathways, 80.6 kcal/mol for i2, and 86 kcal/mol for i3. The i2 isomer of C15H9O dissociates to the i1 isomer of C14H9 with the endothermicity of 44.6 kcal/mol, whereas the i3−i5 isomers of C15H9O decompose to the same C14H9|i2 with endothermicities of 35.7−38.7 kcal/mol. Both C14H9 isomers have two five-member rings next to one another and share one C−C bond; in i1, the five-member rings are located in the top and bottom positions of the molecule, while in i2, one of them is in the bottom and the other in the side location. The i2 isomer of C14H9 is 4.4 kcal/mol more stable than i2. We turn our attention to the oxidation reaction of C15H9|i1 in which the radical site is delocalized on the five-member ring instead of being localized at one of the six-member ring carbon atoms. The reaction mechanism appears to be rather different from those for C15H9|i2−i5 (Figure 3). First, the initial O2 addition is not barrierless and proceeds via a transition state located 1.6 kcal/mol above the reactants. The resulting complex C15H9O2|i1 is located only 17.5 kcal/mol lower in energy than the reactants; therefore, the new C−O bond is ∼30 kcal/mol weaker than those in C15H9O2|i2−i5. However, the strength of the O−O bond in i1, 58.3 kcal/mol, is much higher than those in i2−i5, 27−33 kcal/mol. The O−O bond in C15H9O2|i1 breaks without an exit barrier; a van der Waals productlike complex was also not found. Overall, the C15H9|i1 + O2 → C15H9O|i1 + O reaction was found to be endothermic by 40.8 kcal/mol, in contrast to the exothermic reactions of O2 addition and O elimination involving the other C15H9 isomers where the oxidation occurs

3.1.2. C15H9 + O2. Although the i2−i5 isomers of the C15H9 are not likely to be produced in the pyrenyl + O2 reaction or via subsequent H/H2 rearrangements of the C15H9 radical, for completeness, we briefly consider their reactions with O2. Again, we focus only on the favorable O2 addition−O elimination channel followed by emission of the CO group leading to the conversion of a second six-member ring of the original pyrenyl radical to a five-member ring (Figure 3). The O2 addition reactions to a radical site in six-member rings to form peroxy-type radicals occur without barriers with energy gains of 44−47 kcal/mol. Next, the O−O bonds break, producing oxyradicals via transition states lying slightly below, by 1.6−3.4 kcal/mol, the oxyproducts and product-like weakly bound van der Waals complexes. The reactions C15H9|i2−i5 + O2 → C15H9O|i2−i5 + O are predicted to be exothermic by 12.1−16.9 kcal/mol. The oxyradicals C15H9O|iN, N = 2−5, can undergo unimolecular decomposition by splitting CO, resulting in the formation of C14H9 radical products containing two sixmember rings and two five-member rings. These processes occur by expulsion of the CO group from the six-member ring, contraction of this ring to a five-member cycle, and elimination of CO, and may go via one, two, or three steps: C15H9O|i2 → ts4 → i8 → ts5 → C14H9|i1 + CO, C15H9O|i3 → ts6 → C14H9| i2 + CO, C15H9O|i4 → ts7 → i9 → ts8 → i10 → ts9 → C14H9| i2 + CO, and C15H9O|i5 → ts10 → i11 → ts11 → C14H9|i2 + CO. In the two- and three-step processes, the intermediates, i8− i11, are metastable and reside in shallow wells on the reaction pathways, which are stabilized by low 5−7 kcal/mol barriers at least in one direction. Relative to the C15H9O oxyradicals, the J

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 5. Potential energy diagram for isomerizations of C14H9 radicals. Potential energies calculated at the G3(MP2,CC) level are shown in kilocalories per mole relative to C14H9|i1.

Along the first pathway, the C−C(H)O bond is cleaved first (i4 → i5), then the hydrogen atom of the C(H)O group migrates to the radical site of the six-member ring (i5 → i6); finally, i6 splits the CO2 group, producing the phenanthryl radical with the overall exothermicity of 76.3 kcal/mol. In the second channel, first the H atom of C(H)O in i4 shifts to the out-of-ring oxygen (i4 → i7), then a C−O bond breaks in the C5O ring by rotation of the C(O)OH moiety (i7 → i8); the H atom of the hydroxy group migrates to the radical site of the six-member ring (i8 → i9), and the CO2 group is finally eliminated by a C−C bond rupture. The highest in energy transition states along these two channels, ts6 and ts5, are located 23.9 and 21.1 kcal/mol above i4, respectively; hence, these channels are likely to be competitive. Meanwhile, ts6 and t5 are significantly lower in energy than the initial reactants and C15H9O2|i1 and the highest in energy transition state on the overall reaction pathway is ts2 between i2 and i3, residing 21.8 kcal/mol above the C15H9|i1 + O2 reactants. Thus, the highest barrier from C15H9O2|i1 that needs to be overcome to produce phenanthryl + CO2 is 39.3 kcal/mol; this is 19.0 kcal/mol lower than the energy required to split the O−O bond forming C15H9O|i1. Therefore, the C15H9|i1 + O2 → C15H9O2|i1 → ··· → phenanthryl + CO2 reaction mechanism should be clearly preferable over C15H9|i1 + O2 → C15H9O2|i1 → C15H9O|i1 + O → ··· → phenanthryl + CO + O, both thermodynamically and kinetically. Summarizing, the reaction of the pyrenyl radical or of a pyrenyl radical-like structure in a larger PAH or a graphene sheet with two oxygen molecules is most likely to produce the

on a six-member ring. The further fate of C15H9O|i1 is elimination of the CO group occurring by the following mechanism: hydrogen migration from the C atom linked to O, i1 → ts1 → i6; cleavage of the C(O)−C(H) bond, i6 → ts2 → i7; and splitting of the CO group, i7 → ts3 → C14H9 (phenanthryl). This process is 5.8 kcal/mol endothermic overall, and the highest in energy transition state resides 18.6 kcal/mol above C15H9O|i1. The product, phenanthryl radical, is 11.0 and 6.6 kcal/mol more stable that the i1 and i2 isomers of the C14H9 radical, respectively. However, oxygen atom elimination from C15H9O2|i1 is not the most favorable channel for isomerization−fragmentation of this radical. Figure 4 illustrates the alternative channels first proposed by Raj et al.45 and recalculated here at the higher level of theory. At the first reaction step, the exo-O atom forms a bond with the neighboring carbon via a barrier of 32.4 kcal/ mol, leading to the C15H9O2 intermediate i2, 24.6 kcal/mol above i1; i2 features a four-member C2O2 ring. Next, the O−O and one of the C−C bonds break by rotation of the C(H)O group around the remaining C−C(H)O bonds, overcoming a 14.7 kcal/mol barrier in a highly exothermic process, which results in the intermediate i3 residing 71.1 kcal/mol lower in energy than the C15H9|i1 + O2 reactants. Then, a new C−O bond forms between the former exo-oxygen atom and C of the C(H)O group, producing a six-member C5O ring in i4 via a 32.5 kcal/mol barrier. The i4 intermediate can develop via two alternative channels, both leading to the same phenanthryl + CO2 products. K

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 2. Total Bimolecular Rate Coefficients of Product Formation (cm3 mol−1 s−1) T (K)

C16H9 + O2

1000 1500 2000 2500

1.74 3.27 4.78 6.26

× × × ×

1013 1013 1013 1013

C15H9|i1 + O2 4.10 1.21 2.46 4.18

× × × ×

1011 1012 1012 1012

C15H9|i2 + O2 1.68 3.08 4.46 5.75

× × × ×

1013 1013 1013 1013

phenathryl radical, C16H9 + 2O2 → C14H9 (phenanthryl) + CO + O + CO2, or destroy one six-member ring and form an armchair edge with a radical site on it. A C15H9-like radical with one five-member ring resulting from the reaction with a first O2 molecule would undergo oxidation with a second oxygen molecule attacking the five-member ring, which completely erases the five-member ring, rather than be subjected to an H2/ H assisted isomerization moving the radical site to a sixmember ring. Therefore, two consequent oxidation reactions are not expected to produce relatively unstable structures with two five-member rings next to each other; indeed, such structures are not usually observed in experiment. 3.1.3. Isomerization of C14H9 Radicals. What if the i1 and i2 isomers of C14H9 with adjacent five-member rings are formed somehow? How stable would they be with respect to isomerization to the more stable phenanthryl isomer? The potential energy diagram for this process is shown in Figure 5. One can see that C14H9|i1 can rearrange to i2 via a two-step process (via i4) featuring the highest barrier of 79.9 kcal/mol relative to i1. Alternatively, i1 can lose a hydrogen atom from the CH2 group in the five-member ring, forming a triplet C14H8 molecule without an exit barrier; this fragmentation process requires 84.7 kcal/mol of energy. In turn, C14H9|i2 can rearrange to phenanthryl via a four-step mechanism including a rupture of a C−C bond to the CH2, opening the five-member ring, migration of an H atom from the CH2 to the radical site in the six-member ring, and expansion of the five-member ring to a six-member ring by incorporation of the out-of-ring CH group occurring in two steps via a [5,3]-bycyclic intermediate. The highest in energy transition state on this pathway is on the third step, a three-member ring closure via ts6, and the barrier is 56.2 kcal/mol relative to i2. Alternatively, i2 can simply lose a hydrogen atom from the CH2 group forming a singlet C14H8|i2 molecule via a slightly lower barrier of 53.8 kcal/mol. Another isomerization process pushing two five-member rings in i2 away from each other, to the i3 or i10 intermediates and eventually to the C14H8|i3 singlet product, exhibits much higher barriers of 79−81 kcal/mol and is not likely to be competitive. Thus, the C14H9|i1 → i2 → phenanthryl isomerization is likely to be overcompeted by H losses from i1 and i2 forming the C14H8 i1 and i2 products, respectively, because the activation energies for the H loss are similar or even lower than the highest barriers on the isomerization pathways and the H elimination TSs are loose. On the other hand, i2 and phenanthryl can be produced from the C14H8|i1 + H and C14H8|i2 + H reactions, respectively. However, the calculations show that both i1 and i2 isomers of C14H9 containing two neighboring five-member rings are quite stable kinetically and if formed may survive a rather long time. Therefore, the nonobservation of two five-member rings next to each other in large PAH as well as in soot and graphene particles can be attributed to the fact that the formation channels of such structures are not competitive (as shown for the oxidation reactions in section 3.1.2) rather than to their kinetic instability.

C15H9|i3 + O2 1.66 3.04 4.44 5.76

× × × ×

1013 1013 1013 1013

C15H9|i4 + O2 1.80 3.36 4.92 6.43

× × × ×

1013 1013 1013 1013

C15H9|i5 + O2 1.66 3.09 4.55 5.93

× × × ×

1013 1013 1013 1013

3.2. Reaction Rate Coefficients. 3.2.1. Overall Thermal Rate Coefficients: C16H9 + O2 and C15H9 (i1−i5) + O2. Rate coefficients for the C16H9 + O2 and C15H9 (i1−i5) + O2 reactions computed at the combustion-relevant temperatures of 1000−2500 K and P = 1 atm are collected in Table 2. The data for the other pressures (0.01, 0.1, 10, and 100 atm) as well as the high-pressure-limit rate coefficients of all elementary reactions are given in the Supporting Information, Tables S2−S14. The results demonstrate that ktotal(C16H9 + O2) and ktotal(C15H9|i2−i5 + O2) are in the range of 1.7−6.4 × 1013 cm3 mol−1 s−1 and are close to each other at the same temperatures, with differences not exceeding 10%. In contrast, ktotal(C15H9|i1 + O2) is computed to be over an order of magnitude lower, in the range of 4.1 × 1011 −4.2 × 1012 cm3 mol−1 s−1. Raj et al.27,45 investigated the C16H9 + O2 and C15H9|i1 + O2 reaction pathways using density functional theory at the B3LYP/6-311++G(d,p) level of theory and performed rate calculations using TST. For ambiguous transition states, they conducted intrinsic reaction coordinate calculations. The overall thermal rate coefficient in the high-pressure limit for the C16H9 + O2 reaction was in the range of 5.2−9.5 × 1011 cm3 mol−1 s−1, nearly 2 orders of magnitude lower than the corresponding values computed in the present study. This discrepancy is also seen for reaction C15H9|i1 + O2; the values 1.1 × 1010−2.5 × 1011 reported by Raj et al.27 are an order of magnitude lower than those calculated in the present study. The difference in the rate coefficients can be explained by the differences in the potential energies between C15H9 and O2 fragments and the C15H9O2 adduct calculated by each study. In Raj et al.,27,45 the potential energy difference is 9.8 kcal/mol, while our result is 17.5 kcal/mol. 3.2.2. Product Branching Ratios and Individual Rate Coefficients for Product Formation. The reaction of C16H9 + O2 produces C16H9O + O. The C16H9O can then dissociate into C15H9 + CO. Similar to the phenyl + O2 study,29 collisionally stabilized intermediates were not detected as products in the temperature and pressure ranges studied. Figure 6 exhibits the pressure dependence for the rate coefficient of reaction C16H9 + O2 → C16H9O + O, which becomes significant at temperatures above 1500 K. For instance, at 2500 K and a pressure of 0.01 atm, the rate coefficient of formation of C16H9O + O is over 2 orders of magnitude lower than the high-pressure limit. Figure 7 compares the high-pressure-limit rate coefficient of reaction C16H9 + O2 → C16H9O + O calculated in the present study to the high-pressure-limit rate coefficients of reactions C6H5 + O2 → C6H5O + O, 1-C10H7 + O2 → 1-C10H7O + O, and 2-C10H7 + O2 → 2-C10H7O + O calculated by Kislov et al.29 As can be seen from the displayed results, the reactions forming C6H5O + O and 2-C10H7O + O are 100 and 40% faster, respectively, than the reaction forming C16H9O + O. The rate coefficient of the reaction forming 1-C10H7O + O is within 10% of that forming C16H9O + O. As the number of rings L

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 6. Computed rate coefficient for reaction C16H9 + O2 → C16H9O + O.

Figure 8. Computed rate coefficients for reactions C15H9|i1 + O2 → C14H9+CO2 and C15H9|i1 + O2 → C15H9O+O.

Figure 7. Comparison of high-pressure-limit rate coefficients computed for reactions C6H5 + O2 → C6H5O + O, 1-C10H7 + O2 → 1-C10H7O + O, 2-C10H7 + O2 → 2-C10H7O + O, and C16H9 + O2 → C16H9O + O.

Figure 9. Computed product branching ratios of the C15H9|i1 + O2 reaction.

increases from one to four, the rate of oxidation of the sixmember ring decreases by a factor of 2. A comparison of the high-pressure-limit rate coefficient for the oxidation of phenyl to experimental results can be found in Kislov et al.29 The computed rate coefficients in that study were found to agree with measured rate coefficients from experimental work of Schaugg and co-workers.35 Reaction C15H9|i1 + O2, by which the five-member ring is oxidized, has two main pathways, those producing C15H9O + O and C14H9 + CO2. The C15H9O radical can further dissociate to C14H9 + CO. The computed rate coefficients of these reaction channels exhibited no dependence on pressure and are displayed in Figure 8. As can be seen from these results, the C14H9 + CO2 formation is much faster that the formation of C15H9O + O for all temperatures considered in the present study, with the difference between the two increasing with the decrease in temperature. A closer look at these results, displayed as a branching ratio in Figure 9, indicates that the CO2-producing channel is the dominant one at lower and combustion temperatures. This conclusion is consistent with the isolated pentagon rule30−33 that prohibits formation of two adjacent five-member rings. In Figure 10 we compare our present values with those of Raj et al.27,45 for the high-pressure-limit rate coefficients computed for the oxidation of six- and five-member rings, C16H9 + O2 → C16H9O + O and C15H9|i1 + O2 → C14H9 + CO2, respectively. Raj et al.27,45 published the rate coefficients only for the elementary, well-to-well steps; we used those values to calculate

Figure 10. Comparison of high-pressure-limit rate coefficients computed in the present study to values reported by Raj et al.27

the high-pressure-limit rate coefficients of product formation similarly to our own, as described in section 2.3. Our values are about 2 orders of magnitude higher than those of Raj et al. for the corresponding reactions. Both studies, however, are in agreement that the oxidation of a six-member ring is over 4 orders of magnitude faster than the oxidation of a five-member ring. Further details, including computational results for additional reaction, such as oxidation of C15H9|i2, C15H9|i3, C15H9|i4, and C15H9|i5, are reported in the Supporting Information. 3.3. KMC Simulations. KMC simulations were performed at 1500, 2000, and 2500 K. Two scenarios were simulated: one with pure oxidation and one in which oxidation competes with M

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 11. Substrate size (left-hand panels) and five-member ring fraction (right-hand panels) for the oxidation-and-growth scenario at (a) 1500 K, (b) 2000 K, and (c) 2500 K.

growth. In the latter scenario, the graphene layer was first grown from an initial coronene substrate for 5 ms and then O2 was added to the gaseous environment to simulate oxidation that accompanies growth. The temperature was held constant in each simulation run. The pressure was held constant at 1 atm, and the gas-phase composition was held constant with mole fractions xC2H2 = xH2 = 0.1 and xH = 0.01. After O2 was added, its concentration was also held constant. The pure oxidation scenario was similar to the oxidation-and-growth scenario, except that when O2 was added after the initial 5 ms of growth, acetylene was removed, i.e., xC2H2 was set to 0. The two scenarios were designed to cover conditions encountered in hydrocarbon flames.

3.3.1. Oxidation-and-Growth Scenario. Time evolution of the substrate size and five-member ring fraction for a range of oxygen concentrations are depicted in Figure 11 for each temperature. The results computed for the substrate size, shown in the left-hand panels of Figure 11, demonstrate that as oxygen is added, oxidation begins to compete with growth, increasingly so with an increase in the amount of oxygen added and an increase in temperature. For instance, at the same level of O2 added (e.g., xO2 = 1 × 10−3), the growth of the graphene layer is only inhibited at 1500 and 2000 K, but at 2500 K, the layer decreases in size. The computed five-member ring fractions, depicted in the right-hand panels of Figure 11, reveal that the addition of O2 promotes the formation of five-member rings embedded in the N

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 12. Diagram of two major pathways: thermal decomposition of oxyradicals (green) and regeneration of an aromatic radical site (red).

graphene layer. Analysis of the computed reaction fluxes indicated that the increase in the five-member rings is primarily due to their formation via thermal decomposition of oxyradicals (reactions 51−54). Excluding oxidation of five-member from the model results in faster saturation of the graphene edge with five-member rings, which prevents the graphene edge from further degradation (or growth). The KMC results obtained in this oxidation-and-growth scenario revealed two major, competing pathways for an oxyradical originating in O2 reaction with a surface radical. In the first pathways, the oxyradical undergoes thermal decomposition to form a five-member ring and expel COthe actual oxidation step. In the second pathway, a neighboring-site H adds to the oxyradical forming OH, to which another H atom adds to expel H2O and regenerate a new aromatic radical site a pathway recycling the aromatic-radical site. Figure 12 shows a schematic diagram of the two pathways, and Figure 13 quantifies the competition between these two pathways by comparing numbers of reaction events for the attack of O2 on a radical site (reactions 46−50), thermal decomposition of an oxyradical (reactions 51−54), and regeneration of an aromatic radical site (reactions 85−89). For all three temperatures, the number of oxyradicals formed increases with an increase in oxygen concentration. At 1500 K, the regeneration pathway is dominant over the thermal-decomposition pathway. As the temperature increases to 2000 K, thermal decomposition begins to compete with regeneration, and at 2500 K thermal decomposition is the dominant pathway. This switch is a result of thermal decomposition being a higher activation-energy process than regeneration. Carbon is removed from the graphene edge either by thermal desorption (reactions 5, 7, 17, 25, and 45) or by oxidation (reactions 51−54). For all three temperatures studied, the reaction statistics show that as the concentration of oxygen increases, the frequency of oxidation reactions increases while the frequency of thermal desorption reactions decreases, indicating that competition exists between the two types of reactions. Although the frequency of thermal desorption decreased with an increasing oxygen concentration, the thermal desorption pathway remained dominant over the oxidation pathway for all oxygen concentrations and temperatures tested. This result indicates that when growth reactions (3, 4, 26, 27, 33, 34, 35, 39, 42, and 44) added carbon atoms to the graphene layer, they were more likely to be removed by thermal desorption than oxidation, even for the highest oxygen concentrations. Given a particular interest of the present study in the oxidation of five-member rings, we performed a series of KMC runs to investigate the frequency of reaction 90, C15H9|i1 + O2

Figure 13. Reaction event counts for O2 attack on a radical site (blue), thermal decomposition of an oxyradical (green), and regeneration of an aromatic radical (red) for the oxidation-and-growth scenario.

→ C15H9O|i10 + O. Reaction statistics revealed that reaction 90 occurs rather infrequently. For instance, at 1500 K and xO2 = 0.1, reaction 90 occurred 1.0 times, on average, per simulation, as compared to 23 and 13 counts for reactions 5 and 7, respectively. For lower O2 concentrations, it did not happen. At 2000 K, reaction 90 occurred only 0.1 times per simulation with the highest O2 concentration, and at 2500 K it did not occur at O

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 14. Substrate size (left-hand panels) and five-member ring fraction (right-hand panels) for the oxidation-only scenario at (a) 1500 K, (b) 2000 K, and (c) 2500 K.

decomposition and regeneration, as shown in Figure 12 and exemplified by reaction-event counts in Figure 15, as well as the relative frequency of reaction 90. In the oxidation-only simulations, the frequency of thermaldesorption reactions was over an order of magnitude lower than in the oxidation-and-growth cases because acetylene was not present in the gaseous environment during oxidation; hence, there was a lower occurrence on the graphene edge of lone adsorbates able to desorb. Unlike for the oxidation-andgrowth cases, at 1500 K the oxidation pathway became dominant over thermal desorption for the highest oxygen concentration, xO2 = 0.1. At 2000 K, oxidation dominated thermal desorption for xO2 > 10−4. Still, similar to the oxidationand-growth cases, at 2500 K, the thermal desorption reactions were dominant over the oxidation reactions for all concentrations of oxygen studied.

all. In a repeat of the KMC simulations using the rate coefficients of Raj et al.,27 which are over 2 orders of magnitude lower than ours, reaction 90 did not occur for any of the simulations. 3.3.2. Oxidation-Only Scenario. This set of KMC simulations was performed at exactly the same conditions as the previous one, except that acetylene was removed from the gaseous environment at the time oxygen was added. Figure 14 displays the time evolution of the substrate size and fivemember ring fraction for each temperature and for a range of oxygen concentrations. There is no growth to compete with oxidation in these simulations, so the addition of any amount of oxygen leads to a decrease in the size of the graphene layer. Other than that, the results displayed in Figure 14 are similar to those of the discussed above oxidation-and-growth case. Also similar are the major competition pathways of thermal P

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 16. Representative structures seen in the oxidation-only scenario. The displayed snapshots are from a KMC simulation at 2000 K and xO2 = 0.001: (a) at the end of the growth period and just before the onset of oxidation and (b) after 1.8 ms of the oxidation.

seen in Figure 14, the profiles of the substrate size during the oxidation are nonlinear in time, implying time-varying decayingrates of oxidation. Because the gaseous environment is maintained unchanged during an individual simulation, the only varying property affecting the rate of oxidation is the edge density of reactive sites, which in this case is the fraction of graphene-edge sites able to be removed in reaction with O2. Yet, the oxidation reactions themselves are forming fivemember rings that become embedded in the graphene edge, thereby reducing the edge reactivity and hence the rate of oxidation with reaction time. The variation of the total oxidation rate with time makes its comparison to experimental data on soot and carbon oxidation ambiguous. Furthermore, variation in the gaseous environment may affect the oxidation rate. Thus far, we can say that the initial oxidation rates for the cases presented in Figure 14 are in close agreement with the expression of Nagle−Strickland− Constable67 (NSC) at the lower end of the temperature range and the upper end of the O2 concentrations examined in the present study; at higher temperatures and lower O 2 concentrations, the initial oxidation rates exceed the NSC values. These aspects will be addressed in a future publication.

Figure 15. Reaction event counts for O2 attack on a radical site (blue), thermal decomposition of an oxyradical (green), and regeneration of an aromatic radical site (red) for the oxidation-only scenario.

Similarly to the oxidation-and-growth case, the graphene edge became “non-reactive” quicker when the five-member-ring oxidation was excluded from the simulation. Such “nonreactive” sites are illustrated in Figure 16, which displays two snapshots from an oxidation-only KMC simulation at 2000 K and xO2 = 0.001. The snapshot shown on the left-hand side of the figure is taken at a simulation time of 5 ms, at the instant when C2H2 is removed from and O2 is added to the gaseous environment. The green box exemplifies that before the oxidizer is added, the edge consists of six-member rings with a “free” corner, thus enabling a zipper oxidation.25,26 The righthand side of the figure depicts a structure formed after 1.8 ms of the oxidation period. The red pentagons drawn over the structure highlight some of the five-member rings of the graphene edge. If not removed, the highlighted five-member rings prevent further oxidation from occurring. As seen above, even when the five-member-ring oxidation reaction is included, its rate is relatively low and hence may not prevent the buildup of five-member rings. Accumulation of five-member rings at the graphene edge leads to its reduced reactivity. This is manifested by the decaying rates of graphene-edge oxidation. Indeed, as can be

4. SUMMARY Quantum-chemical calculations were performed for a system of pyrenyl oxidation. On the basis of these results, rate coefficients and branching ratios were calculated for elementary reactions in the system. The latter, in turn, were utilized in modeling of the system kinetics in flame-like environments. The major findings are summarized below, starting with addressing the main questions, posed in the Introduction, that motivated the study. The initial step of pyrenyl oxidation by molecular oxygen

is basically similar, energetically and kinetically, to the oxidation reactions of naphthyl radicals.29 However, a subsequent oxidation of the inner six-member ring of the product does not proceed as depicted in the Introduction by Step IIa, but rather as Step IIb Q

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231. A.M.M. acknowledges the Instructional & Research Computing Center (IRCC, web: http://ircc.fiu. edu) at Florida International University for providing HPC computing resources that have contributed to the research results reported within this paper.

thus conforming to the isolated pentagon rule30−33 that forbids formation of two adjacent pentagon rings next to each other. Oxidation of six-member rings, like in Step I, leads to accumulation of embedded five-member rings in the reacting graphene layer, which causes, in turn, its reduced reactivity. Oxidation of the embedded five-member rings by O2 alleviates this problem to some extent, especially at lower temperatures, but not entirely, because its rate coefficient was computed to be substantially lower than that of six-member rings. The rate coefficients for the oxidation of both six- and fivemember rings by O2 were evaluated to be significantly, roughly by 2 orders of magnitude, faster than those of Raj et al.27 The oxidation system exhibits two basic pathways: thermal decomposition and regeneration of oxyradicals. Their competition is temperature dependent, with the former dominating at higher and the latter at lower temperatures. The overall oxidation of the graphene substrate is computed to be time-dependent, with the initial rates consistent with the NSC67 expression.

■

■

(1) Glassman, I.; Yetter, R. A. Combustion; Academic: San Diego, CA, 2008. (2) Lewtas, J. Air Pollution Combustion Emissions: Characterization of Causative Agents and Mechanisms Associated with Cancer, Reproductive, and Cardiovascular Effects. Mutat. Res., Rev. Mutat. Res. 2007, 636, 95−133. (3) Dedoussi, I. C.; Barrett, S. R. H. Air Pollution and Early Deaths in the United States. Part II: Attribution of PM2.5 Exposure to Emissions Species, Time, Location and Sector. Atmos. Environ. 2014, 99, 610− 617. (4) Shoemaker, J. K.; Schrag, D. P.; Molina, M. J.; Ramanathan, V. What Role for Short-Lived Climate Pollutants in Mitigation Policy? Science 2013, 342, 1323−1324. (5) Pierrehumbert, R. T. Short-Lived Climate Pollution. Annu. Rev. Earth Planet. Sci. 2014, 42, 341−379. (6) Frenklach, M. Reaction Mechanism of Soot Formation in Flames. Phys. Chem. Chem. Phys. 2002, 4, 2028−2037. (7) Wang, H. Formation of Nascent Soot and Other CondensedPhase Materials in Flames. Proc. Combust. Inst. 2011, 33, 41−67. (8) Palmer, H. B.; Cullis, C. F. The Formation of Carbon from Gases. In Chemistry and Physics of Carbon; Walker, P. L., Ed.; Marcel Dekker: New York, 1965; Vol. 1, pp 265−325. (9) Haynes, B. S.; Wagner, H. G. Soot Formation. Prog. Energy Combust. Sci. 1981, 7, 229−273. (10) Kroto, H. W.; Heath, J. R.; Obrien, S. C.; Curl, R. F.; Smalley, R. E. C60: Buckminsterfullerene. Nature 1985, 318, 162−163. (11) Zhang, Q. L.; Obrien, S. C.; Heath, J. R.; Liu, Y.; Curl, R. F.; Kroto, H. W.; Smalley, R. E. Reactivity of Large Carbon Clusters Spheroidal Carbon Shells and Their Possible Relevance to the Formation and Morphology of Soot. J. Phys. Chem. 1986, 90, 525− 528. (12) Frenklach, M.; Ebert, L. B. Comment on the Proposed Role of Spheroidal Carbon Clusters in Soot Formation. J. Phys. Chem. 1988, 92, 561−563. (13) Donnet, J.-B. Fifty Years of Research and Progress on Carbon Black. Carbon 1994, 32, 1305−1310. (14) Frenklach, M. A Unifying Picture of Gas Phase Formation and Growth of PAH, Soot, Diamond and Graphite. In Carbon in the Galaxy: Studies From Earth and Space; Tarter, J. C., Chang, S., DeFrees, D. J., Eds.; NASA Conference Publication 3061: 1990; pp 259−273. (15) Frenklach, M.; Wang, H. Detailed Modeling of Soot Particle Nucleation and Growth. Proc. Combust. Inst. 1991, 23, 1559−1566. (16) Frenklach, M.; Schuetz, C. A.; Ping, J. Migration Mechanism of Aromatic-edge Growth. Proc. Combust. Inst. 2005, 30, 1389−1396. (17) Whitesides, R.; Domin, D.; Salomón-Ferrer, R.; Lester, W. A., Jr.; Frenklach, M. Embedded-Ring Migration on Graphene Zigzag Edge. Proc. Combust. Inst. 2009, 32, 577−583. (18) Whitesides, R.; Domin, D.; Salomón-Ferrer, R.; Lester, W. A., Jr.; Frenklach, M. Graphene Layer Growth Chemistry: Five- and SixMember Ring Flip Reaction. J. Phys. Chem. A 2008, 112, 2125−2130. (19) Whitesides, R.; Kollias, A. C.; Domin, D.; Lester, W. A., Jr.; Frenklach, M. Graphene Layer Growth: Collision of Migrating FiveMember Rings. Proc. Combust. Inst. 2007, 31, 539−546. (20) You, X.; Whitesides, R.; Zubarev, D.; Lester, W. A., Jr.; Frenklach, M. Bay-Capping Reactions: Kinetics and Influence on Graphene-Edge Growth. Proc. Combust. Inst. 2011, 33, 685−692. (21) Violi, A. Cyclodehydrogenation Reactions to Cyclopentafused Polycyclic Aromatic Hydrocarbons. J. Phys. Chem. A 2005, 109, 7781− 7787.

ASSOCIATED CONTENT

S Supporting Information *

Illustration of how a surrogate reaction technique was used to determine the canonical rates variationally for barrierless reactions in Variflex (Figure S1); vibrational frequencies, rotational constants, unsymmetrical hindered rotor parameters, and Cartesian coordinate geometries for the molecular structures of the present study (Table S1); high-pressure rate coefficients of the elementary reactions in the reaction systems discussed in the paper (Tables S2−S7); rate coefficients of the formation of products for each reaction system discussed in the paper (Tables S8−S14 and Figures S2−S5); illustrations of how the high-pressure-limit rate coefficients of product formation were calculated (Figures S6 and S7); an example of the stochastic error for the KMC simulations (Figure S8); average reaction counts for all of the KMC simulations (Tables S15−S32). The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/ acs.jpca.5b00868.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Authors

*Tel: 001-305-348-1945. E-mail: mebela@fiu.edu. *Tel: 001-510-643-1676. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was funded by the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Sciences of the U.S. Department of Energy (Grant DE-FG02-04ER15570 to Florida International University and Contract DE-AC03-76F00098 to Lawrence Berkeley National Laboratory) and by the U.S. Army Corps of Engineers, Humphreys Engineering Center Support Activity (Contract W912HQ-11-C-0035 to University of California at Berkeley). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office R

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (22) Whitesides, R.; Frenklach, M. Detailed Kinetic Monte Carlo Simulations of Graphene-Edge Growth. J. Phys. Chem. A 2010, 114, 689−703. (23) Tokmakov, I. V.; Kim, G.-S.; Kislov, V. V.; Mebel, A. M.; Lin, M. C. The Reaction of Phenyl Radical with Molecular Oxygen: A G2M Study of the Potential Energy Surface. J. Phys. Chem. A 2005, 109, 6114−6127. (24) Zhou, C.-W.; Kislov, V. V.; Mebel, A. M. Reaction Mechanism of Naphthyl Radicals with Molecular Oxygen. 1. Theoretical Study of the Potential Energy Surface. J. Phys. Chem. A 2012, 116, 1571−1585. (25) You, X.; Zubarev, D. Y.; Lester, W. A., Jr.; Frenklach, M. Thermal Decomposition of Pentacene Oxyradicals. J. Phys. Chem. A 2011, 115, 14184−14190. (26) Edwards, D. E.; You, X.; Zubarev, D. Y.; Lester, W. A., Jr.; Frenklach, M. Thermal Decomposition of Graphene Armchair Oxyradicals. Proc. Combust. Inst. 2013, 34, 1759−1766. (27) Raj, A.; da Silva, G. R.; Chung, S. H. Reaction Mechanism for the Free-Edge Oxidation of Soot by O2. Combust. Flame 2012, 159, 3423−3436. (28) Edwards, D. E.; Zubarev, D. Y.; Lester, W. A.; Frenklach, M. Pathways to Soot Oxidation: Reaction of OH with Phenanthrene Radicals. J. Phys. Chem. A 2014, 118, 8606−8613. (29) Kislov, V. V.; Singh, R. I.; Edwards, D. E.; Mebel, A. M.; Frenklach, M. Rate Coefficients and Product Branching Ratios for the Oxidation of Phenyl and Naphthyl Radicals: A Theoretical RRKM-ME Study. Proc. Combust. Inst. 2015, 35, 1861−1869. (30) Kroto, H. W. The Stability of the Fullerenes Cn, with n = 24, 28, 32, 36, 50, 60 and 70. Nature 1987, 329, 529−531. (31) Schmalz, T. G.; Seitz, W. A.; Klein, D. J.; Hite, G. E. Elemental Carbon Cages. J. Am. Chem. Soc. 1988, 110, 1113−1127. (32) Curl, R. F.; Smalley, R. E. Fullerenes. Sci. Am. 1991, 265, 54. (33) Aihara, J. Bond Resonance Energy and Verification of the Isolated Pentagon Rule. J. Am. Chem. Soc. 1995, 117, 4130−4136. (34) Yu, T.; Lin, M. C. Kinetics of the C6H5 + O2 Reaction at Low Temperatures. J. Am. Chem. Soc. 1994, 116, 9571−9576. (35) Schaugg, J.; Tranter, R. S.; Grotheer, H.-H. Rate Coefficients for the Fast Reactions of Phenyl Radicals with NO and O2 from Mass Spectrometric Measurements of Phenyl Decays. In Transport Phenomena in Combustion: Proceedings of the Eighth International Symposium on Transport Phenomena in Combustion, San Francisco, CA, July 16−20, 1995; pp 130−141. (36) da Silva, G. R.; Bozzelli, J. W. Variational Analysis of the Phenyl + O2 and Phenoxy + O Reactions. J. Phys. Chem. A 2008, 112, 3566− 3575. (37) Mebel, A. M.; Lin, M. C. Ab Initio Molecular Orbital Calculations of C6H5O2 Isomers. J. Am. Chem. Soc. 1994, 116, 9577−9584. (38) Carpenter, B. K. Computational Prediction of New Mechanisms for the Reactions of Vinyl and Phenyl Radicals with Molecular Oxygen. J. Am. Chem. Soc. 1993, 115, 9806−9807. (39) Gu, X.; Zhang, F.; Kaiser, R. I. A Crossed Beam Reaction of the Phenyl Radical, (C6H5, X2A1) with Molecular Oxygen (O2): Observation of the Phenoxy Radical, (C6H5O, X2A′). Chem. Phys. Lett. 2007, 448, 7−10. (40) Parker, D. S. N.; Zhang, F.; Kaiser, R. I. Phenoxy Radical (C6H5O) Formation under Single Collision Conditions from Reaction of the Phenyl Radical (C6H5, X2A1) with Molecular Oxygen (O2, X3Σg−): The Final Chapter? J. Phys. Chem. A 2011, 115, 11515− 11518. (41) Albert, D. R.; Davis, H. F. Collision Complex Lifetimes in the Reaction C6H5 + O2 → C6H5O + O. J. Phys. Chem. Lett. 2010, 1, 1107−1111. (42) Parker, D. S. N.; Kaiser, R. I.; Troy, T. P.; Kostko, O.; Ahmed, M.; Mebel, A. M. Toward the Oxidation of the Phenyl Radical and Prevention of PAH Formation in Combustion Systems. J. Phys. Chem. A 2014, DOI: 10.1021/jp509170x. (43) Marinov, N. M.; Pitz, W. J.; Westbrook, C. K.; Vincitore, A. M.; Castaldi, M. J.; Senkan, S. M.; Melius, C. F. Aromatic and Polycyclic

Aromatic Hydrocarbon Formation in a Laminar Premixed n-Butane Flame. Combust. Flame 1998, 114, 192−213. (44) Park, J.; Xu, Z. F.; Lin, M. C. Kinetic Study of the C10H7 + O2 Reaction. J. Phys. Chem. A 2009, 113, 5348−5354. (45) Raj, A.; Yang, S. Y.; Cha, D.; Tayouo, R.; Chung, S. H. Structural Effects on the Oxidation of Soot Particles by O2: Experimental and Theoretical Study. Combust. Flame 2013, 160, 1812−1826. (46) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab-Initio Calculation of Vibrational Absorption and CircularDichroism Spectra Using Density-Functional Force-Fields. J. Phys. Chem. 1994, 98, 11623−11627. (47) Baboul, A. G.; Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-3 Theory Using Density Functional Geometries and ZeroPoint Energies. J. Chem. Phys. 1999, 110, 7650−7657. (48) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Baboul, A. G.; Pople, J. A. Gaussian-3 Theory Using Coupled Cluster Energies. Chem. Phys. Lett. 1999, 314, 101−107. (49) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. Gaussian-3 (G3) Theory for Molecules Containing First and Second-Row Atoms. J. Chem. Phys. 1998, 109, 7764−7776. (50) 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.; et al. Gaussian 09, revision B.01; Gaussian, Inc.: Wallingford CT, 2010. (51) Werner, H.-J.; Knowles, P. J.; Kinizia, G.; Manby, F. R.; Schutz, M.; Celani, P.; Korona, T.; Lindh, R. MOLPRO, 2010.1; Ab Initio Programs: Lexington, MA, 2010. (52) Klippenstein, S. J. Implementation of RRKM Theory for Highly Flexible Transition States with a Bond Length as the Reaction Coordinate. Chem. Phys. Lett. 1990, 170, 71−77. (53) Klippenstein, S. J.; Wagner, A. F.; Dunbar, R. C.; Wardlaw, D. M.; Robertson, S. H. VARIFLEX, 1.00; Argonne National Laboratory: Argonne, IL, 1999. (54) Barker, J. R. Multiple-Well, Multiple-Path Unimolecular Reaction Systems. I. MultiWell Computer Program Suite. Int. J. Chem. Kinet. 2001, 33, 232−245. (55) Barker, J. R. Energy Transfer in Master Equation Simulations: A New Approach. Int. J. Chem. Kinet. 2009, 41, 748−763. (56) Barker, J. R.; Ortiz, N. F.; Preses, J. M.; Lohr, L. L.; Maranzana, A.; Stimac, P. J.; Nguyen, T. L.; Kumar, T. J. D. MultiWell Sof tware, 2014.1b; University of Michigan: Ann Arbor, MI, 2014. (57) Gillespie, D. T. Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem. 1977, 81, 2340−2361. (58) Gillespie, D. T. Markov Processes: An Introduction for Physical Scientists; Academic Press: San Diego, CA, 1992. (59) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell-Scientific: Oxford, U.K., 1990. (60) Wang, H.; Frenklach, M. Transport Properties of Polycyclic Aromatic Hydrocarbons for Flame Modeling. Combust. Flame 1994, 96, 163−170. (61) Whitesides, R.; Frenklach, M. Effect of Reaction Kinetics on Graphene-Edge Morphology and Composition. Z. Phys. Chem. (Muenchen, Ger.) 2015, 229, 597. (62) Frenklach, M. Monte Carlo Simulation of Diamond Growth by Methyl and Acetylene Reactions. J. Chem. Phys. 1992, 97, 5794−5802. (63) Frenklach, M. Monte Carlo Simulation of Hydrogen Reactions with the Diamond Surface. Phys. Rev. B: Condens. Matter Mater. Phys. 1992, 45, 9455−9458. (64) Meana-Paneda, R.; Truhlar, D. G.; Fernandez-Ramos, A. HighLevel Direct-Dynamics Variational Transition State Theory Calculations Including Multidimensional Tunneling of the Thermal Rate Constants, Branching Ratios, and Kinetic Isotope Effects of the Hydrogen Abstraction Reactions from Methanol by Atomic Hydrogen. J. Chem. Phys. 2011, 134. (65) Jodkowski, J. T.; Rayez, M. T.; Rayez, J. C. Theoretical Study of the Kinetics of the Hydrogen Abstraction from Methanol 3. Reaction of Methanol with Hydrogen Atom, Methyl, and Hydroxyl Atoms. J. Phys. Chem. A 1999, 103, 3750−3765. S

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (66) Hidaka, Y.; Oki, T.; Kawano, H.; Higashihara, T. Thermal Decomposition of Methanol in Shock Waves. J. Phys. Chem. 1989, 93, 7134−7139. (67) Nagle, J.; Strickland-Constable, R. F. Oxidation of Carbon between 1000−2000 °C. In Proceedings of the Fifth Carbon Conference; Pergamon: Oxford, 1962; pp 154−164.

T

DOI: 10.1021/acs.jpca.5b00868 J. Phys. Chem. A XXXX, XXX, XXX−XXX