{'fraction_non_alphanumeric': 0.026670483551611442, 'fraction_numerical': 0.00992390085641356, 'mean_word_length': 4.398583386992917, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Collective diffusion of dense adsorbate at surfaces of various geometry

Convenient variational formula for collective diffusion of many particles adsorbed at lattices of arbitrary geometry is formulated. The approach allows to find the expressions for the diffusion coefficient for any value of the system's coverage. It is assumed that particles interact via on-site repulsion excluding double site occupancy. It is shown that the method can be applied to various systems of different geometry. Examples of real systems such as GaAs with specific energetic landscapes are also presented. Diffusion of Ga adatoms on GaAs(001) surface reconstructed in two different symmetries is studied. It is shown how increasing Ga coverage changes the character of diffusion from isotropic two-dimensional into highly anisotropic, almost one-dimensional. It is shown how important the role of the inter-particle correlations is, which influence the value of the collective diffusion coefficient at higher coverages.


I. INTRODUCTION
Collective diffusion is a complex many-body process that plays an important role in many physical applications. During epitaxial crystal growth interplay between adsorption, desorption and surface diffusion is responsible for the character and the quality of the grown layers. Understanding diffusive processes is therefore important to gain a precise control over the crystal growth. Currently of particular interest is synthesis of atomthick layers such as graphene 1,2 , semiconductors like MoS 2 3 or hexagonal-BN 4 , or topological insulators such as Bi 2 Se 3 5,6 or Bi 2 Te 3 6 . The value of the diffusion coefficient (or, in general, the diffusion tensor), which can be either measured in experiments or calculated using theoretical methods, strongly depends on the system's energy landscape, which describes the depth of all the adsorption sites and the jump rates between them. Specifically, both the rate of diffusion and its anisotropy depend on the details of the energetic landscape. Those parameters play an important role during growth of both bulk crystals and various nanostructures at their surfaces.
Even though crystals usually crystallize in highly symmetric structures, the symmetry of their surfaces is often different from that of the bulk. One of the reasons for it is surface reconstruction, which can lead to a change of the height of the energy barriers 7 and the geometry of the surface 8,9 . Other reasons can be presence of steps at the surface 10,11 or defects 12 .
Theoretical calculations of diffusion coefficients require the knowledge of the system's energetic landscape, which is usually determined from ab-initio methods or from experiments. Among common simulation approaches to diffusion problems are Monte Carlo and molecular dynamics, while analytical methods are based for example on Fokker-Planck, Kramers or master equation.  . Many studies have been devoted to the problems of diffusion over concrete crystal surfaces of complex energetic landscape 24,25,31,33,34 . In such case description of singleparticle diffusion already is a difficult task and the pres-ence of many particles increases only the degree of complexity.
Our approach to calculating diffusion coefficients is a variational one and is based on analysis of master equation. It was first proposed in Ref. 26 and has been subsequently developed to deal with various cases of surface diffusion [27][28][29][30][31][32][33][34] . It has been shown that the approach is effective in analyzing diffusion in real systems, such as Ga adatom at GaAs(001) surface 31 or Cu monomer and dimer on Cu(111) and Ag(111) surfaces 33 . Recently using the variational approach we have derived a general formula for the collective diffusion coefficient for interacting particles in an arbitrary one-dimensional potential 34 .
For a small concentration of particles in the system their diffusive paths can be treated as independent and the collective diffusion coefficient can be considered as being equivalent to that of single-particle diffusion. In Refs 24 and 25 analytical formulas for single-particle diffusion coefficients for various lattices of square, rectangular and hexagonal symmetry have been derived. It is known that collective diffusion in an energetic landscape consisting of adsorption sites of equal depth does not depend on the coverage if only site blocking interactions are considered 35 . Therefore, in such cases collective diffusion is equivalent to that of a single particle in the same potential. In Ref. 30 we presented diffusion coefficients for a few lattices of square and hexagonal symmetry with identical adsorption energies calculated using variational method. However, in the case of energetic landscapes with adsorption sites of different depths the collective diffusion coefficient may significantly depend on the system's coverage even if no other interactions than site blocking are present. The degree to which the diffusion is influenced by the inter-particle correlations depends on the difference between the energies of the adsorption sites, the heights of the diffusive barriers and the number of the nearest neighbors, that is the geometry of the system.
In this article we show how to apply the variational approach in calculation of diffusion coefficient of 2D manyparticle systems. We calculate collective diffusion coeffi-arXiv:1801.01374v1 [cond-mat.stat-mech] 4 Jan 2018 cients for several systems. The approach that we show here is general and can be applied in any other case, regardless of the symmetry of the network. Among various examples we study also lattices considered in Refs 24 and 25, for which single-particle diffusion coefficients were calculated. We show density dependence of the diffusion coefficient for all systems studied by taking into account the role of the inter-particle correlations. In the low density limit our results reproduce formulas obtained for single-particle diffusion.
Next we use the presented method to calculate density dependence of the diffusion on GaAs surface, studied before in Ref. 31. We show that the diffusion anisotropy at this surface can be induced by the presence of other particles. Density dependent particle kinetics on lattices of two different symmetries resulting from different surface reconstructions are compared.

A. Definitions
We consider a gas of particles adsorbed at a lattice of arbitrary geometry. The lattice is periodic and consists of elementary cells, each of which contains several adsorption sites. In the system we assume periodic boundary conditions. From the physical point of view adsorption sites are energetic minima of the potential landscape that comes from the underlying surface. We assume that at any time a given adsorption site can be either free or occupied by exactly one particle, that is multiple occupancy is forbidden. If there are no interactions between particles, then the occupancy of a given site is independent of the occupancies of the other sites. In the equilibrium the probability that the site with the adsorption energy E A is occupied is in a grand canonical ensemble, where z = e βµ , µ is the chemical potential and β = 1/k B T is the Boltzmann factor. The mean coverage of the system is defined as where the sum is performed over all the adsorption sites in an elementary cell and n is the number of those sites. Most of the time the adsorbed particles reside at their sites but from time to time they perform thermally activated jumps to neighboring sites. The jumps occur over the energetic barriers, which are the saddle points of the potential landscape between the two involved adsorption sites. The jump rate from the site A to the site B is where E A is the adsorption energy of the site A and E AB is the energy of the saddle point between the sites A and B. The prefactor ν is the attempt frequency, which is the measure of how many attempts per unit time the adsorbed particle performs to jump out of its site. Its value is determined by the vibrations of the atom lattice, that is phonons. Usually we will assume that it is equal to the typical phonon frequency in condensed matter, which is ν = 10 13 /s. At thermal equilibrium detailed balance condition is fulfilled, which means that the probability of every microscopic process is equal to the probability of the reverse process. For our system that condition can be written as It is easy to see that (1) and (3) satisfy the above condition.
Diffusion is the process of random jumps of the particles by which the system approaches the thermal equilibrium. We will consider diffusion that takes place close to the equilibrium, that is when there are only small deviations from that state. Such situation occurs for example during epitaxial crystal growth.

B. Variational approach
In order to calculate collective diffusion coefficients for the lattices considered in the next chapter we use the variational approach that was first proposed in Ref. 26. For one-dimensional collective diffusion it was possible to derive by means of that approach a general formula for the diffusion coefficient expressed in terms of jump rates, interaction constants, temperature and chemical potential 34 . Diffusion in one dimension can be considered as an extreme case of anisotropic surface diffusion. There is always only one channel of diffusion in one dimension since a diffusing particle must overcome all barriers ordered one after another to move to the next elementary cell. However, due to inter-particle correlations the dependence of the diffusion coefficient on the barriers is still not trivial and we showed that its behavior can depend not only on the height of those barriers but also on their order in an elementary cell.
Diffusion in two dimensions is a much more complicated process since there are a lot of different paths along which a particle can move. Because of that the diffusion depends not only on the values of the jump rates but also on the geometry of the lattice. We shall derive a variational formula for a non-interacting case, which depends on a set of variational parameters and can be explicitly written out for a system of arbitrary geometry. Then it can be used for finding the values of the variational parameters and therefore determining the diffusion coefficient.
The starting point in our variational approach is the master equation written for the lattice gas system de-scribed in the previous section which describes the evolution of the probability of finding the system in the microstate {c}. A microstate is defined as (6) which specifies which of the adsorption sites at the lattice are occupied. It distinguishes one of the particles as the reference particle at the position X + l 0 , where X is the position of the origin of its elementary cell and l 0 is its position within that cell. The set of the remaining vectors { m 1 , m 2 , . . . , m N −1 } gives the positions of all the other particles with respect to the reference particle. Since the system has translational symmetry we can define a configuration as which specifies only the position of the reference particle with respect to the origin of any elementary cell and the positions of the remaining particles with respect to the reference particle. All the microstates in which the positions of the particles differ only by a linear combination of the lattice vectors (that is by X) correspond to the same configuration. We can use the definition of configuration to replace W ({c}, {c }) in Eq. (5) by W { m},{ m } . At given time we allow only for one jump of a particle between neighboring sites. For pairs of { m} and { m } which fulfill that condition W { m},{ m } are equal to corresponding jump rates (3) and in the equilibrium P ({c}, t) = P { m} ( X, t) become P eq { m} , which for non-interacting particles is equal to products of θ i and 1−θ i , whose values are calculated by means of (1). For other pairs of configurations W { m},{ m } are equal to zero.
The periodicity of the system allows us also to take a Fourier transform of the master equation (5) to get which is a matrix equation where the position vector X has been replaced by the wave vector k. The components of the vector P( k, t) correspond to the configurations { m} andM ( k) is the rate matrix, whose properties have been described in detail in Refs 26-28, with elements k) equal to 1 except when the reference particle crosses the boundary between elementary cells. In the latter case F { m},{ m } ( k) = exp(±i k A), where A is the lattice vector between those cells and the sign depends on the direction of the jump.
When the system is brought out of the equilibrium state, it will start returning to that state with the rate described by the rate matrix's eigenvalues, which are always real and negative and correspond to the decay of the kth Fourier component of the fluctuation. The components of the leftẽ( k) and the right e( k) eigenvector are related to each other by Close to the equilibrium the eigenvalue λ D ( k) with the lowest absolute value is proportional to k 2 =| k | 2 and gives the diffusion tensor In order to find that eigenvalue we use a variational approach. We denote a trial eigenvector byφ and we obtain its trial eigenvalue λ var which is the variational formula used by us to solve various diffusional problems [26][27][28][29][30][31][32][33][34] . M ( k) and N ( k) have been referred to as the expectation value numerator and the normalization denominator, respectively. The most important step in the variational method is choosing the proper eigenvector. Following our earlier works we assume the form where m i is the position of the ith particle in the { m} configuration, while δ mi and ∆ mi are the variational parameters, called geometrical and correlational phase, respectively, assigned to that particle in that configuration. We insert that trial eigenvector into the variational formula (12), which we then differentiate with respect to the variational parameters in order to find its minimal value. The geometrical phase δ mi depends on the symmetry of the potential with respect to the site m i . By that we understand not only the nearest neighbors of that site but also the global structure of the potential. On the other hand, the correlational phase ∆ mi is in general non-zero when two neighboring sites with different adsorption energies are both occupied. For a given particle the value of δ mi depends only on the position of that particle, while the value of ∆ mi depends also on the occupancy of neighboring sites. Both for the geometrical and the correlational phases one can use symmetry arguments in order to deduce some properties of the values they can assume. For lattices exhibiting some kind of symmetry one can prove that some of the phases are equal to zero, or one can relate the value of one phase to that of another one.
Since N ( k) for our choice of the trial vector does not depend on the variational parameters, we shall consider the factors M ( k) and N ( k) separately. Therefore, in order to minimize the trial eigenvalue given by the formula (12) it is enough to minimize its numerator M ( k). Following our earlier works we write it in the form (14) Then we insert the trial vector (13) into the above expression to find the values of the variational parameters for which the expression has the smallest possible value. In the limit k → 0 and assuming that the configurations { m} and { m } differ only by the position of one particle we can rewrite the numerator as where the first summation is done over all the possible jumps in the system, that is all pairs of nearest neighbors with A denoted as the site from which the jump starts and B as the site to which the particle jumps. The second summation is performed over all possible occupational states of the neighbors of sites A and B. Because of the detailed balance condition (4) it is enough to consider each jump only in one direction, as indicated by the "no rep" comment over the sum symbol.
As mentioned before, the values of the correlational phases depend on the occupancy of the neighboring sites. a AB is the vector connecting the sites A and B.
On the other hand, the normalization denominator can be written as which in the limit k → 0 for particles with only exclusion interactions has a very simple form Once we have the factors (14) and (16) expressed in terms of jump rates and equilibrium occupancies, the diffusion coefficient is the ratio of those factors.

