Effects of Thermal Electronic Excitations on the ... - ACS Publications

Apr 13, 2016 - Copyright © 2016 American Chemical Society. *T. Sun. E-mail: [email protected]., *S. Fabris. E-mail: [email protected]. Cite this:J. Phys...
2 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCA

Effects of Thermal Electronic Excitations on the Diffusion of Oxygen Adatoms on Graphene Tao Sun,*,†,‡ Xinxin Yao,† and Stefano Fabris*,§,∥ †

College of Earth Sciences, University of Chinese Academy of Sciences, Beijing 100049, China Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China § CNR-IOM DEMOCRITOS, Istituto Officina dei Materiali, Consiglio Nazionale delle Ricerche, Via Bonomea 265, 34136 Trieste, Italy ∥ SISSA, Scuola Internazionale Superiore di Studi Avanzati, via Bonomea 265, I-34136 Trieste, Italy ‡

ABSTRACT: We conduct first-principles calculations to study oxygen diffusion on the graphene surface as a function of temperature up to 3000 K. The minimum energy migration path and the corresponding activation energy are determined by the nudged elastic band method with explicit inclusion of thermal electronic excitations. Below 1000 K the activation energy for epoxy oxygen to migrate remains close to its room temperature value (0.72 eV). Above 1000 K the activation energy decreases near linearly with temperature, from 0.70 eV at 1000 K to 0.47 eV at 3000 K. We show that this reduction originates from thermal electronic excitations. In particular, the effect is determined by the large contrasts in the electronic structures of the initial and transition states: the transition state exhibits much larger electronic density of states near the Fermi level and is more susceptible to thermal electronic excitations. The reduction in activation energy leads to appreciable enhancement in the diffusivity of oxygen adatoms. A moderate decrease in the vibrational prefactor, also caused by thermal electronic excitations, does not alter this trend. These findings may facilitate future works to accurately describe the dynamics of O adatoms on graphene at high T, which are critical for determining surface thermodynamic properties such as equilibrium coverage.

1. INTRODUCTION Oxidation of graphitic materials is important in a wide variety of technological fields. Coal and biomass combustion have provided heating and energy for human beings since ancient times.1,2 More recently, oxidation has emerged as an effective route to functionalize and manipulate graphene.3,4 Yet unwanted oxidation can also be hazardous. For instance, erosion at high temperatures by atomic or molecular oxygen may cause structural failures of graphitic components used in spacecrafts5,6 and nuclear reactors.7,8 Because of its great significance, the oxidation process has been subjected to numerous experimental9−13 and theoretical studies.14−19 As oxidation of interest often occurs at high temperatures, including thermal effects in the theoretical modeling is critical.16−19 In principle, there are two types of thermal effects that need to be considered: thermal motions of ions and thermal excitations of electrons. The former have been accounted for in previous theoretical studies using molecular dynamics17,19 or transition state theory via a vibrational prefactor.20−23 The latter, which drive electrons from the electronic ground state to the states of higher energies such that the electron occupancy in thermal equilibrium obeys the Fermi−Dirac distribution, have not received much attention so © XXXX American Chemical Society

far. And their role in the high temperature oxidation process remains unclear and may in principle affect the theoretical predictions. Here we use diffusion of epoxy oxygen on graphene surface as a prototype process to illustrate how thermal excitations of electrons may affect high temperature oxidation of graphitic materials. Such epoxy oxygen atoms are common on the oxidized surface, 22,24 where they bridge between two neighboring C atoms of the hexagonal C network.14,25 They may originate from two sources: direct adsorption of atomic oxygen atoms,13,24 or O spillover after dissociations of O2 molecules at graphene edges.22 Diffusion of epoxy oxygen is believed to be the key process that facilitates the “indirect” CO2 formation during carbon combustion (T ≥ 1000 K).16,18 It is also the rate-limiting step that controls oxidative unzipping and cutting of graphene.25 First-principles calculations based on density functional theory (DFT) have provided valuable insights on the diffusion of epoxy oxygen.21,23,26,27 The transition state was identified to be the configuration where Received: January 14, 2016 Revised: April 12, 2016

