Adsorption of atomic oxygen and nitrogen at beta-cristobalite (100): a density functional theory study.

The adsorption of atomic oxygen and nitrogen on the beta-cristobalite (100) surface is investigated from first principles density functional calculations within the generalized gradient approximation. A periodic SiO2 slab model (6 layers relaxing 4 or 6) ended with a layer of Si or O atoms is employed throughout the study. Several adsorption minima and diffusion transition states have been characterized for the two lowest spin states of both systems. A strong chemisorption is found for either O or N in several sites with both slab endings (e.g., it is found an average adsorption energy of 5.89 eV for O (singlet state) and 4.12 eV for N (doublet state) over the Si face). The approach of O or N on top O gives place to the O2 and NO abstraction reactions without energy barriers. Atomic sticking coefficients and desorption rate constants have been estimated (300-1900 K) by using the standard transition state theory. The high adsorption energies found for O and N over silica point out that the atomic recombination processes (i.e., Eley-Rideal and Langmuir-Hinshelwood mechanisms) will play a more important role in the atomic detachment processes than the thermal desorption processes. Furthermore, the different behavior observed for the O and N thermal desorption processes suggests that the published kinetic models for atomic O and N recombination reactions on SiO2 surfaces, based on low adsorption energies (e.g., 3.5 eV for both O and N), should probably be revised.


I. Introduction
Modeling of the different heterogeneous processes involved in nonequilibrium vibrational kinetics of air hitting over silicon dioxide surfaces is necessary for the study and the design of new thermal protection systems (TPS) to be used in atmospheric terrestrial re-entry vehicles 1 . TPS coatings for space vehicles must have very low catalytic efficiency to reduce the surface heat flux due to atomic recombination effects.
Forthcoming development of reusable launch vehicles (RLV) is a prerequisite for guaranteeing long-term access to space. New materials that allow the reduction of launch mass, improved modeling of aerothermodynamic phenomena, increased performance and reliability of propulsion systems and innovative TPS are required to fulfill RLV objectives. The mastering of the detailed chemical physics phenomena associated with TPS catalysis is therefore a key element in this context. Therefore, theoretical, numerical and experimental developments are necessary in order to model more accurately the corresponding gas-surface interactions.
Currently, silica-based materials (e.g., reaction-cured glass (RCG) with 94 % of SiO 2 , 4% of atomic adsorption-desorption, molecular adsorption-desorption, atom-molecule ER or LH reactions,.. Among the different elementary steps involved in the heterogeneous catalysis: 1) diffusion of reactants to the surface, 2) adsorption of reactants on surface, 3) chemical reactions on surface, 4) desorption of products from surface, and 5) diffusion of products away from the surface, steps 1 and 5 are usually fast, and the rate determining step has to be found among steps 2,3 and 4.
The main goal of our research in progress is the theoretical study from first principles of the steps 2, 3, and 4 with several models of a silica surface in order to provide kinetic and dynamical data useful for Computational Fluid Dynamics (CFD) simulations. Published kinetic models for air recombination over SiO2-based materials incorporate rough estimates of several important physiochemical quantities (sometimes used as adjustable parameters), such as, for example, atomic adsorption energies, atomic sticking coefficients, the number of adsorption sites by unit area and others. We believe that a deep theoretical understanding on these processes will permit the proposal of new less catalytic materials to be used as TPSs.
In a first step of the work we have to study the O and N adsorption processes. Nevertheless, scarce experimental data are available on both processes, which are usually treated in a similar way 3 by means of an average adsorption energy of 3.5 eV 1 regardless of the kind of silica surface (e.g., pyrex, quartz, RCG,..), although several used values are within the range of 2.6-5.5 eV 2-6 . These data are indirectly derived from kinetic modeling of atomic recombination on SiO 2 surfaces to fit some experimental data. Thus, for instance, initial surface sticking coefficients S o of 0.05 e -0.002T (on SiO 2 ) or 1.0 e -0.002T (on RCG) and an adsorption energy of 3.5 eV were reported for oxygen in both surfaces 1 ; for nitrogen on RCG the reported data were S o = 0.95 e -0.002T and an adsorption energy of 3.4 eV 1 .
Among the different forms of silicon dioxide to be studied, either amorphous or crystalline (e.g., quartz, tridymite, cristobalite,..), the b-cristobalite is the most stable polymorph at high temperatures up to the melting point of 1983 K; these elevated temperatures are achieved during the Earth's atmospheric re-entry phase. Moreover, b-cristobalite is the crystalline phase of silica with properties closest to those of amorphous silica (e.g., density, refractive index, band structure,..). Recent experimental measurements 7-8 of the oxygen recombination coefficient gO over silica have shown that this coefficient is four times higher for b-cristobalite than for quartz.
Another related theoretical work 9 has also studied the kinetics and dynamics of the atomic oxygen recombination (ER and LH mechanisms) over b-cristobalite. They used a semiempirical interaction potential for oxygen on b-cristobalite along with a semiclassical collisional model, which allow ascertain the importance of both mechanisms depending on the surface temperature. This potential energy surface, as the authors recognize, seems a crude approximation to the real interaction forces, although theoretical and experimental g values are in reasonable agreement.
In this paper we present a wide theoretical study based on periodic density functional theory (DFT) calculations for the O and N adsorption at b-cristobalite (100), with both O or Si as a first layer. We have considered dehydroxyled b-cristobalite as a first approximation to the real amorphous TPS materials used at high temperatures to simplify the theoretical study and to facilitate its comparison with previous and related theoretical 9 and experimental studies 7-8 , which used a-quartz or b-cristobalite. For instance, the experimental results for gO over RCG are quite similar as for a-quartz 8 . Up to our knowledge, this is the first attempt to quantify from first principles these elementary processes, and it is a first stage in the construction of the ab initio potential energy surfaces needed to study the dynamics of the O and N recombination processes

