Semiclassical second order vibrational perturbation theory for hopping rates of H and D atoms on Pt(111) and H on Ru(0001)

Abstract

Wolynes's theory for tunneling in dissipative systems was constructed for parabolic barriers. Its recent generalization to anharmonic potentials is further developed, most notably by employing the second order vibrational perturbation theory expression for the action, derived by Miller and coworkers. With this construct there is no need to know the full potential energy surface involved in the hopping, it is sufficient to know well and barrier frequencies, barrier heights, friction coefficients and the fourth order derivative of the potential at the barrier top. The resulting theory is applied to model the experimentally measured hopping rates of H and D atoms on a Pt(111) surface and H atom hopping over a barrier on the Ru(0001) surface. In all cases, the results indicate that the barrier frequency is substantially larger than the well frequency. The comparison with experiment sheds light on the information content of the measured data.

---

1 Introduction

In 1940, H. A. Kramers proposed a new theory for the rate of chemical reactions.¹ His approach introduced stochasticity in this field by considering Brownian motion dynamics that take into account the interaction of the reacting system with the surrounding environment. He considered the mechanism of escape processes as a noise-assisted reaction driven by thermal forces related to temperature and friction by the fluctuation–dissipation theorem. The derivation of an energy diffusion equation for the dynamics of reactive particles in the weak-friction case was a pioneering contribution. He also derived analytical expressions for the probability of escape in the moderate-to-large damping limit. Kramers’ result in this damping regime has been applied among others to multidimensional Brownian motion,² non-Markovian processes,³ and dissipative quantum mechanical tunneling systems.⁴

In 1980, Wolynes⁴ generalized Kramers’ spatial diffusion limited rate expression, including in it quantum tunneling effects using a parabolic barrier coupled to a dissipative bath. Wolynes’ expression is valid for temperatures which are above the so called crossover temperature⁵ from deep tunneling to above barrier crossing. Caldeira and Leggett developed a low temperature theory for dissipative tunneling.⁶ Chang and Chakravarty⁷ applied their theory numerically for the entire range of possible friction coefficients with Ohmic dissipation. Grabert, Weiss, and Hänggi,⁸ using the Caldeira–Leggett methodology, derived an expression for the effect of temperature on the tunneling rate, valid also for Ohmic dissipation. This analysis was applied to the experimental results reported by Washburn et al.,⁹ who measured tunneling rates for several Nb edge junctions.

More recently, the helium spin echo (HeSE) surface technique¹⁰ was used to study the diffusion of H and D atoms on metal surfaces over a wide range of surface temperatures. The experimental data was modeled, using Wolynes’ expression for the hopping rate. The analysis suggested that at high temperature, diffusion was dominated by classical activated over-barrier hopping while at lower temperature quantum mechanical tunneling sets in ref. 11. The fit of Wolynes’ theory to the experimental data enabled the extraction of friction coefficients, barrier frequencies and barrier energies. However, the reported friction coefficients for H and D atoms of the order of ns⁻¹ were unreasonably high. Using the theory developed by Grabert and Weiss^12,13 which assumes that the whole process is dominated by tunneling and thus circumvents the crossover divergence problem, Sanz et al.¹⁴ were able to fit the experimental data rather well with friction coefficients that were of the order of a few ps⁻¹.

Most recently, a new approach was suggested whose aim was to overcome the crossover temperature divergence problem which appeared in the Wolynes spatial diffusion expression.^15–17 It is predicated on consistent use of the uniform semiclassical energy-dependent transmission coefficient of Kemble.¹⁸ The resulting theory does not display any singularity at the crossover temperature. This development opens the door to a more precise analysis of experiments that consider dissipative tunneling over a broad temperature range. The central theme of the present work is to develop this “dissipative tunneling theory” to analyze the tunneling diffusion for H and D adsorbates on a Pt(111) surface. For this purpose it is necessary to estimate the tunneling action, which will be done by using the second order vibrational perturbation theory (VPT2) which gives a quadratic energy action relation.^19–21 It also points out that one should not ignore the zero point energy contribution from anharmonic components of the potential. We will show that removing the limitation of tunneling only through a parabolic barrier as in Wolynes’ theory leads to a much more reasonable theoretical fit of the measured diffusion rates than obtained previously.

