1 DFT and kinetics study of O / O 2 mixtures reacting over a graphite ( 0001 ) basal surface

Spin-polarised density functional calculations were used to investigate the interaction of atomic and molecular oxygen on the basal graphite surface at several atomic coverages. Two carbon layers were enough to represent the graphite surface. Oxygen atoms bind mainly over C-C bridge sites forming an epoxide-like structure with a two carbon puckering and with adsorption energies in the 0.95-1.28 eV range, depending on the atomic coverage. Molecular oxygen only shows a very weak physisorption. Atomic adsorption and diffusion along with atomic recombination via Eley-Rideal and Langmuir-Hinshelwood mechanisms were studied. All surfaces processes were activated with energy barriers that decreased for lower atomic coverages. Relaxation effects were non-negligible. A microkinetic model with six elementary surface processes was proposed to see the overall behavior of several initial O/O2 mixtures flowing over a graphite surface at 300-1000 K. Thermal rate constants were derived from Density Functional Theory (DFT) data and standard Transition State Theory. A very low steady-state atomic coverage (qO < 0.5 %) was predicted and also very low atomic recombination coefficients were observed (gO < 5x10-4). The Eley-Rideal together with the adsorption and desorption processes were much more important than the LangmuirHinshelwood reaction.


Introduction
Carbon-based materials have received considerable attention for their use in spacecraft components (e.g., reinforced C is used as a thermal protection system (TPS) in the nose of Space Shuttle) due to their light weight and high strength. In the low Earth orbits (LEO), 160-2000 km above Earth's surface, atomic oxygen is one of the main present species and can collide with the spacecraft surface at very high translational energies (i.e., about 4.5±1.0 eV due to the orbital velocity of around 7.5 km/s) producing an important etching, hence degrading any C-based TPS [1]. Moreover, during hypersonic reentry flights into Earth's atmosphere (e.g., Space Shuttle) not only atomic oxygen but also molecular oxygen processes can be relevant to understand properly the huge heat transfer and the TPS behavior under these extreme conditions [2]. Therefore, a good understanding of all elementary processes involving atomic and molecular oxygen with carbonbased materials is essential for space technology. Thus, in principle, several heterogeneous processes could compete: a) O adsorption, b) O2 dissociative or non-dissociative adsorption, c) O recombination via an Eley-Rideal (ER) mechanism (i.e., O(g) + Oad ® O2 (g)), d) O recombination via a Langmuir-Hinshelwood (LH) mechanism (i.e., Oad + Oad ® O2 (g)) and e) C etching producing CO or CO2, where ad labels adsorbed species and g gas species. Another less important elementary processes could also be added.
There are abundant theoretical studies about the interaction of atomic and molecular oxygen with graphite [3][4][5][6][7], graphene [8][9] and carbon nanotubes [4,[10][11] covering different aspects but without giving a global approach about the O/O2 reactivity over these carbon surfaces. They are mainly concerned with adsorption studies over several type of surfaces (e.g., basal, zig-zag, armchair,..) for different coverages. Thus, there was shown that atomic oxygen becomes mainly chemically adsorbed on a bridge position between two adjacent C atoms of graphite or graphene, forming an epoxy group, with adsorption energies in the range of 1.5-3.2 eV for a basal graphite surface depending of the oxygen coverages and the level of the calculations [3][4][5][6]. The molecular oxygen seems to be only physically adsorbed on basal graphite surfaces with low adsorption energies of about 0.12±0.01 eV at 1 ML coverage as it has been shown using thermal desorption spectroscopy on highly oriented pyrolytic graphite (HOPG) [12]. Common Density Functional Theory (DFT) studies are not accurate enough to provide reliable adsorption energies for O2 over graphite [3,4,7], due to the importance of the missing van der Waals (vdW) contributions.
Nevertheless, a suitable vdW treatment could be done for instance by adding pairwise interatomic C6R -6 terms, where C6 coefficients can be derived from ground-sate electron density for molecules and solids [13]. The inclusion of vdW interactions reproduce quite well the O2 physisorption over graphite [14][15]. The interaction of O2 with defective edge sites on graphite surfaces presents high adsorption energies (2.6-6.2 eV [7]) in a similar way as for atomic oxygen (4.5-5.7 eV [6]).
Moreover, molecular oxygen reactively adsorbs on edge graphite surfaces at 300 °C and 1 atm pressure (the perfect basal surface is hardly active [16]), being the energy barrier for O2 dissociation larger and for CO formation lower when the oxygen surface coverage diminishes [17].
Up to our knowledge there are not theoretical studies about the ER and LH atomic reactions or on the atomic diffusion over graphite surfaces. Characterization of these reaction pathways is missing and it is very important from a kinetic point of view. Only very recently, the migration (diffusion) barrier for O on a graphene sheet has been investigated [8], indicating an energy barrier of 0.58 eV at a coverage of 16.7 %, which would imply a noticeable mobility of isolated adsorbed atoms.
Apart from the mentioned experimental studies about O2 adsorption/desorption on graphite, there are as well a large number of experimental works concerned mainly with the interaction of hyperthermal atomic oxygen on HOPG surfaces, to understand the erosion and degradation of these surfaces under large collision energies [1,[18][19], forming CO and CO2 gas species.
In this work, we study the main elementary processes involving O and O2 over an initially clean graphite basal surface (0001) at the same DFT level, characterizing not only the minima but also the transition states (TS) for the different processes at several oxygen coverages. These data are used in a proposed microkinetic model based on calculated thermal rate constants derived from Transition State Theory, which allows a kinetic study for different temperatures and total and partial pressures (i.e., several O/O2 mixtures). We also determine the atomic oxygen recombination coefficient gO (T, P) (0 ≤ g O ≤ 1), which is defined as the total probability of atomic oxygen recombination per atomic oxygen surface collision; gO values are necessary in standard computer fluid dynamics (CFD) simulations [2] of hypersonic flights aiming to evaluate their aerodynamic heating.

Computational details
The periodic DFT calculations presented in the current work have been carried out using the Vienna ab initio simulation package (VASP) code [20][21][22][23], which uses plane wave basis set. The calculations are mainly based on the generalized gradient approximation (GGA) functional revised Perdew-Burke-Ernzerhof functional (RPBE) [24]. However, some initial tests were also done with Perdew-Burke-Ernzerhof (PBE) [25][26] and Perdew-Wang 91 (PW91) [27][28] functionals. The projector-augmented-wave (PAW) technique within the frozen core approximation has been used to describe the electron-core interaction [29][30]. An energy cutoff of 550 eV has been used in the plane-wave expansion after several tests to confirm the energy convergence of the calculations.
Integration over the Brillouin zone was performed by using a 14x14x1, 11x11x1 and 6x6x1 k-point meshes by means of the Monkhorst-Pack method [31] for slab calculations using several supercells, 1x1, 2x2 and 3x3, respectively. For bulk calculations a 14x14x4 mesh was used. The atomic positions were partially (e.g., for gas/rigid slab system) or fully relaxed by using a quasi-Newton algorithm until the energy decreased below 10 -5 eV, which usually corresponded to forces lower than 10 -2 eV/Å. The energy convergence in the electronic self-consistent cycles was maintained below 10 -6 eV for all geometrical optimizations and vibrational frequency calculations. The vacuum used in the slab model was around 15 Å, which was large enough to prevent significant interactions between periodic images.
Due to the important role that spin plays in both atomic and molecular oxygen (i.e., triplet in their ground states), all calculations were spin-polarised.
Several adsorption sites have been studied for O and O2 species (see next section for details).
Moreover, we have calculated minimum energy paths for the above-mentioned elemental processes (e.g., ER reaction, atomic diffusion,.) by means of the Nudged Elastic Band method (NEB) [32][33] or using more simple reaction coordinates (e.g., ZO distance over a C-C bridge site). Once determined the optimal geometry for each stationary point, we have calculated the Hessian matrix and its corresponding harmonic vibrational frequencies (ni) to verify if they were true minima or transition states. We have also checked the connection among these stationary points. In general we have fixed the optimized slab geometry, although some calculations with full or partial slab relaxation were also performed. Energy barriers (∆E ¹ ) for the elementary processes were mostly computed using the corresponding supermolecule energy, calculated for reactants far enough to give a converged asymptotic energy.
The adsorption energy (Ead) was calculated according to Eq. 1, Ead = (EX + Eslab) -EX/slab (1) where EX is the energy of an isolated O or O2 species at their electronic ground state, Eslab is the energy for the relaxed clean slab (EX+ Eslab calculated as EX+slab) , and EX/slab is the energy of the adsorbed-slab system. Therefore, a positive adsorption energy means that this site is more stable than the X + slab asymptote; on the contrary, a negative adsorption energy corresponds to a less stable one.

Bulk graphite properties
Graphite is formed by a layered structure of six-member carbon rings, stacked in an ABAB sequence with a weak interlayer binding compared to the much stronger binding within the layers.
The experimental [36] lattice parameters for graphite structure are reported in table 1.
We have calculated the graphite lattice parameters using several functionals (PW91, PBE and RPBE). The value of optimum a=b lattice parameter is similar and close to the experimental value [36], although the optimum c lattice parameter is overstimated in all cases. The usual fail of GGA approximation to predict the c-axis lattice parameter and the small interlayer binding energy is well known, along with the better LDA description of c ( where NC is the number of C atoms in the cell, EC(g) is the energy of C (triplet) atom and Ebulk(eq.) is the energy corresponding to the equilibrium bulk structure with optimum a and experimental c lattice parameters.  [40,41]. The correct calculation of the usual bulk modulus (Bo) for highly anisotropic materials imply the use of an uniaxial compression [43] instead of an uniform compression, which could produce values an order of magnitude larger than the experimental ones. Hence, the poor description of the c lattice parameter does not allow a proper calculation of Bo with the present DFT/GGA calculations.
The calculated B'iso are much lower than the experimental values estimated from the initial slopes of elastic constants with pressure at atmospheric pressure and 295 K [40] as shown in table 1.
In general, we conclude that GGA/RPBE calculations give a reasonable good description of bulk graphite properties, although it is necessary to use the experimental c lattice parameter. This will not be a problem for the reliability of the atom or molecule-slab calculations of the next section as it has been also assumed in similar previous studies [6].

Adsorption and diffusion of oxygen
We have studied the adsorption of atomic oxygen on the basal (0001) graphite surface, with a slab formed by one or two carbon layers. The results are almost identical regardless of the number of layers of the slab model, possibly due to the very weak interlayer interaction. This fact was also verified in previous O and K adsorption studies over graphite [3] by using 1, 2 or 3 layers. Hence that another similar O/graphite studies have simulated the graphite as a single sheet of C atoms [4].
Three high-symmetry different sites have been mainly considered ( Fig. 1): on top of a surface C atom, which is either above another C atom (T1) or above an hexagon hollow (T2) from the adjacent plane, and over a bridge between the two nearest-neighbor carbon atoms (B). We checked also that the hollow of the hexagonal ring of the first layer did not present any significant adsorption in agreement with previous studies [5,6]. The atomic coverage effect has been studied considering the adsorption of one O atom over three different unit cells: 1x1, 2x2 and 3x3, namely 50 %, 12.5 % and 5.5 % respect the C atoms (i.e.,T1 and T2 sites). These coverages would be doubled taking into account only the most stable adsorption sites (i.e., true minima), which correspond to the bridges (i.e., B sites); when an O is over a bridge site the closest free bridges become unavailable for another O adsorption because every C atom can make only four bonds. We will indicate this coverage throughout the text giving both values (i.e., % respect C or top sites / % respect bridge sites). T1 and T2 stationary points mainly present two imaginary frequencies, which seem to connect with the adsorption bridge minimum and with a diffusion transition state (see below for more details), although the slab relaxation can affect as well the character of these frequencies (e.g., reducing to one imaginary value for T1 site at 12.5/25% coverage). Table 3 summarizes the effect of slab relaxation in the bridge minimum for several O coverages. This geometry relaxation corresponds to a stretching of C-C bond along with a reduction of both C-O bonds, with also a puckering of the two carbon atoms beneath the adatom. This behavior is in good agreement with a similar DFT study by means of GGA/PW91 with ultrasoft pseudopotentials [4], which shows dO-C = 1.472 Å and dC-C = 1.500 Å, compared to the graphite dC-C = 1.420 Å for a 5.5/11.1 % of coverage. However, a larger adsorption energy (1.91 eV) was calculated in that work. We observe as well that the slab relaxation implies that the adsorption energy tend to be larger for low O coverages (table 3), conversely to the behavior found for rigid slab (table 2). The small changes in Ead (i.e., ± 0.3 eV) produced by coverage and slab relaxation effects seem to be counterbalanced.
The formation of this epoxide-like (COC) ring structure on bridge sites has been observed in theoretical and experimental studies on graphene [9] and graphite [44][45][46]  A Bader charge analysis of the epoxy group shows an important charge transfer between the C atoms of the surface and the adsorbed oxygen (e.g., qO = -0.99, qC = +0.41, +0.58 at 50/100%).
This arrangement of negatively charged oxygen could avoid the nucleophilic attack on C atoms, explaining the relative chemical stability of the epoxide groups in graphite oxides, as it was earlier suggested [44]. Table 4 summarizes the properties of transition states for O adsorption and diffusion processes over the basal (0001) graphite surface at several coverages for rigid slab approximation. Adsorption process is clearly activated as it was also demonstrated for atomic hydrogen over graphene (∆E ¹ ~ 0.14 eV at low coverages [47]) or graphite (∆E ¹ = 0.23 eV at 12.5/25 % coverage [48] The analysis of the imaginary frequency in this transition state confirms that principally the perpendicular movement of oxygen respect the surface (i.e., z O coordinate) corresponds to the adsorption process. This TS has mainly a triplet character due its proximity to the O gas asymptote.
The diffusion of one absorbed O over a bridge to its closest bridge is also an activated process that follows a path roughly shown in Fig. 2, with a TS mainly with a singlet character as in the adsorption minimum. The imaginary frequency involves mostly a parallel movement to the surface as could be expected for this process. The corresponding energy barrier decreases significantly when the coverage is reduced, likewise as for the adsorption process. However, the slab relaxation increases considerably this energy barrier from 0.38 to 0.73 eV at 12.5/25 % (table 4)  We have also studied with the present DFT approach the adsorption of molecular oxygen on the basal surface of graphite. Thus, perpendicular and parallel initial configurations with its center of mass over the different adsorption sites (Fig. 1) were considered, optimizing all degrees of freedom including also the slab relaxation by using a 2x2 unit cell. There were found very low adsorption showed that triplet state was more stable than the singlet one. However, it is clear that is necessary an accurate description of the vdW interactions to describe properly the small experimental physisorption energy (0.12±0.01 eV at 1 ML coverage [12]).

Recombination reactions of atomic oxygen
We have found that the atomic recombination processes (i.e., ER and LH reactions) over the basal (0001) graphite surface are activated processes although their energy barriers considerably decrease with coverage increment, as can be seen in table 5. The ER transition state corresponds to a perpendicular O-O orientation over the bridge site (Fig. 2), being the imaginary frequency also a perpendicular movement that gives place to O2(g). Contrarily, the LH transition state presents a parallel configuration (Fig. 2) with a similar O-O distance as for ER TS (e.g., for a (2x2) unit cell); in this case the imaginary frequency has a mixture of both parallel and perpendicular movements, as could be expected. There is not experimental information to compare with these energy barriers.
However, a similar DFT study with H over graphite (0001) surface for a 12.5/25 % coverage shows that at this low coverage the ER energy barrier would be equally very low ( ≤ 0.05 eV), as could be observed from cross sections vs. incident energy plot of quantum scattering calculations [48]) .
However, the authors assumed a puckering of the C atom below the adsorbed H on top (z = 0.36 Å), which possibly reduces more this ER energy barrier. Conversely, the relaxation of the slab in the present study could increase this energy barrier for O as occurred for the diffusion process ( The LH reaction should be much less important than ER (or hot atom) reaction as can be concluded from the large differences in their energy barriers as will be confirmed in the next section. This fact has also been previously found out for these O recombination reactions over bcristobalite [49].

A microkinetic model for O/O2 mixtures over graphite
We propose a first microkinetic model to estimate the overall effect of all studied heterogeneous processes involving different O/O2 mixtures over a graphite (0001) surface. Six surface elementary processes are included: Chemisorption (k1) and desorption (k-1) of atoms: Langmuir-Hinshelwood reaction (k2) and dissociative molecular adsorption (k-2) Eley-Rideal reaction (k3) and another dissociative molecular adsorption channel (k-3) It is worth noting that a steady-state (ss) is reached rapidly for different initial conditions (i.e., PO, PO2 and T) with a final low atomic coverage in all studied cases. Fig. 5 shows qO at several temperatures for two representative partial pressures although another combinations were also explored. Thus, very low values of a final steady-sate atomic surface coverage (< 10 -3 ) are observed between 300 and 1000 K for both mixtures, with larger values for the higher O partial pressure as could be predictable. This chemisorption of atomic oxygen on the (0001) graphite surface has been also revealed experimentally [16], with a low level of adsorption on the basal surface, mainly enhanced by the concentration of defect sites [17]. Moreover, these low final qO values support as well our use of DFT data for low coverage in the calculated rate constants.
The shape of qO as a function of the temperature shows a maximum at different temperatures (i.e., 400 and 500 K), which depends basically on the initial PO value as could be expected due to at high temperatures the desorption becomes the most important process.
We have also estimated the atomic oxygen recombination coefficient gO from the following equation, which includes only the process that directly involve atomic oxygen, namely adsorption, ER and desorption, once achieved the steady-state ( ), and where ZO is the number of O collisions per unit area and unit time with a planar surface (Hertz-Knudsen relation), where mO is the atomic oxygen mass.
The gO coefficient shown in Fig. 5 almost mimics the shape of the qO plot. An analysis of its three contributions establishes clearly that at low temperatures adsorption and ER reaction are the most important processes with similar contributions, while at high temperatures this occurs for adsorption and desorption processes although in this case their rates are counterbalanced (i.e., r1 = r-1), hence that gO finally decreases as a consequence of the lower values. Although the final steady-state coverage depends as well on LH process (Eq. 7), in the present case its contribution is unimportant due to its low rate constants.
The comparison of the gO(T) shape over graphite respect to other materials (e.g., metals [50], alumina [51], silica [52],..) is similar but with much minor values. Moreover, the gO decline is produced at much lower temperatures for graphite, possibly as a consequence of its lower adsorption energies. The introduction of new rate constants based on DFT data with slab relaxation and /or higher coverages do not produce significant changes respect the present results as we have checked with some additional calculations.
The presence of molecular oxygen (i.e., PO2 ¹ 0) in the initial flux is negligible within the studied temperature range (300-1000 K) for both gO and qO . This is supported by the large energy barriers (Fig. 3) associated to the two dissociative molecular processes (i.e., very small k-2 and k-3) included into the microkinetic model.
The very low gO coefficients and the large melting point of graphite (4800 ±100 K at a pressure of 10-100 MPa [53]) could seemingly indicate that C-based materials would be a good proposal to be used as TPS for atmospheres containing oxygen. Nevertheless, the oxidative etching of graphite, forming CO(g) and also CO2(g), can degrade this material and would prevent its use even at moderate temperatures. In fact, several experimental studies [17,54,55] have shown that O2 (or O) oxidation can produce finally CO(g) in spite of this endothermic process has a high energy barrier around 3-4 eV, arising from the adsorbed O. Although this carbon etching is mainly enhanced by defect sites, at temperatures larger than 1148 K also occurs in basal plane carbon atoms [54]; HOPG samples at 298-493 K exposed to pulsed hyperthermal beams containing atomic and molecular oxygen yield an important and spatially anisotropic etching [19]. For this reason we have only made calculations with the proposed microkinetic model up to 1000 K. Therefore, at higher temperatures these etching processes and another surface processes should be included, as for instance: In spite of this kind of microkinetic models are frequently used to simulate the heterogeneous chemical kinetics of dissociate airflow impinging different kind of surfaces [56,57], with good agreement in comparison to available experimental data (e.g., total heat flux to the surface, g coefficients,..), more sophisticated Monte Carlo simulations (e.g., N over silica [50], H over graphite [58]) can also be performed including for instance different types of Eley-Rideal mechanisms or the coverage effect through the time evolution of the system.

Conclusions
In this work we have made a periodic DFT study (GGA/RPBE) on the main elementary heterogeneous processes involving atomic and molecular oxygen over a graphite (0001) surface to present a comprehensive view on its reactivity.
DFT calculations on bulk graphite (e.g., lattice parameters, cohesive energy, isotropic bulk modulus,..) and on O2 molecule show a quite good agreement with experimental data. The use of one or two layers for the slab model of graphite in the study of the elementary processes seem to be equivalent. The introduction of slab relaxation is relevant for both adsorption energies and the different energy barriers.
Atomic adsorption is an activated process, whose energy barrier decreases for lower coverages and under slab relaxation. The O adsorption is mainly produced on C-C brige sites with a considerable puckering of these carbons. This epoxide-like (COC) ring structure has also been observed in earlier theoretical and experimental studies on graphene and graphite oxides.
Atomic diffusion between bridge sites is also an activated process, with lower energy barriers at low coverages but with higher values when slab relaxation is allowed, due to the larger Oad stabilization.
Atomic recombination processes via ER and LH mechanisms are also activated processes although their energy barriers decrease with coverage increment, especially for ER reaction, which will dominate the atomic recombination reaction from low to high temperatures.
We have proposed a microkinetic model with 6 elementary steps, which confirms the very small contribution of the LH process and predicts a very low atomic coverage (< 0. 5 %), which along with the small rate constants originates very low atomic recombination coefficients (g O  [38] derived from ∆fH o (0 K) for the process: C(graphite) ® C(g,ideal) e Experimental values [40,41] obtained from elastic constants (Cij) through the formula: , [42] f Experimental value [40] obtained from the pressure derivative of the previous equation   a The slab geometry (two C layers) is fixed at its optimum geometry (dC-C = 1.425 Å) derived in absence of O but with the interlayer distance (d12) fixed to its experimental value (3.347 Å) b s or s' correspond to C-C bridge sites c Energy barriers respect reactants (including zero point energy between parentheses) d The transition state for adsorption is practically over the C-C bridge site and for diffusion is between two non-adjacent C atoms (see Fig. 1) e Harmonic vibrational frequencies of the atomic adsorbate with respect to the rigid substrate f In this case the carbon atoms directly bound to oxygen adatom are allowed to relax in all directions (i.e., zC = 0.314 and 0.297 Å for adsorption and diffusion TSs, respectively) g This TS was derived approximately from the crossing between the singlet and triplet curves 22 The slab geometry (two C layers) is fixed at its optimum geometry (dC-C = 1.425 Å) derived in absence of O but with the interlayer distance (d12) fixed to its experimental value (3.347 Å) b s or s' correspond to C-C bridge sites. For LH reaction coverages are necessarily doubled respect ER or adsorption processes, when using the same unit cell c Energy barriers respect reactants (including zero point energy between parentheses) d The transition state geometries for ER and LH reactions are depicted in Fig. 2. <OCCO dihedral is 24.7° at 25/50% and and 16.8° at 11.1/22.2% for LH TS e Harmonic vibrational frequencies of the atomic adsorbate with respect to the rigid substrate. Approximately, the two first correspond to perpendicular and the other to parallel vibrations respect the graphite surface for the ER reaction. In the case of LH reaction the imaginary frequency shows as mixture of parallel and perpendicular movements f This TS was derived approximately from some grids of points Top view of the graphite (0001) slab two-layer model with the different adsorption sites: bridge (B), on top a C with also a C below (T1) and on top a C without a C below (T2). The atoms of the second layer are plotted a bit larger and with light color to be distinguished from first layer atoms.
The transition state for O diffusion between two bridge adsorption sites is also shown. The solid line shows the surface 1x1 unit cell