C. Numerator
Calculating the numerator directly from the Eq. (15) for a specific system is rather tedious as for each jump one has to consider all relevant configurations in the system, that is in general all occupational states of the neighborhood of the sites involved in the jump. However, we shall show now that we can perform the second sum in Eq. (15) in a general case and obtain a formula which can be easily used for calculations. First of all, we assume that the total correlational phase of a particle is the sum of correlational phases related to all of its neighbors Then let us rewrite the second sum in Eq. (15) as where we used the fact that the geometrical phases do not depend on the occupancy of the neighboring sites so we could pull them out before the sum, which is equal to 1 if there are no correlational phases under it. In the above expression we have two sums. The one in the mixed term can be rewritten us where {A(B)} denote the set of all neighbors of the site A(B) excluding the site B(A). On the other hand, the sum with the squared term is where the first sum on the right-hand side can be written as and of course we can get an analogous expression for the second sum in (21). Similarly one can show that the last sum there is Now after using the results above and collecting all terms in the expression (19) we can show that the final formula for the numerator in non-interacting case is After we assign variational parameters to the sites and to the pairs of neighbors in the system this formula is simple to use. In the next chapter we shall show a few examples illustrating how to use that formula to calculate the collective diffusion coefficient.

III. RESULTS
We illustrate this method on several examples. Below we show coverage dependence of collective diffusion coefficients for lattices of complex energy landscapes calculated using variational approach. The formula (24), which is the starting point of our calculations, depends on a few variational parameters, called geometrical and correlational phases, with respect to which it must be minimized. In general, one assigns one geometrical phase to each adsorption site and one correlational phase to each pair of nearest neighbors, however, usually it is possible to reduce the number of independent variational parameters using symmetry arguments.
In order to minimize the variational formula, one has to solve a set of linear equations obtained from differentiating it with respect to the assumed parameters. Such set of equations is always soluble, for example by means of Cramer's rule. Finally, the obtained values of the parameters are inserted into the formula (whose squared terms can be eliminated in the minimum, which simplifies the expression further) to get the diffusion coefficients as functions of jump rates, temperature and the system's coverage.
The above procedure is always the same and does not depend on the lattice symmetry. However, depending on the complexity of the investigated system calculations can be shorter or longer. We start with three examples of square lattices, and then two different hexagonal lattices are analyzed. Finally we show results for lattices that model experimental system of Ga diffusing on GaAs(001) surface. Diffusion coefficients for two different reconstructions of this surface are compared.
A. Square lattice systems 1. Two equivalent sublattices The first system we analyze is shown in Fig. 1. It is a lattice of square symmetry, similar to the checkered lattice considered in Ref. 29, however, in addition to jumps between different sites there are also jumps between identical sites. One elementary cell contains one A and one B site. All jumps are symmetric and there are 8 possible jumps per elementary cell: 2 between A sites, 2 between B sites and 4 between A and B sites.
Because of the central inversion of the lattice in every adsorption site, all the geometrical phases are equal to zero. The only non-zero variational phases will be correlational ones related to the AB bonds. There will be 4 such phases in total, however, due to the lattice symmetry with respect to the rotation by 90 degrees, they will differ only by the direction. Therefore, they can be expressed by one variational parameter, ∆ ±± AB = ∆ AB (±1; ±1), where the index ±± refers to the orientation of the AB bond at the lattice. The numerator of the variational formula is then . After differentiating that formula with respect to ∆ AB and setting the derivative to zero we get Using the condition for the derivative ∂M ∂∆ AB = 0 we can eliminate the square terms in the expression (25) and  Since there is only one adsorption site of each type in an elementary cell, the denominator is The diffusion is isotropic and D x = D y = M/(N (k 2 x + k 2 y )).
The first three terms in Eq. (28), each of which contains one jump constant multiplied by adequate occupation probabilities, can be easily interpreted as contributions to the diffusion from respective paths without considering correlational effects. The term with W AA corresponds to diffusion through A sites, the one with W BB describes the diffusive path through B sites and the last one with W AB accounts for the diagonal path that goes alternately through A and B sites.
On the other hand, the last term in Eq. (28) can be attributed to the inter-particle correlations in the system. It comes from the correlational phase that was initially assumed in the variational formula. That term is negative, which reflects the fact that the variational approach helps to estimate the upper value of the diffusion coefficient. Both in the limit θ → 0 (single particle) and θ → 1 (single hole) the correlation term vanishes. For a single particle the other terms together with the denominator (29) yield the diffusion coefficient which is identical to the expression derived in Ref. 24 for the same lattice. We used the detailed balance conditions for A and B sites while taking the single-particle limit. For a single hole the diffusion coefficient is For W AA = W BB = 0 the diffusion at this lattice is the same as at the checkered lattice from Ref. 29, D p = D h and the maximum of the diffusion coefficient is for θ = 0.5. However, when we add jumps between sites of the same type, represented by W AA and W BB , then the dependence of the diffusion coefficient on the coverage will change significantly. In Fig. 2 we show the dependence of the collective diffusion coefficient on the mean coverage of the system for three different cases. In all of them the depths of the adsorption sites are assumed E A = 0 eV and E B = 0.05 eV and the barrier between them is E AB = 0.07 eV. The barriers for W AA and W BB jumps differ between the curves. The bottom curve represents the diffusion where those barriers are infinite, that is the corresponding jumps are blocked. That case is equivalent to diffusion at the checkered lattice 29 . The middle curve is for E AA = 0.1 eV and jumps between B sites are blocked. We can see that activating a channel between A sites speeds up the diffusion. Because A sites are deeper than B sites, the change in the diffusion coefficient is bigger for lower values of the concentration. Finally, the upper curve shows the diffusion for E AA = 0.1 eV and E BB = 0.12 eV. This time, as expected, the diffusion is more strongly accelerated for high values of θ. Even though E BB > E AA , the effect is stronger than for the channel between A sites. That is because the value that decides about the jump rates is the difference between the barrier and the depth of the adsorption sites that are involved in the jump. In our case E BB − E B < E AA − E A .

