Molecular dynamics study of hydrogen atom recombination over silica, based on a new analytical DFT potential energy surface

A new analytical potential energy surface (PES) based on new density functional theory data is constructed for the interaction of atomic hydrogen with both a clean and an H-preadsorbed β-cristobalite (001) surface. For the atomic interaction, six adsorption sites have been considered, the Si site (T1’) being the most stable one. The PES was developed as a sum of pairwise atom-atom interactions between the gas-phase hydrogen atoms and the Si and O atoms of the βcristobalite surface. A preliminary molecular dynamics semiclassical study of the different heterogeneous processes (e.g., H2 formation via Eley-Rideal reaction, H adsorption) that occur when H collides with an H-preadsorbed β-cristobalite (001) surface was carried out. The calculations were performed for collisional energy in the range (0.06 ≤ Ekin ≤ 3.0 eV), normal incidence and a surface temperature Tsurf = 1000 K. The recombination probability reaches its maximum value of approximately 0.1 for collisional energies in the range 0.3 ≤ Ekin ≤ 0.8 eV. The H2 molecules are formed in medium-lying vibrational levels, while the energy exchanged with the surface in the recombination process is very low.


INTRODUCTION
One of the major technological challenges associated with the access to planetary surfaces is the entry of a space vehicle into these planetary atmospheres at hypersonic speeds [1]. In the case of the Jovian atmosphere, it is necessary to know the interactions between atomic and molecular hydrogen with some materials that could act as a thermal protection system (TPS) of some parts of future vehicles. Currently, silica-based materials are predominantly used in some TPS surfaces, where catalytic atomic recombination (via Eley-Rideal (ER) or Langmuir-Hinshelwood (LH) mechanisms) as well as other important processes (atomic adsorption, reflection, etc.) play an important role in reentry heating at thermal non-equilibrium conditions. Among the different forms of silicon dioxide to be studied, either amorphous or crystalline, the β-cristobalite is the most stable crystalline polymorph at high temperatures up to the melting point of 1983 K, and can be a good material model for the high temperatures expected for future Jovian entries. Scarce theoretical or experimental data are available for the ER reaction with atomic hydrogen over silica materials. Thus, a recent Density Functional Theoretical study (DFT) on H 2 + Si 3 O n clusters reports some information of molecular adsorption energies (E ad = 1.09 -2.29 eV) and about the preference for H 2 dissociation over a single Si surface atom [2]. The atomic recombination coefficient (γ) for H over fused quartz and Pyrex has been formerly measured within the interval 194 -1250 K, showing a small activation energy (E a = 0.42 ± 0.06 eV) at high temperatures but much lower for low temperatures, that would correspond mainly to the Eley-Rideal reaction [3]. In the present study, we have used spin-polarized DFT calculations to investigate the atomic and molecular hydrogen interaction with a clean β-cristobalite surface as well as its atomic ER reaction, which produces H 2 . We have followed a similar theoretical approach as for the previous study on the O/O 2 + β-cristobalite [4,5] system, with major importance for Earth's reentries. From the obtained ab initio results we determined the PES that has been used to perform a preliminary dynamics study on the recombination of H atoms via an ER mechanism. The calculations have been made using the semiclassical collisional method [6], which has been used in previous studies on various atom/surface reactive systems [7,8]. Thus, the main goal of this study is to provide dynamical data for heterogeneous processes involving H/H 2 and β-cristobalite, which can be useful for later computational fluid dynamics (CFD) simulations.

