Jump rate and jump probabilities in the two-dimensional strong-collision model

The diffusion of a particle in a two-dimensional non-separable periodic potential is studied in the framework of the strong-collision model. In this model, the diffusing particle is subjected to well-separated strong collisions with frequency η. After each collision the particle velocity is suddenly thermalized. The model is studied by numerical simulations. The average jump rate, the probability distribution of the jump lengths and the diffusion coefficient D are calculated depending on η and on the x–y coupling term of the potential. For large η, the x–y coupling term decreases the jump rate, while for small η the opposite trend is found. However, the x–y coupling always suppresses a large proportion of long jumps, causing thus a decrease of the diffusion coefficient for all values of η. In the presence of x–y coupling, in the limit η ≪ 1, there is an indication that D behaves as η−α, with α < 1.


Introduction
In the last two decades, several experimental and theoretical works showed that the diffusion of single adatoms on crystal surfaces [1] may proceed quasi-continuously [2,3] or by jumps [4]- [6]. Activated or jump diffusion is recovered at high potential barriers (i.e., at relatively low temperatures), where the adatom spends most part of its time in oscillating at the bottom of the adsorption sites, rarely escaping from them to perform single or multiple jumps. In this case, the diffusion problem is strictly connected to the generalization of the classical Kramers escape problem [7] to periodic potentials. It must be stressed that the most important requirement for a proper definition of the jump rate r j is that the inverse rate is the longest time scale in the system. Numerical calculations [8] have shown that this happens safely for barriers E a 4k B T .
Activated surface diffusion has now been observed in many experiments on different systems [9]- [13], to cite only a few. Theoretically, jump diffusion can be treated by different techniques or models, from stochastic equations with an external periodic field of force, derived from an adiabatic periodic potential, to molecular dynamics, to path integration and to chaotic Hamiltonian dynamics [6], [14]- [20], [22]- [29].
The LE of the classical Brownian motion, where the adatom is dynamically coupled to the thermal bath via a friction force −mηv and a white noise, is derived from the GLE if the velocity v and the position x of the adatom vary on a slower time scale with respect to the correlation time of the fluctuations. The assumptions underlying the LE are valid for the interactions between the adatom and the electrons in the substrate. On the other hand, the time scale separation is more questionable for phonons; in the continuum Debye model [35], the friction force is written as −mηv if the frequencies associated to the motion of the adsorbed atom are well below the highest phonon frequency. It seems, then, reasonable to assume that the LE is valid when the mass of the adatom is large with respect to the mass of the substrate atoms. From a slightly different point of view, x and v are 'slow' variables if the adatom-surface energy exchange satisfies the condition δE k B T , which means 'weak and frequent collisions'. In these conditions, the particle performs a Brownian motion: its velocity is slightly modified by a single collision and a large number of collisions occurs for the thermalization of the adatom in a well. This requirement should again be valid for heavy adatoms. Generally speaking, and especially for light adatoms, the GLE with an appropriate form of the memory kernel should be used [21]; however, the evaluation of the memory kernel, which depends on the physical system, can be cumbersome in practical cases.
A possible way to overcome this difficulty is to try a different approach. In fact, if the particle exchanges energy with the substrate by means of strong and well separated collisions with the phonons, its motion might be described in the framework of the Boltzmann equation [22]. The latter is a very complicated nonlinear equation whose solution is cumbersome. However, if some assumptions are made on the velocity distribution of the particles after a collision, the Boltzmann equation can be noticeably simplified by linearizing the collision kernel. In this framework, the BGK equation is a well-known approximation [31], which assumes that the velocity of the particle is distributed according to the Maxwell law after each collision, due to sudden thermalization. Therefore, this model can be also considered a strong-collision model. This approximation seems to be particularly appropriate when the diffusing particle and the bath particles have the same mass or, in any case, for light adatoms. Both the Langevin and the strong-collision models tend to the same limit for η → ∞, where the parameter η is the collision frequency in the strongcollision case and the friction in the Langevin case [36]. In the one-dimensional (1D) case, strong-collision models have been studied by several authors [37]- [39] also in the case of non-Poissonian statistics of the events. The latter allows the explicit introduction of memory effects in the framework of the strong-collision model.
The above cited stochastic equations have been exactly solved by a computational method, the matrix-continued-fraction method (MCFM) [6,14,15], [17]- [19], due to Risken [36]. The solution has been given in the 1D case, except for the Langevin model which has been solved also in a square 2D periodic potential [15]. Results for hexagonal and honeycomb lattices have been obtained by direct simulation [16]. In the 1D case, both Langevin and strong-collision models present the classical Kramers turnover behaviour depending on η [6,15,17]. Multiple jumps are allowed only below the turnover point. In the case of the GLE with exponential memory kernel, the jump rate again presents a turnover when plotted as a function of the memory decay parameter and once again long jumps are recovered only below this memory turnover [19]. Another important result is that, in the strong-collision model, long jumps are favoured with respect to the Langevin model. In fact, strong but separated collisions are less efficient in thermalizing the velocity than a continuous frictional dissipation [17].
The MCFM for the strong-collision model is a rather cumbersome method already in 1D, requiring considerable computing time at high barriers and at low η. The method becomes impracticable in 2D. In this case, the direct simulation is much more convenient. The latter resembles a classical constant-temperature molecular-dynamics simulation.
In this paper the jump rate, the probability distribution of the jump lengths and the diffusion coefficient are calculated in the 2D strong-collision model by molecular dynamics, in order to study dimensional effects as done in the 2D Langevin case by the MCFM. The same egg-carton periodic potential [15] is employed to model the adiabatic interaction between the adatom and the crystal surface. This type of potential represents the simplest and significant 2D potential to investigate the effects of x-y coupling. The same potential has also been used in theoretical studies about diffusion in superionic conductors [40], about nonlinear conservative dynamics of a classical particle [41] and in chaotic transport driven by ac forces [42].