Lattice of square symmetry -inequivalent sublattices
Now let us consider a lattice shown in Fig. 3. It possesses square symmetry but differs significantly from the square lattice in Fig. 1. In the geometric sense the adsorption sites A and B are not equivalent. Each A site connects to 4 other A sites and to 4 sites B, while each B site connects to 6 sites B and to 2 sites A. Each elementary cell contains one A site and, contrary to the previous examples, two B sites. Those two B sites differ by the angular orientation of the surrounding sites. They have the same adsorption energy but one of them is placed between the two A sites in the horizontal direction and the other in the vertical one. Therefore, the total number of sites per elementary cell is 3. One should also note that there are two types of jumps between the B sites. One of them is a short jump in the diagonal direction and the other one is a long jump in the horizontal or vertical direction. We will denote the corresponding jump rates by W s BB and W l BB , respectively. Again, there will be one correlation phase, which will possess four versions with respect to the angular orien- tation. As in the previous case they can all be expressed by one variational parameter ∆ AB . The numerator is From there we obtain The denominator in this case is

The diffusion coefficient is calculated again by
x + k 2 y )) and its dependence on the mean coverage for three sets of the jump rates is shown in Fig. 4. As can be expected when we decrease energy barrier for long distance-jumps between B and B sites E l BB , diffusion coefficients increase for all surface coverages. However, interestingly the increase of high-coverage diffusion coefficient is much faster. As a result for higher barrier low-coverage diffusion coefficient is higher than high-coverage diffusion coefficient, while for lower barrier high-coverage diffusion becomes higher by three orders of magnitude than the low-coverage one.

