Recombination and chemical energy accommodation coefficients from chemical dynamics simulations : O / O 2 mixtures reacting over a b-cristobalite ( 001 ) surface a

A microkinetic model is developed to study the reactivity of an O/O2 gas mixture over a b-cristobalite (001) surface. The thermal rate constants for the relevant elementary processes are either inferred from quasiclassical trajectory calculations or using some statistical approaches, resting on a recently developed interpolated multidimensional potential energy surface based on density functional theory. The kinetic model predicts a large molecular coverage at temperatures lower than 1000 K, in contrary to a large atomic coverage at higher temperatures. The computed atomic oxygen recombination coefficient, mainly involving atomic adsorption and Eley-Rideal recombination, is small and increases with temperature in the 700-1700 K range (0.01 < gO < 0.02) in good agreement with experiments. In the same temperature range, the estimated chemical energy accommodation coefficient, the main contribution to which is the atomic adsorption process is almost constant and differs from unity (0.75 < bO < 0.80).


Introduction
The reentry of space vehicles into Earth's atmosphere, involves hypersonic speeds. The corresponding hypersonic flows are usually characterized by strong shocks leading to equilibrium or nonequilibrium molecular plasmas [1][2] . Elementary processes involving atomic/molecular oxygen and nitrogen in the boundary layer are catalysed by the materials used as thermal protection systems (TPS) of these vehicles, increasing the total heat flux to their surfaces. The seeming simple heterogeneous processes that happen in these extreme conditions (e.g., O adsorption, O 2 dissociative adsorption, O recombination via an Eley-Rideal mechanism,..) are even nowadays difficult to study experimentally or theoretically [3][4] . Computer fluid dynamics (CFD) simulations are usually carried out to predict the relevant aerothermodynamics 1 . In these simulations two coefficients are required to include appropriately the wall catalycity effects: a) the atomic recombination coefficient g i (also called catalytic efficiency) [3][4] , which is the recombination probability of the i atomic species on the surface and b) the chemical energy accommodation coefficient b i 3,4,5 , defined as the ratio of energy released to the surface per atomic recombination to the maximum energy transferable. Both coefficients are within the range 0 ≤ g i , b i ≤ 1 (for air i will be essentially O and N) and are supposed to depend on temperature and total/partial pressures (i.e., g i (T, P)) although more detailed description would be required for CFD accurate simulations of nonequilibrium flows (e.g., g i (T gas , T surface , P)). Nevertheless, many experiments only measure routinely the effective catalycity (g eff, i = g i ¢ = g i • b i ).
Experimental g i and b i values are characterized by a broad scatter 4 since they depend on the conditions under they were measured. According to this, the use of these coefficients does not properly describe the heat transfer over the whole surface and the entire trajectory of the vehicle. The important effect of g i coefficient on the heat flux has been shown in several simulations 3,6 ; some studies yield up to a factor of two for the wall heat flux of a fully catalytic wall (i.e., g i = 1) wall compared to a noncatalytic one (i.e., g i = 0). The b i coefficient is hardly measured and often it is taken equal to 1 (i.e., fully energy accommodation assumption), although available values are very far from one 4 . Therefore, more experimental or theoretical g i and b i values would be necessary to better design thermal shields for spacecrafts, which would also allow a payload improvement.
Recombination and surface activity data for silica-based materials are important as those are commonly used in TPS (e.g., in Space Shuttle tiles). Measurements of oxygen and nitrogen coefficients for air and O 2 or N 2 pure gases over SiO 2 -based surfaces (i.e., g i , g i ¢, b i on quartz, cristobalite, RCG, vitreosil, pyrex, silicon carbide,.. ) have been recently reviewed 4,7 . g i coefficients show strong temperature dependence, increasing with temperature until a critical value where the thermal desorption becomes significant and hence the recombination coefficient declines. In particular, there are measurements of g O 8 and b O 9 for air at 200 Pa over b-cristobalite at high temperatures (800-1830 K), which were obtained by in situ oxidation of sintered SiC samples.
A large number of kinetic models for dissociated air, O 2 and N 2 gases on silica like surfaces have been developed in the past 6,10,11,12,13,14,15,16,17 . These models take into account the main surface processes: physisorption, chemisorption, dissociative adsorption, thermal desorption, Eley-Rideal (ER) or Langmuir-Hinshelwood (LH) atomic recombinations,… The microscopic parameters used in these models (e.g., desorption energies, activation energies, density of adsorption sites,…) may be either derived from quantum-mechanical calculations or from semiempirical formulas. However, such parameters usually are estimated by fitting the experimental g i (T) coefficients or the heat fluxes 3 .
Thus, g i (T, P) curves are reported although different sets of the kinetic parameters can predict the same heat flux. Apart from macroscopic one-temperature models, more advanced state-to-state models have been also used in the description of the nonequilibrium chemical kinetics of high temperature plasmas 18 . The state-to-state approach describes any internal degree of freedom of the molecules (vibrational, rotational, electronic). Each level is considered as an independent species subject to an appropriate continuity equation, own cross sections and rate coefficients. From the level distributions, the thermal properties of the gas and the global rate coefficients are derived, which may significant differ from those attained by a thermal equilibrium or weak nonequilibrium analysis. Thus, nonequilibrium vibrational kinetics of an O 2 /O mixture 19 or air 20 hitting a catalytic silica surface has been investigated.
In a previous study of the O/SiO 2 system, we carried out Density Functional Theory (DFT) calculations of O adsorption over a b-cristobalite (001) surface, showing a strong chemisorption 21 .
Quasiclassical trajectory studies (QCT) of atomic oxygen on a clean 22 or O-precovered 23,24 bcristobalite (001) surface reveal the importance of the atomic adsorption or penetration into the SiO 2 slab (absorption) along with the formation of adsorbed (O 2(ad) ) and gas phase (O 2(gas) ) molecular oxygen; these latter species mainly were formed through an Eley-Rideal (ER) mechanism. The O 2 scattering on a clean b-cristobalite (001) surface 24

