E. E. Torres-Miyares*ab and
S. Miret-Artés
b
aFundación Humanismo y Ciencia, Guzmán el Bueno 66, 28015 Madrid, Spain. E-mail: elena.torres@iff.csic.es
bInstituto de Física Fundamental, Consejo Superior de Investigaciones Científicas, Serrano 123, 28006 Madrid, Spain. E-mail: s.miret@iff.csic.es
First published on 26th June 2025
In this work, incoherent tunneling surface diffusion is analyzed by means of the stochastic wave function (SWF) method within the Lindblad formalism. Previous calculations based on the transition state theory (TST) were unable to provide reasonable values of friction coefficients for the diffusion motion of H and D adsorbates on a Pt(111) surface. Thermal activation and tunneling regimes were covered by varying the surface temperature, approaching the so-called cross-over temperature. Numerical results for the total hopping/tunneling rates are fairly good when compared to the experimental values by assuming a simple cosine corrugation function. No fitting has been carried out and estimated values of the physical parameters (friction coefficients, well and barrier frequencies and barrier height), extracted from the thermal activation regime, are used. Finally, survival probabilities on the different wells are also reported.
I(K, t) = Be−α(K)t + C, | (1) |
The corresponding measurements indicated that quantum effects are significant at low temperatures. The range of temperatures used went from 250 K up to 80 K, covering thermal activation and tunneling regimes. The crossover temperature was estimated to be 66 K for H and 63 K for D. Furthermore, the diffusion motion seems to correspond to nearest neighbor hopping for a coverage of 0.1 ML, along the [11] direction. Following this CE model, deviations from nearest neighbor random jumps for H and D between fcc hollow sites were reported to be minimal. Diffusion dynamics was also claimed to take place in the moderate-to-high friction regime. The quantum TST of dissipative tunneling for temperatures above the crossover covering tunneling and thermal activation should be applied. From this theory and by assuming a parabolic barrier, one should easily extract information on the corresponding physical parameters. Finally, despite a very good fitting to the experimental jumping/tunneling rates versus the inverse of the surface temperature, the friction values reported were unreasonable large (of the order of 1000 ps−1) and, therefore, the energy barrier can not be the same to the one we have used here. Another attempt to obtain friction values was based on path integrals for a periodic potential.3 The values reported were too small (of the order of a few ps−1) but no information on barriers and frequencies was provided.4 In both cases, the authors obtained very good fittings but with big numerical discrepancies, showing that good fittings are not enough.
The tunneling rates were calculated from closed expressions using Kramers' theory.5 It is generally assumed that tunneling proceeds via a parabolic barrier.6 Within the extended TST, the tunneling rate can be expressed as
![]() | (2) |
![]() | (3) |
Very often, parabolic barriers are not good enough and anharmonic effects are needed in order to extract reliable fitting physical parameters. Moreover, in the vicinity of the crossover temperature, some divergences are found in the theory; the Wolynes factor diverges at the crossover temperature defined as ħβcλ‡ = 2π, where βc = 1/kTc, k and Tc being the Boltzmann constant and the crossover temperature, respectively. This divergence has been recently removed using the uniform semiclassical energy-dependent transmission coefficient. The resulting one-dimensional theory11 has been generalized to dissipative systems12 and may be used to analyze the experimental data. The uniform Wolynes factor is clearly responsible for the change of slope around the crossover surface temperature. However, the bending of the straight line experienced by the tunneling rates around Tc is a real open problem. One of the critical problems in the standard and extended theory used comes from a fundamental assumption. In order to provide closed formulas, and considering Ohmic friction, one has to assume the separable motion between the so-called unstable normal mode and the stable modes along the reaction or symmetry coordinate; in other words, their coupling should be very weak. This assumption could not take place or be in the limit of applicability of the theory in ref. 1. This makes far from trivial the theoretical treatment in order to reach workable expressions for carrying out a fitting procedure. Therefore, in general, several problems should be stressed for this purpose. First, one is faced to fit at the same time at least 4 parameters (friction coefficient, barrier height, well and barrier frequencies and, sometimes, anharmonic corrections to a parabolic barrier). The dimensionality of the space of parameters is thus too high; the fitting procedure is carried out in a global way, not sequential. Second, the different conditions that the theory demands in order to be applicable to reach closed expressions have to be added as constraints; in particular, as mentioned before, if the motion between the stable and unstable modes can be assumed separable. And, third, even if the divergence of the Wolynes factor has been recently removed, it is not clear that a better theory is coming soon; several conditions and/or constraints have also to be fulfilled in the fitting process. In the meantime, one way to skip such difficulties is, in our opinion, the procedure we use in this work. Our starting point is thus different and based on the so-called SWF method.13–15
The ISF, I(K, t), is related to the probability density ρ(R, t) through a space Fourier transforms as follows13–15
![]() | (4) |
![]() | (5) |
A single Lindblad operator Âk is used, given by a linear combination of the position operator and the momentum operator
. In the coordinate representation, this equation leads to an Î to stochastic differential equation which is solved for a high number of stochastic wave functions or realizations. The ISF can then be obtained through an average over stochastic realizations
I(K, t) = 〈e−iK·![]() ![]() | (6) |
A potential with a single Fourier component is assumed: V(x) = V0cos(2πx/a), with an energy barrier V‡ = 2V0 and a lattice length of a = 2.77 Å for the Pt(111) surface. The numerical code is parallelized and the number of realizations Nr = 100
000 is chosen to reach the numerical stability for the mean value expressed in eqn (6). The numerical simulations are run up to 40 ps with a time step of 0.05 ps, by starting with an initial Gaussian wave packet centered at x = 0, with a width of 0.79 Å and momentum distributed according to the Maxwell–Boltzmann distribution
. A one-dimensional space within the interval [−4.2, 4.2] Å is used with Ns = 4096 points, resulting in a space step of Δx = 0.002 Å. At long times, t ≫ γ−1, the ISF is expected to decay as eqn (1) and this allows us to extract the dephasing rate for each value of K. The cut-off time for fitting the tail of the ISF to an exponential function of time is chosen to be 8.2 ps, following the same procedure as in the experimental case; that is, when a certain region of stability of the dephasing rate is observed, which has also to be of the order of the experimental rate value.
Within the CE jump model, if a simple Bravais lattice is assumed as well as instantaneous jumps between different sites, and only transitions between first neighbors are going to be allowed (the diffusion dynamics is assumed to be one-dimensional on a periodic substrate), a Pauli master equation can be written in terms of probabilities as4
Ṗn(t) = Γ+n−1P(t)n−1 + Γ−n+1P(t)n+1 − (Γ+n + Γ−n)Pn(t), | (7) |
Pn(t) = In(Γt)e−Γt, | (8) |
![]() | (9) |
![]() | (10) |
![]() | (11) |
This analytical expression is quite different from the one used in ref. 4, as they arise from different theoretical frameworks.
In Fig. 1, Arrhenius-like plots of the temperature dependence of the tunneling/hopping rates issued from the experimental results1 for H (red circles) and D (red triangles) diffusion and momentum transfer of K = 0.86 Å−1 are plotted. Numerical simulations are also included coming from the SWF method (NS blue points). The experimental total hopping rates have been supposed to have three directions of diffusion but, according to the direction of observation chosen, only two are possible.9 Thus, these ones have been multiplied by a factor 2/3 in order to take into account that in this case N = 2. As can be seen, the agreement is fairly good for all surface temperatures; in particular, at low temperatures close to the cross-over one. The hopping rates for D are consistently smaller than those for H across the entire temperature range, indicating a clear mass effect. For H diffusion, if one assumes that the Arrhenius classical law is due to thermal activation only, one then could compare the corresponding hopping rate to the experimental/numerical one at different surface temperatures. Thus, for 214 K, hops by thermal activation dominate 100–90%; for 80 K, only hops by tunneling predominates, around 90%; and for 121 K, hops for thermal activation and tunneling are around 30% and 70%, respectively.
![]() | ||
Fig. 1 Arrhenius-like plots of the temperature dependence of the tunneling rate issued from the experimental results shown in ref. 1 for H diffusion (red circles) and D (red triangles) diffusion at a momentum transfer of K= 0.86 Å−1. Numerical simulations (NS blue points) using a energy barrier of V = 72 meV and a friction coefficient of 140 ps−1 for H (blue square) and 70 ps−1 for D (blue star) are also plotted and issued from the SWF method (assuming the same friction, barrier height, and vibrational frequencies across all temperatures). |
Finally, another useful information which can be extracted from the experimental results for the total hopping/tunneling rates is the probability to stay in a given surface site (survival probability). This information has not been given so far in this context. From eqn (8), this probability as a function of time is easily calculated once the tunneling/hopping rates are known. In Fig. 2, the time dependence of the probability of staying at the zero, P0(t) (blue curve), first, P1(t) (black curve), and second, P2(t) (red curve), well are plotted for H adsorbates on a Pt(111) surface at three different temperatures, covering the classical and quantum regimes: 214 K, 121 K and 80 K. These calculations are carried out using the experimental tunneling rates shown in Fig. 1. As clearly seen in this figure, at 214 K and 121 K, P0(t) and P1(t) suddenly start with a maximum probability and decreasing with time very rapidly due to the thermal activation (classical regime). Pn(t) functions could also be seen as survival probabilities to remain in a given well. Obviously, P2(t) displays the same pattern but with a smaller maximum value. The third temperature displays a different time behavior since we are in the quantum regime. Here it is supposed that tunneling plays a major role and the sudden increase at very short times is not longer observed. Furthermore, the decreasing of P1(t) and P2(t) with time is not so rapid because the tunneling diffusion takes time since Γ is small. The asymptotic values observed for the survival probabilities indicate us that at a sufficiently long time these probabilities are nearly constant where thermal equilibrium is already reached. At these long times hopping/transition dynamics could repopulate surface sites such as 0, 1, 2, etc.
![]() | ||
Fig. 2 Time dependence of the probability of staying (survival probabilities) at the initial well P0(t) (blue curve), the first P1(t) (black curve), and the second P2(t) (red curves) well from eqn (8) for H adsorbates on a Pt(111) surface is plotted at three different temperatures: 214 K, 121 K and 80 K. These calculations are carried out using the corresponding experimental hopping/tunneling rates obtained for a momentum transfer of K = 0.86 Å−1 along the [11![]() |
Along this work, we have focused on a theoretical challenge, the incoherent tunneling diffusion. The choice of H and D diffusion by a Pt(111) surface in order to illustrate this transition is by no means arbitrary. Unreasonable friction values were reported as well as the standard and extended TST for this case is now questionable. The SWF method within the Lindblad formalism is able to provide reasonable results for the tunneling/hopping rates as a function of the inverse of surface temperature by assuming a simple cosine corrugation function. Physical parameters for friction, well and barrier frequencies as well as barrier height are estimated from the diffusion dynamics at high surface temperatures. Clearly, this numerical method is not convenient at all for a fitting procedure but, waiting for a better analytical theory within the TST formalism, it has shown to be a very good alternative, apart from being an exact theory within the Markovian approach. Obviously, ab initio calculations need to be done for a better description of this tunneling diffusion. ref. 16 should be a good starting point. For other systems, with a different adsorbate or surface symmetry, may require going beyond a 1D model, and the hopping rates cannot be captured by a simple scaling factor. Furthermore, it is also very illustrative to show survival probabilities for different wells which decay exponentially with time, exp(−Γt).
This journal is © the Owner Societies 2025 |