The characteristic function method in surface diffusion

Abstract

In surface diffusion, one of the key observables is the so-called intermediate scattering function which is measured directly from the helium spin echo surface technique. This function is shown to be a characteristic function of probability theory; in fact, it is the central dynamical measurement which clearly emphasizes that the relationship between the structure and dynamical process is a central goal in condensed matter physics. It is also the generating function of the moments and cumulants of the jump distribution function at any order; in particular, the second order which is related to the diffusion coefficient. Furthermore, the frequency moments or sum rules are also calculated from this scattering function and are written in terms of the total jumping rate. The role played by this function, when helium atoms are used for probing adsorbates (here we assume non-interacting adsorbates or very low surface coverages), is essential to clearly show some important differences with respect to neutron scattering. In order to illustrate this new approach, the method of the characteristic function, we have focused on the incoherent tunneling of H and D on a Pt(111) surface. This analysis has led us to some discrepancies with respect to the diffusion coefficient values previously reported. Finally, an extension to jumps beyond nearest neighbors is also considered theoretically.

---

1 Introduction

Diffusion of adsorbates on surfaces sampled by particles such as neutrons or He atoms is a very active field due to the fact that diffusion is one of the primary processes occurring on surfaces. Two well-established surface experimental techniques are widely used, the quasi-elastic He atom scattering (QHAS),^1,2 and neutron scattering (QENS)³ which have been complemented by using the He spin-echo technique (HeSE)^4,5 and neutron spin-echo (NSE).³ Unlike neutron scattering^6–10 where one speaks of incoherent and coherent scattering (depending on if space–time correlation functions are for the same adsorbate or two different adsorbates, respectively), helium scattering is fully coherent.¹¹ However, when non-interacting adsorbates or very low surface coverages are assumed, the difference between coherent and incoherent scattering disappears. An interesting aspect worth mentioning is that, from the very beginning, surface diffusion theory sampled by He atoms has been very much attached to the corresponding neutron scattering theory. The main reason is maybe due to the fact that this scattering was used originally to study diffusion processes and shown to be reduced essentially to a problem of non-equilibrium statistical mechanics¹² (within the linear response framework where the first-order perturbation theory is assumed).

In QHAS experiments, the so-called differential reflection coefficient [scr R, script letter R] (or the probability that the probing particles are scattered when interacting with adsorbates), forming a certain solid angle Ω with an energy exchange ħω and a momentum transfer parallel to surface, K, is proportional to the so-called dynamic structure factor (DSF) or density correlation function, S(K,ω), which is a response function, expressed as

where n_d is the diffusing particle concentration on the surface, F is the atomic form factor which depends on the mutual interaction between probe particles and N interacting adsorbates. Here, capital letters are often used to designate variables on the surface. The DSF also gives the spectrum of spontaneous fluctuations.

On the other hand, in HeSE experiments, the corresponding response function is the intermediate scattering function (ISF), I(K,t), which is the frequency Fourier transform of DSF written as follows

In neutron scattering,^6–10,12,13 theory is discussed in terms of number or particle and current densities. Thus, the ISF is the autocorrelation function of the particle density operator in the reciprocal space with ; R_j(t) being the position operator of the adsorbates as a function of time. Positions at two different times do not commute. Within the Born approximation, I(K,t) can be expressed in terms of the generalized pair-distribution function, also known as the G-function or the van Hove space–time correlation function. This function represents the probability that given a particle at the origin and at time t = 0, any particle, including the same one, can be found at the position R and at time t. The ISF is thus expressed in terms of G as

