Quantum Effects in the Diffusion of Hydrogen on Ru(0001).

An understanding of hydrogen diffusion on metal surfaces is important not only for its role in heterogeneous catalysis and hydrogen fuel cell technology but also because it provides model systems where tunneling can be studied under well-defined conditions. Here we report helium spin-echo measurements of the atomic-scale motion of hydrogen on the Ru(0001) surface between 75 and 250 K. Quantum effects are evident at temperatures as high as 200 K, while below 120 K we observe a tunneling-dominated temperature-independent jump rate of 1.9 × 109 s-1, many orders of magnitude faster than previously seen. Quantum transition-state theory calculations based on ab initio path-integral simulations reproduce the temperature dependence of the rate at higher temperatures and predict a crossover to tunneling-dominated diffusion at low temperatures. However, the tunneling rate is underestimated, highlighting the need for future experimental and theoretical studies of hydrogen diffusion on this and other well-defined surfaces.

Hydrogen (H) atoms can exhibit significant quantum nuclear effects in many and varied materials, such as H-bonded crystals and ferroelectrics,[1−3] high pressure ice,[4] and at the surface and in the bulk of metals,[5−7] where the rate of diffusion is enhanced by quantum tunneling. Studies of H diffusion on metal surfaces provide model systems where tunneling can be studied in great detail under well-defined conditions.

Experimental measurements of H diffusion on metal surfaces are, however, extremely challenging and have at times been contentious. This is illustrated, for example, by the case of H on Ni(111), where initially a sharp classical to quantum crossover was reported along with similar rates for H and D,[8,9] whereas later experiments found only a weak crossover and normal isotope effect.[10] Subsequent calculations from Badescu et al. explained the experimental observations through a mechanism involving tunneling from excited vibrational states.[11] For H on Cu(001), scanning tunneling microscopy (STM) experiments[12] revealed a sharp crossover around 60 K to a nearly temperature-independent tunneling regime, and subsequent density functional theory (DFT) calculations[13−15] of the tunneling rate successfully reproduced the experimental results. Despite these studies, a detailed view of the underlying mechanisms behind the quantum behavior of H and other light adsorbates on surfaces is still lacking, and well-defined comparisons between experiment and theory are uncommon.

We present here a combined experimental and theoretical study of the temperature dependence and mechanism of diffusion of H on Ru(0001), a system well-suited to such an approach: widespread experimental interest exists,[16−24] and previous DFT studies[24−27] found good agreement with experimental adsorption geometries and vibrational frequencies, suggesting that DFT provides an accurate description of the potential energy surface (PES) – a prerequisite for reliable rate calculations.

Experimentally, we use the novel helium spin–echo technique[28] that allows determination of surface dynamics on pico- to nanosecond time scales through measurement of the intermediate scattering function (ISF) I(ΔK,t), a measure of surface correlation after a time t on the direction and length-scale given by ΔK. Our calculations combine a DFT description of the electronic structure with path-integral molecular dynamics (PIMD) to obtain a unified description of jump rates at high and low temperatures with full account of the coupling to surface phonons. The PIMD technique provides a quantum mechanical description of configuration space by representing atomic nuclei as ring polymers of interconnected beads, and the limit of exact quantum mechanical properties can be approached by increasing the number of beads. To obtain jump rates including tunneling contributions, we then combine the PIMD simulations with quantum transition-state theory (QTST).[29,30] Despite a significant computational cost, our calculations demonstrate for the first time the feasibility of investigating tunneling-assisted surface diffusion entirely from first principles using PIMD-based QTST.

We show that the diffusion of H on Ru(0001) exhibits uniquely interesting properties. At the highest temperatures studied overbarrier hopping between nondegenerate fcc and hcp sites is seen, while a constant diffusion rate of ∼3 × 109 s–1 below ∼120 K indicates quantum tunneling at rates many orders of magnitude faster than observed for other systems: 4 decades faster than the highest temperature independent rate for H on Ni(111),[9] 5 decades faster than H on Cu(111),[31] and 12 decades faster than for H on Cu(001).[12] Additionally, for H on Ru(0001), the presence of multiple jumps indicates low adsorbate–substrate friction compared with, for example, the diffusion of H on Pt(111), where HeSE was used to study quantum contributions to the activated motion.[32] We also find evidence of repulsive interadsorbate interactions. The overall temperature dependence of ab initio PIMD rates agrees reasonably with experiment, with good agreement down to around 200 K. At low temperatures the rate is underestimated, and we use a harmonic quantum transition-state theory (HQTST) version of instanton theory[33] to shed light on this discrepancy. This method is also based on Feynman’s path-integral formalism and, in the present implementation, provides the effective lowering of the energy barrier on a rigid PES due to quantum tunneling; it thus does not treat all degrees of freedom explicitly as in PIMD simulations.

