Adatom interaction effects in surface diffusion

Motivated by recent research of Nikitin et al. (J.Phys.D vol. 49, 055301(2009)), we examine the effects of interatomic interactions on adatom surface diffusion. By using a mean-field approach in the random walk problem, we derive a nonlinear diffusion equation and analyze its solutions. The results of our analysis are in good agreement with direct numerical simulations of the corresponding discrete model. It is shown that by analyzing a time dependence of adatom concentration profiles one can estimate the type and strength of interatomic interactions.


I. INTRODUCTION
Diffusion is ubiquitous in Nature. It determines the behavior and controls the efficiency of many biological and technological processes. Examples include the wetting, conductivity of biological membranes, catalysis, growth of crystals, sintering, soldering etc. Surface diffusion is particularly important in nano-technological processes which are aimed to obtain objects of submicron sizes where the surface properties are of the same importance as the bulk ones.
Macroscopic description of diffusion is based on Fick's law, which postulates proportionality between the particle flux and the concentration gradient. Establishing a link between macroscopic laws of diffusion and microscopic non-equilibrium density matrix approach is one of the most challenging and important problems of non-equilibrium statistical mechanics (for reviews see, e.g. [1][2][3][4][5][6]). Surface diffusion is essentially many-particle process. Even at very low coverage when the interaction between adatoms is negligible, the random walk of an isolated adsorbed particle is a collective motion due its interaction with substrate atoms [4,5]. The walker moves in a potential landscape which is changed by the walker [7]. The walk on a deformable medium when the walker leaves behind a trail and due to slow relaxation the trail affects the next walker, is also a collective process [8,9]. At finite coverage particles at surfaces, in addition to the interaction with a substrate, experience lateral interactions of different origin: the attractive van der Waals, oscillating electronic exchange and multipole-multipole electrostatic interactions. The dipole-dipole interaction which is due to a polar (mainly dipolar) character of adsorption bonds is long-ranged (as r −3 ) and generally repulsive since all dipole moments are essentially parallel (see review papers, e.g. [4,5]).
Bowker and King [10] used Monte-Carlo simulations in order to clarify the effect of lateral interactions of adatoms on the shape of evolving concentration profiles in surface diffusion. They showed that the intersection point of the diffusion profiles with the initial stepwise profile lies above θ max /2 in the case of lateral repulsion and below θ max /2 in the case of attraction (θ max is the maximum concentration in the initial step).
In a quite recent paper [11] an approach based on error function expansion was proposed to fit experimental concentration profiles. This algorithm provides a high-accuracy fitting and allows extracting the concentration dependence of diffusivity from experimental data. The goal of our paper is to model and examine the effects of interatomic interactions on adatom surface diffusion. Starting with nonlinear random walk equations where the interatomic interactions are considered in the mean-field approach, we derive a nonlinear diffusion equation and analyze its solutions. The results of our analysis are in good agreement with direct numerical simulations of the corresponding discrete model. It is shown that by analyzing a time dependence of adatom concentration profiles one can estimate the type and strength of interatomic interactions. The paper is organized as follows. In Sec.1 we present the model. In Sec.2 we study both analytically and numerically interaction effects for the case of low adatom coverage. Sec.3 is devoted to analytical treatment of nonlinear diffusion in the case of non-monotonic concentration dependence of the diffusion coefficient. We also compare our results with the results of full scale numerical simulations and results of experimental observation. Sec.4 presents some concluding remarks.