In neutron scattering theory,¹³ it is also customary to analyze the short- and long-time behaviors of time autocorrelation functions; in particular, the ISF. At short-times, the ISF can be Taylor expanded from eqn (2) relating the ISF to the frequency sum rules or frequency moments of S(K,ω) which have a one-to-one correspondence with the time derivatives of the ISF evaluated at t = 0 and are K-dependent quantities. Frequency moments higher than second order depend on the dynamics which implies a previous knowledge of the interaction potential. Moreover, the moments of the jump distribution of the diffusing particles are obtained from the incoherent ISF in a Taylor expansion of the momentum transfer evaluated at its zero value. On the other hand, the long-time behavior is associated with transport coefficients which are a manifestation of dissipation and relaxation processes and where hydrodynamic equations are solved. Fluctuations as well as momentum transfers have to be small enough in order to remain in the linear response framework.

In this work, our aim is to present the role played by the ISF in surface diffusion when He atoms are used as probe particles, being a coherent scattering. Non-interacting adsorbates or very low surface coverages are assumed in this work. Thus, in the next section, the ISF is shown to be a true characteristic function of probability theory^14,15 by analyzing its implications concerning moments and cumulants of the jump distribution as well as frequency sum rules. This approach, the so-called characteristic function method, allows us to exploit the main properties of this function which are well established in probability theory. In Section 3, in order to illustrate the different theoretical steps developed, we have focused on the incoherent tunneling of H and D on a Pt(111) surface where only jumps between nearest neighbor sites have been experimentally reported.¹⁶ By assuming a tight-binding model and solving a Pauli master equation, closed expressions for moments and cumulants at any order are easily extracted; in particular, the second moment which is related to the diffusion coefficient. Some discrepancies are found with respect to the values previously reported.¹⁶ An extension to jumps beyond nearest neighbors is also considered and discussed. Finally, in the last section, some conclusions and future work are presented.

2 Theory

In order to extract relevant physical information about the system of interest, several theoretical methods have been widely used. Namely, molecular dynamics calculations are generally used where a full description of the force fields (adsorbate–adsorbate and adsorbate–substrate interactions) involved is necessary, Langevin calculations by including or not memory effects and reduced density matrix calculations within the Caldeira–Leggett (CL) and Lindblad formalisms. Except in the molecular dynamics calculations, the surface is usually well represented by a thermal bath consisting of an infinite number of harmonic oscillators and, therefore, friction and noise (white or color) appear after integrating over the degrees of freedom of the bath. An equation of motion is then derived for the time evolution of the reduced density matrix which contains both frictional and thermal effects, the so-called CL master equation¹⁷ which is of Markovian type (high surface temperatures and small frictions). The corresponding diagonal matrix elements give the quantum probabilities and the off-diagonal elements, the so-called coherences. The same is found in the Lindblad formalism. The Lindblad master equation^18,19 is usually solved numerically through the so-called stochastic wave function (SWF) method.²⁰ In this context, the ISF can then be expressed as^21–24where ρ(R,t) represents the diagonal elements of the reduced density matrix. We have thus replaced the G-function in eqn (3) by the quantum probability. Interestingly enough, the second part of eqn (4) which expresses the average of the exponential function e^iK·R(t) is indicating that the ISF is a characteristic function (in this case, depending on time) coming from probability theory.^14,15,25 A characteristic function has several important properties: (i) I(0,t) = 1, (ii) I(K,t) ≤ I(0,t) = 1, (iii) I(K,t) is a uniformly continuous function of its arguments for all K values, (iv) all the moments of the position probability distribution are obtained by Taylor expansion of the exponential function and, therefore, one can speak of the momentum generating function, (v) the natural log of the ISF is the corresponding cumulant generating function, and (vi) the Fourier inversion formula exists determining the corresponding probability distribution. Furthermore, the properties of convergence, independent random variables and the sum of independent random variables are also fulfilled. Therefore, the ISF is the central dynamical measurement which clearly shows that the relationship between the structure and dynamical process is also a central goal in condensed matter physics.

On the other hand, from the continuity equation which is written as follows

and by applying the Fourier transform, one can have a very easy way to calculate the current density j(K,t) in terms of the characteristic functionor j(R,t) from the inverse Fourier transform. Furthermore, in the diffusive regime (when time is much greater than the inverse of the friction), the experimental ISF can be well fitted to an exponential function of time according to^4,5,9

