Dynamics of the Oxygen Molecules Scattered from Graphite (0001) Surface and Comparison with Experimental Data

A quasiclassical trajectory dynamics study of molecular oxygen colliding over a free of defects and clean graphite (0001) surface has been performed with a recently published density functional theory based flexible periodic London-Eyring-Polanyi-Sato potential energy surface (PES). Although the PES was mainly constructed for describing accurately the recombination of atomic oxygen over an O-preadsorbed surface, here we show that this PES is also reliable to study the scattering of O 2 over graphite surface. Thus, several initial conditions have been explored: collision energies (0.2 £ E col £ 1.2 eV), incident angles ( q v = 0 ° , 45 ° ) , surface temperatures (100 £ T surf £ 900 K) and some rovibrational O 2 levels (v = 0,1,2 and j = 1,17,25). The calculated polar scattering angular distributions are in good agreement with the experimental ones in a wide range of explored conditions . Moreover, the comparison with hyperthermal experimental data, which was unclear in a previous work, has been finally clarified. The effect of O 2 (v,j) internal state on the scattering is very small.


I. Introduction
Satellites and other orbiting vehicles in low Earth orbits (i.e., at 160 -2000 Km above Earth) must travel at velocities around 7.5 km/s in order to avoid collapsing on Earth's surface due to gravity. At these high velocities the shielding of the vehicle receives a great amount of highly energetic collisions (i.e., about 4.5 ± 1.0 eV 1 ) of the main present species (O/O2). These collisions can produce several heterogeneous processes (etching, adsorption, reflection, recombination, oxidation) and also a flux of heat to the surface of the vehicle. It is worth for space industry to model this kind of elementary processes in order to improve the behaviour of the thermal protection systems of space vehicles at these extreme conditions. Among these processes, the oxidation reaction between O2 and several carbonmade materials such as graphite, graphene and single-wall carbon nanotubes is known to play an important role not only during these reentries but also in many other fields (e.g., heterogeneous oxidation catalysis, solubility and functionalization of single-wall carbon nanotubes). The oxidation of these materials alters their physicochemical properties such as the wettability, adsorption (sticking) of the surface and use to be a drawback for their technological application. It is known that the molecular chemisorption and the dissociative adsorption of O2 molecule over carbon materials depend on the availability of reactive sites (i.e., graphitic edge planes or defects 2 ) and also on the surface temperature. It is important to note that graphite with defect-free basal planes undergoes thermal oxidation at extremely high temperatures whereas at low-temperatures exposure leads to only weakly physisorbed O2. According to this, a pristine (i.e., free of defects) graphite surface would be the more appropriate and desirable for acting as a thermal protection system since the number of surface reactions would be minimized and also the heat transfer to the surface.
Recently, some experiments 3,4 reported the polar scattering angular distributions resulting from a highly collimated and nearly monoenergetic incident beam of O2 scattered from the surface of highly 3 oriented pyrolitic graphite. An important characteristic of the experimental setup is the use of a fixed source-detector geometry in which the angle between the incident beam and the detector direction was held at an angle of 90º. These experiments were carried out at the particular incident energy within the range 0.291 eV £ Ecol £ 0.614 eV and also for a surface temperature in the interval (150 K £ Tsurf £ 500 K). For each collision energy and surface temperature the angular distribution intensities were reported as a function of final polar scattering angle (θv') relative to the surface normal. Moreover, some experiments 5,6,7 were also performed with hyperthermal beams of O( 3 P) atoms with some extent of O2(X 3 Σg -) that impinge over highly oriented graphite. In a previous work 10 , we showed that the quasiclassical trajectory (QCT) polar scattering angular distributions, I(qv'), of the O2 molecules formed by the Eley-Rideal reaction presented a peak centred at qv' = ~40º much lower than the experimental O2 distribution, peaked at qv' = ~60º 5 . This discrepancy is understood since in the experiment, a mixture of O2/O forms the initial beam, so the higher contribution to the O2 polar scattering angle distribution should come from the O2 molecules directly scattered by the surface. Thus, in order to elucidate this point we have completed our previous Eley-Rideal QCT study 10 with the addition of a new O2 QCT scattering analysis for the same initial experimental conditions. This paper presents a QCT dynamics study of O2 molecules colliding with a pristine graphite (0001) surface. Our recent Flexible Periodic London-Eyring-Polanyi-Sato (FPLEPS) potential energy surface (PES) based on density functional theory (DFT) data has been used. This PES was used previously for describing the O collisions with a clean 8,9 and with an O-preadsorbed 10 graphite (0001) surface, obtaining a good agreement with experimental data. As above mentioned, thermal oxidation on free of defects graphite surfaces only occur at extremely high temperatures, so the FPLEPS surface (for a graphite surface without any defects) must be accurate enough for describing the main process (i.e., molecular scattering) occurring during the collision. The use of an analytical expression for the PES allows obtaining easily and efficiently the potential energy and its first order derivatives for performing fast dynamics calculations. Some extra DFT calculations were also performed in order to know the accuracy of the analytical PES in describing the O2 interaction with the graphite surface. This paper is organized as follows: Section 2 provides a brief description about DFT and QCT calculations; Section 3 describes the FPLEPS vs. DFT comparison of the O2 scattering (i.e., reflection) over graphite and Section 4 presents the main dynamical results and discussion, including the comparison with experimental data. Finally, Section 5 gives the summary and conclusions.