II. MODEL AND EQUATIONS OF MOTION
The transport of particles on the surface is described by the set of random walk equations where θ n is the probability for a particle occupying the n-th binding site on the surface (in the literature on surface science this quantity has a meaning of coverage), W n→ n+ ρ gives the rate of the jumps from the binding site n to a neighboring site n + ρ (the vector ρ connects nearest neighbors) . The terms 1 − θ n (t) in Eqs. (1) take into account the fact that there may be only one adatom at a given site or, in other words, a so-called kinematic interaction.
The probability with which the particle jumps from site n to a nearest neighbor n + ρ satisfies the detailed balance condition, where E n is the binding energy of the particle situated at site n , β = 1/k B T, k B is the Boltzmann constant and T is the temperature of the system. For the transition rates we choose which corresponds to setting the activation energy for a jump to the initial binding energy. Here w ρ = ν ρ e −βE b (w ρ = w − ρ ) is the jump rate of an isolated particle with standard notations: ν ρ is a frequency factor and E b − E n is the height of the random walk barrier. Inserting Eq. (3) into Eqs. (1) we obtain that the random walk of particles on the surface is described by a set of equations In the case when the characteristic size of the particle distribution inhomogeneity is much larger than the lattice spacing one can replace θ n and E n by the functions θ( r) and E( r) of the continuous variable r and, by expanding the functions θ( r + ρ) and E( r + ρ) into a Taylor series, obtain from Eqs. (4) that in the continuum approximation the transport of particles on the surface is described by the equation of the form where the notation w = 1 2 ρ ρ 2 w ρ is used.
We will study the particle kinetics in the mean field approach when the binding energy E is assumed to be a functional of the particle density θ( r, t): E( r) = E(θ). In this case Eq. (5) takes a form of nonlinear diffusion equation where is a nonlinear (collective) diffusion coefficient.

III. SURFACE DIFFUSION AT LOW COVERAGE
In what follows we restrict ourselves to studying particle distributions spatially homogeneous along the y coordinate: θ( r, t) ≡ θ(x, t). We assume that initially the particles are step-like distributed where H(x) is the Heaviside step function. By introducing a centered particle density ξ(x, t) = θ(x, t) − 0.5 θ max /θ max , we see that the initial distribution ξ(x, 0) is an odd function of the spatial variable x. It is obvious that in the no-interaction case (E = const), when the diffusion equation (6)  formed in the process of surface diffusion of thorium on tungsten intersected the initial step-like profile at a point lying well above θ max (see [2,12]). Since then, similar behavior has been found for many electropositive adsorbates whose adatoms are known to interact repulsively. A recent example obtained in the case of surface diffusion of Li on the Dy-Mo (112) surface was discussed in [11]. It is worth noticing that the above mentioned behavior was observed even for rather low coverage: θ max < 0.3 (see Fig. 2 in [11]). Therefore to explain such a behavior one may assume that the binding energy E( r) is linearly dependent on the particle concentration: where E 0 is a site energy and V ( r − r ′ ) is an interaction parameter which includes all types of lateral interactions. In this case the nonlinear diffusion coefficient (7) takes a form where is the diffusion coefficient for an isolated particle (or a so-called tracer diffusion coefficient) and the dimensionless parameter characterizes the strength of the lateral interaction.
The aim of this section is to develop an approach which allows estimating the effects of interparticle interactions in the surface diffusion. It is seen from Eqs. (6) that the spatio-temporal behavior of the centered particle density ξ(x, t) is governed by the equation where τ = D 0 t is a rescaled time and the quantity describes the nonlinear properties of the diffusion and vanishes when α → 0. Taking into account that Eq. (11) with the initial condition given by Eq. (8) is invariant under gauge transformations τ → λ 2 τ, x → λ x, θ → θ (λ is an arbitrary number) one can look for a solution of Eq. (11) in terms of the Boltzmann variable z = x 2 √ τ : ξ(x, τ ) = ζ(z) where the function ζ(z) satisfies the equation with the boundary conditions Eqs. (13), (14) can be rewritten in the form of the following integral equation where erf(z) is the error function [19]. It is seen from Eq. (15) that the concentration profiles ξ(x, t) for different time moments intersect at the point 0, ζ(0) with In the weak interaction/low coverage limit when α θ max < 1 one can replace the function ζ(z) in the right-hand-side of Eqs. (15), (16) by its expression obtained in the linear case: ζ 0 (z) = 1 2 erf(z) and obtain approximately that under the step-like initial condition (8) the concentration profiles θ(x, τ ) for different time moments intersect at the point which corresponds to the concentration Thus the concentration value θ 0 at which the concentration profiles intersect changes in the presence of lateral interatomic interactions: θ 0 > θ max /2 (θ 0 < θ max /2) when the interaction is repulsive (attractive). We applied Eq. (17) to analyze the results obtained in [11] for the diffusion of Li on the Dy-Mo (112) surface at low coverage (θ max ≈ 0.33) for which θ 0 ≈ 0.19 and found out that α θ max ≈ 0.97. It is seen that strictly speaking it is not fully legitimate to use our simple analytical perturbation approach (which is valid for α θ max ≪ 1) to analyze the results of experiments [11] but a qualitative agreement takes place. To validate our analytical results we carried out numerical simulations of Eqs. (4), (3) which in the 1D-case for a system with N binding sites have the form where β E n = α θ n . Thus, in our model the total number of particles is a conserved quantity. As an initial state we used a step-like distribution We found out that for θ max = 0.33 the concentration profiles intersect at the point (0, 0.19) ( as it was observed in the experiment [11]) for α ≈ (3 ÷ 3.5) ( see Fig. 1) or α θ max ≈ (1 ÷ 1.5) which is a good agreement with our analytics.
Thus basing on our approach one can conclude that Li adatoms on the Dy-Mo (112) surface mostly repel each other and the intensity of the repulsion is Diffusion of Li adatoms on Dy/Mo(112) was investigated experimentally at T = 600K [11], so the estimated repulsion energy V 0 amounts to 0.16 − 0.18 eV. Let us assess this value in terms of the dipole-dipole interaction. The energy of the repulsive interaction between two dipoles having moments p and situated on the surface at a distance r is The dipole moment can be determined from the work function change ∆ϕ using the Helmholtz formula for the double electric layer: | ∆ϕ |= 4 π n p e , where n is the surface concentration of adatoms and e is the electronic charge. For Li on the Mo(112) surface, the p value at low coverage was found to be 1.4 Debyes [13]. Then, using Eq. (20) we can find that for two such dipoles the interaction energy U dd = 0.16eV can be attained at a distance r ≈ 2.5Å, which is close to the distance between the nearest adsorption sites (2.73Å) within the atomic troughs on Mo(112). This estimation shows that the intensity of the lateral interaction deduced from the diffusion data in the way presented above seems physically reasonable.
Recall, however, that V 0 determines a resultant effect experienced by a jumping particle from all its counterparts which, in the case of heterodiffusion, are non-uniformly distributed over the surface and provide an additional driving force (supplementary to the coverage gradient) that favors a faster diffusion of repulsing particles.
The adatom interaction effects fade away for the late stage of the evolution when the particle density becomes small and the diffusion process transfers into a linear regime ( see Fig. 1 for t = 60000).