Figure 1a shows a typical measurement (an ISF) at 250 K, in which a loss in surface correlation is seen over time as I(ΔK,t) decreases. This is attributed to diffusion of H on the Ru(0001) surface; measurements on a clean Ru(0001) surface show no such decay. By treating each ISF as a combination of exponential decays exp(−αt), dephasing rates (α) can be obtained. Two exponential decays were needed to model the experimental data at 250 K along ⟨11̅00⟩ and ⟨112̅0⟩; this is indicative of hopping between two nondegenerate sites on a non-Bravais lattice, here the fcc and hcp hollow sites are shown schematically in Figure 1b.

Figure 1

(a) Example of an ISF from experiment for T = 250 K, ΔK = 1.76 Å–1 along ⟨112̅0⟩, with the results of the best-fit Monte Carlo simulation (solid red line) and a one exponential fit to the long-time limit data (dashed green line) also shown. (b) Schematic of the Ru(0001) surface showing the ⟨112̅0⟩ and ⟨11̅00⟩ directions. (c) α–ΔK plots along ⟨112̅0⟩ (left) and ⟨11̅00⟩ (right) at 250 K. Data points are shown with black crosses, and the result of the best fit Monte Carlo simulations (see text) is shown with a thick red line.

A Bayesian method was used to fit the entire HeSE data set at 250 K (85 ISFs) with two exponential decays using the analytic model of Tuddenham et al.,[34] with the basic jump rate and energy difference between the fcc and hcp adsorption sites as variable parameters. This gave a peak in the relative probability density at the most probable value of the fcc-hcp site energy difference, 22.2 meV, with a statistical uncertainty of ±0.6 meV. Figure 1c shows the variation of the slow decay constant α with ΔK. The steep approach to the origin at low ΔK is indicative of multiple jumps (low adsorbate–substrate friction).[28] A dip in α at ΔK = ∼1.3 Å–1 is also seen (more obvious along ⟨11̅00⟩ due to the shape of the curves), arising from repulsive interadsorbate interactions, which stabilize ordering when the adsorbates are as far apart as possible. Correlations at values of ΔK corresponding to this separation tend to decay more slowly, giving the dip in the dephasing rate and a corresponding diffuse ring in helium diffraction scans.[35−37]

A Monte Carlo simulation[35] of hopping between fcc and hcp sites for 0.2 ML H coverage provided quantitative insights into the diffusion mechanism and gave a second estimate of the fcc-hcp site energy difference as 18.7 meV, with a statistical uncertainty of ±0.3 meV. Systematic uncertainty from approximations inherent in this and the analytic model[34] exceed the statistical uncertainties, but we can conclude that these independent measures of the site energy difference are consistent with a value of (20 ± 5) meV. The variation of α from the Monte Carlo ISFs with ΔK is shown with the thick line in Figure 1c, which generally fits the data well. The biggest deviation from experiment is at ΔK ≈ 1.3 Å–1, the position of the de-Gennes feature,[36] suggesting that the form of the repulsive interadsorbate interaction is more complex than the simple dipole–dipole repulsion used here.

Figure 2 reports the H diffusion rates measured across a broad range of temperatures, from 250 K down to 75 K. At the highest temperatures a linear dependence of the logarithm of the diffusion rate with 1/T is observed. As the temperature is reduced the diffusion rate levels off, until below ∼120 K the diffusion rate appears to be independent of T, to within the experimental errors. Note that even the lowest rates reported in this Figure are comfortably within the experimental range of the technique, as discussed in ref (35). The linear dependence at high T indicates activated H diffusion, and a fit to all data shown of the form AT exp(−Ea/kBT) + C yields an estimate of the classical activation energy for fcc-hcp hopping Ea = (95 ± 3) meV. The pre-exponential factor is (4.5 ± 0.6) × 109 s–1 K–1, and the temperature-independent tunneling rate C is (1.9 ± 0.1) × 109 s–1.

Figure 2

Experimental rates (black squares), calculated jump rates from PIMD (red line with squares, error bars denoted by thick shaded lines), and instanton calculations on a 1D potential (black line with circles). The dashed line shows instanton results where the fcc-hcp energy difference is reduced to 25 meV. The inset shows high-temperature PIMD rates scaled by 0.25, where PIMD results agree with the slope of the experimental data.