There is an additional aspect to the present analysis. Another system which has thus far evaded a “good” theoretical description is the measured²² thermal hopping rate of H on Ru(0001). In this case, the hopping rate seems to flatten out at low temperatures and attempts to simulate the effect were off by orders of magnitude. Here too, we will show that the present VPT2 uniform semiclassical theory leads to a reasonable description of the experimental results.

This paper is organized as follows. The VPT2 uniform semiclassical theory generalization of Wolynes’ result is presented in Section II. The theory is then applied in Section III to the tunneling of H and D atoms on a Pt(111) surface and H on the Ru(0001) surface. We end in Section IV with some conclusions as well as future perspectives.

2 Theory

2.1 Preliminaries

We assume that the classical motion of an atom of mass M on the surface is well described classically by the generalized Langevin equation (in one dimension)where V(q) is the potential, q is the coordinate, dots denote time derivatives, γ(t) is the time dependent friction function and F(t) the random force related to the friction function through the fluctuation dissipation relation.²³

It is well known²⁴ that the classical Langevin equation may be reformulated in terms of a Hamiltonian in which the system interacts bilinearly with a harmonic bath

where x_j and p_xj are the mass weighted coordinate and momentum of the j-th bath oscillator whose frequency is ω_j which is coupled to the system via the coupling constant c_j. The friction function appearing in the generalized Langevin equation is identified as

Without loss of generality we assume that the potential V(q) has a well located at q = 0 with V(0) = 0 and well frequency ω₀ and a barrier located at q^‡ with barrier energy V^‡ and barrier (imaginary) frequency ω^‡. If the potential is purely quadratic the Hamiltonian may be diagonalized, using a normal mode transformation, the details of which may be found for example in ref. 25. The transformed Hamiltonian expressed in the normal modes around the well is

where y_j, p_yj are the mass weighted phase space variables of the j-th stable bath mode whose frequency is λ_j. Around the barrier the normal mode form of the Hamiltonian iswhere ρ, p_ρ are the unstable normal mode and momentum, y^‡_j and the coordinate and conjugate momentum of the j-th stable normal mode at the barrier. The unstable mode barrier frequency is expressed in terms of the friction function through the Kramers–Grote–Hynes relation^1,26

λ^‡2 + [small gamma, Greek, circumflex](λ^‡)λ^‡ = ω^‡2, (6)

where the ‘hat’ notation denotes the Laplace transform with respect to the argument.

2.2 Wolynes’ expression for the rate

In the spatial diffusion-limited regime, which is relevant to the systems we will study due to the fact that the surface does not have linear channels,²⁷ the rate is written aswhere the classical transition state theory (TST) rate is

The parameter ν takes on the possible values of multiplicity of equivalent paths of escape, for both systems studied here ν = 3.

The Wolynes factor Ξ_W in eqn (7) is the parabolic barrier based ratio of the quantum partition function at the barrier to that of the well expressed in terms of the so-called Matsubara frequencies

ν_n = 2πn/(βħ) = nbω^‡ (9)

as⁴

In the classical limit b → ∞, all the Matsubara frequencies go to infinity and Ξ_W → 1.

2.3 The generalized Wolynes factor

The difficulty with the Wolynes factor is, as already mentioned, that it diverges at the “crossover temperature” ħβ_cλ^‡ = 2π, which is due to the fact that its derivation is based on a parabolic barrier potential. When the friction coefficient γ is not too large γ/ω^‡ ≤ 1 one may assume that the motion along the unstable mode is separable from the other modes. Using quantum transition-state theory,²⁸ the rate can be written aswhere Q₀ is the partition function of the reactants in the well which when assumed to be harmonic isQ^‡ is the partition function of the stable normal modes at the barrierand κ is the thermal quantum transmission factor through the (one dimensional) barrier defined as:where T(E) is the transmission probability through the barrier at energy E.

For a parabolic barrier potential with barrier frequency λ^‡ the parabolic barrier energy dependent transmission factor is

and the parabolic barrier thermal transmission factor (letting the lower limit of the integration in eqn (14) go to −∞) is

The corresponding thermal rate is