B. DFT calculations
Periodic density functional theory calculations have been performed using the Vienna ab initio simulation package (VASP) 11,12,13,14 . All the details about the methodology used (e.g., GGA/RPBE functional, energy cut off, slab model, …) were discussed deeply in previous works 8,9,10 devoted to atomic oxygen interaction with both a clean and an O-preadsorbed graphite surfaces, and are omitted here. New DFT calculations have been performed in order to check the accuracy of the analytical FPLEPS PES in describing the O2/graphite interaction. According to the previous study 9 , it is important to remember that the molecular channel does not present any particular feature; it is essentially repulsive close to the surface and a high endothermicity (5.22 eV) prevents the molecular dissociation at the conditions studied here. Thus, the DFT data used to construct the analytical surface were focused mainly on the Eley-Rideal entrance channel. To check the proper description of the interaction of the incoming molecule with the surface, some additional configurations at several particular f angles ( Figure 1) were calculated at DFT level and compared with the data obtained from the FPLEPS PES, with a good agreement as described below in Section III.

B. QCT calculations
A dynamical study involving the interaction of molecular oxygen with graphite surface was performed by means of the QCT method 15 . Several initial conditions were sampled in order to investigate state-specific and hyperthermal processes. State-specific calculations (v, j, E col , Tsurf, θv) were carried out for fixed initial collision energies in the range of 0.2 eV £ E col £ 1.2 eV for a given surface temperature Tsurf within the 100 -900 K interval and several vibrorotational states (v = 0,1,2 and j = 1,17,25). The incoming velocity angle (θv in Figure 1) is defined with respect to the negative Z axis (i.e., 0° for normal incidence) and its projection onto the X-Y plane (i.e., the azimuthal f v angle) was uniformly sampled within the 0° -360° interval since highly oriented graphite surface consists of many orientations within the topmost layer, being the measurements an average of all azimuthal directions 16 .
The initial collision angles selected for the present study were fixed values at θv = 0° and 45°, and also a random sampling within the 10° £ θv £ 70° range. Initial molecular orientation angles (θ, f in Figure 1) were sampled by the standard Monte Carlo method within the intervals 0° £ θ £ 180° and 0° £ f £ 360°, respectively. The initial r internuclear distance was chosen between the corresponding inner (r-) and outer (r+) turning points of each O2 (v, j) internal state. Only odd j initial levels were chosen for O2 (X 3 Σg) molecule because even rotational levels do not exist due to the nuclear spin of O 16 is zero 17 .
Initial position (Xcm, Ycm) of the molecular centre of mass was randomly selected inside the (1 ´ 1) unit cell while Zcm was set at 7.0 Å, where the interaction with surface is negligible.
A Generalized Langevin Oscillator (GLO) model 18,19,20 was used in order to account for the energy exchange between the molecule and the solid surface, which was considered as a rigid slab (i.e., 6D PES). This model includes the surface and ghost particle motions into the Hamilton equations and uses also the same optimum parameters, which were previously derived and tested for atomic oxygen impinging on both clean 9 and O-preadsorbed graphite surfaces 10 . The oscillator surface frequencies used were wi,x = wi,y =10 -3 and wi,z =3.4 ´10 -4 au for i = 1,2,3, with an effective mass of 60 amu and with friction constants gg,x = gg,y = gg,z = 4.0 ´10 -4 au. The need for an effective mass much larger than that of a single carbon atom has been documented previously, and points out a collective scattering from a large number of surface carbon atoms 3,4 .
The qctsurf code developed in our group was used to calculate the trajectories integrating the Hamilton equations of the system (included the GLO) using the Beeman algorithm. The time step used was 5´10 -17 s, which ensures a total energy conservation along the trajectories lower than 1´10 -4 eV in absence of the thermal bath. The total energy is not constant due to the thermal bath dissipation effect.
Total integration time was set to a maximum value of 2.5 ps, that allows a correct classification of the trajectories as reflected, since no other channels were observed. Anyway, the adsorption classification was set for Zcm values lower than 2.2 Å and for a number of total molecular rebounds with the surface larger than 8, as was also set previously 9,10 for atomic oxygen impinging a clean or an O-preadsorbed graphite surface. On the other hand, O2 was considered reflected if Zcm was higher than 7.1 Å and the direction of the velocity vector pointed to the vacuum.
The number of total trajectories (NT) calculated for each state-specific condition was around 20000.
However, the number of trajectories carried out for comparing with experimental results was increased for some conditions. Thus, for simulating the hyperthermal experimental conditions of Paci et al. 5