Having obtained experimental rates for H diffusion across a broad range of temperatures, we now explore the system with DFT. Our calculations show that the fcc site is the most stable adsorption site, with the hcp site less stable by ∼50 meV, which is ∼30 meV larger than the experimental estimate. The lowest energy diffusion pathway is across the bridge site with an activation barrier of 150 meV, which reduces to 120 meV when zero-point energy (ZPE) effects are taken into account within the harmonic approximation. The experimentally observed rate cannot be described by classical transition-state theory as quantum effects are clearly seen. We therefore turn to ab initio PIMD, an approach that can capture the change of the quantum free-energy barrier for diffusion due to the quantum nature of the proton and thus account for tunneling and ZPE effects beyond the harmonic approximation.

Figure 3a shows quantum free-energy barriers from PIMD at different temperatures and, for comparison, results of classical-nucleus MD simulations at 250 K. The free-energy barrier obtained with classical nuclei is at ∼150 meV very similar to the underlying potential energy barrier. At the same temperature, the free-energy barrier obtained from PIMD is only slightly lower (by ∼15 meV) due mainly to ZPE effects. The fact that the slight lowering of the barrier predominantly comes from ZPE effects and not tunneling can be seen from panel b, where the width of the H probability distribution is plotted as a function of H position along the reaction coordinate. At 250 K the width of the H distribution is unaffected by the (fixed) position of the H centroid along the fcc-hcp path. However, as the temperature is lowered the H probability distribution broadens and the quantum free-energy barrier drops. At 70 K the quantum free energy barrier is 99 meV, 65% of the original classical barrier. This substantial reduction of the quantum free energy barrier arises from tunneling of the H through the potential barrier, as can be seen from the partially bimodal nature of the H probability distribution function (Figure 3b).

Figure 3

(a) Free-energy barriers for H diffusion from fcc to hcp sites from ab initio PIMD compared with the temperature-independent barrier from classical-nucleus ab initio MD. Above, several simulation snapshots are superimposed to show the distribution of beads when the centroid position is constrained at the initial site (fcc, left) and at the transition state (bridge, right) at 70 K. (b) Projected H-atom bead probability distributions at a selected series of fixed centroid positions along the fcc-hcp path. Dashed vertical lines indicate the corresponding fixed centroid positions.

Figure 2 reveals two main effects from the quantum treatment of nuclei in PIMD. First, the crossover to tunneling-dominated diffusion is reflected by the changes in the PIMD rate, with a gradual change of slope in the Arrhenius plot below 100 K. Second, at intermediate to high temperatures the gradient from PIMD agrees well with the experiment due to the decrease in the quantum free-energy barrier and the inclusion of ZPE effects, highlighting the importance of quantum nuclear effects even at high temperatures. Indeed, fitting the same Arrhenius expression to the PIMD results as used for the experimental rates we obtain Ea = 105 meV, close to the experimental value of 95 meV. However, PIMD gives a larger prefactor than the experiment, resulting in rates faster by a factor of ∼4 at high temperatures. This is likely to be due to the use of the classical velocity in the expression for the rate. (See the Experimental Section.) To illustrate the agreement between the experimental and PIMD gradients at high temperatures, we show in the inset of Figure 2 theoretical rates that have been scaled by 0.25.

Despite the overall qualitative agreement with experiment, it is apparent that at low temperatures the PIMD results underestimate the experimental rate by two to three orders of magnitude. This could potentially be explained by the difference between the calculated and experimental estimates of the fcc-hcp energy difference and activation energy. To investigate this discrepancy in the low-temperature rates we use a more flexible and computationally efficient approach than PIMD, specifically instanton theory.[33] For the instanton calculations we used a 1D potential represented by a polynomial fitted to the underlying DFT PES along the fcc-hcp diffusion path. Initially we validated the instanton calculations by comparing the rates it produced with those from PIMD, and as shown in the low-temperature region of Figure 2, the instanton and PIMD results agree quite well. By construction, instantons only provide rates in the tunneling regime. Using the instanton approach we then investigated the sensitivity of the computed rates to the fcc-hcp energy difference. To this end we modified the 1D DFT potential by reducing the activation energy and hence also the fcc-hcp energy difference by 25 meV, bringing the binding site energy difference into close agreement with the experimental result. The dashed lines in Figure 2 show instanton rates for this modified potential. Indeed, significantly better agreement with the low-temperature experimental rates is now seen. The smaller gradient, due to the smaller fcc-hcp energy difference and thus smaller lattice deformation energy required for energy coincidence between fcc and hcp vibrational states, agrees better with the nearly temperature-independent experimental rate, and the magnitude of tunneling is strongly enhanced by the small reduction in the tunneling barrier and the reduction in the energy difference between the initial and final states. These results highlight the sensitivity of the computed jump rates to details of the PES and hence the importance of performing accurate calculations. The rather good agreement between instanton calculations and PIMD rates on the unmodified PES show that polaron effects play a minor role on this surface, in agreement with the experimental observation of low adsorbate–substrate friction.

