ARTICLE pubs.acs.org/JPCA
Molecular-Dynamics-Based Study of the Collisions of Hyperthermal Atomic Oxygen with Graphene Using the ReaxFF Reactive Force Field Sriram Goverapet Srinivasan and Adri C. T. van Duin* Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ABSTRACT: In this work, we have investigated the hyperthermal collisions of atomic oxygens with graphene through molecular dynamics simulations using the ReaxFF reactive force field. First, following Paci et al. (J. Phys. Chem. A 2009, 113, 4677 4685), 5-eV energetic collisions of atomic oxygen with a 24-atom pristine graphene sheet and a sheet with a single vacancy defect, both functionalized with oxygen atoms in the form of epoxides, were studied. We found that the removal of an O2 molecule from the surface of the graphene sheet occurs predominantly through an Eley Rideal-type reaction mechanism. Our results, in terms of the number of occurrences of various reactive events, compared well with those reported by Paci et al. Subsequently, energetic collisions of atomic oxygen with a 25-timesexpanded pristine sheet were investigated. The steady-state oxygen coverage was found to be more than one atom per three surface carbon atoms. Under an oxygen impact, the graphene sheet was always found to buckle along its diagonal. In addition, the larger sheet exhibited trampoline-like behavior, as a result of which we observed a much larger number of inelastic scattering events than those reported by Paci et al. for the smaller system. Removal of O2 from the larger sheet occurred strictly through an Eley Rideal-type reaction. Investigation of the events leading to the breakup of a pristine unfunctionalized graphene sheet and the effects of the presence of a second layer beneath the graphene sheet in an AB arrangement was done through successive impacts with energetic oxygen atoms on the structures. Breakup of a graphene sheet was found to occur in two stages: epoxide formation, followed by the creation and growth of defects. Events leading to the breakup of a two-layer graphene stack included epoxide formation, transformation from an AB to an AA arrangement as a result of interlayer bonding, defect formation and expansion in the top layer, and finally erosion of the bottom layer. We observed that the breakup of the two-layer stack occurred through a sequential, layer-by-layer, erosion process.
1. INTRODUCTION Owing to their lightweight and desirable chemical and mechanical properties, carbon-based materials such as organic thin films, polymers, and carbon-fiber-reinforced composites are commonly used as ablatives and heat shields in spacecraft applications. At low-earth-orbit (LEO) altitudes, ranging from ∼200 to ∼700 km, the rarefied atmosphere consists of various chemical species that include atomic and molecular oxygen, helium, nitrogen, atomic hydrogen, and argon, among others.1 Dominant among these species is atomic oxygen, with a number density of ∼109 cm 3. A spacecraft traveling at an orbital velocity of ∼8 km/s in LEO is subjected to collisions with atomic oxygen in addition to exposure to vacuum ultraviolet light, X-rays, charged species such as O+, and other chemical species.1 These highly energetic collisions result in chemical and physical changes in the material through erosion and oxidation. A combination of the O-atom density and the orbital velocity of the spacecraft yields a flux of approximately 1015 O atoms cm 2 s 1 at a mean collision energy of ∼5 eV.1 To characterize the damage due to these high-energy collisions and predict and improve the long-term behavior of these materials, a fundamental understanding of the chemical reactions and dynamics involved in these processes is essential. r 2011 American Chemical Society
Experimental works done in the past frequently used highly ordered pyrolytic graphite (HOPG) as the model surface to study these highly energetic reactions. Kinoshita et al.2 characterized the surface of HOPG after it had been exposed to atomic O with a translational energy of 5 eV at room temperature using X-ray photoelectron spectroscopy (XPS) and scanning tunneling microscopy (STM). They found from the XPS spectra that the oxygen coverage reached a saturation value of 0.94 with an oxygen fluence of approximately 4 1017 atoms cm 2. Analysis of the C 1s peak by XPS revealed the presence of C—O, CdO, and O—CdO species. STM images showed that the HOPG surface formed hillock-like structures. Ngo et al.3 found that a carbon atom was etched for every eight impinging oxygen atoms from their measurements on an HOPG sample placed aboard the space shuttle Atlantis (mission STS46). On the other hand, the erosion rate of diamond was found to be at least 2 orders of magnitude smaller than that of graphite and also directiondependent, with the (100) direction being much less reactive than the (111) direction. The sp2/sp3 carbon-atom ratio was also Received: July 27, 2011 Revised: September 23, 2011 Published: September 26, 2011 13269
dx.doi.org/10.1021/jp207179x | J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A found to play an important role in the dynamics of the erosion process.4 8 Direct dynamics based on density-functional-theory(DFT-) based tight-binding (TB) methods suggested partial graphitization of the (111) surface and etching of this surface by way of CO2 desorption, whereas exposure to hyperthermal oxygen atoms resulted in full coverage of the (100) surface with ketones. The (100) surface was resistant to etching by hyperthermal oxygen atoms.9 Nicholson et al.10 12 conducted a series of experiments to study the morphological evolution of HOPG upon exposure to hyperthermal O atoms. They found that the erosion rate was strongly dependent on the surface temperature of the HOPG. The erosion rate tripled as the temperature was increased from 298 to 493 K. At 298 K, atomic force microscopy images revealed embedded cylindrical erosion pits whose diameters spanned nanometer through micrometer length scales and depths ranged from a few to tens of nanometers. The pits were found to have convex curvature at the bottom. At 493 K, the surface exhibited large towers and hillock-like features spanning hundreds of nanometers. These effects were attributed to the anisotropic kinetics involved in the lateral versus downward reactivity of HOPG. These experiments, however, did not elucidate the reaction mechanisms. Theoretical work aiming to describe the erosion of graphite has been done by a few researchers. Cohen13 attempted to estimate the loss function of graphite subjected to oxygen collisions using classical trajectories. It was suggested that the erosion of graphite is a two-step process, with the first step involving epoxide functionalization of the surface and the second the removal of this group as CO upon the impact of another O atom. It has been observed in combustion and gasification chemistry that the rate-determining step in the gasification of graphite involves a reaction that produces CO initially. The production of CO is believed to occur through the migration of an epoxide oxygen to defect sites that are functionalized by semiquinones. The semiquinone dissociates to produce CO, and the incoming oxygen forms a new semiquinone or a carbonyl functional group. The epoxide groups neighboring semiquinones are known to weaken carbon carbon bonds and facilitate this reaction.14,15 Hahn16 observed that the thermal oxidation of graphite at 550 950 °C led to the formation through layer-by-layer sequential erosion. The lateral etch rate was found to be 3 orders of magnitude larger than the vertical etch rate. Using ReaxFF-based molecular dynamics (MD) simulations on defect-free graphene oxide (GO), Bagri et al.17 observed that CO and CO2 evolution occurred through the desorption of epoxy groups in close proximity to other saturated sp3 carbon bonds. This led to the creation of vacancies within the basal plane. Defects formed upon high-temperature reduction of GO were always found to be decorated with carbonyl groups. 18 In a recent article, Paci et al.19 reported their findings on the hyperthermal collisions of atomic oxygen with graphite using beam surface scattering experiments and theoretical direct dynamics calculations based on DFT using the Perdew Burke Ernzerhof (PBE) generalized gradient approximation (GGA) functional with a double-ζ plus polarization (DZP) basis set. They observed that the graphite surface was functionalized with epoxides. Incoming oxygen atoms later reacted with the epoxides on the surface to form O2 by way of an Eley Ridealtype reaction. The steady-state concentration of oxygen atoms on the surface was found to be between one and two oxygen atoms per six surface carbon atoms. Theoretical studies were carried out with a periodic 24-atom pristine graphene sheet and a sheet with single vacancy defect. The simulations showed that
ARTICLE
semiquinones can form in functionalized pristine sheets that could later leave the sheet in the form of a CO (or CO2) molecule through direct collision of an incoming oxygen atom with the semiquinone or vibrational excitation due to a nearby oxygenatom collision. They illustrated how the epoxide groups neighboring the semiquinones catalyze the release of CO. Clearly, the fundamental understanding of the reaction mechanisms involved in graphene oxygen interactions is important from a number of different perspectives. Experimental research works2,8 12 have provided very valuable information about the scattering dynamics in such highenergy collision processes and detailed information about the morphology of the eroded graphite specimen. However, the time scales involved in these experiments are on the order of a few seconds at least. Most of the reactions typically have a picosecond time scale.19 As a result, molecular beam scattering experiments are unable to elucidate the reaction mechanisms involved in the process of erosion and oxidation of the HOPG sample. Ab initio calculations are particularly useful in this regard. However, the size of the systems that can be simulated using these methods is small. Use of periodic boundary conditions with a small system can introduce artificially high levels of symmetry in the system. Force-field-based molecular dynamics simulations offer the flexibility to simulate both the small-time-scale reactive events and larger system sizes simultaneously. In the present work, hyperthermal collisions of atomic oxygen with a graphene sheet were studied through molecular dynamics simulations using the ReaxFF reactive force field.20 Three independent sets of simulations were performed in our work. The first set of simulations was aimed at reproducing the results of Paci et al.19 on a 24-atom graphene sheet, whereas the second set of simulations investigated the effect of the size of the graphene sheet on the reactive events. The third set of simulations involved the continuous bombardment of a graphene sheet with hyperthermal atomic oxygen and a study of the effect of the presence of a second layer beneath the graphene sheet. The description of the systems and the details of the molecular dynamics simulations for all of the studies mentioned above are presented in the next section. A discussion of the results of the simulations is presented in section 3, followed by the conclusions in section 4.
2. MOLECULAR DYNAMIC SIMULATIONS Hyperthermal collisions of atomic oxygen with graphene were studied through molecular dynamics simulations using the ReaxFF reactive force field. ReaxFF is a general bondorder-dependent potential that uses relationships between bond distance and bond order and between bond order and bond energy to describe bond dissociation properly. The van der Waals and Coulomb interactions are included from the beginning in the formulation. The dissociation curves are derived from quantum chemical (QC) calculations. This approach retains the central force formalism wherein all atom pairs are assumed to have nonbonded interactions. These nonbonded interactions are described by Morse and Coulomb potentials. The nonbonded interactions are shielded at short range, so that their values become constant as the distance between any two atoms approaches 0. The force-field parameters used in the present simulations were obtained from Chenoweth et al.21 for hydrocarbon combustion. Similar parameters were used by Bagri et al.17 in their ReaxFF-based MD simulations on defect-free GO. 13270
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 1. Model I, sheet with a single vacancy defect functionalized with two oxygen atoms. Figure 3. Model III, sheet with a single vacancy defect functionalized with eight oxygen atoms.
Figure 2. Model II, sheet with a single vacancy defect functionalized with four oxygen atoms.
Three independent sets of simulations were carried out. Graphene sheets were subjected to collisions with 5-eV energetic oxygen atoms. Each oxygen-atom trajectory was simulated for a duration of 1 ps. In all simulations, an MD time step of 0.1 fs was used, and all simulations were carried out in the microcannocial (NVE) ensemble, because a local increase in temperature due to the collision event is essential in driving the reactions. The impact of the oxygen atom with the sheet causes the O graphene system to become vibrationally excited.22 This local vibrational excitation enables the breaking and creation of new bonds. Thus, coupling the system to an external heat bath through a thermostat would be unrealistic. In the first set of simulations, following Paci et al., 19 a 24-atom graphene sheet was considered with three-dimensional periodic boundary conditions. Two different sheets, a pristine graphene sheet and a sheet with a single vacancy defect, were studied in these simulations. A 5-eV oxygen atom is not sufficiently energetic to sputter one or more carbon atoms from the pristine graphene sheet. Thus, the pristine graphene sheet was functionalized with oxygen atoms in the form of epoxides to observe various possible reactions. This gives an oxygen-to-carbon atom ratio of 1:3. Vacancies in the graphene sheet are not expected to remain for long in an atmosphere containing oxygen. Thus, the valences of the dangling bonds at the defect sites were capped
Figure 4. Model IV, pristine sheet functionalized with eight oxygen atoms.
with oxygen atoms. Three different graphene sheets with a vacancy defect functionalized with two, four, and eight oxygen atoms were studied in our simulations. Figures 1 4 show the configurations of the four systems studied in the simulations (models I IV, respectively). The structures were energy-minimized before the MD simulations were performed. The z dimension of the periodic cell (normal to the graphene basal plane) was 30 Å which was sufficiently large for ReaxFF to reduce the interactions between the graphene layers in adjacent periodic cells in the z direction to a minimal value. Only normal incidence of the energetic oxygen atoms was considered. The point of origin (x, y coordinates) of the oxygen-atom trajectory was chosen randomly. The initial separation between the graphene basal plane and the impinging oxygen atom was 5 Å (z separation), although the presence of epoxides on the surface reduced the minimum separation between the functionalized sheet and the oxygen atom to a smaller value. One hundred independent simulations of oxygen-atom impingement were done for each of the models I IV. In the second set of simulations, the effect of the size of the system on the number of occurrences of various reactive events 13271
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 5. (a) Model V, obtained by expanding model IV by 25 times, 5 times along each of the coordinate directions of the basal plane; (b) pristine unfunctionalized graphene sheet; (c) two-layer graphene stack in the AB arrangement, with atoms in the bottom layer colored blue and those in the top layer colored gray.
was investigated for model IV alone, because the results for this case compared most favorably with those of Paci et al.19 As mentioned previously, for periodic systems, size effects can play a very important role in the dynamics. A periodic system that is too small can introduce artificially high amounts of symmetry and increase the stiffness of the system. In addition, small systems also prohibit phonon modes of certain wavelengths in the structure. A combination of these factors could have a profound impact on the reactivity of the system. Then, to study the effect of the system size, the structure from model IV was expanded 25 times, 5 times along each of the coordinate directions (x and y) lying on the basal plane of the graphene sheet. The expanded sheet thus had a total of 800 atoms with the ratio of the number of oxygen to carbon atoms maintained at 1:3 (200 oxygen atoms and 600
carbon atoms). Figure 5a shows the expanded graphene sheet, denoted as model V. The dimension of the periodic cell in the direction perpendicular to the basal plane of the graphene sheet (z direction) was held at 60 Å. The expanded structure was energy-minimized before being used for MD simulations. Again, 100 independent simulations of the impingement of oxygen atoms normal to the graphene basal plane were performed on model V. The origin of the oxygen-atom trajectory was 5 Å above the graphene basal plane, whereas the in-plane coordinates (x and y coordinates) were chosen randomly. The third set of simulations investigated the outcome of continuous bombardment of a pristine unfunctionalized graphene sheet (model VI) with hyperthermal atomic oxygen. Along with the simulations on a single graphene sheet, the effect 13272
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Table 1. Comparison of ReaxFF- and DFT-Based Simulations for the 24-Atom Graphene Sheet model I
model II
DFT
DFT
O1
O3
ReaxFF
ring-opening
36
32
epoxide formation
76
75
epoxide migration
37
carbonyl formation
3
O2 formation CO2 formation CO formation
0
dioxirane formation inelastic O
0 0
sheet damage
0
reaction
model III
model IV
DFT
O1
O3
ReaxFF
16
9
17
94
77
70
21
3
34
2
2
1
1
0
2
0
0
0 0
1
0 0
0 0
2 1
0
1
0
DFT
O1
O3
ReaxFF
O1
26
9
16
13
0
0
0
88
33
29
57
7
17
19
46
0
11
20
17
21
17
25
0
3
0
0
15
0
0
14
17
17
2
55
53
20
79
76
72
1
3
2
2
1
6
0
0
0
6
1
0
0 7
0 2
1 1
1 13
0 11
0 0
0 4
0 0
0
0
0
0
4
3
2
14
of the presence of a second layer beneath the graphene layer in an AB arrangement (model VII) with an interlayer distance of around 3.2 A was also investigated. Parts b and c, respectively, of Figure 5 show these two models. To distinguish between the top and bottom layers in model VII, the bottom layer is shown in a different color (blue). Three-dimensional periodic boundary conditions were imposed on the system. Model VI had a total of 384 carbon atoms, whereas model VII had twice as many. The system was energy-minimized before being used for MD simulations. The dimensions of the unit cell were 30.08 Å 34.76 Å (x and y dimensions) in the basal plane of graphene. The periodic cell was 240 Å long in the direction perpendicular to the graphene basal plane (z direction). Reflective boundary conditions were used in the z direction. These reflective boundary conditions reverse the z component of the velocity of any species that reaches the top (or bottom) face of the periodic cell (it acts like a virtual wall). Thus, the species leaving from the top side of the graphene sheet are confined to the space between the graphene sheet and the top face of the periodic cell. Reflective boundary conditions were used because it would be unrealistic if the species leaving from the top side of the graphene sheet were to enter from the bottom of the periodic cell and participate in the surface chemistry. In addition, the bulk side of a graphite specimen is inaccessible in experiments. Because the oxygen atoms were bombarded onto the graphene structure continuously (unlike in the previous two sets of simulations), the graphene sheet would eventually break up into a number of fragments. To ensure that the leaving species did not interfere with the surface chemistry, the z dimension of the periodic cell was chosen to be 240 Å, and the center of mass of the graphene basal plane was at 30 Å, thereby providing a space of about 210 Å between the graphene surface and the top face of the periodic cell. Continuous bombardment of the graphene sheet with hyperthermal oxygen atoms imparts a significant amount of translational energy to the sheet. Consequently, the sheet is expected to get displaced downward from its position. This would increase the initial vertical separation between the point of origin of the oxygen-atom trajectory and the graphene basal plane. To avoid this scenario, four dummy atoms were placed coplanar with the graphene sheet. A center-of-mass restraint between the atoms in the graphene sheet and the dummy atoms in the z-coordinate direction was imposed with a weak coupling. A weak coupling allows for the buckling of the graphene sheet
O3
ReaxFF
6 0
under the impact of the energetic oxygen atom while holding the z coordinate of the center of mass of the graphene sheet at approximately 30 Å. For model VII, eight dummy atoms were used, four in each of the graphene layers. Two sets of center-ofmass restraints were specified, one between the dummy atoms in the top layer and the carbon atoms in that layer and the other between the dummy atoms in the bottom layer and the carbon atoms in the bottom layer. The initial vertical separation between the graphene basal plane and the point of origin of the oxygenatom trajectory was 7.5 Å. The in-plane coordinates (x and y coordinates) of the oxygen atom were chosen randomly. Every oxygen-atom trajectory was simulated for a duration of 1 ps; that is, one oxygen atom was added to the system from a randomly chosen in-plane location (x and y coordinates) every picosecond.
3. RESULTS AND DISCUSSION In this work, we investigated the hyperthermal collisions of atomic oxygen with a graphene sheet through three independent sets of simulations. In the first set of simulations, we tried to reproduce the results of Paci et al.19 on a 24-atom graphene sheet (models I IV). A comparison of the results in terms of the number of occurrences of various reactive events is presented in Table 1. Most of the reactions observed by Paci et al.19 were observed in our simulations as well, suggesting that ReaxFF is capable of handling the chemistry involved in such highly energetic collisions. It can be seen from Table 1 that, with the present force-field parameters, the agreement between the ReaxFF-based simulations and the DFT-based simulations of Paci et al.19 is very good for model IV (pristine sheet functionalized with epoxides), whereas the agreement for models I III is fair. The discrepancies in the results for models I III (graphene sheet with a single vacancy defect) can be attributed to the nature of the present C/H/O force-field parameters. Although these parameters have been shown to reproduce the chemical reactivity trends of various aromatic and aliphatic hydrocarbons correctly,21 the training set for parameter optimization did not include structures dealing with defective graphene (or GO) species or similar compounds. In all six cases where CO2 formation was observed in our simulations, we observed only incipient CO2; that is, the CO2 molecule was still bound to the sheet. We believe that, if the simulations were to be run for a slightly longer period of time, the molecule would indeed be 13273
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 6. Sequence of events leading to the formation of a CO2 molecule. The pink atom represents the incoming hyperthermal oxygen atom, whereas the blue atom is a surface oxygen atom that forms a CO2 molecule. Carbon atoms are colored gray, and the other surface oxygen atoms are colored red.
released from the surface. Figure 6 shows one such scenario where the CO2 molecule was formed (thereby creating a vacancy defect on the graphene sheet) and yet remained attached to the
sheet at the end of the simulation period. The number of instances where we observed sheet damage and carbonyl formation was more than that observed by Paci et al.19 Apart from these 13274
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
cases, we also observed an instance of CO3 molecule formation. DFT calculations at the B3LYP/6-31G** level of theory indicate that the reaction of a CO2 molecule with atomic oxygen to yield a CO3 molecule is endothermic in nature, with a barrier of around 20 kcal/mol. Although the present C/H/O ReaxFF parameters predicted the reaction barrier precisely, the reaction on the whole was predicted to be exothermic with a heat of reaction of about 50 kcal/mol. Thus, the formation of a CO3 molecule in our simulations was an artifact of the present C/H/O ReaxFF parameters. O2 formation occurred almost always through an Eley Rideal-type reaction. The impinging oxygen atom plucked out an epoxide oxygen from the surface of the graphene sheet to form an O2 molecule. Among the 100 trajectories simulated, we observed only three cases where O2 formation occurred through a Langmuir Hinshelwood-type reaction. Paci et al.19 concluded that the steady-state surface coverage of oxygen was between those in models II and III because, in their DFT studies, the probability of O2 formation was higher in model III whereas the probability of epoxide formation was higher in model II. In contrast, in our ReaxFF-based MD simulations, we observed that the probability of epoxide formation was greater than the probability of O2 formation in all three models with a single vacancy defect (i.e., models I III). In the second set of simulations, we studied hyperthermal oxygen impacts on a version of model IV expanded be 25 times (model V, Figure 5a). We found that the number of trajectories leading to epoxide formation (see Table 2) was slightly more than the number of trajectories leading to the release of an O2 molecule, indicating that the steady-state surface coverage of Table 2. Comparison of ReaxFF- and DFT-Based Simulations on Models IV and V DFT (model IV)
ReaxFF
O1
O3
ring-opening
0
0
0
0
epoxide formation
7
17
19
37
epoxide migration
21
17
25
4
carbonyl formation
0
0
14
20
O2 formation
model IV
model V
79
76
72
30
CO2 formation CO formation
0 0
0
6 0
0 0
dioxirane formation
0
0
0
0
inelastic O
0
4
0
17
sheet damage
3
2
14
5
oxygen is greater than that in model V (or, equivalently, model IV); that is, the steady-state surface coverage is expected to be one oxygen atom for two to three carbon atoms. This indicates that the size of the sheet plays an important role in obtaining the right value for the steady-state coverage of oxygen atoms. A smaller graphene sheet makes the system too stiff and prevents phonon modes of certain wavelengths that can result in the sheet being predicted to be more reactive than it actually is. In all 100 simulations on model V, we observed that the sheet buckled significantly. Interestingly, we also observed that the sheet always buckled along its diagonal. The strain induced in the graphene basal plane by the epoxide functional groups is partially relieved by the formation of such a kink, which makes the configuration energetically more favorable.23,24 Figure 7 shows an instance of sheet buckling. The ability of a larger sheet to deform significantly causes it to act like a trampoline that can scatter the incoming oxygen atom. As a result, we saw an increased number of inelastic scattering events for model V (see Table 2) as compared to the smaller graphene sheets (models I IV). However, this trampoline-like effect of the graphene sheet might be subdued significantly in the presence of a second layer beneath the graphene sheet. This aspect of sheet deformation is yet to be investigated. Watanbe et al.25 report that interactions of hyperthermal Xe atoms with graphite result in only inelastic scattering because of the trampoline-like motion of the graphene sheets. Another noticeable difference between models I IV and model V lies in O2 formation. Whereas a few instances (although a small percentage) of O2 formation by the Langmuir Hinshelwood mechanism were observed in models I IV, O2 formation occurred exclusively by the Eley Rideal mechanism in model V. A comparison between the ReaxFF- and DFT-based simulation results for models IV and V is presented in Table 2. In the third set of simulations, we examined the effect of the periodic bombardment of a graphene sheet and a two-layer graphene stack in an AB arrangement with hyperthermal atomic oxygen. A continuous bombardment of graphene (or a graphene stack) with hyperthermal atomic oxygen eventually led to the breakup of the sheet (or stack) into a number of gaseous species (CO2, CO, C2, etc.). We observed that the breakup of a graphene sheet can be split into distinct regimes. A process timeline of events leading to the breakup of the graphene sheet and the twolayer stack is presented in Figure 8. As mentioned previously, an oxygen-atom trajectory was initiated every picosecond. The sequence of events leading to the breakup of the sheet is shown in Figure 9. The first regime corresponds to epoxide formation (frame 2). In this regime, the impinging oxygen atom binds to the surface in the form of an epoxide. This regime lasts from 0 to 45 ps. In this time period, the temperature and the relative potential
Figure 7. Buckling of model V along the diagonal due to a collision with an energetic oxygen atoms. 13275
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 8. Process timeline for the breakup of a graphene sheet (bottom) and a two-layer graphite stack in the AB arrangement (top) upon successive impacts with hyperthermal oxygen atoms.
Figure 9. Various regimes in the process of breakup of a graphene sheet due to successive bombardments with hyperthermal oxygen atoms. Carbon atoms are colored gray, and oxygen atoms are colored red.
energy (relative to that of the pristine graphene sheet) of the system increase, as can be seen from Figure 10. It must be noted that, in experiments10 12 or during actual flight, the graphitic material dissipates heat. Thus, the material does not experience such high temperatures (∼2000 K) as were encountered in our MD simulations. The first defect in the sheet appears at around 45.5 ps (frame 3) because of the formation of a CO2 molecule. Once a defect is formed, the formation of additional defects (or holes) in the sheet occurs rapidly. The graphene sheet is broken
into fragments by 75 ps. Thus, the time period from 46 to 75 ps corresponds to defect formation and growth and breakup of the sheet. Frame 4 in Figure 9 depicts a snapshot from the simulation results in this time period during which the sheet is on its way to breakup. During the process of defect growth, the temperature of the system remains approximately stable (see Figure 10), indicating that most of the energy input to the system (through the addition of 5-eV oxygen atoms) is consumed in breaking the carbon carbon bonds in the graphene sheet, as evidenced by the 13276
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 10. Variations of system temperature and potential energy with time upon successive impacts of hyperthermal oxygen atoms on a graphene sheet.
Figure 11. Number of C C bonds broken and CO2 molecules formed as a function of time upon successive impacts of hyperthermal oxygen atoms on a graphene sheet.
marked increase in the rate of C C bond breaking (Figure 11). Correspondingly, in this regime, a sharp rise in the potential energy can be seen in Figure 10. Figure 11 shows the number of C C bonds broken and the number of CO2 molecules formed as a function of time. To obtain the C C bond-breaking statistics, any two carbon atoms with an interatomic distance of less than 1.8 Å were considered to be bonded. The leaving species (CO2, etc.) were observed to leave from both the top and bottom sides of the graphene sheet. However, in our simulations on the twolayer graphene stack, we observed that the possibility of leaving
the sheet from the bottom was completely arrested in the presence of the second layer beneath the graphene sheet. The simulations on the two-layer graphene stack revealed some interesting aspects of the breakup process. Figure 12 shows the steps leading to the breakup of the two-layer stack. To begin, successive impingement of hyperthermal oxygen atoms on the graphene stack results in the formation of epoxides (frame 2), similar to that observed in the simulations on single-layer graphene. Soon, however, weak interlayer bonding between the carbon atoms in the top layer and those in the bottom layer starts 13277
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 12. Sequence of events leading to the breakup of a two-layer graphene stack. Carbon atoms in the bottom layer of the stack are colored blue, and those in the top layer are colored gray.
to occur. The first interlayer bond appeared at around 24.25 ps in our simulations (frame 3). Further oxygen-atom impacts, until 52 ps, resulted in the formation of more epoxides on the surface, with only one interlayer bond bridging the two layers. It must be noted, however, that this interlayer bond migrated over different carbon atoms; that is, as the simulation progressed, the carbon atoms in the two layers that were bonded changed. Between 52 and 54 ps, several carbon carbon bonds between the two layers were established, and the graphene stack transformed from an AB arrangement to an AA arrangement (frame 4). In the time period between 54 and 71 ps, impingement of oxygen resulted in the formation of epoxide or the removal of an oxygen molecule from the surface of the graphene stack (top layer) through an Eley Rideal-type direct reaction mechanism. The first defect in the stack appeared at around 71.25 ps as a result of the loss of a carbon atom in the form of a CO2 molecule. Recall that the first defect in the single-layer graphene sheet appeared at around 45.5 ps. Clearly, the presence of a second layer beneath a graphene
sheet delays the onset of sheet breakup. We suspect that this delay is a result of the bottom layer acting as a heat sink that damps the temperature increase in the top layer. Therefore, it requires a longer time to cross the energy barrier for the release of a CO2 molecule. The system temperature of the graphene stack, the bottom-layer temperature of the graphene stack, and the temperature of the single-layer graphene sheet are plotted in Figure 13 for comparison. The defect formed in the top layer (at around 71.25 ps) then expanded upon the release of more carbon atoms from the surface (in the form of CO, CO2, etc.). The top-layer graphene sheet was almost completely destroyed by 120 ps (frame 5). The leaving species from this layer left the surface almost exclusively from the top side of the graphene stack. Thus, the presence of a second layer beneath the graphene sheet prohibited the leaving species from leaving on the bulk side of the graphene stack (or graphite in general). We noticed that the breakup of the second layer statrted when the top layer was almost entirely etched away. 13278
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
Figure 13. Variation of temperature with time for a two-layer graphene stack and a graphene sheet.
As such, the erosion of a graphene stack (or graphite in general) seems to be a layer-by-layer sequential process. The entire graphene stack was broken apart by 150 ps. Frame 6 of Figure 11 shows a snapshot of the system from the simulation a few picoseconds before the stack was destroyed completely.
4. CONCLUSIONS We have investigated the hyperthermal collisions of atomic oxygen with graphene through molecular dynamics simulations using the ReaxFF reactive force field. To begin, we simulated the energetic collisions of atomic oxygen with a 24-atom graphene sheet. Following Paci et al.,19 we considered a pristine sheet and a sheet with a single vacancy defect, functionalized with epoxides. Generally, the results of our simulations in terms of the numbers of occurrences of various reactive events and the reaction mechanisms compared well with those of Paci et al.19 The results of our simulations agreed most favorably with those of Paci et al.19 for the pristine functionalized graphene sheet. The agreement for sheets with a single vacancy defect was only fair, with the largest discrepancy observed for the epoxide formation and migration and ring-opening reactions. This suggests that a retraining of the present C/H/O ReaxFF parameters using suitable ab initio data might be desirable. We observed that O2 formation occurs almost always by the Eley Rideal reaction mechanism, and only a few instances of the Langmuir Hinshelwood reaction mechanism were observed. To investigate the effect of system size, we simulated energetic oxygen-atom collisions on a 25-times-expanded pristine sheet functionalized with epoxides with the same oxygen-to-carbon atom ratio. We found that the impact of an energetic oxygen atom causes the graphene sheet to buckle along its diagonal. The steady-state coverage of oxygen was found to be more than one atom per three carbon atoms, contrary to the conclusions of Paci et al.19 We also observed that the larger sheet can act like a trampoline and scatter the approaching oxygen atom, as evidenced by a larger number of inelastic scattering events in
comparison with the results reported by Paci et al.19 for the smaller system. Thus, the size effect is an important aspect to consider in the dynamics of these collision processes. We also investigated the breakup of a graphene sheet and a twolayer graphene stack subjected to successive hyperthermal atomicoxygen impacts. The breakup of the graphene sheet proceeds through the formation of epoxides and the creation and expansion of vacancy defects. The leaving species can leave from both the top and bottom sides of the graphene sheet. Breakup of the two-layer graphene stack occurs through the formation of epoxides, conversion from an AB to AA arrangement as a result of interlayer bonding, creation and growth of defects in the top layer, and finally erosion of the second layer. The breakup of the graphene stack happens through a layer-wise sequential erosion process. The creation of the first defect in the top layer occurs much later than that in the single graphene sheet because the bottom layer acts as a heat sink that slows the temperature rise due to oxygen-atom collisions.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT We thank Jeffrey Paci for insightful discussions during the course of this work and for his suggestions on the manuscript. This research was supported by the Air Force Office of Scientific Research (AFOSR) under Grant FA9550-10-1-0563. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the AFOSR or the U.S. Government. ’ REFERENCES (1) Dressler, R. Chemical Dynamics in Extreme Environments; World Scientific: Singapore, 2001. 13279
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
The Journal of Physical Chemistry A
ARTICLE
(2) Kinoshita, H.; Umeno, M.; Tagawa, M.; Ohmae, N. Surf. Sci. 1999, 440, 49–59. (3) Ngo, T.; Snyder, E. J.; Tong, W. M.; Williams, R. S.; Anderson, M. S. Surf. Sci. 1994, 314, L817–L822. (4) Joshi, A.; Nimmagadda, R. J. Mater. Res. 1991, 6, 1484–1490. (5) Synowicki, R. A.; Hale, J. S.; Spady, B.; Reiser, M.; Nafis, S.; Woollam, J. A. J. Spacecr. Rockets 1995, 32, 97–102. (6) De Theije, F. K.; van der Lang, N. J.; Plomp, M.; Enckevort, W. J. P. Philos. Mag. A 2000, 80, 725–745. (7) De Theije, F. K.; van Veenendal, E.; Enckevort, W. J. P.; Vileg, E. Surf. Sci. 2001, 492, 91–105. (8) Shpilman, Z.; Gouzman, I.; Grossman, E.; Shen, L.; Minton, T. K.; Hoffman, A. Appl. Phys. Lett. 2009, 95, 174106. (9) Shpilman, Z.; Gouzman, I.; Grossman, E.; Shen, L.; Minton, T. K.; Paci, J. T.; Schatz, G. C.; Akhvlediani, R.; Hoffman, A. J. Phys. Chem. C 2010, 114, 18996–19003. (10) Nicholson, K. T.; Minton, T. K.; Sibener, S. J. J. Phys. Chem. B 2005, 109, 8476–8480. (11) Nicholson, K. T.; Sibener, S. J.; Minton, T. K. High Perf. Polym. 2004, 16, 197–206. (12) Nicholson, K. T.; Minton, T. K.; Sibener, S. J. Prog. Org. Coat. 2003, 47, 443–447. (13) Cohen, L. K. J. Chem. Phys. 1993, 99, 9652–9664. (14) Incze, A.; Pasturel, A.; Chatillon, C. Surf. Sci. 2003, 537, 55–63. (15) Hall, P. J.; Calo, J. M. Energy Fuels 1989, 3, 370–376. (16) Hahn, J. R. Carbon 2005, 43, 1506–1511. (17) Bagri, A.; Mattevi, C.; Acik, M.; Chabal, Y. J.; Chhowalla, M.; Shenoy, V. B. Nat. Chem. 2010, 2, 581–587. (18) Bagri, A.; Grantab, R.; Medhekar, N. V.; Shenoy, V. B. J. Phys. Chem. C 2010, 114, 12053–12061. (19) Paci, J. T.; Upadhyaya, H. P.; Zhang, J.; Schatz, G. C.; Minton, T. K. J. Phys. Chem. A 2009, 113, 4677–4685. (20) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396–9409. (21) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040–1053. (22) Isborn, C. M.; Li, X.; Tully, J. C. J. Chem. Phys. 2007, 126, 134307. (23) Schniepp, H.,C.; Li, J-L; McAllister, M. J.; Sai, H.; HerreraAlonso, M.; Adamson, D. H.; Prud’homme, R. K.; Car, R.; Saville, D. A.; Aksay, I. A. J. Phys. Chem. B 2006, 110, 8535–8539. (24) Li, J-L; Kudin, K. N.; McAllister, M. J.; Prud’homme, R. K.; Aksay, I. A.; Car, R. Phys. Rev. Lett. 2006, 96, 176101. (25) Watanbe, Y.; Yamaguchi, H.; Hashinokuchi, M.; Sawabe, K.; Maruyama, S.; Matsumoto, Y.; Shobatake, K. Chem. Phys. Lett. 2005, 413, 331–334.
13280
dx.doi.org/10.1021/jp207179x |J. Phys. Chem. A 2011, 115, 13269–13280