A

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A the O adatom is on top of a C atom in graphene.21 The activation energy was found to depend on external electric fields,21 neighboring oxygen functional groups,23 curvatures of graphene sheets,26 and topological defects.27 No temperature dependence of the activation energy has been suggested so far, presumably because these DFT calculations were performed with electrons being kept at the electronic ground state. In this study, we take into account thermal electronic excitations explicitly by using the finite temperature electronic structure theory.28 Our results show that strong thermal electronic excitations can alter the profiles of free energy surfaces considerably, leading to T-dependent activation energies. The paper is organized as follows: Simulation methods are reported in section 2. Activation energies as well as vibrational prefactors for oxygen diffusion at different temperatures are shown in section 3. This section also includes a detailed analysis of electronic structures so as to identify the origin of the observed temperature dependency. Section 4 compares the effect of thermal electronic excitations on oxygen diffusion on graphene with that of electron doping. It also discusses how the mobility of O adatoms may affect its equilibrium coverage, and the importance of describing the dynamics of O adatoms accurately. Conclusions are summarized in section 5.

Γ=

⎛ ΔU ⎞ = ν0 exp⎜ − 0 ⎟ ⎝ kBT ⎠

1 2 a Γ 4

where ν0 ≡ is the vibrational prefactor. Equation 5 is known as the Vineyard’s jump rate formula20 widely adopted in first-principles diffusivity calculations. It follows the Arrhenius’ law Γ ∝ exp(−Ea/kBT), with a temperature-independent activation energy Ea equaling ΔU0. For metals or semimetals, besides lattice vibrations thermal electronic excitations also contribute to the system’s free energy. Define the electronic free energy (containing the static lattice energy) at T as Fel(T), the total free energy F = Fel + Fvib. Equation 5 can then be extended by replacing ΔU0 with ΔFel as ⎛ ΔF (T ) ⎞ Γ = ν0 exp⎜ − el ⎟ kBT ⎠ ⎝

(1)

kBT e−FTS/ kBT kT = B e−ΔF / kBT h e−FIS/ kBT h

(2)

where h and kB denote the Planck and Boltzmann constants, respectively, ΔF equals FTS − FIS. For an insulating system, the free energy F consists of the static energy U0 and the vibrational free energy Fvib. Under harmonic approximation, Fvib is expressed in terms of the frequencies of vibrational normal modes νi as Fvib = kBT ln(∏ [2 sinh(hνi /2kBT )]) i

(3)

Equation 2 can then be rewritten as Γ=

3N − 3 ⎛ ΔU ⎞ 2kBT ∏i = 1 sinh(hνi /2kBT ) exp⎜ − 0 ⎟ 3N − 4 h ∏j = 1 sinh(hνj*/2kBT ) ⎝ kBT ⎠

(6)

where ΔFel = − For a given configuration its electronic free energy Fel can be readily determined in finite-temperature DFT simulations by minimizing the Mermin functional.28 If one ignores the temperature dependence of vibrational frequencies induced by thermal electronic excitations, then ΔFel may be assigned as the activation energy at T, in a similar fashion as ΔU0 in eq 5. To summarize, we remark that eq 6 is derived in the framework of transition state theory, which relates the reaction rate to the free energy difference between the initial and transition states. Under the harmonic approximation, the vibrational part of the free energy difference ΔFvib can be factored out from the exponential expression of the rate and expressed as the vibrational prefactor ν0. Thus, entropic contributions from lattice vibrations are included in ν0. The transition state and the associated electronic part of the free energy difference ΔFel are determined by exploring the electronic free energy surface(or potential energy surface for insulators) to identify the minimum energy path. This can be done by conducting nudged elastic band calculations at the corresponding temperature T. Once the transition state is identified and ΔFel is known, ν0 can be evaluated from the vibrational frequencies of the transition and initial states. The atom jump rate as well as diffusion coefficient can then be readily determined from eqs 6 and 1. Besides harmonic transition state theory, diffusion processes can also be studied via molecular dynamics (MD),31 often in combination with accelerated sampling techniques such as metadynamics.32,33 In such MD studies, free energy surfaces are explored directly and explicit decomposition of the total free energy is not needed. A major advantage of MD studies is that lattice vibrations are treated exactly, beyond the harmonic approximation. Moreover, one may consider explicitly the time evolution of electron wave functions and their nonadiabatic coupling with ionic motions.34 Here we limit our scope to evaluating the activation energy for diffusion in the framework of harmonic transition state theory, in the hope that it may serve as a starting point for further investigations. 2.2. Simulation Details. Our density functional theory (DFT) simulations were carried out in the plane-waves pseudopotential framework as implemented in the Pwscfcode of the Quantum ESPRESSO distribution.35,36 Graphene was modeled with a 5 × 5 periodic supercell separated in the FTS el