III. Analytical FPLEPS potential energy surface
The complete expression used for describing the 6D potential V(RA,RB) of the two atoms (A, B) interacting with a rigid graphite surface ( Figure 1) was defined as 9,10 : show that no molecular adsorption is present and that only molecular dissociative adsorption onto no consecutive bridge sites is observed. Thus, for f = 0°, the distance between the two bridge sites where the atomic adsorption occurs is 4.936 Å while for f = 60° the distance is shorter (2.468 Å).
Nevertheless, the process is highly endothermic (5.22 eV). The comparison with the FPLEPS contour plots is also shown in Figure 2, where the agreement is quite good without the necessity of including additional DFT points in the fitting procedure. The root mean square deviation (RMSD) for energy values lower than 1.0 eV (240 points) is 0.3 eV. This value is not too small but the reason is that O2 dissociation energy value used in the FPLEPS PES is the one obtained from more accurate coupledcluster calculations 21 instead of the less accurate DFT ones, hence giving to the FPLEPS PES a good description for both asymptotes. The RMSD value without considering the O2 entrance valley diminishes until 0.12 eV.

A. State-specific initial conditions
When O2 impinges the graphite surface three principal channels are available: molecular adsorption,  (v = 0, j = 1) with normal incidence and at two surface temperatures. In general, collision energy of impinging molecules is lost during the interaction process (Figures 4a and 4b) although at higher temperatures some intensity appears at positive values (i.e., some energy release from the surface to the molecule, Figure 4b). A fraction of the initial collision energy is going to molecular rotational excitation (no vibrational), as can be seen from the positive values of the internal energy exchanged (Figures 4c and 4d). This increment is slightly higher when the surface temperature is augmented. The major part of the initial collision energy lost is transferred to the surface (Figures 4e and 4f). Particular attention has been devoted in simulating the experimental conditions of Oh et al. 3 For each pair of initial collision energy and surface temperature, the polar scattering angular distribution is reported, taken into account that the incident (θv) and the final (θv') velocity angles, both relative here to the surface normal, must accomplish the relationship θv + θv' = 90°. This particular source-detector geometry setup makes that a great amount of QCT trajectories (i.e., 400000) must be integrated, since only those that verify the previous relationship (~1%) match the experimental conditions. Nevertheless, if no constraint is imposed between initial and final angles, for a particular incoming polar velocity angle (e.g., θv = 45°) the final scattering distribution achieved is very similar to those shown in Figures   3d and 3f, depending on the surface temperature considered. Calculations where run at each specific experimental surface temperature (150 K £ Tsurf £ 500 K) and initial collision energy (0.291 eV £ Ecol £ 0.614 eV), choosing randomly the initial polar velocity angle within the range 10° £ θv £ 70°. The angular tolerance in the final scattering angle (θv') was set to 1°. In the experimental conditions the vibrational temperature should be almost the same as the nozzle temperature (i.e., 300 K for 0.291 eV and 700 K for 0.614 eV 16 ) and the rotational temperature is estimated to be very low after a free-jet expansion (i.e., 10 K for 0.291 eV and 55 K for 0.614 eV 16 ). Therefore, we assumed the lowest vibrorotational O2 state (i.e., v = 0 , j = 1) in QCT simulations. From these panels it is important to note that the shape of the QCT distributions are essentially coincident with experimental distributions although the QCT peaks are slightly shifted to lower angles mainly for Ecol £ 0.5 eV. From Table I (Figures 4a and 4b).
In addition, they observed that the collision energy loss becomes smaller at higher surface temperatures. The amount of Ecol loss has been also calculated from QCT results at the peaks of the polar scattering angular distributions ( Figure 5). Results are listed in Table I together

B. Hyperthermal initial conditions
Some simulations were also run at hyperthermal conditions for Tsurf = 503 K and θv = 45° in order to compare with the experimental polar scattering angular distributions of O2 over highly oriented graphite 5 . A large number of trajectories (60000) were calculated simulating exactly the experimental hyperthermal O2 velocity distribution, whose average translational energy corresponds to 10.46 eV, assuming that at these conditions the O2 molecule is at v = 0 and j = 1 vibrorotional state.
The calculated distribution is shown in Figure 6. It presents a maximum centred at θv' = 73°, a bit larger than the experimental one 5 Figure 6, which is peaked at much lower angles (θv' = 37°). However, in spite of the Eley-Rideal probability is high (i.e., 0.489) when an oxygen atom is preadsorbed over the unit cell at these hyperthermal conditions 9 , the necessary previous O adsorption step over a clean surface has a low probability (i.e., 0.0016). Therefore, the intensity of the Eley-Rideal contribution to the total polar scattering angular distribution should be in principle very small compared with that originated only by the scattered O2 molecules from the initial beam. This simple analysis assumes a previous clean and pristine graphite surface. However, in the experiments the steady-state oxygen coverage seems to be very high. Highly oriented graphite surface exposed to a flux of hyperthermal oxygen atoms becomes functionalized with epoxides to a concentration between 1 or 2 oxygen atoms per 6 surface C atoms, as some direct dynamics calculations also support 5 . The high final flux intensity of O2/O ratio observed in the experiments respect to the initial low value of 0.48 is interpreted as the Eley-Rideal mechanism should contribute highly to the final O2 measured distribution. In such case, the addition of this Eley-Rideal contribution to the final QCT total polar scattering angular distribution would produce a curve peaked at lower angles and even much closer to the experimental one.

V. Summary and conclusions
The study of the collisions of molecular oxygen over graphite surface has been carried out over our The short average collision times along with the analysis of several plots support a single collision dynamics mechanism as was concluded in earlier theoretical studies.
The study at hyperthermal conditions shows a good agreement with the experimental polar 16 scattering angular distributions, which would be even better with the inclusion of the calculated O 2 Eley-Rideal contribution to the final total distribution.        25 Figure 6