Eley-Rideal reaction dynamics between O atoms on b-cristobalite (100) surface: a new interpolated potential energy surface and classical trajectory study

We present a theoretical study of the collisions of atomic oxygen with O-precovered b-cristobalite (100) surface. We have constructed a multidimensional potential energy surface for the O2/bcristobalite (100) system based mainly on a dense grid of density functional theory points by using the interpolation corrugation reducing procedure. Classical trajectories have been computed for quasithermal (100-1500 K) and state-specific (e.g., collision energies between 0.01-4 eV) conditions of reactants for different O incident angles (qv). Atomic sticking and O2(adsorbed) formation are the main processes, although atomic reflection and Eley-Rideal (ER) reaction (i.e., O2 gas) are also significant, depending their reaction probabilities on the O incident angle. ER reaction is enhanced by temperature increase, with an activation energy derived from the atomic recombination coefficient (gO(qv=0°,T)) equal to 0.24±0.02 eV within the 500-1500 K range, in close agreement with experimental data. Calculated gO(qv=0°,T) values compare quite well with available experimental gO(T) although a more accurate calculation is proposed. Chemical energy accommodation coefficient bO(T) is also discussed as a function of ER and other competitive contributions.


Introduction
Surface reactions contribute considerably to spacecraft heating during hypersonic re-entry into planetary atmospheres [1]. In the case of Earth atmosphere, atomic recombination processes involving oxygen and nitrogen in the boundary layer are catalyzed by the materials used as thermal protection systems (TPS) of these vehicles. The apparent simple elementary heterogeneous processes that occur in these conditions are difficult to study either experimentally or theoretically, even at present [2]. Along a typical re-entry trajectory the density of air changes by several orders of magnitude. Thus, there is a free molecular to continuum regime change and a laminar to turbulent flow transition. Moreover, different thermal and chemical equilibrium/non-equilibrium situations are observed at several points of the re-entry trajectory. These complex real conditions cannot be well reproduced even in sophisticated wind tunnel facilities (e.g., Von Karman Institute For Fluid Dynamics, Belgium). On the other hand, real flights (e.g., Shuttle or forthcoming EXPERT European Project [3]) can report only limited information. Standard computer fluid dynamics (CFD) simulations [2] carried out for this kind of flights give values of the aerodynamic heating with a factor two or more larger depending if a non-catalytic or a fully-catalytic model is assumed for the principal surface reactions. In these simulations two important type of coefficients are needed: a) the atomic recombination coefficient gi (also called catalycity) and b) the chemical energy accommodation coefficient bi [4]. gi is defined as the probability of atomic recombination per atomic surface collision for each i species and bi represents the amount of energy released to the TPS surface respect the maximum energy transferable per atomic recombination (0 ≤ gi , bi ≤ 1). These coefficients depend on temperature and total and partial pressures (i.e., gi (T, P)). Hence, precise thermal or state-specific theoretical and/or experimental gi and bi coefficients are necessary to simulate these severe equilibrium / non-equilibrium conditions and to improve the design of better TPSs. TPS materials should have high emissivity (to maximize the amount of energy re-radiated), low surface catalycity (to minimize convective heating) and low thermal conductivity (to minimize the mass of material required to insulate the primary structure).
A common material used in Shuttle TPS coatings is the reaction-cured glass (RCG), a borosilicate composed by 94% of SiO2, 4% of B2O3 and 2% of SiB4. Thus, we decided to study several collision processes involving O and O2 over silica, considering a crystalline structure, which simplifies initially the theoretical approach. Moreover, the selected structure, b-cristobalite, is the most stable polymorph at high temperatures up to the melting point of 1983 K, with very similar properties to amorphous silica (e.g., density, refractive index, band structure, etc.). We have carried out a first-3 principles theoretical study on several elementary processes involving these species. In an earlier study [5] we made a full density functional theory study (DFT) of the atomic oxygen and nitrogen interaction with a b-cristobalite (100) surface, because atomic adsorption is the first basic step for a subsequent atomic recombination reaction through an Eley-Rideal (ER), hot atom (HA) or Langmuir-Hinshelwood (LH) mechanism. We observed a large oxygen chemisorption, mainly on top a Si of the Si-ended face, which had a great influence in the atomic oxygen scattering, causing a predominant adsorption or penetration (absorption) into the SiO2 slab [6]. In two more recent studies [7][8] we presented a dynamics study of O2 scattering on a b-cristobalite (100) surface, which showed that reflection and sticking (absorption/adsorption) were the main processes, becoming the dissociative sticking only open at high collision energies (i.e., Ei > 1 eV). We have also recently reported a preliminary study on the oxygen ER reaction over b-cristobalite (100) surface [8] using the same PES and considering several state-specific initial conditions of the reactants.
There are several experimental studies about the ER oxygen recombination reaction over some silica surfaces (quartz, cristobalite, RCG, vitreosil,..) mainly concerned in the determination of g and effective g (i.e., geff = g' = g·b) coefficients (e.g., see a recent review [9]) at high temperatures (400-1800 K) . These coefficients show marked temperature dependence, which increase with temperature until a critical value where the thermal desorption becomes important and then the recombination coefficient decreases. Therefore, the oxygen ER reaction over silica seems to be an activated process (e.g., Ea = 0.29 eV on b-cristobalite for 800 ≤ Ts ≤ 1830 K or Ea = 0.16 eV on quartz for 850 ≤ Ts ≤ 1430 K [10]. The nature of the silica material has as well a great influence on these coefficients. Only a few experimental studies have reported information on the chemical energy accommodation coefficient for oxygen over silica, indicating a strong variation with the temperature (e.g., 0.12 ≤ b ≤ 0.69 within 1770-966 K for b-cristobalite [11]) contrary to the common assumption made in the major part of experimental g' measurements [9], where a constant value equal to 1 is supposed.
There are scarce first principles theoretical studies on oxygen ER reactions over silica surfaces. A previous work [12] presented a semiclassical trajectory study for atomic oxygen recombination over b-cristobalite although they used a tentative empirical PES, not fully justified nor analysed. More recently, the same authors have published a very similar theoretical study over quartz [13], showing similar results as for b-cristobalite, although with a g about 3 times lower.
In the current work, we present some additional DFT calculations for O2/b-cristobalite system and we explain the construction and the properties of the new interpolated PES for this system, than can be used for both the study of molecular oxygen collisions on clean b-cristobalite (100) and for atomic oxygen collisions with O atoms preadsorbed (i.e., O(ad)) onto the same solid surface. The dynamics 4 of the second kind of process is studied for reactants at state-specific or thermal conditions by using classical trajectories and they are compared with available experimental data.
The present results together with our previous data on these elementary processes (i.e., atomic adsorption, molecular dissociative adsorption,..) involving air species and b-cristobalite are planned to be used in some kinetic models and kinetic Monte Carlo simulations to determine appropriately the g and b coefficients as a function of partial and total pressures and temperatures.

Density functional theory calculations
We have performed extensive DFT calculations for the O2/b-cristobalite (100) systems in order to construct a multidimensional PES for the present system, using as well previous DFT data from the O/b-cristobalite (100) [5,6] studies. The energetic calculations have been carried out by using the VASP code [14][15][16][17] to determine the O2 and the O/O(ad) interactions over an Fdd2 b-cristobalite (100) surface [18][19]. We have used a similar procedure as for the previous O/b-cristobalite (100) studies [5,6]. Thus, a lattice parameter of 7.348 Å, determined for the bulk Fdd2 calculation, was used in the present work along with a slab model with 6 layers, ended with a silicon first layer. An extra hydrogen bottom layer was added to saturate the O dangling bonds of the inner layer. The distance between slabs (z = 17-18 Å) was large enough to prevent significant interactions between them. We carried out generalized gradient correction DFT calculations with the Perdew-Wang 91 functional (PW91). The electron-ion interactions were described by using the projector-augmented-wave technique. An energy cut-off of 400 eV for the plane wave basis set was precise enough to obtain converged properties. Integration over the Brillouin zone was made by means of a 3x3x1 k-points mesh by using the Monkhorst-Pack method. Spin-polarized calculations with a fixed magnetization equal to 0, 2 or 4 were carried out to consider the total possible atomic and molecular spin states (i.e., singlet, triplet or quintuplet). The final interpolated PES (see next section) will use the lowest energy for each geometry (basically singlet at short distances and triplet for longer ones), allowing the interpolation method a smoothing of the conflicting crossing regions. Thus, we assume that the evolution of the system is electronically adiabatic, which requires at least one transition between states of different spin, in the same way as for the O/b-cristobalite PES [6].
In order to make easier the search of stationary points we have considered two approaches: a) direct DFT calculations with geometry optimization and normal mode analysis and b) interpolation of 5 a grid of DFT energies calculated for both parallel (q=0º) and perpendicular (q=90º, f=90º) O2 configurations over the central Si site of the unit cell (see next section); q and j are the polar and azimuthal angles of the molecular axis as shown in Fig. 1a. The O2 distance (r) and the Z coordinate of O2 centre of mass (ZCM) were changed for 1 £ ZCM £ 4 Å and 1 £ r £ 3.2 Å (175 points), and 0.5 £ ZCM £ 4 Å and 1 £ r £ 4 Å (216 points), respectively, with the slab geometry fixed at its most stable equilibrium configuration in absence of O2.
We have calculated a minimum energy path for the Eley-Rideal reaction (i.e., from O + Oad(on top Si) to O2(g)) by means of the Nudged Elastic Band method (NEB) [20][21], where the four first layers of the slab have been relaxed. It is found a transition state (TS) geometry, which implies a significant reconstruction of the first layer, with the following oxygen final geometry: XCM = 5.432 Å, YCM = is opposite to the normally assumed high exothermicity in the major part of kinetic models on air over silica materials, because they consider much lower O adsorption energies (e.g., 3-3.5 eV [22][23]) than the present DFT results. The higher oxygen atomic adsorption energy (as it is also observed over metals), the smaller exothermicity of the ER reaction.
We have also found that O2 molecule can become adsorbed on top a Si atom of the first layer, shows very stable minima, also at lower coverages (e.g., using a 2x2 cell). Therefore, the Langmuir-Hinshelwood mechanism will rarely occur because of the very large endothermicity (∆ELH ≥ 6.0 eV). Moreover, we have also shown that the O adatoms migration, even between close sites (a first step necessary for LH processes), will be very difficult due to their high-energy barriers [5].