where a is the distance between two adjacent surface adsorption sites and Γ is the atom jump rate.21 For diffusion proceeds in a NVT ensemble or at ambient pressure, the atom jump rate Γ can be determined from the transition state theory29 via the Helmholtz free energy F of the system in the transition (TS) and initial (IS) states: Γ=

(5)

3N−4 ∏3N−3 i=1 νi/∏j=1 ν* j

2. METHODS 2.1. Diffusion Coefficient from Harmonic Transition State Theory. The diffusivity of a surface adatom is characterized by the diffusion coefficient, expressed in terms of microscopic quantities as D=

⎛ ΔU0 ⎞ exp ⎜− ⎟ 3N − 4 ⎝ kBT ⎠ ∏j = 1 νj* 3N − 3

∏i = 1 νi

(4)

where νi and νj* are the frequencies of the IS and TS, respectively.30 Note that the three translational modes with zero frequencies, as well as the unstable mode in the transition state whose polarization vector corresponds to the migration direction of the adatom, are excluded in the products. In the classical limit (kBT/h → ∞), B

FIS el .

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. Surface structures along the epoxy O diffusion path at T = 300 K: (a) initial state; (b) transition state; (c) final state. The red sphere represents the O adatom. The gray sphere denotes the C atom. The dashed circle marks the position of a second epoxy oxygen at the next-nearest neighbor.

Table 1. Activation Energies (ΔFel) and Geometric Parameters of the Initial, Metastable, and Transition States at Different Temperatures initial state (IS) temp (K)

ΔFel (eV)

dOC (Å)

dCC (Å)

300 1500 2500 3000

0.72 0.65 0.53 0.47

1.47 1.47 1.47 1.47

1.52 1.51 1.51 1.51

metastable state (MS) dOC (Å)

1.41 1.41

dCC (Å)

1.50 1.50

103.2 101.6

transition state (TS)

∠OCC (deg)

dOC (Å)

dCC (Å)

104.1 105.5

1.42 1.41 1.42 1.42

1.49 1.49 1.49 1.49

104.9 105.7

∠OCC (deg) 101.0 101.0 91.5 89.5

101.1 104.1 108.7 109.4

107.5 106.0 108.9 109.6

Figure 2. (a) Activation energies (ΔFel) of epoxy oxygen diffusion at different temperatures. Open squares denote the diffusion barriers of an isolated epoxy oxygen. The solid line connecting the data points serves as a guide of the eye. Filled squares denote the diffusion barriers of an epoxy oxygen in the presence of a second epoxy oxygen at the next-nearest neighbor. (b) Electronic free energies of different configurations along the O diffusion path. Data points correspond to the energies of structural replicas in the NEB calculations. Solid lines are the interpolated minimum energy paths.

and the corresponding Fel was determined by minimizing the Mermin functional.28 The Mermin functional formalism allows accounting for the nonadiabatic effects that may arise at high T in an ensemble-averaged fashion.34 The minimum energy path and saddle point configuration (transition state) were determined by the climbing image nudged elastic band (NEB) method41 employing 8 system replicas. Vibrational frequencies of the initial and transition states were computed from density functional perturbation theory.42 In particular, the effect of thermal electronic excitations on vibrational frequencies were evaluated and the temperature dependence of ν0 was determined. Quantum effects of nuclear motion, as shown in eq 4, were not considered in this study as we are interested in the high T regime where such effects are weak: explicit calculations show that the quantum correction to the classical rate constant decreases with temperature and becomes less than 5% at 1500 K.

perpendicular direction by 14 Å of vacuum. Such supercells have been shown to be sufficiently large to yield converged energetics among various graphene oxide surface structures.37 Because graphene has a tiny negative thermal expansion coefficient (∼−10−6 K−1),38 the effect of thermal expansion can be safely neglected and the in-plane lattice parameter of graphene was fixed to the theoretical static equilibrium value 2.463 Å.12,37 The Perdew−Burke−Ernzerhof generalized gradient approximation (PBE-GGA)39 was used for the electron exchange-correlation energy functional. The energy cutoffs of the plane-wave basis set and the charge density representation were set to 30 and 300 Ry, respectively. A 4 × 4 × 1 Monkhorst−Pack k-point mesh40 was used for Brillouin zone sampling, which yields a convergence of 0.01 eV in relative energy differences between different surface structures. Diffusion coefficients at different temperatures were calculated according to eq 6. At each temperature, electron occupancy was set according to the Fermi−Dirac distribution C

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 3. Partial and total electronic density of states for an isolated oxygen in the 5 × 5 graphene supercell: (a) the initial/final epoxy structure (O at the bridge site); (b) the transition state (O at the top site). Electron occupancy at 300 and 3000 K are shown in dashed lines.