I(K,t) = Be^−α(K)t + C, (7)

where B and C are constants for a given K (magnitude and direction of the momentum transfer), and α(K) is the dephasing rate which is known from experiment.

According to the jump diffusion model described by Chudley–Elliott (CE)⁹ for non-interacting adsorbates (that is, very low surface coverages) and if, for simplicity, a simple Bravais lattice is assumed as well as instantaneous jumps between different sites, the Pauli master equation can be written in terms of the diagonal elements of the reduced density matrix as follows

where the summation runs over all two-dimensional vectors L and Γ_L are the partial jumping rates. Thus, Γ_L⁻¹ represents the average time between successive jumps. It is also assumed that the time for a simple jump is very short compared to the time between successive jumps. The total jump rate is therefore , with Γ_L = Γ_−L. By taking advantage of the linearity property of the Fourier transform and considering the contributions by pairs coming from L and −L, one obtains the following differential equation in terms of the ISFand the solution is simply expressed asThe dephasing rate is thus written aswhere the jump probabilities are given by P_L = Γ_L/Γ. From the experimental dephasing rate obtained after eqn (11), the fitting parameters in this model are the partial jumping rates, Γ_L. Only the jump vectors L not orthogonal to the direction of K contribute to the dephasing rate and ISF.

If K gives us the direction of observation, which is chosen by the experiment, and K the modulus of this vector, the scalar product in eqn (4) can be written, in general, as K·R(t) = KL_‖(t) where L_‖(t) gives the projection of L along K (parallel to the direction of observation) as a function of time. If there are several jump/transition vectors L contributing to the total diffusion process (say m) then . Eqn (4) can then be written as

and where the nth-derivative of the characteristic function or momentum generating function with respect to K, at K = 0, provides us the nth-moment of the jump/transition distribution^14,15,25as well as the nth-cumulant of the same distribution according toIn particular, from the second moment, one can extract the tracer diffusion coefficient when the surface coverage is very low or non-interacting adsorbates are assumed.

Just as a reminder of the relationship between the moments and the cumulants, we list up to the third order for a random variable X(t)¹⁴

and so on. The third moment gives us information about the symmetry of the distribution. The so-called skewness, which is related to the third moment, is a measure of the symmetry of the probability distribution (if negative, the left tail is more pronounced than the right tail; if positive, the other way around).

Following Boon and Yip,¹³ and from eqn (2), the ISF can be written in terms of the frequency moments of S(K,ω) (frequency sum rules) by Taylor expansion with respect to time around t = 0

withand thereforewhich only depends on the projectile along the direction of observation through eqn (12). If the DSF is an even function of frequency, only even frequency moments will survive. In neutron scattering, frequency moments of order higher than three require a previous knowledge of the interaction potential.

3 Results and discussion

As an application of this new approach, we are going to consider the incoherent tunneling of H/D on a Pt(111) surface previously analyzed by means of the HeSE experimental surface technique.¹⁶ Diffusion dynamics were shown to take place in the moderate-to-high friction regime. Following the well-known CE model,⁹ deviations from nearest neighbor random jumps for H and D between fcc sites were reported to be minimal.¹⁶ The experiment covered a range of surface temperatures going from 250 K to 80 K and, therefore, thermal activation and tunneling regimes are coexisting. The crossover temperature was estimated to be 66 K for H and 63 K for D. Diffusion motion was also shown to correspond to nearest neighbor hopping for a low surface coverage of 0.1 ML, along the [112̄] direction with K = 0.86 Å⁻¹ and with a lattice length of a = 2.77 Å. The geometry of the problem has been clearly depicted and discussed in ref. 26 (Fig. 4 and the corresponding discussion). There are 3 high symmetry jump vectors for the Pt(111) surface, two are equivalent, say L₁ and L₂, forming an angle of β = π/6 with respect to [112̄] and, the third one, say L₃, is perpendicular. According to eqn (10) and (11), the last vector does not contribute to the ISF and the dephasing rate; that is, this symmetry direction has no contribution to the total diffusion. Moreover, there is no evidence for a dependency on the azimuthal direction after the experimental work. Our calculated Γ-values issued from the SWF method have been recently reported and showed fairly good agreement with the experimental values.²⁴