Interpolated PES: construction and analysis
We have constructed a multidimensional PES for the O2/b-cristobalite system in a similar way as we carried out for the O/b-cristobalite system [6]. At first, we derived a potential for this system at the slab equilibrium geometry, The first and second terms describe the interactions between each oxygen atom and the atoms of the slab at equilibrium. Both potentials are equivalent and were obtained by using the corrugationreducing procedure (CRP) interpolation method [24] over an extensive grid of DFT points, as was widely explained in a previous work [6]. 7 The third term that contains the O-O interaction, which is modified by the presence of the rigid slab at its equilibrium geometry, is the so-called (within the CRP method) molecular interpolation function [24].
To derive an expression for V OA-OB we have calculated a grid (ZCM, r) of DFT energies (about 700 points) for both perpendicular (q=0º) and parallel (q=90º, f=90º) molecular configurations, with its centre of mass over the T1 site (i.e., XCM=YCM=XSi=YSi, Fig. 1b). The molecular interpolation functions for these two 2D-(ZCM, r) cuts on the T1 site, IT1(q=0º) and IT1(q=90º), were interpolated using cubic splines, and then IT1(ZCM, r,q) was obtained by interpolation over q by using a cos(2q) function, assuming equal weights for both angles: (2) Due to the interaction of O2 with the slab over the hollows (e.g., H1 sites) of the unit cell is almost negligible (like O2 in vacuum), we have assumed that on H1, the molecular interpolation function only depends on the interatomic distance and is equal to the O-O potential in vacuum (i.e., IH1(r)=VO2(r) ).
Then, we have approximated V OA-OB by a 5D potential (fixing f at 90º), which interpolates IT1 and IH1 throughout the surface unit cell using the symmetry-adapted Fourier expansion: To represent the interaction of the O2 molecule with the slab out of its equilibrium configuration, we have used a procedure proposed previously [25]. Therefore, where the pair potentials and (the same used and described in [6]) become zero when (l = i, j), which implies that NSi' < NSi and NO' < NO, being NSi + NO the total number of slab atoms. On the one hand, for the equilibrium configuration of the slab (i.e. {Rij}={Rij}eq), only the first term on the right hand side of eq. (4) [defined in eq. (1)] survives. On the 8 other hand, the variation of the O2-surface interaction with respect to the case of the slab in its equilibrium configuration, is accounted for by the sums of pair potentials and .
Finally, in order to introduce the slab motion during the dynamical simulations, we have added an additional potential (V slab ), which depends of all slab internuclear distances {Rij}. Then, the total potential energy surface V PES describing the interaction of the O2 molecule with the slab (in or out of its equilibrium configuration) is given by: In this work, we use the same empirical potential V slab [26][27] as we did for the O/b-cristobalite system [6], where it was depicted and analysed more extensively. The corresponding O-O and O-Si pair interactions were described by using empirical Born-Mayer-Huggins ionic potentials derived for silica. A possible improvement in V slab silica potential could be done by the introduction of an ab initio (DFT) potential that could predict more accurately several crystalline structures (e.g., cristobalite, quartz, coesite,..) and the liquid silica [28]. The present CRP PES allows the calculation of energies for molecular configurations with ZCM ≥ 0 Å and r ≥ 0.5 Å and for atomic configurations with Z ≥ -1 Å, due to the range of DFT data calculated. Hence, diffusion of O2 into the slab cannot be studied without extra modifications on the 9 current PES. On the other hand this PES has also been used in the dynamics study of O2 (i.e., dissociation and reflection) on b-cristobalite (100) surface [7,8].