COMPUTATIONAL METHOD
We have performed DFT calculations by means of the VASP code [10][11][12][13] , which uses plane wave basis set. The calculations are based on the generalized gradient correction (GGA) functional Perdew-Wang 91 (PW91) [14][15] . The electron-ion interactions were described by using the projectoraugmented-wave (PAW) technique [16][17] , particularly good for transition metals and oxides, and giving similar results as the FLAPW all-electron method. We checked widely the appropriate energy cut-off in several bulk and slab calculations. Thus, our previous bulk calculations to obtain several properties (e.g., bulk modulus, cell lattice parameters, cohesive energy,..) of some silica polymorphs (e.g., a-quartz, b-quartz, b-trydimite or b-cristobalite) showed that an energy cut-off of 400 eV was accurate enough to obtain converged properties. We have also confirmed this value with some extra slab calculations.
Spin-polarized calculations were carried out to check the ground electronic state of the total system, as both atomic adsorbates (i.e., O( 3 P) and N( 4 S)) are open-shell systems. Integration over the Brillouin zone was performed by using a 3x3x1 (in slab calculations) or 9x9x9 (in bulk calculations) k-points mesh by means of the Monkhorst-Pack method 18 . Geometrical optimizations and vibrational frequencies were computed with an energy accuracy of 10 -6 eV.
For adsorption studies we use a 1x1 surface unit cell and several slab models (4, 6 or 8 layers of SiO 2 ) for the b-cristobalite (100) face with either an oxygen or a silicon first layer. An additional hydrogen back layer was added to saturate the Si or O dangling bonds as is commonly accepted for silica and silica-containing oxides [19][20][21] . Thus, to ensure this behavior, the two last layers were always kept fixed while the others were fully or partially relaxed depending on the slab model used. The distance between slabs (ca. z = 17-18 Å) was large enough to prevent significant interactions between them.
In spite of the large lattice parameter of b-cristobalite, we have also carried out some supplementary adsorption calculations with a 2x2 supercell to check that the atomic coverage effect was very small in the calculated properties.
Several adsorption sites were characterized for O and N adatoms (see next section for details). Once determined the optimal geometry for each adsorption site, we calculated the Hessian matrix and its corresponding harmonic vibrational frequencies (ni) for the adatom (i.e., these can be approximately classified as two parallel and one perpendicular movements), keeping fixed the optimized slab geometry. Several calculations introducing a slab relaxation (i.e., 2 or 4 layers) were carried out to verify that the vibrational frequencies were essentially the same. This analysis allowed the complete characterization of minima (i.e., 3 real ni) and surface diffusion transition states (i.e., 1 imaginary frequency (1st-order saddle point) or 2 imaginary frequency (2n-order saddle point)).