Comparing with eqn (7) implies that the Wolynes factor may also be written as

and as shown in ref. 28 this expression is equivalent to Wolynes’ original expression written in the continuum limit as in eqn (10).

2.4 The uniform semiclassical rate expression

As shown in some detail in ref. 15–17 the key to elimination of the divergence at the crossover temperature is the consistent use of Kemble's uniform semiclassical energy-dependent transmission coefficient,¹⁸ formalized by Fröman and Fröman²⁹

Here S(E) is the Euclidean action on the inverted barrier at energy E.

The uniform semiclassical transmission probability is defined by replacing the exact transmission probability with its uniform approximation

Multiplying and dividing the Wolynes factor (eqn (10)) with the parabolic barrier transmission factor (eqn (16)) and then replacing the parabolic transmission factor in the numerator with the uniform semiclassical factor allows us to define a uniform Wolynes factor¹⁶

which does not have a divergence at the crossover temperature.

Using an adiabatic TST implies taking into consideration the net zero point energy which is the difference between the zero point energy at the well and the barrier

The second equality gives a practical algorithm for computing the zero point energy.

Within the adiabatic TST the uniform semiclassical transmission factor is modified to

The working equation for the hopping rate, which shall be used in the next section is then

where the subscript “AD” reminds us that the lower limit for the integration of the transmission factor is the zero point energy as in eqn (23).

2.5 Semiclassical vibrational perturbation theory

It remains to determine the Euclidean action, S(E). For this purpose, we use the second order vibrational perturbation theory (VPT2) approximation. As suggested by Handy, Miller and coworkers^19–21 and further implemented by Stanton, Barker and coworkers,³⁰ second order perturbation theory implies a quadratic relation between the energy and the Euclidean actionwhere the quadratic coupling constant χ depends on the third (V₃) and fourth (V₄) order derivatives of the potential at the barrierand E₀ is the energy shift due to the anharmonicity

The quadratic energy action relation may be inverted, giving an explicit dependence of the action on the energy

In the following section we will employ the VPT2 theory. It is well understood that it typically overestimates the rate in the deep tunneling region.³¹ A more accurate approximation would be to compute the Euclidean action on the potential energy surface.³² However, in practice, this potential energy surface is not available. The strategy we will use is to consider the higher order potential derivatives at the barrier as unknowns to be fitted to the experimental data.

3 Application to H and D atom diffusion

3.1 H and D on Pt(111)

The Pt(111) surface has a three fold symmetry, so the value ν = 3 (see eqn (8)) was used throughout. Within the framework of the theory presented above, we note the following parameters which are needed: the masses, barrier and well frequencies, barrier height, potential parameter V₄ (for simplicity we assume a symmetric potential) and friction coefficient – we assume Ohmic friction so that [small gamma, Greek, circumflex](s) = γ. Not all of them are free parameters. The masses are those of the H and D atoms; the well frequency for the H atom adsorbate has been measured using HREELS³³ and has the value ħω_0H = 30 ± 2 meV, we will use the value of 30 meV throughout. There is also a mass scaling between the various parameters for H and for D. The frequencies scale with mass such that the D frequency is less than the H frequency and the friction coefficients also scale with the mass γ_D = γ_H/2. This implies that to fit the experimental data of both H and D diffusion, we have the following free parameters – the friction coefficient, the barrier frequency for the H atom, the anharmonicity constant χ and the barrier heights for the H and D atoms separately. In principle the potential is invariant to isotopic substitution, however within the vibrationally adiabatic framework, the vibrational modes perpendicular to the reaction coordinate are unknown and affect the net barrier height for H and D differently so that as long as the differences are not too large we may assume somewhat different barrier heights.

Fig. 1 shows an Arrhenius-like plot of the tunneling rates (given in sec⁻¹), comparing experiment with theory. The free parameters used for the H atom are the barrier frequency ħω^‡ = 130 meV, barrier height V^‡ = 129 meV and friction coefficient ħγ = 26 meV. The anharmonicity parameter in both cases was χ = 1.9 a.u., for the D case, the barrier height is V^‡ = 110 meV. The numerical values of the rates are detailed in Table 1.