Square lattice with three different adsorption sites
We analyze diffusion at the lattice shown in Fig. 5. The lattice is of square symmetry with three different adsorption sites. In general they all differ by the adsorption energy. Sites A and B are geometrically equivalent to each other in the sense that they both have four C neighbors, one at each of the four main directions. Sites C, on the other hand, have two A neighbors placed either horizontally or vertically and two B neighbors placed in the perpendicular direction. Each elementary cell has one A site, one B site and two C sites. The two C sites differ by the angular orientation of their neighbors.
There are two types of jumps between adsorption sites: W AC and W BC . Per elementary cell there are four W AC jumps and four W BC jumps. Again, because of the symmetry of the adsorption sites there are no geometrical phases. There are two correlational phases, ∆ AC and ∆ BC , related to the respective pairs of neighbours. The numerator is As before, in order to find the variational parameters we have to differentiate the above expression with respect to them and set the derivatives to zero. That will give us a set of two equations. In the case considered here each equation will contain only one phase so one can easily calculate their values Finally we get the numerator The denominator is From the above expressions one can simply calculate the diffusion coefficient as the ratio D = M/(N (k 2 x +k 2 y )). From Eq. (38) it is seen that there are two main diffusive paths: one through W AC jumps and the other through W BC jumps. Again, the negative terms are due to the inter-particle correlations.