Figure 4. Ratio between the vibrational frequencies at 3000 K and those at 300 K for an isolated oxygen in the 5 × 5 graphene supercell: (a) the initial state (O at the bridge site); (b) the transition state (O at the top site). Three translational modes of zero frequencies and the unstable mode in the transition state are excluded.

3. RESULTS We first consider the diffusion of an isolated epoxy oxygen at different temperatures. Figure 1 shows the calculated diffusion path at 300 K. The transition state (TS) is identified to be the structure in which the oxygen adatom is on top of a C atom while slightly leaning toward the center of the graphene hexagonal cell. The O−C bond length of the transition state is 1.42 Å. The O−C−C bond angles are 101.0°, 101.1°, and 107.5°, respectively, as listed in Table 1. These values are in good agreement with previous studies21 where thermal electronic excitations were not considered. Also, we find the activation energy (ΔFel) at 300 K is 0.72 eV, very close to the previous reported value of 0.73 eV.21,23,25 Apparently at 300 K thermal electronic excitations are weak and their influence to oxygen diffusion is negligible. The activation energy remains close to its room temperature value up to 1000 K. Then it starts to decrease nearly linearly with T, reaching 0.47 eV at 3000 K, as shown in Figure 2a. The reduction of activation energies at high T is accompanied by a flattening of free energy profiles, as

shown Figure 2b. In particular, a small dip emerges in the center region of the free energy profile when T equals 2500 K, indicating the formation of a very shallow metastable state (MS) at high temperatures. The local stability of MS can be verified by structural optimization and vibrational frequencies analysis. Its structure is very close to that of the TS found at 300 K, with a O−C bond length of 1.41 Å and O−C−C bond angles of 103.2°, 104.1°, and 104.9°. This MS easily transforms into the stable epoxy structure by overcoming a negligibly small energy barrier of 0.01 eV. The associated TS has a structure in which the O adatom is more tilted toward the neighboring C atom. Note that although the TSs identified at different temperatures exhibit noticeable structural variations, the initial and final epoxy structures remain nearly identical from 300 to 3000 K. The reduction in activation energies at high T is evident for isolated epoxy oxygen atoms as well as epoxy clusters. With a second epoxy oxygen occupying the next-nearest neighbor site (marked by dashed circles in Figure 1), the activation energy of D

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A epoxy diffusion decreases from 0.66 eV at 300 K to 0.51 eV at 2500 K, shown as filled squares in Figure 2a. To the best of our knowledge, experimental values for the activation energy of O diffusion on graphene are not available, thus a direct comparison of our results with experiment is not possible. There are instead some measurements of O diffusivity on graphite,43 but the proposed diffusion mechanism for that case is different from the one investigated here, as it involves the diffusion of O atoms intercalated between the C sheets. As a result, the activation energy resulting from these measurements at 923 K is quite large (1.5 eV). To understand the origin of this reduction in activation energies, we evaluated the electronic density of states (DOS) of the IS epoxy structure (O at the bridge site) and the TS (O at the top site), as shown in Figure 3. We see that the electronic DOS is strongly affected by the location of the oxygen adatom. When the O adatom is at the bridge site, the energy of its p orbital is more than 2 eV below the Fermi level. When the O adatom is at the top site, its p orbital is just below the Fermi level. Because thermally excited electrons are mostly in the energy range of ∼kBT near the Fermi level, the TS is more susceptible to thermal electronic excitations than the IS epoxy structures. Stronger thermal electronic excitations are associated with larger electronic entropies, resulting in the observed flattening of the free energy profile and reduction in activation energies. Larger electronic entropies also lead to the appearance of the metastable states at T ≥ 2500 K, as seen in Figure 2b. Note the large contrast between the electronic structures of the IS and TS is likely be the key factor here. If the IS and TS have similar electronic structures, then the effects may get canceled, even if thermal electronic excitations are strong for both the IS and TS. Indeed, previous studies44 have included thermal electronic excitations in computing the self-diffusion coefficient of Ni up to 1700 K and found their effects are negligible. We now consider how the vibrational frequencies of the initial state νi, transition state ν*j , and the prefactor ν0 may be affected by thermal electronic excitations. As shown Figure 4, νi is not sensitive to thermal electronic excitations: from 300 to 3000 K νi varies at most 2%. This is not surprising in light of the small electronic DOS near the Fermi level shown in Figure 3a. In contrast, some ν*j can change as much as 20%, even though most of the vibrational frequencies remain nearly identical. The largest frequency shifts are seen in modes whose polarization vectors are predominated by the displacements of the oxygen adatom, indicating that the C−O bonding is affected the most by thermal electronic excitations. The changes in νi and νj* result in a smaller vibrational prefactor ν0: from 2.1 × 1014 Hz at 300 K to 1.5 × 1014 Hz at 3000 K. We have shown thermal electronic excitations reduce the activation energy from 0.72 eV at 300 K to 0.47 eV at 3000 K. If ν0 is assumed to be independent of temperature, the reduction in activation energy would increase the diffusion coefficient at 3000 K by a factor of 2.6. Once the temperature dependences of both ν0 and activation energy are taken into account, the diffusion coefficient at 3000 K gets enhanced by a factor of 1.9. We conclude the net effect of thermal electronic excitations is to enhance the diffusivity of epoxy oxygen.