IV. MEAN-SQUARE DEVIATION
The adatom interactions manifest themselves also in the integral characteristics of kinetics of adatom diffusion. It is well known that in the linear regime the variance behaves (in one-dimensional case) as x 2 = 2τ . Therefore it is naturally to introduce a variance rate whose time dependence provides a useful information about nonlinear effects in the diffusion process.
For this quantity we obtain from Eq. (14) that Assuming that initially the particles are concentrated in a finite domain in a Π-like form: where 2l is the size of the initial domain, for small nonlinearities α and low coverage θ max ≪ 1 we obtain approximately that is the solution of the linear diffusion equation with the initial condition (25). In the limit of small l we obtain that We checked our analytical considerations by carrying out numerical simulations of Eqs. In general, the diffusion coefficient is a non-monotonic function of atomic concentration (see e g [5]). There is a number of physical reasons which can cause a non-monotonic coverage dependence of the diffusion coefficient. In thermodynamic terms, the diffusion flux is proportional to the gradient of chemical potential of adsorbed particles µ which can be written as [14], [15] The first term in this equation is the standard chemical potential of the adsorbate, q(θ) is the differential heat of adsorption and the third term stems from the entropy of mixing of adatoms with the vacant adsorption sites on the substrate. (Note that this simplified expression relates only to the first monolayer and does not take into account the possibility of formation of the second and next monolayers). The diffusion coefficient can be represented as a product where D j is a so-called kinetic factor (or jump diffusion coefficient) [1], [5], [15] and the derivative in the brackets is named the thermodynamic factor. In a simplest case when the cross-correlations between the velocities of diffusing particles are absent, D j coincides with the tracer diffusion coefficient D * given by Eq. (10). Inserting Eq. (29) into Eq. (30), we get It is seen from Eq. (31) that any effect which entails a sharp decrease in the heat of adsorption as a function of coverage will result in a maximum of the diffusion coefficient in this coverage range [16]. For instance, such a situation occurs when all energetically profitable sites at the surface are occupied and adatoms start to fill less favorable sites.
Actually, Bowker and King [10] found in their Monte Carlo simulations that a well-pronounced maximum in the D(θ) dependence observed by Butz and Wagner [20] can be explained by existence of two types of lateral interactions: repulsive one between the nearest neighbors and attractive between next-nearest neighbors. A similar effect is typical for volume diffusion of interstitial atoms in disordered binary alloys having a BCC structure with two nonequivalent interstitial positions [21]. In the framework of local equilibrium statistical operator approach [22] it was shown that the physical reason for a non-monotonic concentration dependence coefficient is a combined action of lateral interaction and adatom density fluctuations [23]. A sharp drop in the heat of adsorption is also observed in the transition from filling the first, strongly bound (chemisorbed) monolayer to filling the second, weakly bound (e.g., physisorbed) monolayer. In such a case, the spreading of the first monolayer proceeds through diffusion in the mobile uppermost (second or next) monolayer (the so called "unrolling carpet" mechanism) [1]. This example shows that a change in the heat of adsorption can be accompanied not only by variation of the diffusion parameters (the activation energy and prefactor D 0 in the Arrhenius equation), but also by a change in the atomistic diffusion mechanism itself.
It is worth noting also that the non-monotonic concentration dependence may be phenomenologically connected with a step-like dependence of the heat of adsorption on the coverage ( see a review paper [24]) . Fig. 4 shows the diffusion coefficient calculated from Eq. (7) by assuming that the on-site adatom energy E(θ) (which in most cases is proportional to the heat adsorption q(θ)) has the form where the parameter θ thr gives the threshold value of the coverage and the parameter κ characterizes the sharpness of the transition to the new state [25].
The aim of this section is to consider the diffusion process with a step-like concentration dependent on-site energy given by Eq. (32) and clarify what kind of new information one can derive by comparing theoretically obtained concentration profiles with experimentally obtained ones.
It is very hard and may be hopeless to solve the equation (6) with the diffusion coefficient given by Eqs. (7) and (32). However, the problem can be solved and some insight into the kinetics can be achieved in the limiting case of very sharp energy concentration dependence: κ → ∞ . In this case the diffusion coefficient (7) and (32) takes the form where a = α e α θ thr (1 − θ thr ), b = e 2 α − 1.
The nonlinear diffusion equation (6) with the diffusion coefficient (33) and the initial condition (8) has a self-similar solution θ(x, τ ) ≡ Θ(z), (z = x/2 √ τ ) which can be presented in the form (see Appendix for a detailed derivation) Here the parameters z 1 and z 2 are determined by the equations e α θ max − θ thr e z 2 2 erfc(z 2 ) = θ thr e z 2 1 erfc(z 1 ) Such a behavior is shown in Fig. 5. The rate with which the length ℓ + p of the plateau increases is determined by the nonlinear parameter α and the threshold coverage θ thr . We carried out also numerical simulations of Eqs. (18) with the the step-like on-site energy E n = α 1 + tanh κ (θ n − θ thr ) and as it is seen from It is worth noting that the diffusion of Dy adatoms absorbed by Mo (1 1 2) for the initial coverage θ( shows a very well pronounced plateau in the concentration profile dependence both on the spatial coordinate for different time moments and on the Boltzmann variable ( see Figs. 7 and 8 in [11]). It means that as it is prescribed by our analytical considerations the length of the plateau increases as t 1/2 . This suggests that our simple analytical model may be a useful tool in analyzing an experimentally observed concentration behavior.

VI. CONCLUSIONS AND DISCUSSION
In this paper, we have investigated the role of interactions between adatoms in surface diffusion. The problem was considered analytically in the mean-field approach. By analyzing discrete nonlinear random walk equations and  From Eqs. (A1) and (33) we see that the function where z(Θ) is an inverse function with respect to Θ(z), satisfies the equation and the continuity condition y(θ thr + 0) = y(θ thr − 0) = y(θ c ) .