Model and method
In the strong-collision model, a particle suffers strong and isolated collisions with frequency η. On each collision, the velocity of the particle is changed, and extracted from a Maxwell-Boltzmann distribution at the given temperature T . In our case, the particle is also subjected to a deterministic force F (x, y) which varies periodically in space and is given by a is the spacing of the square lattice. The amplitudes A 0 and A 1 are both positive and A 0 A 1 . When A 1 = 0, the x and y coordinates are decoupled, and the system behaves as two independent 1D systems. The energy barrier between neighbour cells is given by E a = 2(A 0 − A 1 ). When A 1 > 0, the saddles are narrower than the minima. This can be seen by calculating the curvatures at the minimum ω m and at the saddle point ω s in the direction across the diffusion path: Note that at the minimum, the oscillation frequencies along and across the diffusion path coincide.
To explore the effects of the x-y coupling, it is then useful to introduce the parameter ξ, defined as the ratio of ω m and ω s In the following, we consider a fixed E a /k B T = 6, and vary ξ and the normalized collision frequency γ In our simulations, the equations of motion for the diffusing particle are solved by the velocity Verlet algorithm [43]; the time step δt is chosen as being much smaller than the smalloscillation period at the well bottom and than 1/η. The particle must suffer collisions with frequency η. To this end, at each step we calculate the quantity η δt 1 and extract a random number r (0 r 1). If r is smaller than η δt the particle suffers a collision and its velocity is extracted from the Maxwell distribution at the desired temperature T . Otherwise, the particle continues its deterministic motion following Newton's equations.
Our results are restricted to the analysis of the jump-diffusion regime. As discussed in section 1, this regime takes place when the activation barrier is sufficiently high, since jumps can be defined only when a clear time-scale separation between the oscillation period in a well and the average residence time in a given cell is achieved. In order to quantify jump diffusion, a criterion for counting jumps, distinguishing single from long jumps, is necessary. In a long jump, a particle starts from a given cell and thermalizes finally in a distant cell, without stopping in the intermediate cells. How do we determine quantitatively whether a particle has stopped in an intermediate cell or not? This is a point which deserves some attention. We have counted that a particle is thermalized in a given cell if it has spent there a time interval (without interruptions) which is longer than any characteristic time of the system. In this case, the characteristic times include the oscillation period, the inverse collision frequency and the prefactor of Kramers' highfriction formula. The latter amounts to 2πη/ω 2 m , and needs to be considered in order to count the correct number of jumps at high η. For each value of the triplet (E a /(k B T ), γ, ξ), we calculate the diffusion coefficients D xx and D yy as follows. We calculate the rates r x j and r y j , which are associated to jumps with displacements along x and y, respectively. In this way, a jump with displacement along both x and y contributes to both r x j and r y j ; however, jumps of this kind are very rare in square geometry, since most jumps are either along x only or along y only. The total jump rate is thus given by r j = r x j + r y j . After having associated a length to each jump, we calculate the probability distributions of the jump lengths p n x and p n y . The diffusion coefficient D xx is then evaluated as where n 2 x is the mean square jump length. The expression for D yy is analogous. In this geometry, D xx = D yy = D, and thus D is numerically evaluated as D = (D xx + D yy )/2.