electron density of −3.82 × 1013 cm−2 gives an activation energy of 0.45 eV.21 Also, the variations of free energy profiles with respect to electron density closely resemble those shown in Figure 2b. This resemblance may originate from the fact that thermal electronic excitations prompt electrons to occupy conduction bands, which has similar effects on the activation energy as electron doping. Still, there is a big difference. Electron doping takes place at room temperature. The thermal energy kBT equals 0.026 eV, considerably smaller than the activation energy for diffusion. As a result, decreasing the activation energy from 0.73 to 0.45 eV enhances the diffusion rate at room temperature by a factor of 5.2 × 104. In contrast, significant reduction in the activation energy due to thermal electronic excitations occurs at high T, where the thermal energy is about 1/3 to 1/2 of the activation energy. Accordingly, the enhancement in the diffusion rate is less prominent compared to the room temperature electron doping case. From the previous section we see oxygen adatoms are quite mobile on graphene at high T. The diffusion coefficient at 1500 K is on the order of 10−5 cm2/s, comparable to that of H2O molecule in liquid water at ambient conditions (2.4 × 10−5 cm2/s). In contrast, at 300 K the diffusion coefficient is on the order of 10−16 cm2/s. The increased mobility of oxygen adatoms at high T may have interesting thermodynamic implications. In the following we present a preliminary analysis on how the equilibrium coverage of oxygen adatoms may be affected by the mobility of O adatoms. The equilibrium coverage Θ is a basic surface thermodynamic property that relates the amount of surface adsorbates with that of free molecules in the gas phase (quantified by partial pressure). For localized adsorbates Θ is determined by the Langmuir formula45

4. DISCUSSION In some respects, the effect of thermal electronic excitations on epoxy diffusion is similar to that of electron doping. As pointed out by Suarez et al.,21 the activation energy of epoxy diffusion decreases near linearly as the electron density increases. An

where A is the effective surface area corresponding to one adsorption site.45 For an O adatom, the adsorption energy on graphene is 1.9 eV, which we find is not affected by thermal electronic excitations. As each graphene unit cell contains three bridge

Θ(P ,T ) =

P P + P0(T )

(7)

where P is partial pressure of the adsorbate molecules in the gas phase. P0 is the pressure at which the coverage equals 1/2, determined as ⎛ E ⎞ ⎛ 2πmkBT ⎞3/2 kBT ⎟ exp⎜ − ad ⎟ P0(T ) = ⎜ 2 2 ⎝ h ⎠ zp zv ⎝ kBT ⎠

(8)

In eq 8, zp and zv are the vibrational partition functions in the planar and vertical directions, respectively. They are expressed in terms of vibrational frequency ν as 1/z = 2 sinh(hν/2kBT). m is the mass of the adsorbate molecule. Ead is the adsorption energy. For mobile adsorbates the situation becomes much more complex, as the adsorbates may migrate in various preferred directions with a wide range of possible rates. The simplest approximation is to model the planar movements of adatoms as a 2D ideal gas and treat the vertical movements as vibrations. As such, an effective equilibrium coverage may be evaluated as Θ(P ,T ) =

E

⎛E ⎞ PA h2 z v exp⎜ ad ⎟ kBT ⎝ kBT ⎠ 2πmkBT