Classical trajectory method
We have carried out a dynamics study of atomic gas oxygen colliding with O-precovered bcristobalite (100) surface by means of the classical trajectory method [29]. The silica slab includes a total of 104 atoms distributed into 9-layers for a 2x2 square surface unit cell with a cell parameter of 7.348 Å. The slab temperature (TS = 100 -1500 K) was controlled by means of a Generalized Langevin equation (GLE) approach [29,30]. The corresponding random and friction forces were applied only to the atoms of the bottom Si layer, which are also forced to oscillate around their bulk equilibrium positions (by using extra "springs") to avoid a global translation of the slab. We checked that the empirical slab vibrated with an average slab energy corresponding to the selected temperature (<Eslab> = 3NkbTS), hence the energy will not be constant. Generalized friction and fluctuating forces balance, according the fluctuation-dissipation theorem, so that the proper temperature is reasonably maintained in the primary zone (i.e., the lattice atoms that can interact directly with colliding O atom).
The colliding oxygen atom is located initially at a distance of 4.5-6.0 Å from the surface and its aiming point (XO,YO) is determined by a random uniform sampling only on an 1x1 unit cell.
Preadsorbed O atom was located over the central Si atom (T1 site) at an initial Z of around 1.5 Å.
This atom was quickly thermalized at the selected slab temperature before the arrival of the incoming O atom. We checked as well that only at high temperatures (e.g., 1000 K) this adatom in T1 could evolve to T1' (penetration into the slab) for collision times longer than about 2 ps. Therefore, this initial geometry for O(ad) seems to be appropriate for the study of all processes derived from O(g) + O(ad) collisions within the time scale considered.
Kinetic energy (Ei) of the colliding atom was within the range 0.05 -4.0 eV. The corresponding velocity angles qv and fv (Fig. 1a) were fixed (0º and 45º) and sampled by an uniform random sampling in the interval (0-360º), respectively.
Hamilton's equations were integrated by using a modified Beeman algorithm with a fixed step of 1x10 -4 ps, which permits an energy conservation of 10 -4 -10 -5 eV in absence of the GLE bath. Two kind of calculations were considered: quasithermal (T = TO(g) = TS, qv) and state-specific (Ei, TS, qv).
Batches of 3000-25000 trajectories were integrated for the different initial conditions, which ensured statistical errors below 1% for the reaction probabilities of the main channels. Trajectories were integrated up to a maximum time of 2 ps. Some tests at longer times (e.g., 5 ps) were also made to check the convergence of reaction probabilities. Nevertheless, some trajectories were ended up at shorter collision times due to the before mentioned PES limitations, although we think that the final probabilities are even so quite well converged.
Several main processes can occur: atomic reflection, atomic sticking, and O2 formation (i.e., Eley-Rideal reaction). We have also differentiated the molecules that become retained by the slab (O2(ad)) at least during the integration time. In spite of these adsorbed molecules could finally desorb or produce two adsorbed or gas atoms at a much longer time scale, the PES topology seems to prevent a possible significant increase of the ER probability for longer collision times. In the present study we do not separate between atomic (or molecular) adsorption and absorption (penetration) processes as in our earlier paper [6], but this is practically irrelevant for the main goal of the present work, the ER reaction dynamics study.
We have classified the reaction channels analysing r, ZO, ZO' and ZCM values and the corresponding z velocity components of both atoms [6] and molecule. For instance, the ER channel was achieved when ZCM was much larger the 5-6 Å, the r distance was not far from the O2 equilibrium distance (i.e., r < 1.9 Å) and vz CM > 0. We verified that once the integration was finished (usually in 2 ps) and the trajectory classified, this was equally classified after an extra number of steps. In case of discrepancy the trajectory was integrated more steps until a converged classification. This criterion was especially useful for channels involving finally adsorbed species.
The final internal state (v', j') for O2(g) was assigned rounding up or down to the nearest integer its rotational angular momentum j' (quasiclassical approach) and by means of the radial action variable semi-classical method for the vibrational quantum number (v') [31]. In spite of O2 molecule presents only odd rotational levels at its fundamental electronic state, this was not taken into account for final j' truncation due to the high closeness of the rotational levels and because it was necessary to use bins of at least 5j' consecutive values to reduce the statistical errors in the final distributions. formation and ER reaction for normal incidence and a surface temperature of 1100 K. Atomic sticking is the principal process and it increases with translational energy increment, in agreement with the O scattering over a clean b-cristobalite surface [6]. O2(ad) formation is the second channel in importance and decreases with the augment of translational energy. Atomic reflection is lower than 10 % and ER reaction is almost negligible in comparison with the other channels. Atomic desorption channels (i.e., involving O(ad)) are practically forbidden in these conditions.

