Ab initio molecular dynamics study of the Eley-Rideal reaction of H + Cl-Au(111) → HCl + Au(111): Impact of energy dissipation to surface phonons and electron-hole pairs.

The reaction between an impinging H atom and a Cl atom adsorbed on Au(111), which is a prototype for the Eley-Rideal mechanism, is investigated using ab initio molecular dynamics at different incidence angles. The reaction yielding gaseous HCl with large internal excitation proceeds via both direct and hot-atom mechanisms. Significant energy exchange with both surface phonons and electron-hole pairs has been observed. However, their impact on the reactivity and final state distributions was found to be limited, thanks to the large exothermicity and small barrier of the reaction.


I. INTRODUCTION
Fundamental molecule-surface encounters are the basis for understanding many interfacial problems such as aerosol chemistry, corrosion, materials processing, heterogeneous catalysis, and plasma chemistry. All chemical processes at gassolid interfaces are initiated by collisions of gas-phase species with solid surfaces, but the subsequent bond formation and rupture may proceed through different mechanisms. An idealized limit is that of Langmuir and Hinshelwood (LH), in which reactions take place between adsorbed reactants that are fully thermalized on the surface. 1 The other limiting case is the Eley-Rideal (ER) mechanism, 2 which is nowadays ascribed to reactions between gas phase species and adsorbed ones through direct collision. As the impinging reactant is not equilibrated with the surface, an ER reaction is decidedly nonthermal. Between these two limits, an intermediate case is the hot-atom (HA) mechanism, in which the impinging atom (or molecule) does not collide with the adsorbate directly as in the ER scenario but remains trapped on the surface for some time before finding its reaction partner. 3 The HA mechanism differs from LH in that the hot atom is not completely thermalized on the surface, but it is also difficult to distinguish this mechanism from ER experimentally. While an LH reaction can be characterized with canonical transition-state theory, the fast and non-thermal ER and HA reactions require a dynamical treatment.
Recently, there has been increasing interest in understanding the dynamics of ER and HA reactions. High vacuum studies have been reported for several systems a) Authors to whom correspondence should be addressed: bjiangch@ ustc.edu.cn and hguo@unm.edu including H/D + D*/H*, [4][5][6][7][8] H + Cl*, [9][10][11] and N/O + N*, 12,13 where * denotes an adsorbed species. Fast products with internal excitation are found with narrow angular distributions, due apparently to the direct nature of these exothermic reactions. These pioneering experiments have stimulated many theoretical studies.  While earlier models were reduceddimensional, more recent studies have included all molecular degrees of freedom. The potential energy surfaces (PESs) used in these simulations have also metamorphized from empirical ones 14,16,21 to ones parameterized by density functional theory (DFT) calculations. 18,20,22,28,31,34,35 The dynamics on these PESs have been simulated using both quantum mechanical (QM) 14,15,21 and quasi-classical trajectory (QCT) methods. 16,18,30,31 Since the ER/HA dynamics are often barrierless and exothermic, QCT yields good results as quantum effects are often small. Indeed, most recent dynamics calculations in high dimensions have been exclusively performed using QCT. These theoretical studies have succeeded in capturing the essential features of ER/HA reactions, such as product vibrational excitation and near specular scattering. Importantly, theoretical simulations have revealed that the HA mechanism typically plays a dominant role in such processes. [18][19][20]23,27,30,35 A less-investigated aspect in ER/HA dynamics is the involvement of the substrate surface and the extent of energy exchange. There are two possible energy dissipation pathways on metal surfaces. One is the commonly known surface phonons, namely, the lattice vibration, which are particularly important for molecules containing heavy atoms. The less investigated is the electron-hole pairs (EHPs), due to the scattering of metal electrons on the moving atom. Such a process is electronically nonadiabatic because of the involvement of electronic excitation in the metal. The J. Chem. Phys. 148, 014702 (2018) involvement of EHPs in surface dynamics has been shown to be omnipresent on metals. [40][41][42][43] Recent studies have, for example, demonstrated that impinging H atoms can quickly dissipate their energy through coupling to low-energy EHPs in the metal surface. [44][45][46] More recently, it was predicted that EHP-induced dissipation can greatly reduce the H + H* HA recombination reactivity on tungsten surfaces. 37,38 Thus, it is imperative to include both dissipation channels in dynamical studies, particularly, if they involve light atoms or adsorbates. 36,47 Many earlier theoretical models for ER/HA dynamics have relied on analytical PESs assuming a rigid surface. In more recent studies, energy exchange with surface phonons has been approximated either with simple surface oscillators 24,29 or a generalized Langevin oscillator model. 32,33,36 On the other hand, the effects of EHPs have been explored [36][37][38][39] using the local density friction approximation (LDFA), 48 which expresses the non-adiabatic interaction with the surface EHPs as a friction term in the dynamical equation. 49 For the current system, the importance of EHPs has been suggested previously, but only a semi-empirical correction was considered. 27 Here, we investigate the dynamics of a prototypical ER reaction, in which a gas phase H atom reacts with a Cl atom adsorbed on the Au(111) surface to form desorbed HCl, using ab initio molecular dynamics (AIMD). 50 The AIMD approach differs from previous theoretical studies on the recombination reactions 15,[25][26][27] in two aspects. First, it avoids the need to build a high-dimensional PES as the forces are calculated on the fly, and thus eliminating possible fitting errors. Second, the surface phonons are explicitly included in the simulation. The AIMD approach has very recently been used successfully in studying HCl scattering and dissociation on Au(111). 51,52 In addition to the AIMD studies, we have also carried out AIMD with electronic friction (AIMDEF) 47,53,54 calculations in which the possible involvement of EHPs in the energy dissipation of hot H atoms is explored with the LDFA approach. 48 The AIMDEF approach eliminates the uncertainties in the previous empirical treatment of the EHP effects. 27 Thus, the results presented here are expected to provide a much more quantitative and complete picture of the ER/HA dynamics of the title reaction and the role of energy dissipation.