(9)

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 5. Equilibrium coverage of O adatoms on the graphene surface predicted by treating O adatoms as localized (red) or mobile (blue) adsorbates: (a) PO equals 1 Pa = 1 × 10−5 atm; (b) PO2 equals 1 atm.

sites, the effective surface area of one adsorption site A is 1.75 Å2. The partition functions are evaluated approximately by assuming a frequency of 10 THz for both planar and vertical vibrations. For the partial pressure PO, we consider two cases: (i) PO is fixed at 1 Pa = 1 × 10−5 atm, and (ii) the atomic oxygen is in thermal equilibrium with O2, whose partial pressure PO2 is fixed at 1 atm. In the latter case, PO is determined at each temperature according to the Gibbs energy change associated with the reaction O2(g) = 2O(g).46 The equilibrium coverages corresponding to these two cases are shown in Figure 5. We see that, for case i, Θ decreases rapidly with temperature, as higher percentages of adatoms have sufficient kinetic energy to leave the surface. In contrast, for case ii, Θ increases with temperature. This is because more O2 molecules get dissociated into atomic oxygen at high T. The increase in the total number of O atoms has a dominant effect on Θ with respect to the rise in kinetic energies. In both cases, treating O atoms as mobile adsorbates leads to a coverage 10 times higher than the one predicted for localized adsorbates. We remark that these results are based on simplified models for the dynamics of O adatoms, such as localized adsorbates or 2D ideal gas. To evaluate the equilibrium coverage more accurately, one needs to go beyond these models and acquire a precise description of the dynamics of oxygen adatoms. The quantitative analysis of the thermal electronic effects presented in this study may facilitate such endeavors in the future.

particular interest is mineralogy, where extensive theoretical modelings are now underway to elucidate the diffusion processes in mantle minerals and silicate melts.47,48 We anticipate thermal electronic effects can be significant in some of these mineralogical processes, which take place under the high temperature and pressure conditions of the Earth’s deep interior.

■

AUTHOR INFORMATION

Corresponding Authors

*T. Sun. E-mail: [email protected]. *S. Fabris. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS S.F. is thankful for support from CAS President’s International Fellowship Initiative (PIFI) for visiting scientists. T.S. is supported by Ministry of Science and Technology of China grant No. 2014CB845905 and National Science Foundation of China grant No. 41474069. S.F. is also supported by Regione FVG with the project NANOCAT and by the MIUR PRIN project 20105ZZTSE. We thank Changru Ma and Dong-Bo Zhang for useful discussion. Some of the calculations were performed on the Tianhe-1A supercomputer at National Supercomputer Center of China (NSCC) in Tianjin.

■

5. CONCLUSIONS We have investigated theoretically the diffusion of epoxy oxygen on the graphene surface from 300 to 3000 K. The activation energy exhibits little temperature dependence at T ≤ 1000 K. It then starts to decrease nearly linearly with T because of the large contrasts in the electronic structures of the initial and transition states. Compared to the initial states, the transition states are more susceptible to thermal electronic excitations and therefore have larger electronic entropies. The reduction in activation energy leads to appreciable enhancement in the diffusivity of epoxy oxygen. A decrease in the vibrational prefactor, also caused by thermal electronic excitations, does not alter this trend. These findings highlight the importance of thermal electronic effects in the process of oxygen diffusion and may facilitate future works on accurately describing the dynamic of O adatoms at high T. Such electronic effects may also be relevant to other high temperature dynamical processes, insofar as the process is accompanied by drastic changes in the electronic structure. One area of

REFERENCES