Fig. 1 Comparison of experimental³⁴ and theoretical hopping rates for H and D on Pt(111). Data for H are colored blue, for D red. The solid lines are the theoretical result based on eqn (24), the dotted lines are obtained with the parabolic barrier theory (eqn (7)).

Table 1 Experimental tunneling rates (Exp) extracted from ref. 34, the corresponding numerical values (VPT2) from the second order vibrational perturbation theory (eqn (24)), the percent relative error for each temperature (%-error), and the average error for H and D. The fitted parameter values are given in the text

Tunneling rate × 10^11 [s^−1]
1/T	Exp-H	VPT2-H	% error
4.679	1.364 ± 0.09	1.271	7.3
4.937	1.02 ± 0.06	1.037	1.6
5.194	0.949 ± 0.05	0.852	11.4
5.503	0.723 ± 0.04	0.678	6.6
5.828	0.592 ± 0.03	0.538	10.0
6.618	0.337 ± 0.02	0.319	5.6
7.115	0.266 ± 0.02	0.235	13.2
8.298	0.126 ± 0.011	0.122	3.3
9.498	0.0572 ± 0.003	0.0680	18.9
9.978	0.057 ± 0.003	0.0550	3.6
10.51	0.043 ± 0.005	0.0439	2.1
11.07	0.0383 ± 0.003	0.0351	9.1
11.09	0.0321 ± 0.002	0.0348	8.4
11.74	0.0251 ± 0.003	0.0272	8.4
12.49	0.0154 ± 0.005	0.0208	35.1
Average error:		9.6

---

1/T	Exp-D	VPT2-D	% error
4.344	1.449 ± 0.12	1.487	2.6
4.535	1.232 ± 0.09	1.258	2.1
4.764	1.031 ± 0.06	1.033	0.2
4.994	0.758 ± 0.07	0.851	12.3
5.261	0.655 ± 0.05	0.683	4.3
5.548	0.522 ± 0.04	0.542	3.8
5.873	0.39 ± 0.03	0.420	7.7
6.255	0.326 ± 0.02	0.315	3.5
6.656	0.225 ± 0.02	0.235	4.4
8.338	0.081 ± 0.004	0.0784	3.3
9.522	0.039 ± 0.003	0.0404	3.6
10.02	0.032 ± 0.005	0.0314	1.9
10.54	0.022 ± 0.004	0.0245	11.4
11.11	0.017 ± 0.003	0.0189	11.2
11.78	0.014 ± 0.005	0.0143	2.1
Average error:		5.0

---

The resulting fits are not overly sensitive to the parameters used. One may lower the barrier height, this would lead to increased tunneling and so it becomes necessary to lower the barrier frequency and increase the friction coefficient. For example by using a barrier height of 104 meV for H and 93 meV for D, with a barrier frequency of 120 meV, a friction coefficient of 41 meV and χ = 1.4 a.u. one still gets an “acceptable” fit to the experimental data. At this point, one cannot pinpoint the results any further, for this one needs some quantum chemistry computations.

To get a feeling for the magnitude of V₄ we note that the zero point energy lowering of the barrier (eqn (27)) is E₀ = −9.5 meV for H and −5.2 meV for D. These values are small compared to the full barrier height, pointing out that the fourth order derivative is not excessively large. It is also of interest to note the ground state energy at the well in the presence of dissipation (eqn (22)). One finds that ε₀ = 18 meV for H and 12 meV for D. Summing the two results, implies that the adiabatic barrier for the D atom (not including perpendicular modes) is ∼93 meV, sufficiently close to the adiabatic value of 101 meV for H.

We note the good agreement between the theoretical fits and the experimental data. Secondly, one sees that the parabolic barrier based results are “good” in the high temperature limit, implying that the assumption¹⁴ of tunneling for all measured temperatures is questionable. It is though necessary to introduce the uniform semiclassical approach to allow for extension of the theory to temperatures below the divergence temperature of the parabolic barrier approximation. Thirdly, the reduced friction coefficient is which is sufficiently lower than unity to justify the one dimensional separable model used for the theory.