II. COMPUTATIONAL DETAILS
AIMDEF 46,53 is an integration of AIMD 50 with LDFA. 48 The latter is an effective way to treat energy dissipation due to EHPs on metal surfaces, in which the nonadiabatic effect on the nuclei dynamics is approximated as a friction force. The atomic friction constants are approximately determined by the friction coefficients that those atoms would have in a homogeneous free-electron gas. 55,56 In case of a molecular projectile, the independent atom approximation (IAA) is often employed, in which the friction coefficient of each atom in the molecule is calculated independently to be that of this atom embedded in the homogeneous free-electron gas with the electron density equivalent to that on the bare surface. 48 Thus, this approximation neglects the contributions of other atoms in the adsorbate in the embedding electron density. In AIMD simulations where the surface atoms are allowed to move, the total electron density varies with the substrate structure, making it difficult to find the embedding electron density for each atom. To circumvent this shortcoming, the Hirshfeld partitioning scheme 57 is used to decompose the overall charge to each atom and then remove its contribution from the electron density of the entire system obtained at each time step, yielding the embedding bare-surface electronic density on the fly. 47,54 This AIMDEF treatment has been successfully applied to previous studies on the relaxation of hot atoms and molecules 47 and vibrational lifetimes of molecular adsorbates. 58 The AIMD and AIMDEF calculations were performed using a modified version of the Vienna Ab Initio Simulation Package (VASP), 59,60 which includes the original LDFA package of Refs. 47 and 54. Spin-polarized DFT with projector augmented wave (PAW) 61 potentials were used for electronion interactions, and the electron exchange correlation effects were treated within the generalized gradient approximation (GGA), using the revised Perdew-Burke-Ernzerhof (RPBE) functional. 62 The Cl-covered Au(111) surface was modeled with a 4-layer slab of a 2 × 2 surface unit cell (1/4-ML coverage) with a vacuum region of 16 Å between the neighboring slabs. The Monkhorst-Pack k-points grid mesh 63 is 4 × 4 × 1 and the planewave expansion is cut at a kinetic energy of 350 eV. The optimized lattice constant for bulk Au is 4.199 Å, which was multiplied by an expansion coefficient of 1.0032 (4.212 Å) in order to mimic the small thermal expansion observed in experiments. 52 To simulate a Cl-pre-covered thermalized surface, the Cl atom was first relaxed on one of the four fcc hollow sites of the simulated Au(111) surface together with the three top layers, and then the initial velocities for the Cl and surface atoms were taken from the Maxwell-Boltzmann distribution with the initial temperature set as twice of the target temperature, followed by an equilibration run in the microcanonical (NVE) ensemble for 3 ps with a 1 fs time step. An extra 1 ps propagation was performed to take 1000 subsequent configurations for serving as possible initial configurations of the thermalized slab in the corresponding AIMD and AIMDEF simulations.
The initial position of the H atom was placed 5.5 Å above the Au(111) surface and randomly distributed in the unit cell. Here, we consider an incidence kinetic energy of 0.3 eV with an incidence polar angle (θ i ) of 60.0 • from the surface normal, as used in the experiments. 10,11 For comparative purpose, we have also performed simulations for H impinging at normal incidence (θ i = 0 • ) with the same incidence energy. In both cases, the azimuthal angle (φ i ) is randomly sampled. The AIMD trajectories were propagated using the leap-frog algorithm implemented in VASP, with a maximum propagation time of 1 ps and time step of 0.5 fs. The total energy was well conserved within ∼60 meV for most AIMD trajectories. On the other hand, the AIMDEF trajectories were propagated using the Beeman propagator 53 with a time step of 0.5 fs. In the AIMDEF calculations, the surface temperature enters the generalized Langevin equation as a random force.
A trajectory is considered reactive if the adsorbed Cl is abstracted by the impinging H, resulting in HCl ejection from the surface. In practice, this is defined by the HCl center of mass at 5.5 Å above the surface with the corresponding velocity pointing away from the surface. An abstraction is considered as ER when the H atom and Cl atom directly recombined at 2.0 Å or more above the Au(111) surface. Otherwise, the abstraction is counted as HA. In both cases, the vibrational action of the HCl product was determined by the EinsteinBrillouinKeller (EBK) method, and the traditional histogram binning (HB) method was used to determine the vibrational and rotational quantum numbers. A one-dimensional potential curve for HCl was calculated ahead of time using the same DFT method and used in the EBK calculation. Two other possible outcomes include exchange (EX) if the Cl atom desorbs and the incident H atom remains bound to the surface and reflection (RF) if the incident H atom is scattered back to vacuum. Finally, the event was considered as trapping if the H atom stays on the surface after 1 ps.

