A molecular dynamics simulation of hydrogen atoms collisions on an H-preadsorbed silica surface

The interaction of hydrogen atoms and molecules with a silica surface is relevant for many research and technological areas. Here, the dynamics of hydrogen atoms colliding with an H-preadsorbed β-cristobalite (0 0 1) surface has been studied using a semiclassical collisional method in conjunction with a recently developed analytical potential energy surface based on density functional theory (DFT) calculations. The atomic recombination probability via an Eley–Rideal (E–R) mechanism, as well as the probabilities for other competitive surface processes, have been determined in a broad range of collision energies (0.04–3.0 eV) for off-normal (θv = 45°) and normal (θv = 0°) incidence and for two different surface temperatures (TS = 300 and 1000 K). H2,gas molecules form in roto-vibrational excited levels while the energy transferred to the solid surface is below 10% for all simulated conditions. Finally, the global atomic recombination coefficient (γE–R) and vibrational state resolved recombination coefficients (γ(v)) were calculated and compared with the available experimental values. The calculated collisional data are of interest in chemical kinetics studies and fluid dynamics simulations of silica surface processes in H-based low-temperature, low-pressure plasmas.


Introduction
The interaction of plasmas with surfaces can lead to different surface processes active in different collisional energy regimes according to the behavior of the gas-phase species and substrate involved. In particular, atom recombination at surfaces can be an effective source of roto-vibrationally excited molecules and, at the same time, an effective process for surface atom abstraction and atom removal from the plasma region [1]. Such processes have a strong impact on the chemistry of the bulk region and at the plasma-surface interface. It is, in fact, well-established that the reactivity of molecular plasmas under low-pressure, low-temperature conditions depends on, and is often controlled by, the formation of energetically vibrationally activated molecules [2][3][4][5]. Therefore, it is important to understand the processes by which one can store vibrational quanta in the vibrational manifold of the active molecules in the plasma or gaseous media. There are many such processes and they are quite often difficult to describe and identify both from a theoretical and experimental point of view; although the identification and description of these processes at a molecular dynamical level has long been the objective of many research groups active in the wide area of plasma-chemistry and plasmaphysics [6][7][8].
This paper is focused on the molecular surface processes due to the interaction of a flux of H atoms colliding with a silica surface at low-collisional energies.
Hydrogen plasmas are currently of great interest in different research areas of fundamental and technological interest such as micro-electronics, modern solar-cells, fusion reactors, astrophysics and aerothermodynamics. Fluxes of hydrogen atoms on silicon or silicon dioxide surfaces are used in the deposition of thin solid films, material growing and etching under low-pressure plasma conditions where physical and structural properties must be well defined [9]. Hydrogen atom interaction with interstellar solid grains (graphite-or silicate-based surfaces) is of interest in astrophysics for the study of the chemistry of the early universe, and, in particular, to explain the great abundance of molecular hydrogen observed in the interstellar media [10]. Hydrogen is also the main constituent of Jovian atmosphere so that, in view of a possible future Jovian exploration, one of the major technological challenges is to plan the safe entry of spacecraft into Jovian's atmosphere by protecting the walls of the spacecraft from the heat fluxes produced by the heterogeneous recombination reactions of H atoms.
There aren't many theoretical calculations and experimental determinations for H atom recombination on silicon and silicon dioxide substrates, so that basic collisional data related to this heterogeneous process, such as the recombination probability, the related global recombination coefficient, and the roto-vibrational distribution of the H 2 molecules formed at surfaces, are almost unknown. On the other hand, these and other collisional data are also of fundamental importance for kinetic modeling and computer fluid dynamics (CFD) simulations of the plasma-chemistry occurring at the plasma-surface interface. The important effect of the recombination coefficient (γ E-R ) and of vibrational state resolved recombination γ (v) coefficients on the heat flux has been already shown in several simulations [20][21][22][23][24][25][26]. The lack of collisional data for the system under study, and in general for heterogeneous reactions in plasmas, is mainly due to the highly specificity of the wall processes. In fact, reaction probabilities depend on properties of incident flux (kinetic energy, approaching geometry, etc.) and on the surface properties (temperature, morphology, coverage, etc.). However, some caution must be taken in the comparison of collisional data obtained via molecular dynamics (MD) and those determined in bulk experiments due to the fact that under plasma conditions different processes can act simultaneously and influence each other [1].
This work is devoted to MD simulation of hydrogen molecules' formation by H atom recombination with an H atom preadsorbed on a silica surface assuming the Eley-Rideal (E-R) recombination reaction mechanism. The silica polymorph used is β-cristobalite, since it is the most stable silica crystalline polymorph at high temperature up to the melting point (i.e., 1983 K).
The E-R recombination mechanism is a two-step process. First, an H atom is adsorbed on the silica surface in a specific active surface site: then a hydrogen atom approaching the surface from the gasphase reacts with the preadsorbed atom in a single or multiple collision(s): thus forming a hydrogen molecule in a specific roto-vibrational state. E (i) is the exothermicity of the process (i). In our study, the rate-determining step (2) is assumed in the MD simulation. It is worth noting that for the E-R reaction to occur the atomic adsorption process (i.e., equation (1)) is a mandatory condition and this fact should be taken also into account depending on the experimental setup for a good comparison between the theoretical and the experimental recombination coefficients.
The energetics of the recombination reaction (2) is of special interest. The energy released in this process can be shared among the internal degrees of freedom of the newly formed H 2 molecules (i.e., electronic transitions, molecular rotations, vibrations and translations) and the degrees of freedom of the silica substrate (phonons and/or electrons). It is of great importance to know, on the one hand, the extent to which the desorbing H 2 molecules can be formed in vibrationally and rotationally excited states and, also, the heat flux transferred to the silica surface. Indeed, this issue, i.e. energy sharing mechanism, is among the most complex and intriguing problems relevant to recombination reactions on surfaces.
Experimentally, the global recombination coefficient γ for H reacting over fused quartz and Pyrex has been formerly determined within the interval 194-1250 K [14] showing poor Arrhenius behavior with a small activation energy (E a = 0.42 ± 0.06 eV) that would correspond mainly to the E-R reaction (γ E-R ). The recombination coefficient at room temperature for hydrogen atoms in low-temperature plasma on fused silica was also determined by a pulsed plasma excitation technique in [27].
In this work we focus mainly on the study of reaction (2) in order to obtain γ E-R . Nevertheless, all the elementary processes originated by the collisions were also analyzed and discussed. The energetics of H 2,gas (v, j ) molecules formed in reaction (2) is solved by determining the energy transferred in the internal degrees of freedom of the molecule. The rotovibrational distributions as well as the average vibrational quantum number are determined together with the vibrational state-to-state recombination coefficient γ (v).
The reaction dynamics has been determined using a semiclassical dynamics method [28]. The procedure followed can be summarized in three main steps: (1) building-up of a 3D model of silica surface for which the phonon structure is determined by performing a normal mode analysis of the silica lattice atom vibrational motions; (2) analytical representation of the Potential Energy Surface (PES) obtained from electronic structure DFT calculations; (3) application of the semiclassical method to describe the dynamics of hydrogen atoms (one in gas phase and the other one preadsorbed) interacting with the silica surface.
We provide only a few details of steps (1) and (2) because these two points have already been extensively described in previous papers. The surface model used (step 1) was already considered in [16] while the details of the DFT calculations as well as the analytical PES construction were reported in [29]. Thus, the present study is mainly focused on a detailed account and discussion of the main dynamical results (step 3). This paper has the following structure: in section 2 some minor details about the recently published analytical PES for the interaction of H gas and H 2,gas with a β-cristobalite surface are reported; in section 3 a brief overview of the dynamics method used is also given; in section 4 MD calculations are described, whereas in section 5 the principal results obtained for the E-R reaction together with other competitive processes for off-normal and normal incidence and at different surface temperatures (i.e., T s = 1000 and 300 K) are presented.
The last section is devoted to the summary of the main results.