The parameters used are not unique. One obtains similar results when changing the barrier frequency to 120 or 140 meV with small changes in the other parameters. However, it is clear that for agreement to be reasonable the barrier frequency must be four to five times as large as the well frequency. Importantly, the friction coefficient is now much smaller than the value of ħγ = 431 meV used in ref. 34 in the context of the parabolic barrier theory. It was only through such a large value of the friction that one could obtain a “crossover temperature” below the experimentally measured range. However, previous studies indicate that it was too large,³⁵ the present value seems much more “reasonable”.

3.2 H on Ru(0001)

The experimental hopping rates for this system were reported in ref. 22. Here too there is a three fold symmetry as may be seen from Fig. S3 in the SI of ref. 22. The fit to the experimental data is shown in Fig. 2. The frequency in the well – 84 meV – was measured experimentally.³⁶ The fitted parameters are ħω^‡ = 130 meV, V₄ = 1.93 a.u., V^‡ = 231 meV, ħγ = 2.5 meV. These parameters imply that the zero point energy shift at the barrier is E₀ = −9.5 meV and the zero point energy at the well is ε₀ = 19 meV so that the adiabatic barrier height, not including perpendicular modes is ∼202 meV. This is somewhat higher than the value reported in ref. 22 for the free energy barrier of 135 meV.

Fig. 2 Comparison of experimental²² and theoretical hopping rates for H on Ru(0001). The solid line is the theoretical result based on eqn (24), the dotted line shows the parabolic barrier theory result (eqn (7)).

Inspection of Fig. 2 shows good agreement between the experimental and theoretical hopping rates at high temperatures and somewhat less so at the lowest temperatures. However, the agreement is much better than that reported between path integral molecular dynamics (PIMD) and transition state theory rate estimates, which are two orders of magnitude lower than the experimental. In this context, it is interesting to note that considering the barrier zero point energy E₀ would raise the rate by a factor of ∼5 at T = 70 K. This effect is not included in the PIMD computations.

4 Discussion

The theory and computations presented in this paper stress the importance of using the uniform semiclassical rate theory which includes in it zero point energy corrections due to the anharmonicity of the potential. In both cases, whether motion on Pt(111) or Ru(0001), the fitted barrier frequency is substantially larger than the well frequency and this accounts for the extensive tunneling at low temperatures.

The theory used here may be improved. The use of the VPT2 action at low temperatures is questionable. It would be much better to have a global potential energy surface and compute instantons on it within the framework of the uniform semiclassical theory as described in ref. 32. Yet, even this is not necessarily sufficient as one should include multidimensional effects, which especially influence the activation energy. The main difficulty is though obtaining a reliable force field which includes both the motion of the adatom, its interaction with the surface modes and the surface modes themselves. This is a tall order, forcing us at this point to rely on semiempirical methods, as also described in this paper.

There are additional points to be considered. We assumed harmonic motion in the well. In principle, using perturbation theory one may add third and fourth order vibrational derivatives around the well bottom. But this would lead to additional fitting parameters, we wanted to reduce the fitting to the bare minimum. We assumed symmetric motion about the barrier, that is we set V₃ = 0. Since the potentials in both systems studied have three fold symmetry, this too is not precise. Including V₃ would also imply an additional empirical parameter. The theory as used in this paper disregards the existence of a vibrational band structure, assuming a continuous energy from the ground state energy and upwards. To account for the band structure one would need further information, which is not at our disposal.

To summarize, the theory and results presented in this paper do indicate that the hopping of H and D on Pt(111) and H on Ru(0001) is not “strange”. Good agreement may be achieved between experiment and theory. This paper does challenge theory to show or to negate the observation that the barrier frequency seems to be much higher than the well frequency. The bottom line, the motion of H and D atoms on metallic surfaces will continue to intrigue us for a while to come.

Author contributions

Not relevant, the author is the sole author.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data used in this paper appears in the paper itself and is readily generated. The necessary Maple files are with the author and he would readily make them available to anyone who requests them.

Acknowledgements

I thank Professor Salvador Miret-Artés and Dr Elena Torres for many discussions and the Ben May center at the Weizmann Institute for its gracious support.

Notes and references