System of two hexagonal sublattices
Next we analyze a lattice of hexagonal symmetry shown in Fig. 6. Here each site contains 6 connections to sites of the same type and 3 connections to sites of the other type.
Similarly to the previous cases, using symmetry arguments one can reduce the number of independent variational phases to one correlation phase between A and B sites. The numerator of the diffusion coefficient is From there one can find the correlational phase which gives the final value of the numerator The denominator is the same as for the case of two square sublattices, that is N = θ A (1−θ A )+θ B (1−θ B ) and again D x = D y = M/(N (k 2 x + k 2 y )). Expression (42) has in fact the same structure as the numerator (28) for the square lattice. The only difference is in the coefficients that stand by various terms. That is not surprising since both lattices contain two different adsorption sites per elementary cell and they are connected with each other similarly at both lattices. They differ by the number and the geometry of those connections but at both lattices the three main diffusive paths are along the W AA , W BB and W AB jumps.
We show the collective diffusion coefficient for the hexagonal lattice in Fig. 7 for the same set of parameters as for the square lattice from Figs 1 and 2. The behavior of the diffusion is qualitatively the same in both cases but we can see that here the diffusion is generally lower and the curves are flatter than in Fig. 2, especially the one for E AA = E BB = ∞.

Lattice of hexagonal symmetry -three sublattices
In Fig. 8 we show another lattice with three different adsorption sites. All of them have the same geometry, each site of a given type has six connections to sites of two other types in that way that A site connects to 3 B and 3 C sites and so on.
There are three possible types of jumps in the system: W AB , W BC and W CA and per elementary cell each of Unlike previous examples, here mixed terms (for example ∆AB∆BC) enter the expression for the numerator. That means that the derivatives will contain also other phases than that with respect to which the expression is differentiated. We shall write the set of linear equations for the correlational phases in the matrix form where As we can see this case is more complex than the previous ones and also it is more difficult to obtain the solution. However, one can solve the set of equations using Cramer's rule Again, we can simplify the expression (43) by using the condition of vanishing derivatives to get Then we put the phases (46) into the above expression and after dividing it by the denominator we get the collective diffusion coefficient D.