Results
As written before, we consider a fixed barrier E a = 6k B T , which is sufficient to ensure that the jump diffusion regime is fully achieved, and study the behaviour of the jump rate r j , of the probability of single jumps p 1 and of the diffusion coefficient D as functions of the normalized collision frequency γ, for four different values (1, 0.5, 0.25 and 0.1) of the ratio ξ between the width of the saddle and the width of the minimum.
In figure 1, the results for the jump rate r j are reported. r j shows the usual turnover behaviour, resembling the case of the LE, that was already found in the 1D case [17]. However, the important point in the 2D case is the effect of the x-y coupling on r j . For high γ, the decrease of ξ (corresponding to the narrowing of the saddle with respect to the well bottom) causes the decrease of r j . Considering that for large η, the Langevin and the strong-collision models tend to be equivalent (even though some numerical differences are present at any finite value of η [38]), this result agrees qualitatively with the prediction of the multi-dimensional high-friction Kramers' rate theory [30]. We have indeed verified that r j for large η is well approximated by the following Kramers' expression in the absence of x-y coupling where we have taken into account the existence of four equivalent saddles to escape from a given well. When the x-y coupling is present, Kramers' expression predicts that the rate is reduced, being simply multiplied by ξ. The rate r j is indeed reduced, but less than the value predicted by Kramers' formula. This could be due to the fact that Kramers' formula does not incorporate anharmonic corrections, and also to the fact that at any finite η, the Langevin and the strongcollision models do not coincide. For small collision frequency, the effect of narrowing the saddle has the opposite effect, since it causes an increase of the jump rate. This result shows a further analogy between the strong-collision and the Langevin models, even though in this limit the two models are not perfectly equivalent as in the case η → ∞. In fact, Hershkovitz and Pollak [44], showed that in the Langevin model, an increase of the coupling between degrees of freedom enhances the escape rate at low friction. The presence of chaotic motion in the system enhances the rate of energy diffusion (energy diffusion dominates the escape rate when friction is low). Therefore, Hershkovitz and Pollak concluded that more the chaos, the faster the energy diffusion and the larger the rate. An increase of escape rate at increasing x-y coupling was found also by Vega et al [45]. However, their case was qualitatively different from ours, because they considered a potential where saddle points become wider at increasing coupling, so that the crossing of the saddle point becomes easier for any value of η. The effects of the x-y coupling on the jump rate are thus quite evident. However, from the quantitative point of view, the effects on the probability of making either single or long jumps are much stronger, as can be seen in figures 2 and 3. In figure 2, the probability of making single jumps p 1 is reported as a function of γ for fixed E a /(k B T ) and for different values of ξ. A decrease of ξ from 1 to 0.5, corresponding to a narrowing of the saddle by a factor of two, strongly enhances p 1 : for γ = 0.1, p 1 changes from 0.2 to 0.4, and for γ = 0.01, the change is from 0.015 to 0.25. The probability of making long jumps p lj = 1 − p 1 is correspondingly reduced, and the effect is more and more evident for decreasing collision frequency. This can be better seen from figure 3, where the jump length probability distribution is reported for the  case of γ = 0.01. In the absence of x-y coupling, the probability p n of making a jump of n cells decays exponentially with n, even in the case of small γ. This result was already found in [17], where an accurate analytical approximation to p n was also derived. The introduction of the x-y coupling changes the behaviour of p n in a qualitative way. In fact, the decay is faster for small n, and becomes exponential only in the large n limit. In this respect, the strong-collision and the Langevin models behave in a different way, because deviations form the exponential decay of p n are present in the Langevin model for low friction also in the absence of x-y coupling [6]. These deviations are of the same kind as those found now in the strong-collision model with coupling. The considerable suppression of long jumps has clear consequences on the diffusion coefficient, as reported in figure 4. Since the diffusion coefficient is proportional to the mean square jump length n 2 x , and the latter is dominated by very long jumps, we expect that the x-y coupling strongly reduces the diffusion coefficient when γ is small. Let us consider the small-γ case of figure 3. For ξ = 0.1 (stars in the figure) jumps of more than 20 cells are practically negligible. In contrast, for ξ = 1 (circles in the figure) jumps much longer than 30 cells (up to more than 100 cells indeed) are still likely. Correspondingly, the diffusion coefficient is larger by a factor 10 2 for ξ = 1 than for ξ = 0.1. However, in analogy with the case of the Langevin model [15,16,46,47], the effect of the x-y coupling on D is not only quantitative, but also qualitative. From the lower panel of figure 4, it turns out that in the absence of coupling, D diverges as γ −1 in the limit γ → 0. The introduction of the x-y coupling changes this behaviour (at least in the range down to γ = 0.01) to γ −α , with α < 1. From the data in figure 4, one has α = 0.74, 0.51 and 0.32 for ξ = 0.5, 0.25 and 0.1, respectively.