1. H. A. Kramers, Physica, 1940, 7, 284–304.
2. J. S. Langer, Ann. Phys., 1969, 54, 258–275.
3. D. Chandler, J. Chem. Phys., 1978, 68, 2959–2970.
4. P. Wolynes, Phys. Rev. Lett., 1981, 47, 968–971.
5. I. Affleck, Phys. Rev. Lett., 1981, 46, 388–391.
6. A. Caldeira and A. Leggett, Phys. A, 1983, 121, 374.
7. L. Chang and S. Chakravarty, Phys. Rev. B:Condens. Matter Mater. Phys., 1984, 29, 130.
8. H. Grabert, U. Weiss and P. Hänggi, Phys. Rev. Lett., 1984, 52, 2193.
9. S. Washburn, R. A. Webb, R. Voss and S. Faris, Phys. Rev. Lett., 1985, 54, 2712.
10. A. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison and J. Ellis, Prog. Surf. Sci., 2009, 84, 323–379.
11. A. P. Jardine, E. Y. M. Lee, D. J. Ward, G. Alexandrowicz, H. Hedgeland, W. Allison, J. Ellis and E. Pollak, Phys. Rev. Lett., 2010, 105, 136101.
12. H. Grabert and U. Weiss, Phys. Rev. Lett., 1985, 54, 1605–1608.
13. U. Weiss and H. Grabert, Phys. Lett. A, 1985, 108, 63–67.
14. A. Sanz, R. Martínez-Casado and S. Miret-Artés, Surf. Sci., 2013, 617, 229–232.
15. S. Upadhyayula and E. Pollak, J. Phys. Chem. Lett., 2023, 14, 9892–9899.
16. E. Pollak, J. Chem. Phys., 2023, 159, 224107.
17. E. Pollak, Phys. Rev. A, 2024, 109, 022242.
18. E. C. Kemble, Phys. Rev., 1935, 48, 549–561.
19. W. H. Miller, R. Hernandez, N. C. Handy, D. Jayatilaka and A. Willets, Chem. Phys. Lett., 1990, 172, 62–68.
20. M. J. Cohen, N. C. Handy, R. Hernandez and W. H. Miller, Chem. Phys. Lett., 1992, 192, 407–416.
21. R. Hernandez and W. H. Miller, Chem. Phys. Lett., 1993, 214, 129–136.
22. E. McIntosh, K. T. Wikfeldt, J. Ellis, A. Michaelidis and W. Allison, J. Phys. Chem. Lett., 2013, 4, 1565–1569.
23. C. Gardiner, Handbook of Stochastic Methods, Springer-Verlag, Berlin, 1999.
24. U. Weiss, Quantum dissipative systems, World Scientific Publishing, Singapore, 2nd edn, 1999.
25. J.-L. Liao and E. Pollak, Chem. Phys., 2001, 268, 295–313.
26. R. F. Grote and J. T. Hynes, J. Chem. Phys., 1980, 73, 2715–2732.
27. E. Hershkovitz and E. Pollak, J. Chem. Phys., 1997, 106, 7678–7699.
28. E. Pollak, Chem. Phys. Lett., 1986, 127, 178–182.
29. N. Fromän and P. Fromän, JWKB approximation, Contributions to the Theory, North-Holland, Amsterdam Amsterdam, 1965.
30. T. L. Nguyen, J. F. Stanton and J. R. Barker, Chem. Phys. Lett., 2010, 499, 9–15.
31. P. Goel and J. F. Stanton, J. Chem. Phys., 2018, 149, 134109.
32. S. Upadhyayula and E. Pollak, J. Phys. Chem. A, 2024, 128, 3434–3448.
33. S. Badescu, P. Salo, T. Ala-Nissila, S. Ying, K. Jacobi, Y. Wang, K. Bedürftig and G. Ertl, Phys. Rev. Lett., 2002, 88, 136101.
34. A. Jardine, E. Lee, D. Ward, G. Alexandrowicz, H. Hedgeland, W. Allison, J. Ellis and E. Pollak, Phys. Rev. Lett., 2010, 105, 136101.
35. Y. Georgievskii, M. A. Kozhushner and E. Pollak, J. Chem. Phys., 1995, 102, 6908–6918.
36. K. Kostov, W. Widdra and D. Menzel, Surf. Sci., 2004, 560, 130–144.

---