C. Diffusion on GaAs(001)surface
The variational method is not limited to lattices whose adsorption sites possess symmetry. It can be just as well used to analyze diffusion in systems where there is an asymmetry of the jumps from a given site. In such cases one should assume a geometrical phase related to that site. In general, one should look not only at the symmetry of the jumps from that site but also at the symmetry of the entire system. Even when the jumps are symmetric, the geometrical phase can be non-zero when the whole potential is not symmetric with respect to the site considered. The whole procedure of finding the diffusion coefficients is the same as for symmetric landscapes. The variational formula is minimized with respect to the parameters by simple differentiation, which gives a set of linear equations.
In this section we shall consider two examples of lattices which contain at least one site with respect to which the energetic landscape is asymmetric and consequently its geometrical phase cannot be set to zero. These examples model Ga adatoms diffusion on GaAs(001) surfaces in two different reconstructions. We show how diffusion anisotropy changes with the system coverage in both cases.

Ga adatoms on GaAs(001)-c(4x4)
In Fig. 9 a simplified lattice for Ga adatom at reconstructed GaAs(001)-c(4x4) surface is shown. The scheme of jumps is based on ab-initio calculations for that system in Ref. 36. Based on those data diffusion coefficients for a single Ga atom were calculated in Ref. 25. A more detailed landscape was found in Ref. 37, which contained a few more adsorption sites and energetic barriers. Single-particle diffusion coefficients based on those new data were found in Ref. 31 by means of our variational approach.
Here we shall calculate collective diffusion coefficients for the lattice shown in Fig. 9. Unlike earlier approaches, the result will depend on the coverage of the system as our method takes into account inter-particle correlations.
As usual, we start by writing out the variational formula for the lattice studied. As shown in Fig. 9 one elementary cell contains 1 A site, 2 B sites and 2 C sites. There are two types of jumps: W AB and W BC . Per elementary cell there are 2 W AB jumps and 4 W BC jumps. Because of that we will have two correlational phases, each assigned to one pair of neighboring sites: ∆ AB and ∆ BC . Depending on the orientation of a given pair in the lattice they will assume the values ∆ ±x AB = (±∆ AB , 0) and ∆ ±x±y BC = (±∆ x BC , ±∆ y BC ). Since sites A and C are inversion centers of the lattice, the geometrical phases at those sites are equal to zero. On the other hand, site B does not have such symmetry, however, it has symmetry with respect to the x axis and two inequivalent B sites in an elementary cell are related to each other by the symmetry with respect to the y axis. Because of that the y component or its geometrical phases will vanish, while the x component will differ by the sign between the B sites neighboring to the left and to the right of the A site, that is δ ± B = (±δ B , 0). Taking all those facts into account we can write the numerator as We can see that the diffusion is anisotropic here, however, its main directions are always along x and y axes. The xcomponent depends on three parameters: δ B , ∆ AB and ∆ x BC , while the y-component depends only on one parameter, ∆ y BC . One has to minimize both components with respect to those parameters.
For the x-direction we get a set of three equations, which can be again written in the matrix form similar to (44) and again we can obtain the phases The simplified formula (50) for the x-direction is to which we insert the phases (53) and obtain the diffusion coefficient where Calculating diffusion coefficient in the y-direction at this lattice is much easier as it requires to minimise the terms in the formula (50) which stand by k 2 y with respect to only one variational parameter, namely ∆ y BC . After which inserted into simplified formula (50) for the ydirection gives and In Fig. 10 we show the collective diffusion coefficients at that lattice as functions of the mean coverage. We take the values of the adsorption energies and the energetic barriers from Ref. 36. We see that both coefficients start at around 2 · 10 −8 cm 2 /s for θ = 0, which is the value of the same order of magnitude as obtained in the Monte Carlo simulations in Ref. 36. For low coverages the diffusion is almost isotropic. However, starting from θ = 0.2 it becomes faster in the y-direction, while in the x-direction it remains constant. At θ = 0.6 D y reaches its maximal value while D x starts to decrease rapidly. Diffusion in the y-direction becomes two orders of magnitude faster than in the x-direction for the single hole limit (θ = 1). Such anisotropy should be easy to verify in STM measurements or other experimental techniques.