The adsorption energy (Ead) was defined as
where E(atom) is the total calculated ground state energy of the atom in gas phase (i.e., O( 3 P) or N( 4 S)), E(slab) is the value for the relaxed clean slab in its singlet state and E(atom+slab) is the value of the full relaxed slab containing the adatom. Therefore, when we report a positive adsorption energy for a particular site (e.g., in tables) it means that this site is a minimum (stationary point), which is more stable than the atom + slab asymptote; on the contrary, a negative adsorption energy corresponds to another stationary point (usually a transition state) less stable than the atom + slab asymptote. Thus, minima should be considered as the true adsorption sites.
The atomic energies were calculated with the atom inside large broken symmetry boxes (i.e., 8x7.5x7.6 Å 3 ); we used some higher boxes (e.g., 16x16.5x16.6 Å 3 ) to be sure that we obtain the same converged energies. that the pseudocubic structure is the most stable ( (or ) < Fd3m < P213) even though the energy difference for the two first is very small (i.e., ˜ 0.05 eV). Moreover, they considered a bigger unit cell ( ) instead of the equivalent one.
We have also studied the properties of these structures to check the reliability of the present DFT method to reproduce b-cristobalite in bulk and slab models. Thus, Table 1 presents the calculated internal parameters and the cohesive energy of these bulk structures in comparison with available experimental data for the ordered structures. The calculated structural parameters compare well with earlier DFT results [30][31] and with experimental data. The major discrepancy is observed for dSiO in the Fd3m symmetry, where the calculated value is somewhat higher (i.e., 0.064 Å) than the experimental one. Respect the energy stability, we found that P213 is the less stable structure in agreement with some preceding results 31 , although on the contrary we find that Fd3m is a bit more stable that the structure. Nevertheless, the differences are within the expected error of the present DFT method, as it was previously discussed in section II (e.g., due to the O atom energy). In addition, we have performed a Fd3m structure optimization using the PBE functional obtaining that the structural parameters (e.g., a=b=c= 7. We have also studied the lowest triplet state to confirm that singlet state was the ground state in the different slab models. Thus, the calculated excitation energies (i.e., singlet-triplet splitting energy defined as Eexc=Etriplet-Esinglet) for fixed structures involving 4, 6 or 8 layers were 2.97, 2.99 and 3.01 eV, respectively, compared with 6.42 eV for the bulk calculation. These results corroborate that the singlet state is the ground state for silica. On the other hand, it is observed that the slab models underestimate this energy splitting. This fact agrees with the conclusions derived in an earlier study 34 on self-trapped excitons in a-quartz, which were triplet excited states that distort the crystal locally. The authors concluded that the PW91 functional underestimates the S-T splitting in cases where the triplet state is delocalized, as would correspond to our calculation.

Atomic oxygen adsorption on b-cristobalite
Firstly, we have studied the atomic oxygen adsorption on the ideal b-cristobalite (100) surface, with a first layer formed by Si atoms. Four different sites have been considered: on top Si (T1), a hollow (H) and two bridges (B1 and B2), as shown in Figure 1. Table 2  full relaxed structures in models with 6-8 layers. This means that the errors due to finite slab thickness will probably not be relevant. Hence, we have finally chosen the 6-layer slab model (relaxing the first-four layers) to be used in the rest of the work unless in one case shown later. Table 3 summarizes the results for the 4 sites. Two strong chemisorptions are found for T1 and B1 sites. In both cases, singlet and triplet states give place to large adsorption energies.
Moreover, the singlet state is much more stable than the triplet state in T1 site, whereas a slight opposite behavior is found in B1 site, though with much closer energies in this latter case. This change on the ground spin-state of the system depending on the adsorption site has been also reported for atomic oxygen adsorption over oxides (e.g., MgO(001) 35 ), metals (e.g., Pt(111) 36 ) or alloys (e.g., Pt/Ni(111) 37 ).
T1 and B1 structures are truly minima as we have verified by the analysis of harmonic vibrational frequencies ( Table 3). The different values observed for their vibrational frequencies show the different curvatures of the potential energy surface at these geometries, and can be explained analyzing the position of the adatom in each site and its interactions when the adatom is vibrating along the corresponding normal mode. Thus, for O adsorption on T1 site it is observed that the perpendicular frequency is much higher than the two parallel ones (Table 3)  The d Oad-Si values for T1 and B1 minima (Table 3) are quite close to the bulk b-cristobalite values (Table 1). However, only a slight relaxation of the silica slab is produced in T1, but an important reconstruction is observed in B1. In this latter case, the Si-Si distance decreases from 5.196 Å (clean slab) to 3.416 Å (singlet state) or 3.179 Å (triplet state), and the d12 interlayerdistance changes from 1.007 Å (clean slab) to 0.708 Å (singlet state) or 0.516 Å (triplet state).
The calculated adsorption energies for T1 and B1 sites (in this latter case its half value can be taken as the binding energy due to the formation of two SiO bonds) are slightly lower than the usual SiO bond dissociation enthalpies in gas phase (e.g., 6.50 eV in SiO2 or 8.29 eV in SiO at 25 ºC 32 ) or in b-cristobalite (e.g., 9.61 eV derived from its cohesive energy (Table 1)). In spite of the absence of direct experimental adsorption (or desorption) energies for the present system 38   The B2 and H structures correspond to very small adsorptions, mainly in the triplet state, even though these stationary points are not truly minima but transition states (i.e., first-order saddle points). In these structures there is no any significant geometrical relaxation with respect the clean slab geometry. We have also checked that the H structure seems to be a diffusion transition state that connects two equivalent B1 minima. Nevertheless, as surface reconstruction is important in B1 minimum, the analysis of the vibrational frequencies in H structure has to be taken as approximate due to the vibrational frequencies were calculated keeping fixed the slab geometry in both stationary points. In B2 structure this problem was not found and it was much easier to establish that it corresponds to a diffusion transition state for the T1 interconversion. atomic collisions over silica the initial direct adsorption in B1 and T1 minima will be more relevant because the initial slight adsorption in H or B2 sites will evolve to the T1 and B1 minima.
In previous section we discussed the bulk and slab models for clean b-cristobalite, and we proposed a slab model (relaxed) derived from a Fd3m bulk structure. To see the possible influence of the SiOSi angles of the slab in the O adsorption energies we have also made calculations by using a slab model derived from a Fdd2 bulk structure (Table 1), which essentially means than the SiOSi angles are close to 147º instead of 180º. The results are presented in Table 4. The comparison with the results shown in Table 3 indicates that the adsorption energies and the geometries are very similar for all structures, and both slabs seem to be almost equivalent for the present study.
We have also studied the possible adsorption of oxygen on the b-cristobalite (100) face terminated with O atoms. Three several sites were finally considered: on top O (T2) and two bridges (B3 and B4), as shown in Figure 2. In this case we had to use the 6-slab model with full relaxation to obtain converged results, which are shown in Table 5. In the first site (T2) we did not find an atomic adsorption but the abstraction reaction to produce the O2 molecule without energy barrier. In both bridge sites a chemisorption was obtained either in the singlet or in the triplet state, with similar adsorption energies and very small geometrical relaxation in comparison with the clean slab geometry. The stability of both B3 structures (singlet and triplet) is coherent with the most stable Lewis structures, with two electron lone pairs in each O atom for the singlet B3 structure and with two unpaired electrons in two of the three oxygen atoms for the triplet.
Nevertheless, B4 has two imaginary frequencies in both states, and it could be expected that a firstorder saddle point should be close to the reported geometries. We have also calculated the vibrational frequencies for B4 structure (triplet state), relaxing the two first layers of the slab. All  Table 5 (i.e., 547.7, 136.7i, 206.3i) and show that the existence of two imaginary frequencies is possibly due to the true transition state is slightly shifted from the reported geometry.
Therefore, the rigid slab approximation for the vibrational frequency calculations seem to be appropriate enough to characterize the stationary points.
The number and the energy of the stationary points found points out that the Si face on bcristobalite (100) will be more important than the O face for oxygen adsorption processes and probably also for the successive Eley-Rideal or Langmuir-Hinshelwood reactions, such as it has been assumed up to date in all kinetic modeling studies on air chemical processes over silica materials 1 .

Atomic nitrogen adsorption on b-cristobalite
We have carried out a similar study as for oxygen using now atomic nitrogen. The same sites (Figures 1 and 2) were initially considered for this system in the two lowest spin states derived from the N( 4 S) (i.e., doublet and quartet states).  Table 6, are within the range of the usual N-Si distances in silicon nitrides in gas phase (e.g., 1.57 Å in SiN ( 2 S + ) or 1.70 Å in Si 2 N ( 2 P g ) 32 ) or in nitride films on silicon (e.g.,  43 on the N adsorption on SixOy (silica clusters) reports an adsorption energy for the T1 site (i.e., 2.84 eV) and a N-Si distance (i.e., 1.7 Å) quite similar to our present values (Table 6), which is nonetheless based on periodic (slab) DFT-GGA/PW91 calculations.
The difference between doublet and quartet SiN binding energies in T1 structures (Table 6) could be justified in a similar way as it was made for O adsorption structures. Thus, the doublet state allows the possibility of a Si=N double bound as N can keep one unpaired electron; conversely, the quartet state makes even difficult the formation of a single bond due to the difficulty for N or Si (Si is also bonded to two O atoms of the slab) to have three unpaired electrons and simultaneously to make a bond, which is consistent with the unimportant bonding in quartet state for T1 structure. Similar arguments can be valid for B1 structures.
The stationary point corresponding to N adsorption in H site is a diffusion transition state that connects both equivalent B1 minima (i.e., B1 ® H ¹ ® B1) with high energy barriers of 3.60 eV (doublet state) and 6.61 eV (quartet state). The B2 structure turns out to be a diffusion transition state connecting both equivalent T1 minima (i.e., T1 ® B2 ¹ ® T1) with energy barriers of 3.96 eV and 0.04 eV for doublet and quartet states, respectively. This behavior is similar as we observed for O adsorption, and does not exclude the existence of even minor diffusion energy barriers for another directions (e.g., line T1-B1 in Figure 1).
The adsorption of N into the first oxygen face on the b-cristobalite (100) shows as well similar sites ( Figure 2) and behavior as for O adsorption, and its results are indicated in Table 7.
Thus, it is also observed a molecular formation (i.e., NO) in T2 site for both spin states instead of the N adsorption process, with energy barrierless in both cases. Moreover, an important chemisorption is obtained in B3 site (true minimum) for both spin states. Their corresponding NO binding energies (6.10/2 eV and 5.66/2 eV in doublet and quartet states, respectively) are within the ordinary range of the NO bond dissociation enthalpies in gas phase (e.g., 3.18 eV in NO2( 2 A1) or 6.55 eV in NO( 2 P) at 25 ºC 32 ). The Nad-O distances (1.336 Å and 1.442Å eV in doublet and quartet states, respectively) are rather longer than the NO distances in small gas phase molecules (e.g., 1.197 Å in NO2( 2 A1) or 1.151 Å in NO( 2 P) 32 ); only a very small slab relaxation is produced as in the case of O adsorption. In B3 minima, N is bonded to two oxygen atoms, which are also bonded to the same tetra-coordinated Si atom (i.e., SiO4 group in Figure 2). The slight difference in their binding energies (derived form Ead values shown in Table 7) can be explained again taking into account the major facility of the N atom to keep one unpaired electron (and a lone pair) when is also simultaneously making two NO bonds in the doublet state, which seems a bit more difficult for the quartet state.
The B4 stationary point corresponds to a diffusion transition state for the B3 interconversion (i.e., for B3 ® B4 ¹ ® B3), with a very small energy barrier of 0.18 eV for the quartet state.
However, this transition state has a slightly lower energy than B3 minima for the doublet state; even though B3 minimum should be deeper than B4 transition state, the closeness of their energies could point out that B3 is a local minimum and the most stable one should be located not very far.

Atomic sticking coefficients and desorption rate constants
The and (3) where mA is the A atom mass, n½½ and n^ are the two parallel and one perpendicular vibrational frequencies of the adatom A, respectively, and the rest of symbols have their usual meanings. In both equations we do not introduce the catchall factor (i.e., Pdes = 1), which often is taken as an adjustable parameter in some kinetic models 1 . Our calculated TST values agree with the low So values derived in both kinetic models.
Moreover, its temperature dependence is much closer to model II shape, which is in between T1 and B1 curves for both systems.
In Figure 5 45 , which would be much closer to our current DFT results.
These high adsorption energies imply that the chemisorbed atoms will be desorbed only at relatively high temperatures (T > 900 K for N and T > 1800 K for O from Figure 5), so that these atoms will be held irreversibly on the surface and they will be removed only by atomic recombination (e.g., ER or LH processes) at lower temperatures. Experimental studies by using single-crystal adsorption calorimeter 38,46 or thermal desorption spectroscopy 46 would be necessary to validate entirely our theoretical results. The calculated So(A) and kdes (A) for N and O would tend to increase the corresponding gA coefficients if the kinetic data 6 for ER and LH were kept constant.
However, a DFT study on the minimum energy paths and the transition states involved in ER and LH reactions (in progress in our group) will be necessary to obtain an accurate calculation of gA(T) coefficients.
Finally, despite the assumptions made in the So(A) and kdes (A) calculations: TST equations, single-crystal, Si face, small effect of the coverage in DFT adsorption energies,.., we believe that our values are more suitable to be used in the mentioned kinetic models than the usual adjustable (semiempirical) parameters not derived from first principles, which differ considerably among the numerous empirical kinetic studies made up to date.

Summary and conclusions
In this work we have made a periodic DFT study (GGA/PW91) on the adsorption of atomic oxygen and nitrogen over b-cristobalite (100) considering the two lowest spin states The approach of either O or N on top oxygen of b-cristobalite (100) gives place to the abstraction reactions, which produce O 2 or NO molecules, respectively, without energy barrier.
Several diffusion transition states have also been characterized, which connect some of the adsorption minima with rather high energy barriers; this does not exclude the existence of much lower energy barriers in another directions of the cell (e.g., in the line between B1 and T1 minima).
Atomic sticking coefficients and desorption rate constants have been estimated within the temperature range 300 -1900 K by means of the usual transition state equations and the standard approximations. The high adsorption energies found for O and N over silica imply that the atomic recombination processes (i.e., Eley-Rideal and Langmuir-Hinshelwood mechanisms) will have a more important role in the atomic detachment processes than the thermal desorption processes at temperatures lower than 1800 K for O and 900 K for N. Further experimental data should confirm these conclusions along with DFT calculations on these reactions, in progress in our research group.
Finally, the magnitude of the DFT adsorption energies also suggests that the published kinetic models of atomic O and N recombination on SiO 2 surfaces, which make use of low adsorption energies (e.g., 3.5 eV equal for O and N 1 ), should be revised, particularly for oxygen processes. It can be expected that the structural differences in silica based-materials (e.g., pyrex, cristobalite, RCG,..) will not give a large divergence in adsorption energies and hence that the present values could be still suitable for them.