In a tight-binding description, we are dealing with the probability of an adsorbate to stay in a given potential well or surface site n at time t, P_n(t), and therefore eqn (4) should be written as a discrete Fourier series. Within the formalism of the master equation (the CE model can be considered as a special case),^15,25 if a simple Bravais lattice is assumed as well as instantaneous jumps between adjacent sites/wells along one symmetry direction (L_‖,1 or L_‖,2), on a periodic substrate, a Pauli master equation can be written in terms of such probabilities as^24,27,28

Ṗ_n(t) = Γ⁺_n−1P_n−1(t) + Γ⁻_n+1P_n+1(t) − (Γ⁺_n + Γ⁻_n)P_n(t), (19)

with Γ^±_n±1 being the tunneling/hopping transition rates from the (n ∓ 1)-th well to the n-th well and ± denotes if diffusion goes to the right (+) or to the left (−). Usually, these transition rates are fitting parameters or can be calculated from the transition state theory. If the initial condition is such that P_n(0) = δ_n0 and Γ = Γ⁺ + Γ⁻ describes the total rate (with Γ⁺_n = Γ⁻_n, Γ⁺_n = Γ⁺, and Γ⁻_n = Γ⁻), one obtains^24,28

P_n(t) = I_n(Γt)e^−Γt, (20)

where I_n(x) is the modified Bessel function of integer order n. Thus, P_n(t) gives us the probability to stay in the nth-well at time t. Then, the ISF can be written aswhere the 2-factor in the exponent comes from the two equivalent symmetric contributions (do not confuse the ISF which depends on K and t from the Bessel function of integer order with argument Γt). Now, after eqn (7) and (21), one has that the total transition rate can be expressed asThis rate is a function of the projection along the direction of K and also implicitly of surface temperature and friction coefficient. The total jumping rate, Γ, is usually plotted versus the inverse of the temperature in an Arrhenius-like plot.^16,24

It is now quite straightforward to calculate the first two moments and cumulants. Thus, one has that

where the diffusion coefficient D_K is expressed as

D_K = 2Γa² cos² β, (24)

with β = π/6. When the first moment is equal to zero, the three first moments and cumulants are equal. From the fourth order, they are no longer equal. Eqn (24) is different from the existing expressions in the literature; in particular, it differs from Weiss and Grabert²⁷ (D = a²Γ) by the 2-factor and the projection along the direction of observation due to their assumption of considering only a one-dimensional diffusion model. It also differs from the one used in the experimental work¹⁶ (D = a²Γ/4) by two important aspects, namely, the 1/4-factor and again the lack of the projection mentioned above. The experimental paper seems to have considered the three symmetry directions equivalent, disregarding the corresponding projections. Thus, their total and theoretical jumping rates Γ are multiplied by a factor of 3 (not 2). In order to consistently compare the diffusion coefficients issued from this work with the corresponding ones obtained by the experimental work, the experimental jumping rates have been chosen in each case and multiplied by a factor of 2/3.

Thus, our diffusion coefficient value along the direction of observation for the highest surface temperature of 230 K is 5.62 × 10⁻⁵ cm² s⁻¹ whereas the value obtained by using D = a²Γ/4 is 1.88 × 10⁻⁵ cm² s⁻¹. At this temperature, only activated diffusion takes place. On the contrary, at 140 K where both activated and tunneling diffusion regimes coexist,²⁴ our value is 1.02 × 10⁻⁵ cm² s⁻¹ versus the value of 3.40 × 10⁻⁶ cm² s⁻¹ reported previously.¹⁶ On the other hand, a previous experimental work by Graham et al. by using the helium atom scattering technique (without spin echo)²⁶ provides a value of around 5.0 × 10⁻⁶ cm² s⁻¹ at 140 K and around 5.0 × 10⁻⁵ cm² s⁻¹ at 250 K. This discrepancy with the older experimental values, apart from a slightly different surface temperature, is attributed to jumps beyond nearest neighbors which were considered. Moreover, these experimental values are issued from fitting a Lorentzian line shape to the quasi-elastic peak given by the DSF and extracting its full width at half maximum which is proportional to the diffusion coefficient.²⁶