Ga adatoms on GaAs(001)-(2x4) β2
We consider a different reconstruction of the GaAs(001) lattice, namely (2x4) β2. We show the scheme of the lattice in Fig. 11, which was based on ab-initio calculations in Ref. 38.
The lattice has two types of adsorption sites, to which we assign geometrical phases δ A = (0, δ A ) and δ ± B = (±δ x B , δ y B ), where the sign ± refers to the B site lying closer either to the left or to the right of the A site. Each cell has one A site and two B sites. There is also one correlational phase ∆ ± AB = (±∆ x AB , ∆ y AB ) assigned to the AB bond, where ± again refers to the position of the B site to the left or to the right with respect to the A site. The jumps per elementary cell are: one W AA , two W AB , one W h BB (horizontal), two W v BB (vertical) and two W d BB (diagonal). We write the numerator as Again one has to minimize the above expression separately for the x-and the y-direction with respect to the variational parameters. For the x-direction we have The corresponding diffusion coefficient is In the above we used the denominator N = θ A (1 − θ A ) + 2θ B (1 − θ B ) and the variational parameters should be taken from Eq. (60). The obtained formula is valid for any value of the coverage θ. For θ → 0 we have ∆ x AB = 0 but δ x B = 0 and which is the same expression as the one obtained in Ref. 25. Now let us consider diffusion in the y-direction. There are three variational phases related to that direction: δ A , δ y B and ∆ y AB . However, they enter the numerator (59) only as the difference δ y B −δ A . Therefore, we can use that difference as a new parameter δ y AB reducing the number of independent variational phases to two. After differentiating we get  Fig. 11. a=1.13 nm, EA=0 eV, EB=0.3 eV, EAB=0.8 eV, E h BB =0.7 eV and E v BB = E d BB =0.9 eV. The bottom curve is for the x-direction and the upper one for the y-direction. Temperature T=473K and attempt frequency ν = 2 · 10 13 /s. Even though we allowed for a non-zero value of ∆ y AB it turned out to be zero. It reflects the fact that the two channels of diffusion in the y-direction, one through A sites and the other through B sites, are independent from each other. The final result for the diffusion coefficient is Since the contributions to the diffusion in the y-direction from the W AB jumps cancel out, the above expression does not depend on that jump rate. It is easy to see that in the θ → 0 limit one again gets the result obtained in Ref. 25. Again we illustrate the results by plotting the dependence of the diffusion coefficients as functions of the coverage in Fig. 12 with the values of the barriers taken from Ref. 38. We can see that there is a high anisotropy in the whole range of the coverage. The behaviors of diffusion in the two directions are completely distinct. Diffusion along x-direction is mostly through three adsorption sites: one deeper (A) and two shallow ones (B). Therefore, it is similar to the one-dimensional case discussed in Ref. 34. There is a characteristic peak at θ = 1/3 related to ordering of the system at that coverage, in which sites A are completely filled while sites B are free. On the other hand, diffusion along y-direction exhibits a behavior of switching between the deep and the shallow diffusive channel. Such behavior was earlier observed in the model of hills and valleys 32 , where the number of sites in the deep channel and the shallow one was equal and the transition between them occurred at θ = 1/2 instead of θ = 1/3.