Classical trajectory results
The comparison with Ts = 300 K results [8] seemingly shows small differences, mainly at low collision energies (approximately for Ei ≤ 0.5 eV) for O sticking and O2(ad) formation probabilities.
Nevertheless, there is a relevant difference for the ER process, which is completely absent at 300 K but its contribution is significant at 1100 K, becoming more important for the lowest collision energies (see inset in Fig.3a).
Off-normal incidence enhances strongly the atomic sticking (Fig. 3b)  The atomic reflection for off-normal incidence is now larger for all collision energies when compared with normal incidence. Fig. 4 shows more clearly the effect of the incident angle on the different probabilities. This reveals that the comparison between experimental and theoretical results should take into account an averaging over this angle (qv), specially for O sticking and O2(ad) formation; for ER reaction or O reflection processes this effect is smaller although it is observed that O collisions with small qv angles enhance the ER process (e.g., probability changes from 1.33x10 -3 at 0° to 0 at 60° for 1 eV and 1100 K) and decrease the reflection more significantly as shows Fig. 4. The formation of O2 molecules, which remain strongly adsorbed, is much wider than the ER reaction (i.e., O2(g) formation) in good agreement with an earlier study over this system, but based on an empirical PES [12]. A semilogarithmic plot of the quasithermal ER probability (gO(qv=0°,T)) is shown in Fig. 6 to appreciate much better its behaviour compared with some experimental gO(T) values for different silica materials. In fact, we have assumed a kind of effusive atomic beam (i.e., thermal O distribution) colliding perpendicularly over the solid surface, which possibly does not match exactly the non-fully characterized experimental conditions [10], but which is supposed as well in some earlier similar studies [12,13,34]. The averaging over the incident angle (qv) would decrease the presented g(T) values (possibly in a factor that could be up to 10) , as could be deduced comparing fig. 5a and 5b.
Due to the LH process practically will not be open at these temperatures because it is highly endothermic and that the molecular dissociative adsorption is also negligible in these conditions, as we demonstrated in some previous O2/b-cristobalite dynamical studies [7][8], we have made the assumption that gO » gER for the classical trajectory values (calculated as gER = 2PER) shown in Fig.   6. This approximation can be supported by some previous theoretical and experimental studies for O recombination over b-cristobalite [12,34] or quartz [13], that conclude that gER > gLH and consequently gO » gER, particularly at high temperatures.
However, a more accurate gO determination is planned to be done in our research group by means of a kinetic Monte Carlo technique [37], which can elucidate the pressure dependence of gO (i.e., surface coverage dependence), the importance of O2(ad) formation (possibly increasing g) and also its values for non-steady conditions.  Fig. 6 correspond to normal incidence, the good agreement between classical trajectory and experimental g values and also for the activation energies should be taken with caution, as qv averaging would decrease these coefficients.
Moreover, an accurate direct comparison cannot be done without using an expression like (6).
First, gO for air and pure oxygen are similar but not exactly equal as it has been stated in a recent review of catalycity values for silicon based materials [38]. Second, the type of silica structure (different crystalline or amorphous phases) shows as well a marked effect on these coefficients [38], as it can be easily observed in Fig. 6. In fact, the experimental values from ref. [10] were obtained forming in situ SiO2 b-cristobalite by oxidizing sintered SiC. Therefore the sample was not a perfect crystalline structure. Surface roughness tends to increase gO respect its value (called intrinsic or absolute g) at a microscopic surface (e.g., a 2.4 roughness factor was observed for O or N over silica [23]). Third, catalycity usually increases with surface pressure [38] at moderate temperatures (experimental values for b-cristobalite and quartz of Fig. 6 are measured at a total air pressure of 200 Pa [10]). Fourth, some trajectories classified as O2(ad) formation should be finished at collision times shorter than 2 ps (e.g., a 36 % at 1000 K and 0º) due to the interpolated PES is not defined for molecules penetrating into the slab (ZCM < 0 ); thus, the integration at longer collision times could increase the number of O2 molecules finally desorbed (specially at high temperatures), increasing accordingly the ER probability, although this PES limitation will be less important for another reaction channels. Therefore, the inclusion of all these four aspects into the gO calculation would increase a bit more its values (although qv averaging would tend to decrease), expecting a possibly even better agreement with the experimental values. Nevertheless, additional measurements (e.g., with arc-jet facilities or inductively coupled plasma generators) should be also necessary to confirm the validity of experimental data [10], due to the usual wide dispersion observed in catalycity coefficients over several thermal protection materials [38].
Present results show some significant differences respect reaction probabilities and gO(T) or gO(Ts, Ei) coefficients when compared with some previous classical trajectory calculations based on an empirical PES for the O2/ b-cristobalite system [12,34]. Thus, these differences are originated not only by the topology of the PES but also by the selection method used for the initial adsorbed atom positions (O(ad)). Three cases were considered in the first study [12]: O adsorbed on top a Si, over a bridge between two Si atoms or in a random position over the unit cell, for the Si ended surface. This mentioned study showed important differences depending on the type of choice. Thus, for instance, the values of PER(Ts =1000 K , Ei = 0.1 eV) for normal incidence were 0.225, 0.030 and 0.025, respectively, compared with our classical trajectory value of (0.44±0.09)x10 -2 for O(ad) on top selection. However, the authors finally chose the O(ad) random sampling over the cell instead of locating O(ad) on top Si to calculate the catalycity coefficient at 1000 K (gO = 2.9 x10 -2 ), with possibly a fortuitous agreement with the experimental value (2.7 x10 -2 [10]) because their coefficients calculated for O(ad) on top Si would be much higher, with a much worse agreement with the experiment. The same authors in a recent study of N atom recombination over a silica surface [13] place more appropriately the N(ad) on top a Si atom of the silica surface. Their PES for O2/bcristobalite neither does support any temperature effect on the ER probability. Such effect was only found for the LH probability.
A noteworthy conclusion together from the present and previous classical trajectory studies for O reaction with O-precovered b-cristobalite is that O2(ad) formation will be important at low collision energies, competing with the ER reaction.
We have also analysed the energy distributions of the O2 molecules produced in ER reaction.
Nevertheless, its low reaction probabilities imply that the statistical uncertainties will be larger in comparison with another reaction channels (e.g., O2(ad) formation). Thus, Fig. 7 shows the translational, vibrational and rotational distributions for normal incidence and reactants at 1000 K.
Molecules are translationally hot, with an average Ei' equal to 0.84 eV, much larger than if products remained at the same temperature as the reactants (at 1000 K the average energy would be 0.13 eV).
The new O2 molecules are also internally excited (see Fig. 7b  These results are quite different from the ones obtained with the before mentioned empirical PES [12], which showed a very small translational excitation for low O collision energies, Ts = 1000 K and O(ad) located randomly over the unit cell or on a Si bridge (results for O(ad) on top Si, like in the present study, were not reported to compare with).
Due the interest in the determination of the b coefficient for pure oxygen/b-cristobalite system, we have calculated the fraction of energy transferred with the slab (fslab) caused by the ER reaction for 1000 K and normal incidence, by means of the equation, The values obtained are negative, with an average value of <fslab> = -0.31 ± 0.09, showing that in these conditions the slab becomes colder during the ER reaction as a result of the energy transfer from the slab to the new formed O2 molecule. Nevertheless, the absence of energy conservation due to the thermal bath applied by means of the GLE approach could generate some significant uncertainties in this fraction. The energy transferred with the solid can be related with the b coefficient, which is usually defined as [4,38], However, some authors have made the assumption that b » bER, but defining this coefficient in a different way to (7); for instance, for the oxygen/b-cristobalite system it was calculated by using, where Ej indicates the total energy for each j species. A value of 0.127 for reactants at 1000 K was calculated in this way [34]. Nevertheless, this coefficient should not be straightforwardly compared with the experimental b coefficient because the clear different definition. Up to our knowledge, there exists only a work reporting some direct experimental measurements of bO for air and b-cristobalite [11], whose values are in the range 0.69-0.12 for 1000-1800 K at a total pressure of 200 Pa; similar values for quartz were also measured [39]. This behaviour was also observed for collisions with the clean cell [6] and was justified by the strong oxygen adsorption and penetration (i.e., absorption), which is much more hindered along the other diagonal (Fig. 1b). The formation of both gas and adsorbed molecular oxygen is mainly produced for O collisions on a region of the cell very close to the initial adsorbed atomic oxygen, with a clear Eley-Rideal instead of a hot-atom mechanism. The O reflection processes are mainly produced over the corners of the cell, far from the central O(ad). Thus, it seems that the zones of the cell tend to present a different selectivity to produce all possible elementary gassurface processes. Moreover, it is observed that the slab is not a simple spectator for these processes.
Its geometry is changed along the trajectories not only by the thermal effects but also by the chemical processes. For instance, the position of the Si atom below the O(ad) is shifted from its initial position in all three directions, with a different extension depending on the kind of process. Nevertheless, an accurate analysis of these slab changes would require the integration again of a large number of trajectories at shorter collision times (e.g., t < 2 ps) to be able to distinguish suitably between thermal and reaction channel induced geometry distortions.

Conclusions
The main aim of this work has been to study the interaction of atomic oxygen with the O- Atomic sticking is the principal process, which increases with translational energy increment.
O2(ad) formation is the second channel in importance and decreases with the augment of translational energy. Atomic reflection and ER reaction are much less important although also significant. The new O2 molecules produced by the ER reaction at high temperatures become internally and translationally excited, while the slab becomes colder. These results along with the positive experimental bO(T) coefficient could indicate that another processes (as for g calculation) should be also considered to be able to balance this effect in order to produce a total b positive value.
Several regions of the cell tend to present a different selectivity for the possible elementary gassurface processes. Moreover, it is observed that the slab is not a simple spectator for these processes.
It seems that its geometry is changed along the trajectories not only by the thermal movements but also by the type of the chemical processes. Thus, a rigid slab approximation for the study of these processes possibly would not be too appropriate.
The present PES model seems to be more realistic than the previous empirical one and allows the kinetics and dynamics study of O or O2 collisions over b-cristobalite in quite good agreement with the limited experimental data. Nevertheless, some improvements (e.g., additional O-b-cristobalite DFT data, an accurate DFT V slab potential,…) should be necessary to permit as well the study of some subsurface situations or processes, than could have some influence on the studied reaction probabilities. Tables   Table 1 Relative       Calculated gO(qv=0°,T) (solid circles) and experimental gO(T) atomic oxygen recombination coefficients versus the reciprocal temperature for several silica materials. The solid lines show Arrhenius fits of some of these data. Experimental data correspond to: air with b-cristobalite (solid squares) [10], air with a-quartz (solid triangles) [10] and oxygen with fused quartz "vitreosil" (solid diamonds) [33]. Experimental bar errors are ± 30% [10] and classical trajectory values represent one standard deviation.