(1) Smith, I. W. The intrinsic reactivity of carbons to oxygen. Fuel 1978, 57, 409−414. (2) Leistner, K.; Nicolle, A.; Berthout, D.; da Costa, P. Kinetic modelling of the oxidation of a wide range of carbon materials. Combust. Flame 2012, 159, 64−76. (3) McAllister, M. J.; Li, J.-L.; Adamson, D. H.; Schniepp, H. C.; Abdala, A. A.; Liu, J.; Herrera-Alonso, M.; Milius, D. L.; Car, R.; Prud’homme, R. K.; et al. Single sheet functionalized graphene by oxidation and thermal expansion of graphite. Chem. Mater. 2007, 19, 4396−4404. (4) Dreyer, D. R.; Park, S.; Bielawski, C. W.; Ruoff, R. S. The chemistry of graphene oxide. Chem. Soc. Rev. 2010, 39, 228. (5) Ngo, T.; Snyder, E. J.; Tong, W. M.; Williams, R. S.; Anderson, M. S. O atom etching of graphite in low earth orbit. Surf. Sci. 1994, 314, L817−L822. (6) Bianchi, D.; Nasuti, F.; Martelli, E. Coupled analysis of flow and surface ablation in carbon-carbon rocket nozzles. J. Spacecr. Rockets 2009, 46, 492−500. (7) Luo, X.; Jean-Charles, R.; Yu, S. Effect of temperature on graphite oxidation behavior. Nucl. Eng. Des. 2004, 227, 273−280. F

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Could subsurface hydrogen play an important role? J. Chem. Phys. 2006, 124, 044706. (31) Kumar, P. V.; Bardhan, N. M.; Chen, G.-Y.; Li, Z.; Belcher, A. M.; Grossman, J. C. New insights into the thermal reduction of graphene oxide: Impact of oxygen clustering. Carbon 2016, 100, 90− 98. (32) Koizumi, K.; Boero, M.; Shigeta, Y.; Oshiyama, A. Atom-scale reaction pathways and free-energy landscapes in oxygen plasma etching of graphene. J. Phys. Chem. Lett. 2013, 4, 1592−1596. (33) Negreiros, F. R.; Camellone, M. F.; Fabris, S. Effects of thermal fluctuations on the hydroxylation and reduction of ceria surfaces by molecular H2. J. Phys. Chem. C 2015, 119, 21567−21573. (34) Doltsinis, N.; Marx, D. First principles molecular dynamics involving excited states and nonadiabatic transitions. J. Theor. Comput. Chem. 2002, 1, 319−349. (35) Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G. L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (36) Ultrasoft pseudopotentials C.pbe-rrkjus.UPF and O.pberrkjus.UPF, from the QUANTUM ESPRESSO distribution were used to describe the electron-ion interaction. (37) Sun, T.; Fabris, S.; Baroni, S. Surface precursors and reaction mechanisms for the thermal reduction of graphene basal surfaces oxidized by atomic oxygen. J. Phys. Chem. C 2011, 115, 4730−4737. (38) Pozzo, M.; Alfè, D.; Lacovig, P.; Hofmann, P.; Lizzit, S.; Baraldi, A. Thermal expansion of supported and freestanding graphene: lattice constant versus interatomic distance. Phys. Rev. Lett. 2011, 106, 135501. (39) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865. (40) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188−5192. (41) Jonsson, H.; Mills, G.; Jacobsen, K. W. In Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds.; World Scientific: Hoboken, NJ, 1998. (42) Baroni, S.; de Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Phonons and related crystal properties from density-functional perturbation theory. Rev. Mod. Phys. 2001, 73, 515−562. (43) Yang, R. T.; Wong, C. Kinetics and mechanism of oxidation of basal plane on graphite. J. Chem. Phys. 1981, 75, 4471. (44) Hargather, C. Z.; Shang, S.-L.; Liu, Z.-K.; Du, Y. A firstprinciples study of self-diffusion coefficients of fcc Ni. Comput. Mater. Sci. 2014, 86, 17−23. (45) Hill, T. L. An Introduction to Statistical Thermodynamics; Addison-Wesley: Reading, MA, 1960. (46) Barin, I. Thermochemical Data of Pure Substances, 3rd ed.; VCH: Weinheim, 1995. (47) Ammann, M. W.; Brodholt, J. P.; Dobson, D. P. Simulating diffusion. Rev. Mineral. Geochem. 2010, 71, 201−224. (48) Dominguez, G.; Wilkins, G.; Thiemens, M. H. The Soret effect and isotopic fractionation in high-temperature silicate melts. Nature 2011, 473, 70−73.