Quasiclassical trajectory studies
We have performed a QCT dynamics study 31 of atomic and molecular oxygen colliding with a clean b-cristobalite (001) surface. We have also studied the O collisions with an O-precovered bcristobalite (001) surface to investigate the ER reaction. These processes were previously studied at several state-specific conditions (e.g., initial collision energy (E col ), internal O 2 (v,j) states,..) and only in few thermal (or quasithermal) conditions 22,23,24,25 . Here we increase the temperature range (300 -1600 K) and improve the statistics of the calculated reaction probabilities. The present study makes use of the previous developed PES for the global O 2 /b-cristobalite system 23 , based on DFT data interpolated via by the corrugation reduction procedure 32 . In an earlier work 21 we extensively showed that spin-polarized DFT-GGA-PW91 calculations seemed to be accurate enough to describe not only the bulk properties of silica but also the interaction of oxygen and nitrogen over several b-cristobalite faces. The interpolated PES was based on a dense grid of DFT-GGA-PW91 points calculated for several O and O-O configurations over the 1x1 surface unit cell 23 , taking the lowest energy state (essentially singlet at short surface distances and triplet for longer ones). An empirical silica slab was also used to introduce the slab motion. Therefore, the present multidimensional PES allows both the oxygen and b-cristobalite movements. A wide description of this adiabatic PES (e.g., expressions, parameters, plots,,..) are reported in previous papers 22,23 . We have used the same PES (available upon request) in the present study although we have added a small correction in two subroutines to make sure that the Cartesian gradients were fully periodic for atom/slab calculations with cells larger than the 1x1 one. Previous studies should not be affected by this fact as we used also the 1x1 unit cell, mainly for normal incidence calculations.
For O collisions with the clean surface we keep the second O atom at a large enough distance to the surface to avoid any undesirable interaction.
The QCT description has been widely presented elsewhere [22][23] . Briefly, the b-cristobalite slab used is formed by 104 atoms distributed into 9-layers for a 2x2 surface unit cell with a lattice parameter of 7.348 Å. The slab temperature (T S ) was controlled by means of a Generalized Langevin equation approach 31,33 . The initial position of the incoming atomic or molecular oxygen was sampled over an 1x1 surface unit cell, with an initial distance higher than 6 Å to avoid significant surface interaction.
For the ER study the adsorbed O atom (O ad ) was located over the central Si atom of the first layer (T1/T1" sites, Fig. 1) at an initial z distance of around 1.5 Å, being also thermalized together with the slab. Since the experimental measurements of g O on b-cristobalite 8 correspond approximately to a flux of air colliding normally to the surface, we have restricted our simulation to normal incidence (q v = 0º) for all atomic and molecular collisions. Thus, as we are interested in quasithermal conditions (i.e., T = Batches of 5000-15000 trajectories were integrated at each temperature to ensure a reasonable good statistics for the reaction probabilities of the main channels. Different processes could be observed: atomic and molecular reflection, atomic and molecular sticking (adsorption + absorption), O 2(ad) formation and O 2(g) formation (i.e., ER reaction). Classification criteria for these events were broadly discussed in earlier works 22,23 . These were based mainly on the analysis of Z O , Z O' , Z CM (molecular centre of mass Z coordinate) and r (O-O distance) variables, and the corresponding z velocity components of the atoms or the molecule (v z CM ) after a minimum collision time (e.g., t = 2 ps). Thus, molecular oxygen (gas) was formed by the ER reaction when Z CM was much larger than 6 Å, the r distance was shorter than the double of the equilibrium distance (i.e., r < 2.5 Å) and v z CM > 0, whereas O 2(ad) was classified when Z CM < 3.0 Å with r < 2.5 Å. For this latter case, the integration was really finished if after some extra steps (e.g., 5-10) the classification was preserved. This criterion was especially useful for channels involving finally adsorbed/absorbed species. Nevertheless, some trajectories were ended up at shorter collision times (t < 2 ps) due to some PES limitations for extrapolation (e.g., the PES is not defined for Z CM < 0, Z O < -1.0 A or r < 0.5 Å), although we think that the final thermal probabilities are even pretty well converged and accurate enough for rate constant calculations. as new processes such as O 2 formation appear. In fact, we distinguish between the formation of the O 2(g) via the ER reaction and the production of molecular oxygen that is strongly retained on the surface (O 2(ad) ). This second reaction channel is much more probable than the first one, which is shown in a semilogarithmic scale in Fig. 2d. This trend was also observed in similar studies of atomic oxygen recombination over quartz 26 and in previous molecular dynamic simulations for O/O 2 over b-quartz 29 .
In this latter study, O 2(ad) formation was established if adsorption was longer than 5 ps, although simulations at very long times (i.e., 1.5 ns) showed important O 2 desorption, which increased the final ER probabilities and hence their g coefficients, although possibly this should be classified instead as molecular desorption. The comparison of the present probabilities ( Fig. 2c and 2d) with our previous ones 23 show small differences (i.e., lower O 2 sticking and consequently larger O sticking for O collisions over an O-precovered surface, because the sum of probabilities is the unity) mainly originated by the introduction in our QCT code of a different classification of O 2(ad) events (e.g., Z CM < 3.0 instead of 1.9 Å, which was used in the previous study 23 ). The use of longer integration times in our trajectories (i.e., t > 2.5 ps) could also increase a bit the final O 2 desorption although the mentioned PES limitations (i.e., energy extrapolation at very short O or O 2 distances to the surface) avoid to describe more accurately such process.

Surface processes
We propose a first microkinetic model to ascertain the global effect of the main heterogeneous processes involving different O/O 2 mixtures over a b-cristobalite (001) surface at quasithermal conditions. Ten surface elementary processes are included, where i label will be used for the direct processes and -i for their reverse ones (i = 1,..,5). These are: chemisorption (k 1 ) and desorption (k -1 ) of oxygen atoms, molecular nondissociative adsorption (k2) and molecular desorption (k-2), O2(ad) formation by atomic recombination (k3) and molecular dissociative desorption (k-3), Eley-Rideal reaction (k4) and molecular dissociative adsorption (k-4), (4) and Langmuir-Hinshelwood reaction (k5) and molecular dissociative adsorption (k-5), where s indicates a free site on the surface. DFT studies show that molecular oxygen is adsorbed in parallel over a Si atom, using two free adjacent sites of the four available in each (1x1) unit cell 23

Thermal rate constants
We have calculated the thermal rate constants for the elementary processes (1), (2), (3) and (4) from the corresponding thermal QCT reaction probabilities (P i (T)) reported in section 2 (Fig. 2) by using conveniently the following standard expressions 6,12 , where k B is the Boltzmann constant, T is the temperature, m l is the atomic mass (molecular mass for k 2 ) and [s] 0 is the initial surface density of free sites (i.e., 4 sites per unit cell or 7.41x10 18 sites •m -2 for a b-cristobalite (001) surface). The i values in the eqn (6) and (7) show the elementary processes for which these rate constants are calculated.
We have used fittings of the QCT probabilities by using some polynomials for (1), (2) or (4) processes or an Arrhenius expression for ER reaction (3), which are also plotted in Fig. 2. As the LH reaction was not studied by the QCT method, we have estimated its rate constant from the following equation where ∆ is the mean distance between active sites, the square root term is the mean velocity of O reacting adatoms, assuming a two-dimensional gas and is the LH energy barrier including zero point energy, derived from our PES data (i.e., 3.29 eV). The rate constants for reverse processes (k -i ) were obtained by means of the principle of detailed balance, using the DFT reaction exo-or endothermicities of each process ( , including zero point energies). We have not introduced spin statistical factors into the rate constants. This assumption could be acceptable as we use the same adiabatic PES for all O 2 /b-cristobalite processes and because these factors should produce only small and similar changes into the calculate rate constants. Fig. 3 shows the calculated QCT rate constants within the 400-2000 K range of temperature. The great differences in these rate constants will have an important influence on the global recombination kinetics of the system. It is for instance expected that LH reaction contributes much less than ER reaction to the global recombination within this model.
The rate constants are finally fitted by using an Arrhenius formula within the 700-1700 K range, whose parameters (A i , E a,i ) are indicated in Table 1 together with values. This more restricted temperature range is selected to obtain better Arrhenius parameters to be used for the following g O and b O calculations; the initial T interval is too wide (400-2000) so that significant curvature is observed at lower temperatures.

Surface coverage and g O recombination coefficient
The time evolution of the surface concentration of both atomic and molecular adsorbed species resulting from the microkinetic model is given by the solution of the coupled equations (9) and (10) where as each O 2(ad) uses two sites, the density of free adsorption sites at any time will be (11) Fig. 1 neighbours (i.e., two O on the same Si) will be effective for several reactions (e.g., two adjacent sites inside the unit cell, which means one of the three pairs of sites, are necessary for reaction (2)). In fact, this small correction introduces some correlation between the locations of the reactants on the surface that a phenomenological kinetics approach does not take into account as compared with more accurate kinetic Monte Carlo simulations 34 .
The numerical integration of these differential equations ( however, the formation of O 2(ad) by the competitive processes (2) and (3), which have similar rate constants (Fig. 3), increases mainly with the raise of P O2 and P O , respectively, hence alternating their influence for different initial mixtures. For pure O mostly atomic coverage is observed.
The difference between both behaviours are mainly originated from the temperature dependence of both molecular desorption processes (-2) and (-3). At very high temperatures their rate constants (k -2 , k -3 ) are very large (being k -2 >> k -3 ) and only a final atomic adsorption is observed (e.g., qO = 0.99 at 1500 K, Fig.4b), being negligible the final O 2(ad) formation by (2) or (3) processes. At lower temperatures atomic oxygen is initially adsorbed until a suitable atomic coverage (maximum in Fig.   4a), when reaction (3) start to produce O 2(ad) , which is not desorbed enough because both k -2 and k -3 (also with k -2 >> k -3 ) are now much lower (e.g., qO = 0.11 at 900 K, Fig. 4a). The calculated atomic and molecular coverages are much larger that the ones observed for graphite (0001) (qO < 10 -3 and qO2 » 0) 30 , which can be justified by the very strong atomic and molecular chemisorptions on b-cristobalite.
We have also calculated the atomic oxygen recombination coefficient gO as the ratio of the total flux of recombining oxygen atoms (O(g) consumption and formation) to the initial flux of impinging oxygen atoms (ZO) over the surface once achieved the steady-state concentrations, using the following equation,  Small recombination coefficients are observed (0.01 < gO < 0.02) in the 700-1700 K range, which are essentially originated from atomic adsorption and ER reaction with a similar proportion. No LH contribution is observed at all. In spite of the reaction (3) competes with ER reaction (4) with much higher reaction probabilities ( Fig. 2c and d), its lower rate constants (i.e., k3 < k4; see Fig. 3a and c) together with the small amount of free surface sites (see [s]ss in Fig. 4), can justify its negligible contribution to gO.
These two contributions in gO are in agreement with previous assumptions made in similar kinetic models 5, 37 . At very large pressures the LH reaction seems to gain importance for b-quartz 29 , leading to larger gO values (0.01 ≤ gO ≤0. 35). In a previous work 23 we directly compared the experimental gO values with QCT ER reaction probabilities, obtaining similar values although with a slope closer to the experimental one than that of the present work. The seemingly better previous results could be fortuitous because another processes (e.g., specially the atomic adsorption) should be included in a more reliable calculation of this coefficient. On the other hand, as the calculated gO coefficients are based on QCT reaction probabilities, some improvements on those could increase slightly the gO values. Thus, the inclusion of physisorbed O atoms (not only chemisorbed on T1/T1" sites) in the ER study would probably increase the O2 formation. Moreover, the integration at longer times (e.g., 2.5-5 ps) for reaction (3), which produces adsorbed O2, could increase the molecular desorption with following final augment of the ER contribution. For instance, simulations for b-quartz 29 were ran at least for t < 5 ps, showing also a significantly increasing of gO but possibly at too long times (i.e., 1.5 ns).
Finally, although microkinetic models are commonly used to simulate the heterogeneous chemical kinetics of dissociated airflows impinging different types of surfaces 16,36 with good agreement in comparison to available experimental data (e.g., total heat flux to the surface, g coefficients,..), more developed Monte Carlo simulations (e.g., N over silica 41 , H over graphite 42 ) could also be needed to provide a deeper microscopic understanding of these processes.

b O chemical energy accommodation coefficients
Following a similar approach as for the calculation of g coefficient, we can calculate the b coefficient once achieved the final steady-or quasisteady-state from the ratio of the total chemical energy that can be transferred to the surface to the energy that would be released if all consumed O atoms gave place to O2(g), using the expression (14) where DO2 is the O2 dissociation energy, Fi represents the atomic or molecular fluxes at the steadystate (shown in Table 1) and Qi are the energies per O atom required or released in each process, which are derived from the expressions and activation energies (Table 1). We assumed that adsorbed species (e.g., process (1)) release all of their available energy to the surface and that gas products (e.g., process (4)) will escape with all of their available energy (i.e., assuming a very short interaction with the surface), which would give a minimum b. This latter point is supported by previous QCT study of the ER reaction on b-cristobalite 23 , which showed that formed O2 molecules become translationally and internally excited, taking rather some energy from the surface. Similar expressions were used in earlier studies 5 contributions to bO is shown in Fig. 8, which displays the total and the partial energy fluxes originated by the main processes. Clearly, the atomic adsorption process is largely the one producing almost the total energy released to the surface.
Up to our knowledge there is only a communication 43

Summary and conclusions
In this work we present a theoretical study of a gas mixture of atomic and molecular oxygen reacting over a b-cristobalite (001) surface, which is relevant for the simulation and the understanding of the aerodynamic heating produced during the reentry of spacecrafts into Earth's atmosphere.
We use a previous multidimensional potential energy surface based on DFT data to determine thermal rate constants for atomic or molecular oxygen colliding with a clean surface and for atomic oxygen colliding with an O-precovered surface, by means of quasiclassical trajectories.
The present dynamical studies confirm earlier conclusions: atomic sticking over the clean surface is the predominant process but decreases for an O-precovered surface due to the appearance of two competitive processes that produce O2(ad) and O2(g). Moreover, molecular nondissociative sticking is also observed.
We propose a microkinetic model with ten heterogeneous processes to study this system for An estimation of a minimum chemical energy accommodation coefficient has been also presented, obtaining an almost constant value (0.75 < bO < 0.80) within 700-1700 K; therefore the common assumption that bO = 1 and hence gO' = gO seems inaccurate for O recombination over b-cristobalite.
The atomic adsorption process is mostly the one producing almost the total energy released to the      Pa) at (a) 900 K and (b) 1500 K.             Table 1 Arrhenius parameters derived from the calculated thermal rate constants using analytical P i (T) QCT probabilities within 700-1700 K. are the exo-or endothermicities energies. Q i energies and F i fluxes are used later for b calculations.
reaction process A a E a,i (eV) (eV) b Q i c F i (m -2 s -1 ) d  c Energies by O atom required (< 0) or released (>0) to the surface for each i process.
d Partial atomic or molecular fluxes for each reaction process.