3D solid surface model and analytical PES
The recently published analytical PES [29] is based on periodic DFT-GGA calculations performed by means of the Vienna ab initio simulation package (VASP) [30][31][32][33][34]. A (1 × 1) unit cell was used for the adsorption studies along with a slab model formed by six layers of SiO 2 with an outermost silicon layer. An additional hydrogen back layer was added for saturating the dangling bonds of terminal oxygen atoms. Details of the calculations can be found in [29]. The calculations were based on the generalized gradient correction (GGA) functional Perdew-Burke-Ernzerhof (PBE) [35]. The previous DFT study [29] reports adsorption energies for several adsorption sites (i.e., T1, T1 , B1, S2, B2, H1 drawn in figure 1, obtaining results that the most stable adsorption occurs on T1 site (E ad = 2.10 eV), where the H-Si bond is tilted to some extent with respect to the top T1 site (α = 55.5 • , φ = 45.0 • ) [29], being d(H ad − Si) = 1.51 Å. The adsorption on the other sites is less favorable (E ad = 1.08, 0.50, 0.28, 0.01 and 0.009 eV for T1, S2, B1, H1 and B2 sites, respectively).
The analytical PES [29] was made up of sums of pairwise atom-atom interactions between each of the gas phase hydrogen atoms (H gas (1) and H gas (2) , hereafter H (1) and H (2) ) and each of the Si and O atoms that form the β-cristobalite 3D model. The β-cristobalite (2 × 2) 3D slab model consisted of 149 atoms distributed in ten layers according to the atom placement in the crystallographic Fd3m unit cell (nine of these layers are shown in figure 1).
The phonon dynamics of the surface atoms were then developed by numerical diagonalization of the 3D-dynamical matrix of the force constants obtained assuming the Born-Mayer-Huggins as the interaction potential of the Si and O lattice atoms [36]. The phonons' eigenvalues and eigenvectors that enter into the scattering equations of motion are then obtained while the calculated phonon frequency spectrum [16] is compared with the available literature data [37], leading to a good agreement.
The complete expression used for describing the 6D potential of the two atoms (i.e., 1, 2) interacting with a rigid β-cristobalite surface was defined as: This equation is formed by the atom-surface terms (V 3D ), one for each hydrogen atom, the H-H gas phase interaction term (V 1D ) and the term that accounts for the two hydrogen atoms interacting with the surface (V 6D ). The V 3D term is split in two gas-solid pair-wise atomic interactions (i.e., H-Si and H-O): corresponds to the i-th hydrogen atom Cartesian coordinates and R ij stands for the distances between the i-th gas phase hydrogen atom and the j -th Si or O atoms that belong to the 3D slab model. This term fitted a total of 160 DFT data points for H adsorption on the abovementioned six selected high-symmetry sites (i.e., T1, T1 , B1, S2, B2, H1 in figure 1) by means of modified Morse-like functions not shown here for purposes of brevity [29]. This kind of function is appropriate since the adsorption over all these sites was not activated. During the fitting procedure, much more attention was paid to satisfactorily describing the energetics of minimum adsorption sites. The fitting reports the following adsorption energies: 2.06, 0.75, 0.14, 0.16, 0.004 and 0.006 eV for T1 , T1, S2, B1, H1 and B2 sites. A good agreement with DFT data was achieved in the final fit, mainly for the most important site (T1 ). The global fit was good enough to reproduce the stability of the T1 site (where the hydrogen adatom is preadsorbed during the dynamical treatment) whereas it also reproduced qualitatively the adsorption on the other possible sites around the surface.
The V 1D term corresponds to the fitting of DFT data for the H 2,gas curve by means of a Morse-like potential, taking as zero of energies the asymptotic value: where r is the internuclear distance between the two hydrogen atoms. The final fitted parameters reported in [29] lead to dissociation energy and equilibrium internuclear distance values that compare quite well with the experimental ones [38]. For the V 6D term 2650 DFT data points were fitted for different approaches of the second H (2) atom keeping the H (1) adatom fixed on the T1 site, the most stable adsorption site. Thus, V 6D term was expressed as: where the analytical expressions for the interactions were also modified Morse-like functions and the switching functions (i.e., f 1 ,f 1 and f 1 ) are responsible for delimiting the strong interaction area [18].
The energy values associated with atomic adsorption (equation (1)) and E-R (equation (2)) processes (i.e., E (1) and E (2) , respectively) derived from the global PES (equation (3)) are −2.4 eV and −2.1 eV, respectively. Figure 2(a) shows the reference frame used during the study of H interaction with an H-preadsorbed surface (see next paragraph).
The internal coordinates used are the internuclear H (1) -H (2) distance (r), and the two angles (θ and φ) that define the orientation of the approaching. The angle that forms the internuclear distance and the negative z axis corresponds to θ and the azimuthal angle φ corresponds to the angle formed by the impact parameter, that is the orientation of projected internuclear distance onto the xy plane, with respect x axis. The atomic incoming velocity angle (θ v ) is defined with respect to the negative z axis (e.g., 0 • for normal incidence and 45 • for off-normal incidence) and its projection onto the xy plane (i.e., the azimuthal φ v angle) was uniformly sampled within the 0-360 • interval.
With respect to the E-R channel, and as discussed in our previous paper [29] focused on the construction of the analytical PES, the barrier for this reaction depends strongly on the direction of the incoming H (2) atom. The resulting analytical surface (equation (3)) shows that the most favorable incoming angle corresponds to φ = 55 • (figure 2(b)) and for an impact parameter value of 0.5 Å the barrier is very small (i.e., 0.02 eV with respect reactants), but when the φ angle changes the barrier becomes higher. Thus, for φ = 45 • in the B1 direction, the barrier respect reactants increases to a value of 0.16 eV. These values are in agreement with DFT-based results reported in [29] that predict a barrier height close to 0.08 eV for the preferred approaching angle of φ = 55 • .

The semiclassical dynamics method
The method assumed to describe the reaction dynamics of the two hydrogen atoms interacting with the β-cristobalite surface is the semiclassical method described in great detail in G Billing's book [28]. The method has been used to describe the dynamics of various atom recombination and dissociative chemisorption processes catalyzed by surfaces (see for example [16,17]) and takes into account several chemicalphysics effects which are typical of the heterogeneous character of the system under study, e.g., the crystallographic structure of the silica surface, the silica phonon structure and dynamics and the interaction between the silica phonons and the H atoms hitting the surface.
According to the semiclassical approach [28], the dynamics of the both gas phase and adsorbed hydrogen atoms is described classically, whereas the lattice atom vibrations are quantized.
The equations of motion of the two hydrogen atoms are solved using an effective potential defined by: where R is the two-body distance vector between the interacting H atoms and the lattice atoms. V int (R) is given by equation (3) and ph is the wave function for the quantum phonons. Thus, the motion of the atoms and molecules in the gas phase or at the surface is obtained by numerical solution of the Hamilton's equations of motion using the effective Hamiltonian H eff (t, T s ) given by, where the first term is the kinetic energy of the gas phase atoms and the last term is the energy exchanged between the hydrogen atoms and the silica substrate [28]. This method is expected to be valid within the limits imposed by a classical treatment of the translational and rotovibrational motion of the gas phase particles interacting with the surface. Those limits are mainly related to quantum tunneling effects and to the non-conservation of the zero-point energy.
Because of the large energy delivered in the recombination process and with such an insignificant barrier, the semiclassical approach is expected to be in substantial agreement with quantum results as well as quasiclassical ones. This fact has been confirmed previously, in the study of the isotope effect on H recombination on a rigid Cu (1 1 1) surface and also on graphite [39][40][41].
The potential energy surface for H 2 formation through H recombination on silica displays a very low energetic barrier to reaction (see the previous section). This point makes a classical dynamics treatment feasible even for the lightest nuclei in a wide collisional regime (from sub-thermal energies up to several eV).

Molecular dynamic calculations
Reaction (2) has been simulated by considering collisions of H gas with an H-preadsorbed β-cristobalite surface. The adatom is located on the T1 site at a variable distance around the chemisorption equilibrium site and it is assumed to be in thermal equilibrium with the silica surface so that the normal component of initial momentum is set equal to k B · T s , (where k B is the Boltzmann constant). The azimuthal angle (φ v ) has been sampled by randomly selecting the initial (x, y) position coordinates of the incoming H gas atom within the (1 × 1) unit cell, while its collision energy (E kin ) is varied in the range (0.04-3.0) eV. The effect of the incidence angle of the atomic flux hitting the surface and of the surface temperature on the reaction dynamics has also been explored. Then, to investigate the incidence angle effect two sets of initial conditions are assumed for H gas : in one case the gas phase H atom strikes the silica surface in an off-normal direction (θ v = 45 • ), while in the second simulation H gas approaches the surface in the normal direction (θ v = 0 • ). To explore the surface temperature effect two surface temperatures have been selected, T s = 1000 and 300 K. For each initial condition and collision energy at least 15 000 trajectories (collisional events) have been propagated. For the complete set of values of E kin the number of trajectories run was enough to assure a convergence of the calculated collisional observables (i.e., probabilities of the various surface processes, vibrational distributions in the final molecular states, etc.). As mentioned in the Introduction the simulated conditions are only an approximation of those more complex ones that exist in plasma conditions. Just to mention a few differences: here a perfect surface is considered while real surfaces can have defects and chemical contaminations; besides, since the method is single-collision, the initial coverage (1.97 × 10 14 atoms cm −2 , corresponding to one adsorbed atom in the crystallographic unit cell) is preserved in the trajectory evolution, while surface coverage is continually evolving if a flux is supposed.
The final state analysis of the trajectories shows that the following surface processes are energetically open to the collision dynamics: The classification of a trajectory is done according to both geometrical and/or energetical requirements. Thus, a H 2,gas molecule is considered as formed when the interatomic H-H distance r is less than 2.5 Å and the distance from the surface to the centre of mass of the molecule is larger than 7 Å. The rotational and vibrational states of hydrogen molecules, assuming the Morse oscillator given in equation (5), are analyzed in terms of the action-angle variables using the semiclassical quantization rules [42]. A hydrogen atom is considered adsorbed on the surface when the available energy to escape from the chemisorption potential well ( E esc ) is lower than the effective potential (V eff ) defined in equation (7). In particular, the E esc term is defined as E tot − E ph − E 0 where E tot is the total energy for a specific collision event, necessarily conserved in the course of the trajectory, E 0 is the zero point energy of the chemisorbed atom and E ph was already defined as the energy exchanged between both H atoms and the surface. Finally, the atoms are considered in gas phase when their distances from the surface are larger than 7 Å and the direction of their velocity vectors points to the vacuum. According to these criteria, the probability associated with all the processes detailed in equations (9)(a-f) (i.e., E-R recombination, atomic reflection, atomic adsorption + desorption, atomic adsorption, molecular formation + adsorption, and atomic desorption, respectively) is determined as a function of E kin of the incoming H gas atom. For the case of E-R reaction (equation (9)(a)) the energetic of the newly formed H 2,gas molecules is also resolved by determining the average energy fraction transferred to their translational, vibrational and rotational degrees of freedom as well as the average energy fraction interchanged with the surface. The global recombination coefficient γ E-R (T s , T gas ) is obtained by averaging the atomic recombination probabilities over a Boltzmann kinetic energy distribution of the incident atomic hydrogen flux hitting the silica surface [17]: where T gas is the gas temperature and E kin is the collision energy of the H gas atom impinging the H-preadsorbed one on the β-cristobalite surface. Likewise, the state-to-state recombination coefficients γ (v) (v = 1-12) can be obtained by replacing P E-R (E kin , T s ) with P E-R (v, E kin , T s ), that is the probability that for a given E kin and T s , atoms recombine to form a molecule in the vibrational quantum number v, whatever its rotational state.

Results and discussion
In figure 3 probabilities for the complete set of elementary surface processes, described in equations (9)(a-f), for both off-normal and normal incidence of H gas and T s = 1000 K are reported as a function of collision energy. The probabilities for H 2,gas formation through E-R reaction (equation (9)(a) and figure 3(a)) present a similar behavior for the two approaching angles. Thus, they rise sharply up to a maximum value for E kin ∼ 0.4 eV and then they decrease until reaching a constant value for E kin higher than 1.5 eV. Nevertheless, the off-normal incidence favors the E-R reaction, whose maximum value of probability is a factor of two greater than that for normal incidence. It is noticeable that the principal process at off-normal incidence corresponds to the reflection of the incoming atom, equation (9)(b) and figure 3(b) for all E kin values explored, whereas in the case of normal incidence the adsorption of the second hydrogen atom, equation (9)(d) (figure 3(d)) is the main process within the range E kin = (0.4-0.8) eV. This latter process, instead, has a probability significantly lower in the case of off-normal incidence, mainly in the energy range (0.4-1.5) eV. Finally, no significant differences are observed with regard to the effect of the impact angle of incoming H gas , on the probability of the other surface processes shown in figures 3(c) and (f ), i.e., adsorption + desorption and desorption, respectively.
The effect of the incident angle is stronger when considering the energy partitioning in the E-R recombination process. In particular, figure 4 reports the average energy fractions distributed into the final states of the newly formed H 2,gas molecules for the case of off-normal and normal incidence. Although all the degrees of freedom of the newly formed H 2,gas molecules are involved in the dynamics, in the case of off-normal incidence ( figure 4(a)) and E kin < 0.5 eV there is a clear close correlation between the vibrational and translational motions since a gain/loss of vibrational energy corresponds to an almost identical loss/gain of translational energy. The average energy transferred to the surface ( E ph ) is very low for both incidence angles and even for the higher E kin explored it remains below 4 and 8% for off-normal and normal incidence, respectively. The different interplay between the rotational and vibrational energy of formed H 2,gas (v, j ) molecules can be inferred comparing the roto-vibrational distributions for E kin = 0.06 eV for off-normal and normal incidence (figures 5(a) and (b), respectively).
The vibrational distributions peak at v = 5 and v = 6 for off-normal and normal incidence, respectively. The small difference obtained in vibrational excitation, for the two impinging geometries, can be ascribed to the higher rotational excitation obtained in the case of off-normal incidence at this E kin . In fact, in the case of off-normal incidence the average rotational quantum number is higher than that obtained in the case of normal incidence: j = 12 and j = 10 respectively. For off-normal incidence the corresponding roto-vibrational distributions appear broader and more densely populated in the medium-high levels.
2D contour plot for vibrational distributions, P (v, E kin ), for H 2,gas are reported in figures 6(a) and (b) for off-normal and normal incidence, respectively. It appears that, for both incidence angles considered in the calculations, and for low E kin , the vibrational distributions peak at medium vibrational levels. When collision energy is increased the distribution becomes wider and the peak moves towards lower vibrational levels, as can be seen by considering the average vibrational quantum number v . v has a step-like trend, being v = 6 for 0.04 < E kin < 0.1, v = 5 for 0.2 < E kin < 2.5 and v = 4 for E kin = 3.0. Looking to the vibrational distributions, it can also be inferred that recombination reaction of hydrogen atoms on silica surface produces H 2 molecules with non-thermal distributions.
These non-thermal distributions are in agreement with the non-early energy barrier character of the E-R channel ( figure 2(b)) that enables the activation/de-activation of the formed molecules.
In addition, the highest populated vibrational level v = 12 corresponds to an internal energy of E(v = 12, j = 0) =∼ 2.44 eV that is in complete agreement with the exothermicity of E-R reaction (2) in the PES (i.e., The non-thermal roto-vibrational distributions obtained in our study are in agreement with the results reported in [43] where high rotational levels up to j = 17 are observed for the vibrational state v = 2 and confirm the fact that the atom recombination is the main source of roto-vibrational excited molecules for hydrogen plasmas. Obviously, the rotovibrational energy content of formed H 2,gas (v, j ) molecules is strictly related to the specific reaction exothermicity; besides, the vibrational distributions detected in plasmas might not necessarily be the nascent distributions, as determined in the calculations, but also the results of multiple gas-phase and gas-surface collisions at the plasma-surface boundary layer. Production of highly vibrationally excited H 2 molecules has been also observed in experiments over metal surfaces (e.g., Cu and W) [5,8,44] and confirmed by dynamical calculations [45,46].
From the analysis of the trajectories it is also possible to ascertain the mechanism followed by E-R recombination as a function of the incoming collision energy range. The mechanism is determined by the adsorbed atom. Thus, at low collision energy (i.e., 0.06 eV), when the incoming atom approaches the surface (approximately 600 fs after the beginning of trajectory), H ad is already out of the minimum potential well geometry so that the molecule is formed far from the surface. At higher kinetic energies (i.e., 0.6 eV) when the incoming atom reaches the surface (approximately 200 fs after the beginning of trajectory), H ad is still trapped on the surface by the minimum potential well so that the recombination process takes place closer to the surface. This behavior can be inferred from figure 7 where the z coordinate for the two hydrogen atoms along with the interatomic distance r as a function of collision time for two typical trajectories at E kin = 0.06 eV and 0.6 eV are reported as an example of the general trends observed examining the complete set of reactive trajectories.
The reaction mechanism followed in the case of offnormal incidence is similar to that observed in the case of normal incidence.
The surface temperature effect has been also investigated only in the case of normal incidence. Figure 8 shows that at T s = 300 K the probabilities for an E-R reaction (equation (9)(a)), inelastic scattering of H gas (equation (9)(b)) as well as for the atomic adsorption + desorption process (equation (9)(c)) increase slightly with respect to the values at T s = 1000 K in the complete range of collisional energies explored. This increase is concomitant with a decrease of the probabilities for the atomic hydrogen and molecular hydrogen adsorption processes (equations (9)(d) and (e), respectively).
The surface temperature effect is difficult to predict due to the different factors that may contribute to it, but from an accurate analysis of the trajectories it can be pointed out that the different surface temperature effects observed on the processes described in equations ((9)(a,b,d and e)) can be ascribed to the greater mobility of H ad on the surface at T s = 1000 K. In fact, in the case of T s = 1000 K the preadsorbed atom moves over the surface changing the interaction site and obviously the final reaction channel. However, it should be noted that the observed surface temperature effect is due only to phonon dynamics. It is also interesting to point out that the same surface temperature effect on the recombination probability was found for H recombination on graphite (0 0 0 1) surface [41].
From the analysis of trajectories it emerges that the E-R recombination reaction at T s = 300 K follows a mechanism similar to that found at T s = 1000 K and already described. Regarding the average energy partitioning among the degrees of freedom of the formed H 2,gas molecules at T s = 300 K, a coupling between the vibrational and rotational energy fractions is observed mainly for E kin lower than 0.5 eV. Moreover, the average energy exchanged with the solid surface at T s = 300 K is lower than the value obtained at T s = 1000 K, as shown in figure 9. Once again, this effect is strictly related to the greater mobility of the preadsorbed atom as the surface temperature is increased.
An example of this occurrence is shown in figure 10 where it is reported the gas-solid interaction potential of the same reactive trajectory at two different surface temperatures as a function of the collision time. H 2,gas molecule is formed at t = 22.0 fs and t = 22.3 fs at T s = 300 K and T s = 1000 K, respectively. It is worth noting that the interaction potential experienced by the H atoms in the trajectory, before the reaction, is different and in particular at that time when the gas phase molecule is formed the interaction potential on the hottest surface (red line in figure 10) is higher. The different interaction potential is the result of different paths followed on the cold and hot surface by the preadsorbed H ad atom.
The average vibrational quantum number of the newly formed molecules for T s = 300 K is similar to that for T s = 1000 K, as a consequence, the determined vibrational distributions have, on the whole, the same behavior for both surface temperatures.
The different surface processes can be subdivided into those that deplete (equations (9)(a) and (f)) and those that populate (equations (9)(d) and (e)) the surface of H atoms. Processes of equations (9)(b) and (c) can be substantially neglected leaving the number of both adsorbed and gas-phase atoms unchanged. This analysis is made in figure 11 where, for all the simulated conditions, the probabilities of those processes that deplete and populate the surface are reported. Figure 11 shows that for the case of off-normal incidence (θ v = 45 • ) at T s = 1000 K and for normal incidence at T s = 300 K (except for E kin = 0.5 eV for which there is a small increase in the population probability) there is a high prevalence of depopulation processes. Instead in the case of normal incidence and T s = 1000 K, the hydrogen atoms adsorbed on the surface increase for kinetic energy values in the range (0.3-0.6) eV, outside of this range the amount of hydrogen atoms adsorbed diminish significantly.
This circumstance can have an impact on the material technology. In fact, more than the recombination processes, the prevalence of H atoms adsorption on the silica surface  for well-defined collision conditions can change the chemical composition of the gas-phase near the solid substrate and can also be responsible for possible thermal damage to the surface itself because of the large amount of energy transferred to the surface as heat flux. The global recombination coefficient γ E-R calculated from equation (10), using the probability obtained for normal incidence and for T s = 300 K and T s = 1000 K is reported in figure 12 as a function of gas temperature (T gas ) within the range 300-3000 K. It is observed that γ E-R increases as T gas does but decreases when T s becomes higher. However, if thermal conditions are imposed (i.e., T = T s = T gas , hereafter T) in the γ E-R calculation, the values obtained are 0.063 and 0.048 at T = 300 and 1000 K (highlighted with an orange box in figure 12), respectively. This trend is apparently opposite to the experimental behavior observed for silica between T = 300-900 K [47], for fused quartz between 194 and 1250 K (γ = 1.1 × 10 −3 for T = 1000 K) [14], for Pyrex within the range 275-450 K (γ = 1.2×10 −3 for T = 300 K) [14] and for fused silica at T = 300 K γ = 2.3×10 −3 (near post-discharge) to γ = 2.1 × 10 −4 (late post-discharge) [27]. Moreover, the calculated γ E-R values are one order of magnitude greater than the experimental γ values but an important point has to be taken into account: the experimental samples of silica did not present any preadsorption of H atoms. Thus, for an accurate comparison with experiments, it is important to note that the experimental γ coefficient is the result of at least a two-step process: atomic adsorption (equation (1)) and atomic E-R recombination (equation (2)), so the global probability of the recombination process should be given as a product of the two corresponding probabilities (i.e., P = P ad · P E-R ). This effect was already discussed in the recombination study of O/O 2 over β-cristobalite through a microkinetic model [26] leading to a good comparison with experiments. Since the calculations were done for an H-preadsorbed surface, adsorption probability (P ad ) can be approximated as the probability value obtained from equation (9)(d) (i.e., the probability that the incoming atom gets adsorbed when there is another one already adsorbed). Obviously, the correct value of P ad should be an average value obtained for the adsorption over a clean surface and also over a surface with several H atoms preadsorbed. Thus, the value considered here should be enough as a first approximation, taking into account that in the experiments no data about coverage are reported. Thus, the global recombination coefficient can be recalculated by replacing in equation (10) P E-R (E kin , T s ) by P (E kin , T s ) = P E-R (E kin , T s ) · P ad (E kin , T s ). By doing this, two approximate values of γ = 3.720×10 −3 and 5.185×10 −3 for T = 300 and 1000 K, respectively are obtained. These two values better reproduce the T-dependence and also the order of magnitude of the experimental γ values [14,27,47].  (9)(a-f)) for normal incidence and T s = 300 K (color lines and solid symbols). For comparison the probability for the same process and normal incidence but for T s = 1000 K are reported (black lines and open symbols). Note the different Y-scale in panel (b) for the case of H gas reflection (equation (9)(b)). The slight differences between the results reported here and the experimental ones [14,27,47] could be ascribed to the fact that in the present study the adsorbed atom is placed in a well-defined site on the surface, while in the experiment the atom can also be adsorbed in other surface sites, so that the measured recombination coefficient is an average value over all active sites on the surface. Another point concerns the silica structure used. Thus, in the present work a silica β-cristobalite surface is used while in [14,27,47] quartz, Pyrex or fused-silica surfaces were studied. Indeed, different silica polymorphs exhibit a different behavior in the same conditions, both for recombination and for adsorption processes [14,48,49]. Furthermore, well-defined values  of both gas and surface temperatures, together with fixed impinging angles, were assumed in the present study while in the experimental setup they were not.
As aforementioned, the recombination coefficients are very useful in kinetic modeling and CFD calculations for hydrogen plasmas kinetics and even more the state-to-state recombination coefficients γ (v) reported in table 1 for T s = 1000 K, with the greatest values within the range 5 < v < 7.

Concluding remarks
The recombination of hydrogen atoms on a β-cristobalite surface via an E-R mechanism was studied together with the other competitive elementary processes. In order to determine the influence of the incident angle of the atomic flux hitting the surface on the dynamics, two sets of initial conditions have been assumed: normal (θ v = 0 • ) and off-normal (θ v = 45 • ) incidence. At high surface temperature (T s = 1000 K) the main process at all collision energies (E kin ) explored for offnormal and normal incidence corresponds to the reflection of the incoming atom. Regarding the adsorption of the second hydrogen atom (as separated atoms or as a newly formed and adsorbed molecule), it becomes the main process for θ v = 0 • at E kin ∼ 0.5 eV. E-R reaction is favored for off-normal incidence whereas the other processes benefit from normal incidence. The shape of the E-R probabilities are similar for the two incidence angles showing a maximum value at E kin = 0.4 eV, which is a factor of two greater in the case of off-normal incidence. The atomic recombination dynamics has been investigated in detail.
Thus, for the H 2,gas formed molecules the energetics, as well as their roto-vibrational distributions have been determined, showing a different coupling between the degrees of freedom for the two different incident conditions. Particularly, these molecules are formed in the medium-high lying vibrational levels and rotationally excited. The energy transferred to the surface is very small, below 4 and 8% for off-normal and normal incidence, respectively. The E-R mechanism is similar in the case of both incidences and it corresponds to a direct mechanism.
The T s effect was also investigated at normal incidence. At T s = 300 K, the probability of E-R reaction and the inelastic scattering of H gas increase slightly in the complete range of E kin considered. In addition, a decrease of the probabilities for the atomic hydrogen and molecular hydrogen adsorption processes is also observed. The different T s effect on the processes is attributed to the greater mobility of H ad on the surface at T s = 1000 K. For E-R recombination, the mechanism is similar at the two temperatures. The T s effect observed is the same reported for H recombination on graphite (0 0 0 1) [41]. At T s = 300 K the energy exchanged with the surface is also lower than that at T s = 1000 K, while the average vibrational quantum number of the newly formed molecules is similar for both T s .
Two atomic recombination coefficients were calculated (γ E-R and γ ). In the first one only the E-R probability (P E-R ) was taken into account, while in the second one the probability of the preadsorption of hydrogen atom was also considered (i.e., P = P ad · P E-R ), since it is a necessary step to produce the E-R reaction. A more accurate comparison with experiments was obtained when the preadsorption process of H atoms was included. The calculated γ values reproduce the T dependence and also the order of magnitude of the experimental values. Nevertheless, slight differences may be attributed to the polycrystalline nature of the surfaces used in experiments and also to the fact that trajectories were run at precise values of gas and surface temperatures and for selected incidence angles.
The obtained results confirm H atom recombination at β-cristobalite as a source of roto-vibrationally excited H 2 molecules in low-temperature hydrogen plasmas also in the case of silica substrate as previously observed for metal surfaces.