(8) El-Genk, M. S.; Tournier, J.-M. P. Development and validation of a model for the chemical kinetics of graphite oxidation. J. Nucl. Mater. 2011, 411, 193−207. (9) Stevens, F.; Kolodny, L. A.; Beebe, T. P. kinetics of graphite oxidation: monolayer and multilayer etch pits in HOPG studied by STM. J. Phys. Chem. B 1998, 102, 10799−10804. (10) Sima-Ella, E.; Yuan, G.; Mays, T. A simple kinetic analysis to determine the intrinsic reactivity of coal chars. Fuel 2005, 84, 1920− 1925. (11) Nicholson, K. T.; Minton, T. K.; Sibener, S. J. Spatially anisotropic etching of graphite by hyperthermal atomic oxygen. J. Phys. Chem. B 2005, 109, 8476−8480. (12) Larciprete, R.; Fabris, S.; Sun, T.; Lacovig, P.; Baraldi, A.; Lizzit, S. Dual path mechanism in the thermal reduction of graphene oxide. J. Am. Chem. Soc. 2011, 133, 17315−17321. (13) Larciprete, R.; Lacovig, P.; Gardonio, S.; Baraldi, A.; Lizzit, S. Atomic oxygen on graphite: chemical characterization and thermal reduction. J. Phys. Chem. C 2012, 116, 9900−9908. (14) Li, J.-L.; Kudin, K. N.; McAllister, M. J.; Prud’homme, R. K.; Aksay, I. A.; Car, R. Oxygen-driven unzipping of graphitic materials. Phys. Rev. Lett. 2006, 96, 176101. (15) Carlsson, J. M.; Hanke, F.; Linic, S.; Scheffler, M. Two-Step mechanism for low-temperature oxidation of vacancies in graphene. Phys. Rev. Lett. 2009, 102, 166104. (16) Orrego, J. F.; Zapata, F.; Truong, T. N.; Mondragón, F. Heterogeneous CO2 evolution from oxidation of aromatic carbonbased materials. J. Phys. Chem. A 2009, 113, 8415−8420. (17) Paci, J. T.; Upadhyaya, H. P.; Zhang, J.; Schatz, G. C.; Minton, T. K. Theoretical and experimental studies of the reactions between hyperthermal O(3P) and graphite: graphene-based direct dynamics and beam-surface scattering approaches. J. Phys. Chem. A 2009, 113, 4677−4685. (18) Radovic, L. R. Active sites in graphene and the mechanism of CO2 formation in carbon oxidation. J. Am. Chem. Soc. 2009, 131, 17166−17175. (19) Xu, S. C.; Chen, H.-L.; Lin, M. C. Quantum chemical prediction of reaction pathways and rate constants for the reactions of Ox(x = 1 and 2) with pristine and defective graphite (0001) surfaces. J. Phys. Chem. C 2012, 116, 1841−1849. (20) Vineyard, G. H. Frequency factors and isotope effects in solid state rate processes. J. Phys. Chem. Solids 1957, 3, 121−127. (21) Suarez, A. M.; Radovic, L. R.; Bar-Ziv, E.; Sofo, J. O. GateVoltage control of oxygen giffusion on graphene. Phys. Rev. Lett. 2011, 106, 146802. (22) Radovic, L. R.; Silva-Tapia, A. B.; Vallejos-Burgos, F. Oxygen migration on the graphene surface. 1. Origin of epoxide groups. Carbon 2011, 49, 4218−4225. (23) Radovic, L. R.; Suarez, A.; Vallejos-Burgos, F.; Sofo, J. O. Oxygen migration on the graphene surface. 2. Thermochemistry of basal-plane diffusion (hopping). Carbon 2011, 49, 4226−4238. (24) Barinov, A.; Malcioǧlu, O. B.; Fabris, S.; Sun, T.; Gregoratti, L.; Dalmiglio, M.; Kiskinova, M. Initial stages of oxidation on graphitic surfaces: photoemission study and density functional theory calculations. J. Phys. Chem. C 2009, 113, 9009−9013. (25) Sun, T.; Fabris, S. Mechanisms for oxidative unzipping and cutting of graphene. Nano Lett. 2012, 12, 17−21. (26) Kawai, T.; Miyamoto, Y. Simulations for the formation dynamics and electronic states of carbon nano materials: Diffusion and alignment of oxygen atoms on graphene. Curr. Appl. Phys. 2011, 11, S50−S54. (27) Mehmood, F.; Pachter, R.; Lu, W.; Boeckl, J. J. Adsorption and diffusion of oxygen on single-layer graphene with topological defects. J. Phys. Chem. C 2013, 117, 10366−10374. (28) Mermin, N. D. Thermal properties of the inhomogeneous electron gas. Phys. Rev. 1965, 137, A1441. (29) Levine, I. N. Physical Chemistry, 6th ed.; McGraw Hill: New York, 2009. (30) Henkelman, G.; Arnaldsson, A.; Jonsson, H. Theoretical calculations of CH4 and H2 associative desorption from Ni(111): G

DOI: 10.1021/acs.jpca.6b00423 J. Phys. Chem. A XXXX, XXX, XXX−XXX