Conclusions
In summary, the analysis of the strong-collision model in a 2D periodic potential has shown that the presence of an x-y coupling term has strong effects on the diffusive motion of the particle. In our analysis, we have considered a fixed activation barrier E a = 6k B T , which has been chosen in order to obtain a jump-diffusion regime, and several values of the ratio ξ between the width of the saddle points and the width of the minima. For each value of ξ and E a , we have varied the collision frequency η in a wide range of values. Our results have shown that the x-y coupling reduces the jump rate at large η, because the saddles become narrower than the minima. In contrast, for small η, the jump rate is increased by the x-y coupling term, because the latter allows the energy exchange between the two degrees of freedom. However, the x-y coupling term always suppresses a large proportion of long jumps. This effect is again due to the narrowing of the saddles. For small η, the decrease of the long-jump probability overcompensates the increase of the jump rate, so that the overall effect is a strong decrease of the diffusion coefficient. In this limit, there is an indication of an anomalous dependence of the diffusion coefficient on the jump frequency η. In fact, while in 1D D diverges as η −1 , in coupled 2D the behaviour is of the kind η −α , with α < 1. The results of the strong collision model show several analogies with those of the Langevin model, especially when the coupling is present. In fact, the main qualitative difference in the results of the two models in 1D is eliminated by the x-y coupling. In 1D, the decay of the jump probabilities with distance is perfectly exponential for the strong-collision model, while the Langevin model gives an exponential decay only asymptotically for large jump lengths. In coupled 2D, the behaviour becomes generally non-exponential for the strong-collision model also, and the deviations from the exponential decay are qualitatively of the same kind as in the Langevin model.