IV. CONCLUSIONS
We have formulated in a simple way a variational method of analyzing collective diffusion of particles adsorbed at lattices of arbitrary geometry, interacting by infinite on-site repulsion, excluding double site occupancy. We have derived a convenient variational formula which is expressed by all jump rates and equilibrium occupancies of the adsorption sites in a given system, that is its energetic landscape. Finding diffusion coefficients from that formula is a simple and straightforward procedure. It has been shown how the method can be applied to various systems, including real ones with specific values of the energetic barriers.
We showed explicit calculations of the diffusion coefficient based on our variational approach for a few lattices taken from Refs 24 and 25. In each case our results reproduce the expressions for the single-particle diffusion coefficient (taking the limit θ → 0) calculated in those publications. Our approach allows to find the expressions for the diffusion coefficient for any value of the system's coverage θ. For the values of θ between 0 and 1 the effects of inter-particle correlations become important and the collective diffusion coefficient is different from the one obtained without assuming correlations in the system. Only in the case with identical adsorption energies of all sites the diffusion coefficient does not depend on the coverage 30,35 . However, in real systems that condition is rarely fulfilled.
The presented method was applied to two different re-constructions of GaAs(001) surface and we can study the coverage dependence of the diffusion coefficient. Increasing Ga coverage not only changes value of the diffusion coefficient, but also dramatically influences the ratio of diffusion in two main directions. At the surface of c(4x4) reconstruction diffusion changes from isotropic for low Ga coverages into strongly one-dimensional for high Ga coverages, for which diffusion coefficients in two directions differ by three orders of magnitude. At the surface of (2x4) β2 reconstruction symmetry diffusion of Ga adatoms is faster in one direction for all coverages, but ratio between D x and D y changes from three orders of magnitude, then decreases to one for coverage 0.35 and for higher coverages increases again. Diffusion is never isotropic for this reconstruction of GaAs(001) surface.
It is important to note that the approach presented in this article is not limited to the examples shown in the section with the results. The variational method is not dependent on the symmetry of the lattice, it is the variational parameters whose values adjust to the geometry of the energetic landscape and the correlations between the diffusing particles to give a proper expression for the diffusion coefficient. Our approach is also not limited to two dimensions, without any modifications in the formalism one can apply it also to analyze bulk diffusion.