If hops are not restricted to only nearest neighbors, for a simple Bravais lattice and one of the symmetry directions (say L₁ or L₂), the corresponding master equation is expressed now as

The ISF can then be written as a product of backward and forward diffusionwithwhere again a and β are the unit cell length along the symmetry direction considered and the angle formed by this direction and K, respectively. An alternative expression for the ISF can also be written as (from the CE model)where and P̄_n = Γ_n/Γ and a factor of 2 has been added in order to take into account the two equivalent directions. Notice the presence of the 2-factor with respect to eqn (21). From this expression, the moments and cumulants of the jump distribution are again straightforward calculated. Thus, we have thatwhere the diffusion coefficient D_K is expressed asUsually, the jump probabilities decrease exponentially with n, P̄_n = exp(−bn), and the summation over n could be replaced by an integral from 1 to infinity (resulting expression written in the second line). It should be stressed here the appearance of an extra 2-factor and the mean 〈n²〉 in eqn (30) with respect to eqn (24).

In older works,^29,30 the diffusion coefficient is written as

where . Remember that the Γ value and diffusion coefficient for this expression are valid for only one of the symmetry jump directions allowed by the corresponding projection on K; similar values are obtained for other symmetry directions but with different β-values.

Our next step is to calculate the first frequency moments (or frequency rules) when the ISF is given by eqn (21). Due to the fact that this function is an exponential function of time, the DSF is a Lorentzian function which is an even function of frequency and, therefore, only even frequency moments are allowed. Thus, after eqn (16), the ISF is expressed now as

withInterestingly enough, these frequency moments also depend on the total jump rate (a dynamical property). Similar expressions are easily obtained when the ISF is given by eqn (28). In this sense, surface diffusion analyzed by He atom scattering differs from other scattering experimental techniques such as neutrons and light. For example, for neutron scattering, the zeroth moment provides the static structure factor and, the second moment, the square of the product of the momentum transfer and thermal velocity.¹³

4 Conclusions

Through this work, we have shown that the ISF is a characteristic function of probability theory. In this new approach, the so-called characteristic function method has important consequences since by means of the moment and cumulant generating functions, one can access in a simple analytical way to a complete information for the jump distribution function of the adsorbate with time. In particular, the (tracer) diffusion coefficient is easily extracted from the second order moment and cumulant since both quantities coincide when the first moment is zero. Frequency moments or sum rule are also calculated from the ISF, being expressed in terms of the total jumping rates. In particular, the D_K-values along the K-direction are obtained from the second moment of the ISF instead of the full width at half maximum of the Lorentzian line shape for the quasi-elastic peak (the DSF). This analysis has also allowed us to provide values of the diffusion coefficient which are greater than previously reported.

The extension of this approach to surface diffusion at higher coverages should be straightforward. An important issue is also how to evaluate the so-called collective diffusion coefficient which describes the density fluctuation relaxation of the adsorbates on the surface.³¹ In this context, the ISF can be calculated from molecular dynamics or by solving the CL or Lindblad master equations. The second moment should be now related to the collective diffusion constant. In the diffusive regime, the ISF still displays an exponential function of time, eqn (7), where the dephasing rate has also to be expressed in terms of the surface coverage but it should be considered as a parameter of the theory. From this knowledge, de Gennes narrowing has been reported in surface diffusion analyzed by He atom scattering.^4,5 Work in this direction is now in progress.

Author contributions

All authors have contributed equally to this work.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data supporting the findings of this study are included within the article.

