- Email: [email protected]
Embedded-ring migration on graphene zigzag edge
Available online at www.sciencedirect.com
Proceedings of the
Proceedings of the Combustion Institute 32 (2009) 577–583
Combustion Institute www.elsevier.com/locate/proci
Embedded-ring migration on graphene zigzag edge Russell Whitesides a,b,*, Dominik Domin c, Romelia Salomo´n-Ferrer c,d, William A. Lester Jr. c,d, Michael Frenklach a,b a Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740, USA Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA c Kenneth S. Pitzer Center for Theoretical Chemistry, Department of Chemistry, University of California, Berkeley, CA 94720-1460, USA d Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
b
Abstract Reaction pathways are presented for hydrogen-mediated isomerization of a zigzag graphene edge containing a five-member carbon ring surrounded by six-member rings. A new reaction sequence in which this embedded five-member ring moves, or migrates, through the edge has been identified. The elementary steps of the pathways were analyzed using density functional theory (DFT). Rate coefficients were obtained by classical transition state theory utilizing the DFT energies, frequencies, and geometries. The results indicate that this new reaction sequence is competitive with the other important zigzag edge reactions allowing embedded five-member rings to move freely within a zigzag edge. The embedded rings have slight thermodynamic preference for the interior of the edge over the corner for large substrates. Published by Elsevier Inc. on behalf of The Combustion Institute. Keywords: Soot; PAH; Surface growth; Reaction mechanisms; HACA
1. Introduction The mechanisms of soot formation have been studied for quite some time, and there is general consensus that polycyclic aromatic hydrocarbons (PAH) are key intermediates: their growth leads to nucleation of particles and the latter continue to add mass via surface growth [1–5]. Detailed understanding of the underlying phenomena often invokes HACA [5], which describes the growth of both gaseous PAH molecules and graphitic edges *
Corresponding author. Address: Department of Mechanical Engineering, 244 Hesse Hall No. 1740, University of California, Berkeley, CA 94720-1740, USA. Fax: +1 510 643 5599. E-mail address: [email protected] (R. Whitesides).
of soot particles as a repetitive reaction sequence of substrate activation by gaseous hydrogen atoms followed by addition of carbon (usually acetylene) to the activated site thereby leading to the formation of the aromatic rings. In the initial formulation [6,7], the HACA sequence was presumed to take place on the arm-chair edge of graphene sheets H, C 2 H 2
;
ð1Þ
as the geometry of the arm-chair edge leads to energy-favorable ring closure. Recent studies have focused on zigzag graphene edges. On a flat zigzag edge, the HACA sequence forms a five-member ring, which can move (or migrate) along the zigzag edge,
1540-7489/$ - see front matter Published by Elsevier Inc. on behalf of The Combustion Institute. doi:10.1016/j.proci.2008.06.096
578
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
H
;
ð2Þ
and transform into a six-member ring at corners of the edge [8–10]. The mobility of surface rings was suggested to be of critical importance to the growth rate and evolving morphology [9], thus pointing to the need for understanding the ring migration chemistry. The ring collision reaction, H
;
ð3Þ
was found to be one of the important reaction steps [11]. The collision reaction allows a new layer to be initiated at the center of a layer instead of at the corners, implying that five-member rings could be included in the layer building process instead of being transformed into six-member rings. Most recently, the product of the ring collision has been examined and a pathway in which it isomerizes to reverse its orientation, or flips, H
;
ð4Þ
has been reported [12]. The rate of the flip reaction was found to be of the same order of magnitude as those of collision and migration, thus further implicating five-member rings in zigzag edge processes. One of the intriguing questions that reaction (4) brings is whether a five-member ring will be embedded in the growing layer or whether it will be ‘‘pushed out” to the corner of the edge via reactions like H
:
ð5Þ
This constitutes the subject of this study. First, we compute the minimum energy path of reaction (5) on the smallest substrate, determine its rates, and compare them with previously reported zigzag edge reactions. We then examine the reaction in systems of increasing size. 2. Computational methods 2.1. Energetics Density functional theory (DFT) was used to calculate the molecular and energetic parameters of all stable species and transition states. The molecules studied range from seven to nine aromatic rings with substrates modeled by the zigzag edges of tetracene and pentacene. The smallest system was the one on which the reaction can occur and the largest was that feasible for the quantum chemical method employed. Geometry optimizations were performed with the B3LYP hybrid functional [13] and the
6-311G(d,p) basis set. Previous studies have shown energetic predictions of B3LYP calculations at the 6-311G(d,p) level to be in good agreement with experimental and high-level ab initio results for stable species [14–16]. The energies of transition states predicted by this method, however, are often underestimated by about 5 kcal mol1 [17,18]. This shortcoming lessens the accuracy of rate constants derived from the calculated energetics yet allows for an order-of-magnitude analysis. There has been some discussion and debate in the literature about the nature of the ground state of large aromatic molecules. For large acenes, a triplet ground state was originally predicted by extrapolation from experimental data [19] and DFT calculations [20]. This conclusion has recently been put in question by broken symmetry open-shell DFT [21] and density matrix renormalization group [22] calculations, which show a difference between closed-shell and open-shell calculations and predict singlet polyradical character for the ground state of large acenes. To test the behavior of species in this study, we performed broken symmetry open-shell singlet and triplet calculations for species with an even number of electrons at the B3LYP/3-21G level. We found the open-shell singlet and triplet energies to be within a few kcal mol1 of each other and lower than the corresponding closed-shell singlet energy for most of the species. However, the open-shell singlet calculations suffer from significant spin contamination (hS2i 1). In addition, the triplet calculations are not reliable because they depend on the use of an inherently ground state method [23] to calculate excited state energies. We therefore treat species in this study with an even number of electrons as closed-shell singlets and those with an uneven number as doublets. The differences between energies calculated with open-shell and closed-shell treatments are similar to the expected uncertainty in barrier heights at the current level of theory. While the true electronic configuration of these species is of interest in a purely chemical sense and may aid in understanding the nature of observed paramagnetism of young soot particles [1], use of either the triplet or open-shell singlets in the rate coefficient calculations would not affect the major conclusions of this paper. Force calculations were performed at each predicted stationary point to confirm the geometry to be an energetic minimum (no imaginary frequencies) or a saddle point (one imaginary frequency). Transition states were confirmed to connect the reactant and product stable species by visual inspection of normal modes corresponding to the imaginary frequencies calculated at the B3LYP/6-311G(d,p) level and by intrinsic reaction coordinate calculations at the B3LYP/ 3-21G level. Zero-point energies were determined from the force calculations and scaled by a factor
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
of 0.9668 [24]. All calculations were performed using the Gaussian 03 suite of codes [25] on an Intel Xeon cluster. 2.2. Kinetics The kinetics of the reaction pathways were examined using classical transition state theory (TST) using version 2.08 of the MultiWell suite of codes [26,27]. The key inputs—reaction barriers, frequencies, and moments of inertia—were assigned from the DFT calculations at the B3LYP/6-311G(d,p) level of this study. Both bimolecular and unimolecular reactions were treated by transition state theory. We have found [9,11,12] that systems of this size are at the high pressure limit in the temperature range of interest and so TST and RRKM treatments are essentially identical. To evaluate the effect of tunneling at the conditions of our interest, we calculated the Wigner tunneling correction [28]. Even at the low end of the studied temperature range, the tunneling correction for hydrogen abstraction from a substrate molecule is about 3%, which is insignificant in the context of the current analysis and hence tunneling corrections were not applied. 3. Results and discussion 3.1. Reaction mechanism The minimum energy path (MEP) for the embedded-ring migration reaction on the smallest system examined is shown in Fig. 1. The overall reaction consists of a three-step mechanism initiated by hydrogen abstraction, followed by the unimolecular reaction migration step, and finished through hydrogen addition. The three systems examined are shown in Fig. 2. The energies of the individual species in the systems are compared in Table 1. Detailed information (energies, geom-
579
etries, etc.) about the molecular species involved can be found in Tables S1 and S2 of the Supplemental Material. The MEPs for all three systems are similar when comparing overall reaction energy, the isomerization barrier, as well as the energy of the individual species. One notable exception to this similarity is the energy of species 2 which varies in the three systems by as much as 10 kcal mol1 with respect to the reactant. The MEPs are also similar to that of the flip reaction initiated by hydrogen abstraction presented previously [12]. The latter reaction was found to dominate the flip reaction initiated by hydrogen addition presented in that work, and hence the embedded-ring migration analogous to the flip by addition has been disregarded in the current study. The energetic similarity of the flip reaction initiated by hydrogen abstraction and the embedded-ring migration by itself indicates that the embedded-ring migration reaction is a viable competitor with other zigzag edge reactions and the reaction rate analysis of these systems bears this out. The classical transition state theory rate coefficients computed for the individual steps in each system are given in Tables 2–4. The concentrations of H and H2 define the rate determining step for the overall embedded-ring migration reaction. For [H]/[H2] < 0.1, in the temperature and pressure range of interest, the rate determining step is the initiating abstraction reaction with a rate coefficient on the order of 1013 [H] cm3 mol1 s1, which is comparable to that found for the flip reaction [12] and two orders of magnitude larger than the rate of lone-ring migration, reaction (2) [9,11]. We have examined the hydrogen abstraction reaction (1 ? 2) in further detail, because it is the rate-limiting step in the embedded-ring migration and because of the differences in barrier heights between systems a, b, and c. While this methodology is limited in accuracy for this reaction, we can compare our results to higher-level
30 25
-1
E (kcal mol )
20 +H2
15 10
+H
5
+H2 2
+H
0 1
-5 3
-10
4
Fig. 1. Minimum energy path for embedded-ring migration reaction on system a.
580
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
a +H
+H
+H
+H
+H
+H
b
c
Fig. 2. Reacting systems analyzed.
Table 1 Energies of species in reaction systems (kcal mol1) from B3LYP/6-311G(d,p) including scaled ZPE Species 1+H 2 + H2 3 + H2 4+H TS 1–2 TS 2–3 + H2 TS 3–4
Reaction system a
b
c
0.0 4.5 5.4 8.1 8.8 24.6 3.7
0.0 5.7 5.7 0.0 0.0 26.2 0.0
0.0 4.1 9.5 6.5 2.5 20.0 4.8
calculations and experimental data for the case of hydrogen abstraction from benzene and we have done so in our previous study [12], where we found reasonable agreement among rate coefficients obtained experimentally and at different levels of theory. In comparing the reaction barriers, we note a lower barrier at the B3LYP level (10 kcal mol1) than at higher-level theories or inferred from experiment (16 kcal mol1) [29,30]
as expected from the documented performance of B3LYP discussed in Section 2.1. However, this does not explain the substantially lower barriers of systems b and c compared with that of system a. These lower barriers could be the result of the more constrained ring geometry or the inherent inaccuracy of the quantum chemical methodology. A higher barrier will lead to a slower reaction rate. In the context of this study, the slower the rate of abstraction, the stronger the argument that H-abstraction is the rate-limiting step of the embedded-ring migration reaction sequence. As a result, a higher barrier for H-abstraction will not alter our mechanistic conclusions. If one were to take the lowest reported rate for H-abstraction, the embedded-ring reaction may play a less prominent role in graphene chemistry, however, we still expect it to have a rate comparable to those of other important graphene edge reactions. Understanding the nature of competition between these reactions requires modeling of the surface kinetics including steric effects, which will be the subject of our future study.
Table 2 TST rate coefficients for elementary reactions of system a Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
4.88 1012 9.78 1012 1.70 1013 2.67 1013 3.92 1013
8.23 1011 1.39 1012 2.18 1012 3.21 1012 4.51 1012
1.24 1010 3.59 1010 8.00 1010 1.49 1011 2.46 1011
2.13 109 9.68 109 3.02 1010 7.31 1010 1.48 1011
2.90 1011 5.79 1011 1.02 1012 1.64 1012 2.48 1012
1.63 1012 3.61 1012 6.75 1012 1.12 1013 1.72 1013
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
581
Table 3 TST rate coefficients for elementary reactions of system b Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
3.39 1013 4.41 1013 5.53 1013 6.77 1013 8.10 1013
1.12 1012 1.98 1012 3.19 1012 4.79 1012 6.85 1012
7.05 108 3.49 109 1.16 1010 2.95 1010 6.23 1010
7.05 108 3.49 109 1.16 1010 2.95 1010 6.23 1010
1.12 1012 1.98 1012 3.19 1012 4.79 1012 6.85 1012
3.39 1013 4.41 1013 5.53 1013 6.77 1013 8.10 1013
Table 4 TST rate coefficients for elementary reactions of system c Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
1.46 1013 2.15 1013 2.95 1013 3.88 1013 4.92 1013
6.93 1011 1.27 1012 2.09 1012 3.21 1012 4.65 1012
9.76 109 3.33 1010 8.36 1010 1.71 1011 3.03 1011
3.08 109 1.37 1010 4.17 1010 9.94 1010 1.99 1011
6.68 1010 1.74 1011 3.73 1011 7.00 1011 1.20 1012
2.25 1012 4.94 1012 9.16 1012 1.52 1013 2.31 1013
3.2. Migration length In the interior of the layer, the embedded ring will not have a significant preference for the travel direction and will move through the edge as a onedimensional random walk. This diffusion-like process can also be compared to migration of adsorbed species on a diamond surface [31]. Just as in that study, we can evaluate the average number of migration steps, N, as the ratio of the embedded-ring migration rate to the rate of migration terminating processes, N ¼ Rmigration =Rtermination : ð6Þ In the case of embedded-ring migration, we assume the likely terminating process to be H, C 2 H 2
: ð7Þ
We estimate the rate of reaction 7 from the rate of analogous acetylene adsorption and cyclization reactions on the zigzag edge, reaction sequence S1 from ref. [9]. The average length of a ring’s movement within the layer, L, can be estimated from random walk theory [32] as L ¼ L0 N 1=2 ;
ð8Þ
where L0 is the distance travelled in a single migration step, reaction (5). We define the unitless average migration length, K, as L/L0 to simplify the notation in the following discussion. For mole fractions xH = 0.01 and xC2 H2 ¼ xH2 ¼ 0:1, typical of laminar premixed flame simu-
lations of past studies [9,10], K for the embeddedring migration varies from hundreds at 1500 K down to 10s at 2500 K. The high average migration lengths indicate that an embedded five-member ring will freely move through the entire edge available to it. In comparison, K for lone-ring migration, reaction (4), is on the order of 1–0.1 in the same temperature range. In the case of lonering migration, the rate of terminating processes, Rtermination, is dominated by ring desorption, H (-C 2H 2)
, ð9Þ
which is much faster than ring adsorption at high temperature [9], and which is not a path available to embedded rings. Because adsorption (via reaction (7)) is the dominant terminating process for embedded-ring migration, the average migration length, Kembedded, has a much stronger dependence on the relative concentrations of acetylene and hydrogen atom than the average length for the lone-ring migration, Klone. If [C2H2]/[H] is increased from its value of 10 in the example above to 106, Kembedded decreases by an order of magnitude. Klone is virtually unchanged for the same change in [C2H2]/[H]. For values of [C2H2]/[H] greater than 106, Kembedded varies only weakly. 3.3. Equilibrium The embedded-ring migration reaction has an equilibrium constant, Keq, equal to unity for reactions occurring in the interior of the layer,
582
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
7 6 5
K eq
4 3 2 1 0 1500
1750
2000 2250 Temperature (K)
2500
Fig. 3. Equilibrium constants as a function of temperature for systems a (closed circles), c (open squares), d (closed diamonds), e (closed triangles), and f (open circles). Systems a and c are those shown in Fig. 2, while systems d through f are the analogous systems with hexacene, heptacene, and octacene substrates, respectively.
H ,
ð10Þ because reactant and product are essentially identical. At the corner of the zigzag edge, the reaction no longer has symmetry, and so either the product, with the five-member ring on the corner, or the reactant, with the five-member ring in the interior, will be thermodynamically favored. We examined equilibrium constants for embedded-ring migration reactions at the corner of the edge (e.g. a1 ? a4 and c1 ? c4), for a series of systems of increasing size up to eight substrate rings (octacene) using DFT calculations at the B3LYP/3-21G level. Calculations at this level of theory were found to be in good agreement (to 0.2 kcal mol1) with those at the B3LYP/6-311G(d,p) level for the singlet molecules of reactant and product. The computed equilibrium constants are plotted in Fig. 3 as a function of temperature for all systems examined. The equilibrium constants are all order 1 in the temperature range examined. The system with Keq farthest from unity is a, due largely to two-fold rotational symmetry in the reactant of this system. For the three smallest systems, Keq is greater than one and decreases as temperature increases. Keq for the two largest systems varies only slightly with temperature and is approximately 0.5 over the entire temperature range studied. The leveling off is indicative of approaching arbitrarily-large system size with respect to equilibrium.
Fig. 4. Combined sequence of embedded- and lone-ring migration and collision reactions minimizing inclusion of five-member rings in a growing graphene layer.
Therfore, for large graphene edges, the less than unity equilibrium constant indicates a weak thermodynamic driving force causing embedded rings to be found more often away from the corners. Nonetheless, the embedded rings will still have ample access to the edge corner (e.g., with the probability of 1/3 for Keq = 0.5) where they may interact with migrating rings on the neighboring edge or with gas phase species. In other words, embedded rings are not static as previously suggested [9,10]. Instead, embedded rings appear to be highly mobile and will lead to annealing of surface defects. Figure 4 shows one possible sequence in which the presence of five-member rings can be minimized as a layer grows.
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
4. Conclusions A new reaction pathway has been identified in which an embedded five-member ring migrates through the zigzag edge of a graphene layer. Rate coefficients were found to be comparable to the previously investigated flip reaction and other competitive zigzag edge reactions. The fast kinetics indicates that the embedded ring moves essentially freely within the zigzag edge. On larger substrates, the reaction thermodynamically slightly favors the configuration with the fivemember ring in the interior of the edge as opposed to at the corner, causing embedded rings to be found more often away from the corner of zigzag edges. In spite of this slight thermodynamic preference, the occurrence of the embedded-ring migration reaction gives embedded rings ample access to the edge corner where they may interact with migrating rings or with gas phase species. The high mobility of embedded rings enables the layer to minimize the inclusion of five-member rings, and thus should contribute significantly to annealing and smoothing of growing surfaces. Acknowledgments Russell Whitesides, Romelia Salomo´n-Ferrer, William A. Lester Jr., and Michael Frenklach were supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division of the US Department of Energy, under Contract No. DE-AC03-76F00098. Dominik Domin was supported by the CREST Program of the National Science Foundation under Grant No. HRD-0318519. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.proci.2008.06.096. References [1] B.S. Haynes, H.G. Wagner, Prog. Energ. Combust. Sci. 7 (1981) 229–273. [2] K.H. Homann, Proc. Combust. Inst. 20 (1985) 857– 870. [3] I. Glassman, Proc. Combust. Inst. 22 (1989) 295– 311. [4] H. Richter, J.B. Howard, Prog. Energ. Combust. Sci. 26 (2000) 565–608.
583
[5] M. Frenklach, Phys. Chem. Chem. Phys. 4 (2002) 2028–2037. [6] M. Frenklach, D.W. Clary, W.C. Gardiner Jr., S.E. Stein, Proc. Combust. Inst. 20 (1985) 887–901. [7] M. Frenklach, H. Wang, Proc. Combust. Inst. 23 (1991) 1559–1566. [8] M. Frenklach, N.W. Moriarty, N.J. Brown, Proc. Combust. Inst. 27 (1998) 1655–1661. [9] M. Frenklach, C.A. Schuetz, J. Ping, Proc. Combust. Inst. 30 (2005) 1389–1396. [10] M. Frenklach, Proc. Combust. Inst. 26 (1996) 2285– 2293. [11] R. Whitesides, A.C. Kollias, D. Domin, W.A. Lester Jr., M. Frenklach, Proc. Combust. Inst. 31 (2007) 539–546. [12] R. Whitesides, D. Domin, R. Salomo´n-Ferrer, W.A. Lester Jr., M. Frenklach, J. Phys. Chem. A 112 (2008) 2125–2130. [13] A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. [14] M.R. Nimlos, J. Filley, J.T. McKinnon, J. Phys. Chem. A 109 (2005) 9896–9903. [15] N.W. Moriarty, N.J. Brown, M. Frenklach, J. Phys. Chem. A 103 (1999) 7127–7135. [16] J. Cioslowski, P. Piskorz, D. Moncrieff, J. Am. Chem. Soc. 120 (1998) 1695–1700. [17] J.P.A. Heuts, R.G. Gilbert, L. Radom, J. Phys. Chem. 100 (1996) 18997–19006. [18] J.L. Durant, Chem. Phys. Lett. 256 (1996) 595–602. [19] H. Angliker, E. Rommel, J. Wirz, Chem. Phys. Lett. 87 (1982) 208–212. [20] K.N. Houk, P.S. Lee, M. Nendel, J. Org. Chem. 66 (2001) 5517–5521. [21] M. Bendikov, H.M. Duong, K. Starkey, K.N. Houk, E.A. Carter, F. Wudl, J. Am. Chem. Soc. 126 (2004) 7416–7417. [22] J. Hachmann, J.J. Dorando, M. Aviles, G.K.L. Chan, J. Chem. Phys. 127 (2007) 134309– 134317. [23] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [24] R.D. Johnson, III (Ed.) NIST Computational Chemistry Comparison and Benchmark Database. National Institute of Standards and Technology, Aug 2005, available at http://srdata.nist.gov/cccbdb. [25] M.J. Frisch, et al., Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT, 2004. [26] J.R. Barker, N.F. Ortiz, J.M. Preses, L.L. Lohr, A. Maranzana, P.J. Stimac, MultiWell, 2.08, 2007, avilable at http://aoss.engin.umich.edu/multiwell/. [27] J.R. Barker, Int. J. Chem. Kinet. 33 (2001) 232–245. [28] E.P. Wigner, Z. Phys. Chem. Abt. B 19 (1932) 203. [29] J.H. Kiefer, L.J. Mizerka, M.R. Patel, H.C. Wei, J. Phys. Chem. 89 (1985) 2013–2019. [30] A.M. Mebel, M.C. Lin, T. Yu, K. Morokuma, J. Phys. Chem. A 101 (1997) 3189–3196. [31] M. Frenklach, S. Skokov, J. Phys. Chem. 101 (1997) 3025–3036. [32] D.T. Gillespie, Markov Processes: An Introduction for Physical Scientists, Academic, San Diego, CA, 1992.
Proceedings of the
Proceedings of the Combustion Institute 32 (2009) 577–583
Combustion Institute www.elsevier.com/locate/proci
Embedded-ring migration on graphene zigzag edge Russell Whitesides a,b,*, Dominik Domin c, Romelia Salomo´n-Ferrer c,d, William A. Lester Jr. c,d, Michael Frenklach a,b a Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740, USA Environmental Energy Technologies Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA c Kenneth S. Pitzer Center for Theoretical Chemistry, Department of Chemistry, University of California, Berkeley, CA 94720-1460, USA d Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
b
Abstract Reaction pathways are presented for hydrogen-mediated isomerization of a zigzag graphene edge containing a five-member carbon ring surrounded by six-member rings. A new reaction sequence in which this embedded five-member ring moves, or migrates, through the edge has been identified. The elementary steps of the pathways were analyzed using density functional theory (DFT). Rate coefficients were obtained by classical transition state theory utilizing the DFT energies, frequencies, and geometries. The results indicate that this new reaction sequence is competitive with the other important zigzag edge reactions allowing embedded five-member rings to move freely within a zigzag edge. The embedded rings have slight thermodynamic preference for the interior of the edge over the corner for large substrates. Published by Elsevier Inc. on behalf of The Combustion Institute. Keywords: Soot; PAH; Surface growth; Reaction mechanisms; HACA
1. Introduction The mechanisms of soot formation have been studied for quite some time, and there is general consensus that polycyclic aromatic hydrocarbons (PAH) are key intermediates: their growth leads to nucleation of particles and the latter continue to add mass via surface growth [1–5]. Detailed understanding of the underlying phenomena often invokes HACA [5], which describes the growth of both gaseous PAH molecules and graphitic edges *
Corresponding author. Address: Department of Mechanical Engineering, 244 Hesse Hall No. 1740, University of California, Berkeley, CA 94720-1740, USA. Fax: +1 510 643 5599. E-mail address: [email protected] (R. Whitesides).
of soot particles as a repetitive reaction sequence of substrate activation by gaseous hydrogen atoms followed by addition of carbon (usually acetylene) to the activated site thereby leading to the formation of the aromatic rings. In the initial formulation [6,7], the HACA sequence was presumed to take place on the arm-chair edge of graphene sheets H, C 2 H 2
;
ð1Þ
as the geometry of the arm-chair edge leads to energy-favorable ring closure. Recent studies have focused on zigzag graphene edges. On a flat zigzag edge, the HACA sequence forms a five-member ring, which can move (or migrate) along the zigzag edge,
1540-7489/$ - see front matter Published by Elsevier Inc. on behalf of The Combustion Institute. doi:10.1016/j.proci.2008.06.096
578
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
H
;
ð2Þ
and transform into a six-member ring at corners of the edge [8–10]. The mobility of surface rings was suggested to be of critical importance to the growth rate and evolving morphology [9], thus pointing to the need for understanding the ring migration chemistry. The ring collision reaction, H
;
ð3Þ
was found to be one of the important reaction steps [11]. The collision reaction allows a new layer to be initiated at the center of a layer instead of at the corners, implying that five-member rings could be included in the layer building process instead of being transformed into six-member rings. Most recently, the product of the ring collision has been examined and a pathway in which it isomerizes to reverse its orientation, or flips, H
;
ð4Þ
has been reported [12]. The rate of the flip reaction was found to be of the same order of magnitude as those of collision and migration, thus further implicating five-member rings in zigzag edge processes. One of the intriguing questions that reaction (4) brings is whether a five-member ring will be embedded in the growing layer or whether it will be ‘‘pushed out” to the corner of the edge via reactions like H
:
ð5Þ
This constitutes the subject of this study. First, we compute the minimum energy path of reaction (5) on the smallest substrate, determine its rates, and compare them with previously reported zigzag edge reactions. We then examine the reaction in systems of increasing size. 2. Computational methods 2.1. Energetics Density functional theory (DFT) was used to calculate the molecular and energetic parameters of all stable species and transition states. The molecules studied range from seven to nine aromatic rings with substrates modeled by the zigzag edges of tetracene and pentacene. The smallest system was the one on which the reaction can occur and the largest was that feasible for the quantum chemical method employed. Geometry optimizations were performed with the B3LYP hybrid functional [13] and the
6-311G(d,p) basis set. Previous studies have shown energetic predictions of B3LYP calculations at the 6-311G(d,p) level to be in good agreement with experimental and high-level ab initio results for stable species [14–16]. The energies of transition states predicted by this method, however, are often underestimated by about 5 kcal mol1 [17,18]. This shortcoming lessens the accuracy of rate constants derived from the calculated energetics yet allows for an order-of-magnitude analysis. There has been some discussion and debate in the literature about the nature of the ground state of large aromatic molecules. For large acenes, a triplet ground state was originally predicted by extrapolation from experimental data [19] and DFT calculations [20]. This conclusion has recently been put in question by broken symmetry open-shell DFT [21] and density matrix renormalization group [22] calculations, which show a difference between closed-shell and open-shell calculations and predict singlet polyradical character for the ground state of large acenes. To test the behavior of species in this study, we performed broken symmetry open-shell singlet and triplet calculations for species with an even number of electrons at the B3LYP/3-21G level. We found the open-shell singlet and triplet energies to be within a few kcal mol1 of each other and lower than the corresponding closed-shell singlet energy for most of the species. However, the open-shell singlet calculations suffer from significant spin contamination (hS2i 1). In addition, the triplet calculations are not reliable because they depend on the use of an inherently ground state method [23] to calculate excited state energies. We therefore treat species in this study with an even number of electrons as closed-shell singlets and those with an uneven number as doublets. The differences between energies calculated with open-shell and closed-shell treatments are similar to the expected uncertainty in barrier heights at the current level of theory. While the true electronic configuration of these species is of interest in a purely chemical sense and may aid in understanding the nature of observed paramagnetism of young soot particles [1], use of either the triplet or open-shell singlets in the rate coefficient calculations would not affect the major conclusions of this paper. Force calculations were performed at each predicted stationary point to confirm the geometry to be an energetic minimum (no imaginary frequencies) or a saddle point (one imaginary frequency). Transition states were confirmed to connect the reactant and product stable species by visual inspection of normal modes corresponding to the imaginary frequencies calculated at the B3LYP/6-311G(d,p) level and by intrinsic reaction coordinate calculations at the B3LYP/ 3-21G level. Zero-point energies were determined from the force calculations and scaled by a factor
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
of 0.9668 [24]. All calculations were performed using the Gaussian 03 suite of codes [25] on an Intel Xeon cluster. 2.2. Kinetics The kinetics of the reaction pathways were examined using classical transition state theory (TST) using version 2.08 of the MultiWell suite of codes [26,27]. The key inputs—reaction barriers, frequencies, and moments of inertia—were assigned from the DFT calculations at the B3LYP/6-311G(d,p) level of this study. Both bimolecular and unimolecular reactions were treated by transition state theory. We have found [9,11,12] that systems of this size are at the high pressure limit in the temperature range of interest and so TST and RRKM treatments are essentially identical. To evaluate the effect of tunneling at the conditions of our interest, we calculated the Wigner tunneling correction [28]. Even at the low end of the studied temperature range, the tunneling correction for hydrogen abstraction from a substrate molecule is about 3%, which is insignificant in the context of the current analysis and hence tunneling corrections were not applied. 3. Results and discussion 3.1. Reaction mechanism The minimum energy path (MEP) for the embedded-ring migration reaction on the smallest system examined is shown in Fig. 1. The overall reaction consists of a three-step mechanism initiated by hydrogen abstraction, followed by the unimolecular reaction migration step, and finished through hydrogen addition. The three systems examined are shown in Fig. 2. The energies of the individual species in the systems are compared in Table 1. Detailed information (energies, geom-
579
etries, etc.) about the molecular species involved can be found in Tables S1 and S2 of the Supplemental Material. The MEPs for all three systems are similar when comparing overall reaction energy, the isomerization barrier, as well as the energy of the individual species. One notable exception to this similarity is the energy of species 2 which varies in the three systems by as much as 10 kcal mol1 with respect to the reactant. The MEPs are also similar to that of the flip reaction initiated by hydrogen abstraction presented previously [12]. The latter reaction was found to dominate the flip reaction initiated by hydrogen addition presented in that work, and hence the embedded-ring migration analogous to the flip by addition has been disregarded in the current study. The energetic similarity of the flip reaction initiated by hydrogen abstraction and the embedded-ring migration by itself indicates that the embedded-ring migration reaction is a viable competitor with other zigzag edge reactions and the reaction rate analysis of these systems bears this out. The classical transition state theory rate coefficients computed for the individual steps in each system are given in Tables 2–4. The concentrations of H and H2 define the rate determining step for the overall embedded-ring migration reaction. For [H]/[H2] < 0.1, in the temperature and pressure range of interest, the rate determining step is the initiating abstraction reaction with a rate coefficient on the order of 1013 [H] cm3 mol1 s1, which is comparable to that found for the flip reaction [12] and two orders of magnitude larger than the rate of lone-ring migration, reaction (2) [9,11]. We have examined the hydrogen abstraction reaction (1 ? 2) in further detail, because it is the rate-limiting step in the embedded-ring migration and because of the differences in barrier heights between systems a, b, and c. While this methodology is limited in accuracy for this reaction, we can compare our results to higher-level
30 25
-1
E (kcal mol )
20 +H2
15 10
+H
5
+H2 2
+H
0 1
-5 3
-10
4
Fig. 1. Minimum energy path for embedded-ring migration reaction on system a.
580
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
a +H
+H
+H
+H
+H
+H
b
c
Fig. 2. Reacting systems analyzed.
Table 1 Energies of species in reaction systems (kcal mol1) from B3LYP/6-311G(d,p) including scaled ZPE Species 1+H 2 + H2 3 + H2 4+H TS 1–2 TS 2–3 + H2 TS 3–4
Reaction system a
b
c
0.0 4.5 5.4 8.1 8.8 24.6 3.7
0.0 5.7 5.7 0.0 0.0 26.2 0.0
0.0 4.1 9.5 6.5 2.5 20.0 4.8
calculations and experimental data for the case of hydrogen abstraction from benzene and we have done so in our previous study [12], where we found reasonable agreement among rate coefficients obtained experimentally and at different levels of theory. In comparing the reaction barriers, we note a lower barrier at the B3LYP level (10 kcal mol1) than at higher-level theories or inferred from experiment (16 kcal mol1) [29,30]
as expected from the documented performance of B3LYP discussed in Section 2.1. However, this does not explain the substantially lower barriers of systems b and c compared with that of system a. These lower barriers could be the result of the more constrained ring geometry or the inherent inaccuracy of the quantum chemical methodology. A higher barrier will lead to a slower reaction rate. In the context of this study, the slower the rate of abstraction, the stronger the argument that H-abstraction is the rate-limiting step of the embedded-ring migration reaction sequence. As a result, a higher barrier for H-abstraction will not alter our mechanistic conclusions. If one were to take the lowest reported rate for H-abstraction, the embedded-ring reaction may play a less prominent role in graphene chemistry, however, we still expect it to have a rate comparable to those of other important graphene edge reactions. Understanding the nature of competition between these reactions requires modeling of the surface kinetics including steric effects, which will be the subject of our future study.
Table 2 TST rate coefficients for elementary reactions of system a Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
4.88 1012 9.78 1012 1.70 1013 2.67 1013 3.92 1013
8.23 1011 1.39 1012 2.18 1012 3.21 1012 4.51 1012
1.24 1010 3.59 1010 8.00 1010 1.49 1011 2.46 1011
2.13 109 9.68 109 3.02 1010 7.31 1010 1.48 1011
2.90 1011 5.79 1011 1.02 1012 1.64 1012 2.48 1012
1.63 1012 3.61 1012 6.75 1012 1.12 1013 1.72 1013
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
581
Table 3 TST rate coefficients for elementary reactions of system b Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
3.39 1013 4.41 1013 5.53 1013 6.77 1013 8.10 1013
1.12 1012 1.98 1012 3.19 1012 4.79 1012 6.85 1012
7.05 108 3.49 109 1.16 1010 2.95 1010 6.23 1010
7.05 108 3.49 109 1.16 1010 2.95 1010 6.23 1010
1.12 1012 1.98 1012 3.19 1012 4.79 1012 6.85 1012
3.39 1013 4.41 1013 5.53 1013 6.77 1013 8.10 1013
Table 4 TST rate coefficients for elementary reactions of system c Temp. (K)
1500 1750 2000 2250 2500
Reaction 1?2 (cm3 mol1 s1)
2?1 (cm3 mol1 s1)
2?3 (s1)
3?2 (s1)
3?4 (cm3 mol1 s1)
4?3 (cm3 mol1 s1)
1.46 1013 2.15 1013 2.95 1013 3.88 1013 4.92 1013
6.93 1011 1.27 1012 2.09 1012 3.21 1012 4.65 1012
9.76 109 3.33 1010 8.36 1010 1.71 1011 3.03 1011
3.08 109 1.37 1010 4.17 1010 9.94 1010 1.99 1011
6.68 1010 1.74 1011 3.73 1011 7.00 1011 1.20 1012
2.25 1012 4.94 1012 9.16 1012 1.52 1013 2.31 1013
3.2. Migration length In the interior of the layer, the embedded ring will not have a significant preference for the travel direction and will move through the edge as a onedimensional random walk. This diffusion-like process can also be compared to migration of adsorbed species on a diamond surface [31]. Just as in that study, we can evaluate the average number of migration steps, N, as the ratio of the embedded-ring migration rate to the rate of migration terminating processes, N ¼ Rmigration =Rtermination : ð6Þ In the case of embedded-ring migration, we assume the likely terminating process to be H, C 2 H 2
: ð7Þ
We estimate the rate of reaction 7 from the rate of analogous acetylene adsorption and cyclization reactions on the zigzag edge, reaction sequence S1 from ref. [9]. The average length of a ring’s movement within the layer, L, can be estimated from random walk theory [32] as L ¼ L0 N 1=2 ;
ð8Þ
where L0 is the distance travelled in a single migration step, reaction (5). We define the unitless average migration length, K, as L/L0 to simplify the notation in the following discussion. For mole fractions xH = 0.01 and xC2 H2 ¼ xH2 ¼ 0:1, typical of laminar premixed flame simu-
lations of past studies [9,10], K for the embeddedring migration varies from hundreds at 1500 K down to 10s at 2500 K. The high average migration lengths indicate that an embedded five-member ring will freely move through the entire edge available to it. In comparison, K for lone-ring migration, reaction (4), is on the order of 1–0.1 in the same temperature range. In the case of lonering migration, the rate of terminating processes, Rtermination, is dominated by ring desorption, H (-C 2H 2)
, ð9Þ
which is much faster than ring adsorption at high temperature [9], and which is not a path available to embedded rings. Because adsorption (via reaction (7)) is the dominant terminating process for embedded-ring migration, the average migration length, Kembedded, has a much stronger dependence on the relative concentrations of acetylene and hydrogen atom than the average length for the lone-ring migration, Klone. If [C2H2]/[H] is increased from its value of 10 in the example above to 106, Kembedded decreases by an order of magnitude. Klone is virtually unchanged for the same change in [C2H2]/[H]. For values of [C2H2]/[H] greater than 106, Kembedded varies only weakly. 3.3. Equilibrium The embedded-ring migration reaction has an equilibrium constant, Keq, equal to unity for reactions occurring in the interior of the layer,
582
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
7 6 5
K eq
4 3 2 1 0 1500
1750
2000 2250 Temperature (K)
2500
Fig. 3. Equilibrium constants as a function of temperature for systems a (closed circles), c (open squares), d (closed diamonds), e (closed triangles), and f (open circles). Systems a and c are those shown in Fig. 2, while systems d through f are the analogous systems with hexacene, heptacene, and octacene substrates, respectively.
H ,
ð10Þ because reactant and product are essentially identical. At the corner of the zigzag edge, the reaction no longer has symmetry, and so either the product, with the five-member ring on the corner, or the reactant, with the five-member ring in the interior, will be thermodynamically favored. We examined equilibrium constants for embedded-ring migration reactions at the corner of the edge (e.g. a1 ? a4 and c1 ? c4), for a series of systems of increasing size up to eight substrate rings (octacene) using DFT calculations at the B3LYP/3-21G level. Calculations at this level of theory were found to be in good agreement (to 0.2 kcal mol1) with those at the B3LYP/6-311G(d,p) level for the singlet molecules of reactant and product. The computed equilibrium constants are plotted in Fig. 3 as a function of temperature for all systems examined. The equilibrium constants are all order 1 in the temperature range examined. The system with Keq farthest from unity is a, due largely to two-fold rotational symmetry in the reactant of this system. For the three smallest systems, Keq is greater than one and decreases as temperature increases. Keq for the two largest systems varies only slightly with temperature and is approximately 0.5 over the entire temperature range studied. The leveling off is indicative of approaching arbitrarily-large system size with respect to equilibrium.
Fig. 4. Combined sequence of embedded- and lone-ring migration and collision reactions minimizing inclusion of five-member rings in a growing graphene layer.
Therfore, for large graphene edges, the less than unity equilibrium constant indicates a weak thermodynamic driving force causing embedded rings to be found more often away from the corners. Nonetheless, the embedded rings will still have ample access to the edge corner (e.g., with the probability of 1/3 for Keq = 0.5) where they may interact with migrating rings on the neighboring edge or with gas phase species. In other words, embedded rings are not static as previously suggested [9,10]. Instead, embedded rings appear to be highly mobile and will lead to annealing of surface defects. Figure 4 shows one possible sequence in which the presence of five-member rings can be minimized as a layer grows.
R. Whitesides et al. / Proceedings of the Combustion Institute 32 (2009) 577–583
4. Conclusions A new reaction pathway has been identified in which an embedded five-member ring migrates through the zigzag edge of a graphene layer. Rate coefficients were found to be comparable to the previously investigated flip reaction and other competitive zigzag edge reactions. The fast kinetics indicates that the embedded ring moves essentially freely within the zigzag edge. On larger substrates, the reaction thermodynamically slightly favors the configuration with the fivemember ring in the interior of the edge as opposed to at the corner, causing embedded rings to be found more often away from the corner of zigzag edges. In spite of this slight thermodynamic preference, the occurrence of the embedded-ring migration reaction gives embedded rings ample access to the edge corner where they may interact with migrating rings or with gas phase species. The high mobility of embedded rings enables the layer to minimize the inclusion of five-member rings, and thus should contribute significantly to annealing and smoothing of growing surfaces. Acknowledgments Russell Whitesides, Romelia Salomo´n-Ferrer, William A. Lester Jr., and Michael Frenklach were supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences, Geosciences and Biosciences Division of the US Department of Energy, under Contract No. DE-AC03-76F00098. Dominik Domin was supported by the CREST Program of the National Science Foundation under Grant No. HRD-0318519. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.proci.2008.06.096. References [1] B.S. Haynes, H.G. Wagner, Prog. Energ. Combust. Sci. 7 (1981) 229–273. [2] K.H. Homann, Proc. Combust. Inst. 20 (1985) 857– 870. [3] I. Glassman, Proc. Combust. Inst. 22 (1989) 295– 311. [4] H. Richter, J.B. Howard, Prog. Energ. Combust. Sci. 26 (2000) 565–608.
583
[5] M. Frenklach, Phys. Chem. Chem. Phys. 4 (2002) 2028–2037. [6] M. Frenklach, D.W. Clary, W.C. Gardiner Jr., S.E. Stein, Proc. Combust. Inst. 20 (1985) 887–901. [7] M. Frenklach, H. Wang, Proc. Combust. Inst. 23 (1991) 1559–1566. [8] M. Frenklach, N.W. Moriarty, N.J. Brown, Proc. Combust. Inst. 27 (1998) 1655–1661. [9] M. Frenklach, C.A. Schuetz, J. Ping, Proc. Combust. Inst. 30 (2005) 1389–1396. [10] M. Frenklach, Proc. Combust. Inst. 26 (1996) 2285– 2293. [11] R. Whitesides, A.C. Kollias, D. Domin, W.A. Lester Jr., M. Frenklach, Proc. Combust. Inst. 31 (2007) 539–546. [12] R. Whitesides, D. Domin, R. Salomo´n-Ferrer, W.A. Lester Jr., M. Frenklach, J. Phys. Chem. A 112 (2008) 2125–2130. [13] A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. [14] M.R. Nimlos, J. Filley, J.T. McKinnon, J. Phys. Chem. A 109 (2005) 9896–9903. [15] N.W. Moriarty, N.J. Brown, M. Frenklach, J. Phys. Chem. A 103 (1999) 7127–7135. [16] J. Cioslowski, P. Piskorz, D. Moncrieff, J. Am. Chem. Soc. 120 (1998) 1695–1700. [17] J.P.A. Heuts, R.G. Gilbert, L. Radom, J. Phys. Chem. 100 (1996) 18997–19006. [18] J.L. Durant, Chem. Phys. Lett. 256 (1996) 595–602. [19] H. Angliker, E. Rommel, J. Wirz, Chem. Phys. Lett. 87 (1982) 208–212. [20] K.N. Houk, P.S. Lee, M. Nendel, J. Org. Chem. 66 (2001) 5517–5521. [21] M. Bendikov, H.M. Duong, K. Starkey, K.N. Houk, E.A. Carter, F. Wudl, J. Am. Chem. Soc. 126 (2004) 7416–7417. [22] J. Hachmann, J.J. Dorando, M. Aviles, G.K.L. Chan, J. Chem. Phys. 127 (2007) 134309– 134317. [23] P. Hohenberg, W. Kohn, Phys. Rev. 136 (1964) B864. [24] R.D. Johnson, III (Ed.) NIST Computational Chemistry Comparison and Benchmark Database. National Institute of Standards and Technology, Aug 2005, available at http://srdata.nist.gov/cccbdb. [25] M.J. Frisch, et al., Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford, CT, 2004. [26] J.R. Barker, N.F. Ortiz, J.M. Preses, L.L. Lohr, A. Maranzana, P.J. Stimac, MultiWell, 2.08, 2007, avilable at http://aoss.engin.umich.edu/multiwell/. [27] J.R. Barker, Int. J. Chem. Kinet. 33 (2001) 232–245. [28] E.P. Wigner, Z. Phys. Chem. Abt. B 19 (1932) 203. [29] J.H. Kiefer, L.J. Mizerka, M.R. Patel, H.C. Wei, J. Phys. Chem. 89 (1985) 2013–2019. [30] A.M. Mebel, M.C. Lin, T. Yu, K. Morokuma, J. Phys. Chem. A 101 (1997) 3189–3196. [31] M. Frenklach, S. Skokov, J. Phys. Chem. 101 (1997) 3025–3036. [32] D.T. Gillespie, Markov Processes: An Introduction for Physical Scientists, Academic, San Diego, CA, 1992.