III. RESULTS AND DISCUSSION
The preferred binding site of both the H and Cl atoms was found to be on fcc hollow sites at 0. 8  The experimental values are generally larger than the calculated ones, presumably due to differences in coverage. 26 The climbing-image nudged elastic band (CI-NEB) method 67 was used to determine the transition state (TS) for the dissociative chemisorption (and the reverse recombinative desorption) of HCl. As shown in Fig. 1, the barrier height is 1.02 (0.34) eV for the HCl dissociation (reverse recombinative desorption) on Au(111). The dissociation barrier agrees well with the recent theoretical results of 1.08 68 and 1.05 eV 52 using the RPBE functional. The optimized TS geometry has bond length r H-Cl = 1.94 Å, polar angle θ = 44.7 • , and Z = 2.48 Å, which is also similar to those previous RPBE results. 52,68 On the other hand, the recombination barrier has not been reported in previous theoretical studies. 25,27 As shown in Fig. 1, the ER reaction (H + Cl* → HCl) is highly exothermic (∆E = 2.68 eV), although the exothermicity is smaller than that reported earlier (3.0 eV), 26,27 and barrierless as indicated by the dotted line in Fig. 1. In addition, the LH reaction (H* + Cl* → HCl) is still exothermic but with smaller exothermicity and a small barrier lower than the energy of the reactant asymptote.
Our AIMD and AIMDEF simulations show that the recombination reaction (ER + HA) is extremely efficient at this coverage, with large probabilities of about 77%-85%. The probabilities and corresponding cross sections for the different channels observed in our simulations at normal and off-normal incidences are listed in Table I together with the total number of trajectories used in each simulation type (last column). In all cases the HA mechanism dominates over the ER one, although the latter increases considerably at the off-normal incidence. We interpret such an increase in terms of a larger geometric cross section of the adsorbed Cl* atom for the H atom impinging with an oblique angle. In addition to recombination, there is significant trapping at normal incidence, but the percentage of this channel is much smaller at the off-normal incidence. It is expected that the trapped trajectories will eventually lead to LH reactions, but they are not followed in this work beyond the 1 ps time limit. Finally, note that there are only small fractions of trajectories undergoing reflection and exchange.
This overall picture from our AIMD simulations is consistent with previous theoretical results of Quattrucci and Jackson using an analytical PES. 27 Interestingly, the inclusion of EHPs has a relatively small impact on the final outcome of the H+Cl*/Au(111) interaction. This point will be further addressed later. Below, we will focus on the ER and HA dynamics, ignoring the other channels. It should be noted that due to the relatively small numbers of AIMD trajectories, the statistical error can be quite significant, resulting in sometimes rugged distributions. In Fig. 2, the vibrational and rotational state distributions of the HCl product obtained from both the AIMD and AIMDEF simulations are displayed for the two incidence angles. It is clear that the ER pathway produces high vibrational excitations with a relatively narrow range peaking at v = 7, due apparently to the direct nature of the collision. The HA channel, on the other hand, yields a broader distribution (0 ≤ v ≤ 9) with less vibrational excitation on average. The corresponding ER rotational excitations are however much lower and narrower than the HA ones. These observations, which are also consistent with the previous theoretical work of Quattrucci and Jackson, 27 can be interpreted as follows. In ER, the direct hit of Cl* by the impinging H atom imparts a large amount of energy into the H-Cl relative motion, thanks to the strong H-Cl attraction, which leads to a vibrational excitation of HCl. 27 Its rotational excitation is minor because of the small impact parameter in the direct ER mechanism. In contrast, the HA mechanism involves the initial trapping of H on the surface, which approaches Cl* with potentially large impact parameters. Such a scenario thus results in much higher and broad rotational state distributions because H* readily enters an orbital motion around Cl. 27 Again, the inclusion of EHPs does not change the overall conclusion regarding the rovibrational state distributions of the formed HCl. Unfortunately, a direct comparison with experiment is not possible because the measured rovibrational distributions of the HCl product are not complete. [9][10][11] The existence of measured velocity and angular distributions of the HCl product allows us to make comparisons with experimental results in Fig. 3. Since the experimental data were measured within the scattering plane, which according to our results constitutes ∼45% of the total reaction flux, the calculated results in Fig. 3 are thus obtained from outgoing HCl trajectories within the scattering plane defined by a small aperture of ±3 • in the azimuthal angle. In addition, note that the simulation temperature of T s = 300 K is higher than in experiment (T s = 100 K), but it is not expected that a lower surface temperature will significantly alter the qualitative outcome of our simulations of the ER dynamics, 11,32 although the temperature might significantly alter the percentage of the thermal channel. 11 The agreement for the velocity distributions is quite reasonable, which confirms indirectly the high internal excitation of the HCl product. The angular distributions of the outgoing HCl also reproduced the overall envelop of the experimental distributions, both with a slight bias near the specular angle, further confirming the direct nature of the reaction. The dip in the angular distributions near θ f = 0 • is due to the quasi-classical nature of the Cl* vibration prior to the reaction, as explained earlier by Quattrucci and Jackson. 27 While the overall results presented above are mostly consistent with both available experiments 9-11 and previous  theoretical studies, 27 the use of AIMD and AIMDEF allows us to provide more concrete details of the energy dissipation into the surface during the recombination process. The energyloss distributions obtained from these simulations are shown in Fig. 4. Here, the total energy exchanged with the surface is defined as the exothermicity (2.68 eV) plus the initial incident energy of the H atom (0.3 eV), plus the zero point energy of Cl* on the Au(111) surface (0.02 eV), and finally subtracted by the total energy of the HCl product (internal and translational energies). Thus, a positive (negative) energy exchange corresponds to the energy transfer from (to) the surface to (from) the HCl product. The mean energy loss (eV) in AIMD and AIMDEF calculations is listed in Table II. It is clear from the AIMD results that the HCl formation is accompanied by important energy exchange with the substrate phonons. The majority of the energy exchange results in energy loss to the surface phonons, but there are significant energy gains as well. The average energy loss is 0.13 and 0.02 eV for AIMD calculations at normal and off-normal incidences, respectively. In the same figure, the energy dissipation to surface phonons is decomposed for the ER and HA mechanisms. It is clear that the amount of energy exchanged with the surface in the ER channel is relatively small, and the corresponding distributions are roughly symmetric around zero (i.e., the energy gain and loss is quite similar). By contrast, there is much more energy loss in the HA channel, presumably due to the collisions that the hot H atoms experience with the surface atoms prior to recombining with the adsorbed Cl. This is consistent with earlier findings concerning the HA dynamics. 24 The AIMDEF results indicate additional energy loss due to surface EHPs. Interestingly, there is a significant difference in the energy dissipated into EHPs between the normal and off-normal incidences, as evidenced by the average energy loss of 0.44 and 0.09 eV, respectively. At the off-normal incidence, the energy loss to EHPs, which can roughly be estimated as the AIMDEF energy loss minus the AIMD energy loss, is relatively small, with an average ∆E of ∼0.07 eV. However, the dissipation to EHPs is much larger (∆E ∼ 0.31 eV) at normal incidence. As recognized in many previous studies, 46,54,69,70 the amount of energy dissipation to surface EHPs depends on the access to high electron density regions and on the time of such interaction. As shown in Fig. 5, the distribution of the residence time of H atoms on the surface, which is defined as the time spent by these atoms at vertical positions below 2 Å, is quite different in the two cases. The H atoms with a 60 • incidence angle have relatively short residence times. The short residence times (∼50 fs) deprive the H atoms the opportunity to efficiently interact with the surface electrons and create EHPs, leading to a relatively small energy loss in this channel. On the other hand, the normal incidence affords the H atoms large kinetic energies perpendicular to the surface, which allows these atoms to access regions with high electron densities. Indeed, a significant portion (∼15%) is found to penetrate the second layer of the metal surface. As a result, the averaged residence time of the H atoms (∼110 fs) is much longer, and the energy dissipation into EHPs is accordingly much larger. This observation is a reminiscence of the recent experimental observation of H-atom slowing down on a metal surface and its nonadiabatic origin suggested by theory. [44][45][46] It is also similar to the large EHP-induced dissipation on the H + H* recombination on tungsten surfaces. 37,38 In that case, however, the energy loss was predicted to be sufficient to reduce the reactivity substantially because the ER exothermicity of the H + H* → H 2 reaction is significantly smaller (1.5 eV) than in the current case (2.68 eV), and the corresponding LH reaction is strongly endothermic (∆E = 1.5 eV). Thus, slower H atoms may not be able to overcome the endothermicity and thus get trapped.
Interestingly, the larger energy dissipation to EHPs at normal incidence has little impact on the HA reactivity for the current system. This can be understood as the barrier for the LH recombination of H* and Cl* is merely 0.34 eV. Even if the energy of H* is reduced by 1.0 eV, it still has sufficient energy to overcome such a low LH barrier.
It should perhaps be noted that the IAA in which the LDFA is based does not account for the tensorial nature of the friction. As a result, the simplification of the friction tensor by atomic Finally, let us discuss the reaction cross section, which is defined as the sampling unit cell area multiplied by the reaction probability. The values are listed in Table I for the HA and ER channels at both incidence angles. For normal incidence, the HA cross section obtained from AIMD simulation is 20.6 Å 2 , which is much larger than that of the ER one (3.1 Å 2 ). The total non-LH cross section is thus 23.7 Å 2 . The corresponding values for off-normal incidence (16.9, 8.4, and 25.3 Å 2 ) are quite similar. As discussed above, the inclusion of EHPs does not qualitatively change this picture. These theoretical cross sections are much larger than the experimentally reported value of 3 ± 1 Å 2 . 10,11 It should be noted, however, the experimental cross section was estimated assuming that all non-LH reactivities derive from direct ERtype collisions, which is now known to be unrealistic. Our calculated cross sections are quite close to what Quattrucci and Jackson reported in their earlier work when the EHP effect was not considered, 27 although their semi-empirical damping model seems to overestimate the EHPs effect when comparing with the more quantitative results reported here. Interestingly, Rettner stated that in experiment, "the total probability for Cl removal is >0.5 at 600 K," 11 which is not inconsistent with theory (75%-80%) although the experiment number includes all (ER, HA, and LH) mechanisms. Another factor is the Cl coverage. At a smaller Cl coverage, the cross section is expected to be smaller. However, it is difficult to simulate low coverages with AIMD because of the large unit cell required. We also note in passing that the setup used in this work might not correctly characterize the so-called secondary HA events, which involve collisions with more than one adsorbate, because the Cl atoms are not independent because of the periodic boundary conditions.