Acknowledgements

E. E. T.-M. and S. M.-A. acknowledge support of a grant from the Ministry of Science, Innovation and Universities with Ref. PID2023-149406NB-I00.

Notes and references

1. F. Hofmann and J. P. Toennies, Chem. Rev., 1996, 96, 1307–1326.
2. A. P. Graham, Surf. Sci. Rep., 2003, 49, 115–168.
3. I. Calvo-Almazán and P. Fouquet, Eur. Phys. J.-Spec. Top., 2012, 213, 149–163.
4. A. P. Jardine, G. Alexandrowicz, H. Hedgeland, W. Allison and J. Ellis, Phys. Chem. Chem. Phys., 2009, 11, 3355–3374.
5. A. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison and J. Ellis, Prog. Surf. Sci., 2009, 84, 323–379.
6. L. V. Hove, Phys. Rev., 1954, 95, 249.
7. G. Vineyard, Phys. Rev., 1958, 110, 999.
8. P. G. De Gennes, Physica, 1959, 25, 825–839.
9. C. T. Chudley and R. J. Elliott, Proc. Phys. Soc., 1961, 77, 353.
10. G. R. Kneller, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 9450–9455.
11. D. J. Ward, A. Raghavan, A. Tamtögl, A. P. Jardine, E. Bahn, J. Ellis, S. Miret-Artès and W. Allison, Phys. Chem. Chem. Phys., 2021, 23, 7799–7805.
12. S. Lovesey, Theory of neutron scattering from condensed matter, Clarendon Press, 1984, vol. 2.
13. J. P. Boon and S. Yip, Molecular hydrodynamics, Courier Corporation, 1991.
14. C. W. Gardiner et al., Handbook of stochastic methods, Springer, Berlin, 2004, vol. 3.
15. G. Haag, Modelling with the Master equation, Springer, 2017.
16. A. Jardine, E. Lee, D. Ward, G. Alexandrowicz, H. Hedgeland, W. Allison, J. Ellis and E. Pollak, Phys. Rev. Lett., 2010, 105, 136101.
17. A. O. Caldeira and A. J. Leggett, Phys. A, 1983, 121, 587.
18. G. Lindblad, Commun. Math. Phys., 1976, 48, 119–130.
19. A. Kossakowski, Rep. Math. Phys., 1972, 3, 247.
20. H. Breuer and F. Petruccione, The theory of open quantum systems, Oxford University Press on Demand, 2002.
21. E. Torres-Miyares, G. Rojas-Lorenzo, J. Rubayo-Soneira and S. Miret-Artés, Mathematics, 2021, 9, 362.
22. E. Torres-Miyares, G. Rojas-Lorenzo, J. Rubayo-Soneira and S. Miret-Artés, Phys. Chem. Chem. Phys., 2022, 24, 15871–15890.
23. E. Torres-Miyares, D. Ward, G. Rojas-Lorenzo, J. Rubayo-Soneira, W. Allison and S. Miret-Artés, Phys. Chem. Chem. Phys., 2023, 25, 6225–6231.
24. E. Torres-Miyares and S. Miret-Artés, Phys. Chem. Chem. Phys., 2025, 27, 14739–14743.
25. U. Weiss, Quantum Dissipative Systems, World Scientific Publishing, Singapore, 1999.
26. A. P. Graham, A. Menzel and J. P. Toennies, J. Chem. Phys., 1999, 111, 1676–1685.
27. U. Weiss and H. Grabert, Phys. Lett. A, 1985, 108, 63–67.
28. A. Sanz, R. Martnez-Casado and S. Miret-Artés, Surf. Sci., 2013, 617, 229–232.
29. J. Frenken and B. Hinch, Helium Atom Scattering from Surfaces, Springer, 1992, pp. 287–313.
30. A. Graham, A. Menzel and J. Toennies, J. Chem. Phys., 1999, 111, 1676–1685.
31. R. Gomer, Rep. Prog. Phys., 1990, 53, 917.

---