DFT STUDY
We have performed periodic DFT calculations by means of the VASP code [9,10,11], which uses a plane wave basis set. The calculations are based on the generalized gradient correction (GGA) functional Perdew-Burke-Ernzerhof (PBE [12]). The electron-ion interactions were described by using the projector-augmented-wave (PAW) [13] technique. An energy cutoff of 400 eV was high enough to obtain converged properties. Spin-polarized calculations were carried out to check the ground electronic state of the total system, as the atomic adsorbate is an open-shell system (i.e., H( 2 S)). Integration over the Brillouin zone was done by using a 3×3×1 (in slab calculations) or 3×3×3 (in bulk calculations) k-point mesh by means of the Monkhorst-Pack method [14]. Geometrical optimizations and vibrational frequencies were computed with an energy accuracy of 10 -6 eV. For adsorption studies we use a (1×1) unit cell and a slab model formed by 6 layers (Si-O-Si-O-Si-O) of SiO 2 with an outermost silicon layer. An additional hydrogen back layer was added for saturating the dangling bonds of terminal oxygen atoms. A vacuum of 11.6 Å was enough to avoid significant interactions between slabs in the z direction.
The 6-layer slab model described above and used throughout all study is showed in figure 1, where the six sites considered are labeled as T1', T1, S2, B1, H1 and B2, following the notation used in our previous works [4,5]. The analysis of the stationary points over the surface for a single hydrogen atom adsorption (coverage, θ = 0.25) have been obtained for a rigid slab model as well as for a four-layer relaxed slab, although with similar results. The geometrical and energetic data are listed in Table 1. The most stable adsorption site corresponds to the tilted T1' site (a true minimum) located over the positive diagonal of the cell, just in the opposite direction of the oxygen atoms belonging to the second layer (figure 1). The other sites explored correspond to transition states, all of them related with the diffusion of the hydrogen atom over the surface because the imaginary frequencies correspond to parallel movements over the surface. From the adsorption energy values, it is worth noting that the H atom needs to overcome at least 1 eV of potential energy for discussion to occur. This value is practically a half of the one obtained for the diffusion of atomic oxygen over βcristobalite [4], being the diffusion of hydrogen more important that for the case of oxygen. We have also studied the coverage effect ( Table 2). As mentioned above, the most stable adsorption site for a single H atom was on the Si atom (T1' site, E ad = 2.10 eV/atom) whereas the adsorption of a second H atom depends on the new T1' position with respect to the preadsorbed one. Thus, a second adsorption over the same Si atom (geminal adsorption) generates a more stable configuration (E ad = 3.02 eV/atom) than the case of two H atoms adsorbed over different Si atoms (vicinal adsorption) for both rigid and 4-layer relaxed slab model. Therefore, a stronger adsorption is obtained for higher coverage where all SiH bonds are geminal. This extra-stabilization obtained for geminal adsorptions is due to the saturation of the two dangling bonds of each silicon atom. Once the main geometrical details are considered, we focused our attention on calculating the approaching of a single H atom over each of the sites showed in figure 1 in order to obtain a correct description of the adsorption process throughout the surface. The DFT data obtained for each site are showed in figure 2. Afterwards, all these data were fitted over each site with the analytical expression described in the next section. Figure 2 shows that all the adsorption sites studied are not activated. Moreover, we paid more interest in describing satisfactorily the minimum adsorption sites and their neighborhood (T1' and T1 sites). A reasonable agreement is achieved but DFT data allows some penetration mainly in T1' and in greater extent in B1 and H1 sites. For the case of T1' the fitted data do not describe this penetration but in the case of B1 and H1 sites it is allowed in some extent. The effect of trying to obtain the best fit for all the sites considered makes that the width of T1' DFT curve was reduced in the final fitted analytical PES. For the DFT study of the ER reaction (ΔE = -2.37 eV), we have located one hydrogen atom (H (1) ) preadsorbed over the most stable minimum site (T1') whereas the incoming one (H (2) ) samples all the φ angles at a particular impact parameter b for different z distances to the surface of the two hydrogen atoms (figure 3). The minimum of the potential energy at φ = 55° in figure 3b indicates that the incoming H (2) atom prefers the approaching in the direction of B1, that corresponds to an open area of β-cristobalite. According to this, we have fitted the main characteristics of the DFT data set at this particular φ angle. The DFT ER barrier is small with a value of ΔE ≠ = 0.22 and 0.08 eV for rigid and 4-layers relaxed slab models, respectively. The LH reaction is endothermic (∆E = 1.63 eV) with a large energy barrier of ΔE ≠ = 4.65 eV for rigid slab model with two geminal adatoms. This fact makes LH reaction clearly less favorable than the ER reaction and it is the reason for centering our attention in the ER recombination process. (b) φ angle dependence for H (2) approach to H (1) for b = 1 Å and z (2) = z (1) = 0.74 Å.
The details of transition states for ER and LH reactions are listed in Table 3. The main result comes from the fact of obtaining two imaginary frequencies but it is a logical behavior since the two reactions are closely related. Thus, for the ER transition state two possible reactions can occur at this point, the first one if the H-Si distance becomes greater giving to the H 2 gas phase (namely the ER reaction, achieved by means of the 243.4i (z) imaginary frequency) and the second one if the HHSi bonds bend allowing the outermost H atom to become adsorbed in the vicinal adjacent T1' site (through the 168.3i (x,y) parallel normal mode). For the case of LH transition state the two imaginary normal modes are also related with the two possible channels. The highest one (2612.9i (x,y) parallel mode) drives to the dissociative adsorption of H 2 molecule on the two geminal T1' adsorption sites. The other one (340.4i (z)) is related with the detachment of the H 2 molecule that finally becomes separated from the surface (namely the LH reaction).

ANALYTICAL POTENTIAL ENERGY SURFACE
In consonance with the semiclassical dynamics method, the interaction potential is given as sum of pairwise atomatom interactions between the gas-phase hydrogen atoms (H (1) and H (2) ) and the Si and O atoms of the β-cristobalite surface. The global expression for the PES potential (1) is formed by the atomic-surface terms (V 3D ), one for every hydrogen atom, the gas phase hydrogen atoms interaction term (V 1D ) and the term that accounts for the two hydrogen atoms interacting with the surface (V 6D ). (1) ,y H (1) ,z H (1) ,x H ( 2 ) ,y H ( 2 ) ,z H ( 2 ) ( )

V x H (1) ,y H (1) ,z H (1) ,x H ( 2 ) ,y H ( 2 ) ,z H ( 2 )
The β-cristobalite (2×2) crystal model was built up, consisting of 149 atoms disposed on ten layers according to the atom placement in the Fd3m crystallographic unit cell. The phonon dynamics of the surface atoms was then developed by numerical diagonalization of the 3D-dynamical matrix of the force constants obtained assuming the BMH (Born-Mayer-Huggins) as interaction potential of the Si and O lattice atoms [15]. The phonons eigenvalues and eigenvectors that enter the scattering equations of motion are then obtained. The calculated phonon distribution function was previously reported [17], and it was in qualitative agreement with the experimental phonon spectrum [18].
The V 3D term is split into the gas-solid pairwise atomic interactions, ) are the i-th hydrogen atom Cartesian coordinates. Moreover, R ij stands for the distances between the gas phase i-th hydrogen atom and the j-th Si or O lattice atoms. This term fitted 160 DFT points ( figure  2) for H adsorption on the six selected high-symmetry sites (T1, T1', B1, S2, B2, H1 drawn in figure 1) by means of modified Morse-like functions (3) and (4), with, α 3 = c 9 ⋅ R ij 2 + c 10 ⋅ R ij + c 11 , α 4 = c 12 ⋅ R ij 2 + c 13 ⋅ R ij + c 14 (6) T ij x H ( i ) ,y H ( i ) The optimal parameters are listed in Table 4 and they are able to describe accurately the DFT data. The DFT calculations for H 2(g) curve were fitted by means of a Morse like potential, taking as zero of energies the asymptotic value (2H (g) in their lowest doublet states), where r corresponds to the internuclear distance between the two hydrogen atoms. The final fitted parameters are 4.4700 eV, 2.1787 Å -1 and 0.7500 Å for D H 2 , α and r eq , respectively. The dissociation energy and equilibrium internuclear distance compare quite well with the experimental values ( D H 2 = 4.4780 eV and r eq = 0.7414 Å [19]).
For the V 6D term we have fitted a total of 2650 DFT points for different approaches of the second H (2) atom, keeping fixed the X and Y coordinates of H (1) adatom over T1' site. Thus, V 6D term was expressed as, where the analytical expressions for the H (i) -Si/O (j) interactions are, ( 2 ) ,y H ( 2 ) ,z H ( 2 ) ( (12) and the switching functions ( f i and f i for i = 1, 2 and f 1 ), which delimits the strong interaction area are, b being the impact parameter between the two hydrogen atoms. The optimal parameters are given in Table 4.

SEMICLASSICAL DYNAMICS RESULTS
We simulate the rate-determining step of the ER recombination mechanism: H (g) + H (ad) / SiO 2(s) → H 2(g) (v, j) + SiO 2(s) taking the adsorbed hydrogen atom (H (ad) ) on the T1' site at a distance variable around the chemisorption distance and in thermal equilibrium with the surface at T surf = 1000 K, while the gas-phase hydrogen (H (g) ) atom strikes the silica surface in the normal direction with respect to the surface plane. The kinetic energy (E kin ) of the impinging atom is in the range 0.06 -3 eV. For each collisional energy we propagate a number of trajectories enough to assure a convergence for the calculated probabilities. In Figure 5a we report the ER recombination probability (i.e., P ER =N ER /N TOTAL ) as a function of collision energy. This rises sharply to its maximum value of nearly 0.1 for collisional energies around E kin = 0.4 eV, then decreases for collisional energy higher than 0.6 eV, reaching an asymptotic and constant value for energies higher than 1.5 eV.
We can calculate the γ(T surf , T gas ) recombination coefficient (figure 5b) by averaging P ER (T surf , E kin ) over a Boltzmann distribution of the initial kinetic energies of the atomic hydrogen flux at a given T gas [20]. FIGURE 5. ER recombination probability for H 2(g) formation as a function of the atomic hydrogen E kin (a) or its gas temperature T gas (b) over a silica (β-cristobalite) surface, with T surf = 1000 K.
When we compare the results obtained with the experimental ones for T surf = T gas =1000 K [3], we found that the recombination coefficient calculated by us is nearly two order of magnitude larger with respect to the extrapolated values for high temperatures [3]. This discrepancy can be ascribed to the fact that in our case the adsorbed atom is placed in a well-defined site on the surface, while in the experiment the atom can be adsorbed also in other surface sites, so that the measured recombination coefficient is an average over all active sites on the surface. Another point is that in our case we consider a silica β-cristobalite surface, while in Ref. [3] it is considered a generic silica surface. Indeed, different silica polymorphs exhibit a different behaviour under the same conditions, both for recombination and for adsorption processes [20,21].
The recombination reaction is exothermic so that it is very important to evaluate the part of this exothermicity that goes to the surface, as phonons excitation/de-excitation, and to the molecular internal states of the formed molecules, as vibrational (E vib ), rotational (E rot ) and translational (E tr ) motions. The final energy fractions for molecule and solid are reported as a function of E kin in figure 6a. We can observe that for lower collisional energies molecules are formed vibrationally excited. For E kin higher than 0.5 eV, the molecules are more rotationally excited while the rotational and vibrational reductions correspond to an increment of the translational motion, which for the highest collision energies is dominant. The energy transferred to the surface by the ER process is very low and also for the higher collisional energies remains below 8%. Therefore the main contribution to the energy transfer toward the surface arises from the atomic adsorption and the H 2(ad) formation, mainly from this latter process having a probability twice that for H 2(g) production.
Looking at the final H 2(g) vibrational distributions (figure 6b) we can infer that they are inverted and peaked around v = 6 for the lowest collisional energies, while the maximum shifts towards lower vibrational levels when increasing the energy of the impinging H atom.

CONCLUDING REMARKS
A DFT study on the interaction of atomic H over a β-cristobalite (001) surface reveals the existence of a strong chemisorption mainly over the Si atoms without an energy barrier. The adsorption of two H atoms over the same Si atom (geminal adsorption) generates a more stable configuration than the case of two H atoms adsorbed over different Si atoms (vicinal adsorption).
A new analytical potential energy surface was constructed based on present DFT data, which describes the interaction of atomic hydrogen with an H-preadsorbed β-cristobalite (001) surface.
A semiclassical dynamics study of atomic H colliding with a preadsorbed β-cristobalite (001) surface has been performed with this PES. Different competitive processes were observed (e.g., H adsorption, H 2(ad) formation) although the attention was focused here on the H atomic recombination to form H 2 gas molecules (i.e., ER reaction). Its probability was found to be small (< 10% for 0.3 ≤ E kin ≤ 0.8 eV and T surf = 1000 K) and the new H 2 molecules become vibrationally excited, while the energy released to the solid surface is very low.