In summary, we have presented a combined experimental and theoretical study of the diffusion of H on Ru(0001), showing a transition from overbarrier hopping between fcc and hcp adsorption sites at high temperatures to quantum tunneling at low temperatures. The experimental results using HeSE give evidence of low adsorbate–substrate friction and repulsive interadsorbate interactions, and our experimental tunneling rates for H are much higher than those seen previously for other surfaces such as Ni(111),[8−10] Cu(111),[31] and Cu(001).[12] Tunneling from excited vibrational states was found to play a key role for H diffusion on Ni(111).[11] For H on Ru(0001), however, the low-temperature rate in Figure 2 is effectively nonactivated (the data are fitted by an activation energy of (− 2 ± 2) meV[35]) and so tunneling from the excited states (the lowest of which has been determined from HREELS to be 84 meV[22]) does not contribute significantly at low temperatures. The fcc-hcp energy barrier is 150 meV for Ru(0001) but closer to 200 meV for Ni(111),[10] which would lead to faster tunneling for the case of H on Ru(0001) in the ground state. Recent STM results for H diffusing on Cu(111) at 5 K[31] (comprising streaks on STM images) were interpreted to give an H jump rate of 30 Hz. However, in addition to the possible presence of tip induced effects, the imaging rate was only 5.5 lines s–1 and so would be insensitive to diffusion as fast as reported here.

Ab initio PIMD rates show reasonable agreement with the temperature dependence of the experimental results at high temperatures, while the tunneling rates at low temperatures are significantly underestimated, similar to previous studies on other surfaces employing empirical potentials.[38−40] This discrepancy is highly interesting, and by performing instanton calculations on a modified PES we showed that small changes in activation energy and energy difference between fcc and hcp adsorption sites may play a role. However, other features of this system may also be important in driving the faster experimental than theoretical rates, which highlights the need for future experiments and theoretical calculations of H diffusion on this and other well-defined metal surfaces.

Experimental Section

Experimental measurements were carried out using the Cambridge HeSE spectrometer,[28] which measures the ISF I(ΔK,t), related to the dynamic structure factor S(ΔK,Δω) and the pair correlation function G(R,t) by temporal and spatial Fourier transforms, respectively.[28] As such, HeSE is a uniquely noninvasive technique that can be used to study the detailed mechanism of individual atomic jumps on well-defined surfaces. The single crystal Ru(0001) sample (Surface Prep. Lab., The Netherlands) was cleaned by repeated cycles of argon sputtering, flash annealing, oxidation, and reduction,[41] and the surface quality was monitored using helium reflectivity measurements. Molecular H (Air Liquide, 99.999% purity) was dosed by backfilling the scattering chamber, and surface uptake of atomic H was followed using the specularly scattered helium beam to achieve a coverage of 0.2 monolayer (ML), determined from the known dose of H and sticking probability.[16]

All DFT calculations were performed within a periodic plane-wave framework using projector-augmented wave potentials in the VASP code[42,43] modified for PIMD simulations.[44] The PIMD simulations were performed with 16 replicas (i.e., beads of the ring-polymers). Tests with 32 and 64 beads showed very good agreement with the 16 bead simulations. (See the Supporting Information for computational details.[35])

PIMD-based QTST[29,30] gives the jump rate over the reaction coordinate ξ (the fcc-hcp diffusion path) aswhere v̅ is the thermal velocity: v̅ = (2kBT/πm)1/2. P(ξ‡,T) and P(ξ0,T) are the centroid (center of mass of the ring polymer) probability densities at the transition state ξ‡ and initial (fcc) state ξ0, respectively, and ΔF is the corresponding free-energy difference obtained from a series of constrained-centroid PIMD simulations.[35] Although representing the reaction coordinate in terms of centroid positions becomes inaccurate for highly asymmetric potentials at temperatures far below the crossover to quantum tunneling,[45−47] it should provide an accurate description for the moderately asymmetric potential just below the crossover temperature as studied here.