IV. CONCLUDING REMARKS
In summary, we present in this work extensive AIMD simulations of a prototypical Eley-Rideal type reaction between impinging H atoms and chlorine atoms adsorbed on Au(111). Such simulations allow a quantitative analysis of the energy dissipation into both the surface phonons and EHPs, and their impact on the reaction dynamics. Our results confirm those found in previous theoretical simulations, namely, the domination of the indirect HA mechanism over the direct ER counterpart and the high internal excitation of the HCl product. In addition, the angular and velocity distributions are in reasonably good agreement with available experimental data.
The most significant aspect of the current work is the quantitative elucidation of the role of surface phonons and EHPs in energy dissipation, and their impact on the reaction dynamics. To this end, the simulation results indicate significant energy exchange with both the surface phonons and EHPs. Specifically, the nonadiabatic EHP-induced dissipation depends strongly on the incidence angle. For normal incidence, the H atoms access high electron-density regions in both the first and second layers of the metal surface, resulting in much stronger dissipation than that in the off-normal incidence. However, the energy dissipation is insufficient to alter the reactivity and internal state distribution of the HCl product.
While many insights have been gained from the AIMD simulations, the number of trajectories and propagation time are quite limited. Further investigations of more trajectories will help one to establish more quantitative conclusions and to understand the dependence of dynamics on different experimental conditions, such as the incidence angle and energy as well as surface structure and coverage. This is now achievable by high-dimensional potential energy surfaces that include the surface degrees of freedom. 74,75 Work in that direction is underway.
Furthermore, the description of the EHP-induced dissipation can be further improved. It has recently been shown that the corresponding friction can be strongly directional, 71,72 and a tensorial treatment of the friction might be needed to capture the inherent complexity of the possibly mode dependence of energy dissipation that is largely ignored in the LDFA approach.