In the growth process of various crystals, the diffusion of clusters plays an important role in the evolution of the morphology of the crystal to its final equilibrium shape. Therefore the kinetics processes which control the mobility of clusters are of interest both experimentally and theoretically. The mobility of clusters has been of interest to researchers for a long time. ^{[ 1 , 2 ]} One of the early studies in this area concentrated on the diffusion of small rhodium clusters on rhodium surface and found that the diffusion is predominant due to the atoms that move along the perimeter of the cluster. ^{[ 1 ]} In the same study the dependence of diffusion coefficient D on the number of particles N in the cluster is in the form of D ∝ N ^α where α is found to be approximately –1.76. In a Monte Carlo simulation of cluster diffusion on a triangular lattice gas model, the dependence of cluster diffusion coefficient D on the attachment/detachment, vacancy motion and particle motion along the cluster edge processes are investigated among others. ^{[ 2 ]} The main process by which the cluster moves is found to be by the motion of particles along the cluster edge. The oscillatory nature of small clusters’ (2–12 atoms) diffusion coefficient is investigated in Ref. [ 3 ] by field ion microscopy. Diffusion of large Ag clusters on Ag is studied by scanning tunneling microscopy and the main contribution to the diffusion coefficient is found to be the evaporation/condensation mechanism of particles from/to a cluster. ^{[ 4 ]} They also found that the diffusion coefficient changes very little with the size of the cluster. In another experimental study which is coupled with a theoretical model suggestion, the mobility of monoatomic deep vacancy islands was studied. ^{[ 5 ]} The motion of vacancy island is observed by scanning tunneling microscope and its diffusion coefficient is studied as a function its size. The main contribution to diffusion coefficient is found to be the motion of atoms across the vacancy. In a kinetic Monte Carlo study, the diffusion of clusters due to the motion of the particles along the cluster edge only is studied in Ref. [ 6 ] where they found that the diffusion coefficient depends on the size of the cluster and it is fitted to a power law best when small size and large size clusters are fitted separately. The size dependence of cluster mobility is studied in Ref. [ 7 ] by using a lattice model. For small clusters, they found that the diffusion of particles along the cluster edge is the dominant mechanism. The defects on the surface also play a role in the mobility of clusters, they capture small clusters but above a critical size the cluster becomes mobile. ^{[ 12 ]} The behavior of clusters in two and higher dimensional lattices is investigated from a different perspective in Refs. [ 13 ] and [ 14 ] where they study the scaling properties of cluster heterogeneity defined as the number of distinct cluster sizes.

In the present work we consider clusters of various sizes on a crystal surface and study their approach to equilibrium and their mobility under equilibrium conditions as a function of temperature. We study the motion of clusters using kinetic Monte Carlo (KMC) method. ^{[ 10 ]} We take into account particle attachment/detachment to/from clusters, diffusion of particles along cluster edges and free diffusion of particles on the surface. The cluster is assumed to reside on a square lattice grid and the motion of particles is determined by probabilities assigned for each particle along the possible four directions on the square grid. Only the motion of a single cluster is studied at a time. A cluster of a predefined size is created on the surface (see the inset of Fig.  1 on the left) and the system is simulated until the cluster reaches an equilibrium state (see the inset of Fig.  1 on the right) where the number of particles in the cluster on average is a constant. By calculating the mean square displacement of center of mass of the cluster, the diffusion constant of the cluster is obtained. The dependence of diffusion constant on temperature, cluster size, and diffusion barrier along the cluster edge is also considered. The paper is organized as follows: Section 2 gives a brief account of simulation method and surface processes considered in this study. Section 3 provides the main results of our simulations and their discussion. A conclusion is provided in Section 4.

Fig. 1.

Figure Option View Download New Window

Fig. 1. The evolution of a cluster to its equilibrium shape where the number of particles in the cluster fluctuates around an average value. A typical initial cluster shape is shown in the inset on the left, and a typical equilibrated cluster shape is shown in the inset on the right where some freely diffusing particles on the surface are visible. The cluster resides on a square grid of lattice constant a of dimensions 350×350 in all simulations presented in this work. The energy parameters used in the simulation (unless otherwise indicated) are as follows (all in units of meV): E b = 125, E d = 150, and E ds = 100. The grey color scale in the insets represents the height of the single atomic layer of the cluster, the lowest level being the darkest in the figure.

As we have already indicated in the Introduction, we consider the mobility of a cluster that resides on a square grid of lattice constant a . The grid has a dimension of 350×350 in units of a in all simulations presented here. A typical initial cluster shape on the surface is shown on the left inset of Fig.  1 . The cluster consists of a single layer of particles. We consider the transfer of particles by surface diffusion only and desorption of particles from the surface is not allowed. We also assume that there is no flux of particles to the surface. The single layer cluster ejects particles from its perimeter to the surface or captures particles that are moving freely on the surface. The shape of the cluster also changes due to the motion of particles along the cluster edge. We employ the KMC method to simulate the motion of the particles and as a result the change of the center of mass of the cluster. In KMC method the energy barriers the particles may encounter before a move are required. The kinetic processes we consider are as follows: i) free diffusion of particles on the surface, ii) attachment/detachment of particles to/from cluster edges, and iii) diffusion of particles along a cluster edge. We also consider Ehrlich–Schwoebel barrier, ^{[ 15 , 16 ]} which takes into account the difference between a particle’s up-cluster and down-cluster motion. Since we do not consider any flux to the surface this process has no effect on the mobility of the clusters. Only the last two processes just mentioned, namely the attachment and detachment of particles to/from cluster edge and the motion of particles along the cluster edge lead to the motion of the center of mass of the clusters.

The energy barrier a particle encounters before it moves to a vacant nearest neighbor position on the surface is assumed to be in the form

In our simulations, first a circular cluster of a predetermined size is created on the surface (Fig.  1 ). Then this surface is allowed to equilibrate by considering the kinetic processes mentioned above in the KMC method. The cluster exchanges particles with the surface by the attachment and/or detachment of particles from the cluster edge. After a while as can be seen in Fig.  1 , the cluster and the free particles on the surface come to an equilibrium situation where the number of particles in the cluster begin to fluctuate around an average value. We call the cluster obtained in this manner the equilibrated cluster. After that the equilibrated cluster and those particles on the surface that are diffusing freely are taken together as the initial surface. That initial surface is simulated using the KMC method and the center of mass of the cluster is recorded at equal MC steps. Using the same initial equilibrated cluster and the surface atoms, the simulation experiment is repeated at least five times with different random number sequences and their average is taken.

The diffusion coefficient of a cluster can be obtained from the Einstein relation

The dependence of diffusion coefficient D of the cluster on temperature and number of atoms N in the cluster is predicted to be in the form ^{[ 2 , 4 , 6 , 8 , 9 ]}

We consider the motion of an equilibrated cluster as defined in Section 2 using the kinetic processes mentioned in that section. The motion is studied as a function of temperature. The kinetic parameters of the problem are E _b , E _d , and E _ds defined in the previous section and unless we indicate otherwise we use the following values for these parameters (all in meV’s): E _b = 125, E _d = 150, and E _ds = 100. Rather than the precise values of these parameters, their relative values are important for the purposes of our simulations. The cluster is located on a square grid of size 350×350 in units of lattice spacing a in all simulations presented. The number of particles in a cluster varies from 150 to 6000. Periodic boundary conditions are used at the edges of the square grid. The approach of a cluster to its equilibrium shape where the number of particles in the cluster fluctuate around an average value is shown in Fig.  1 . The initial cluster shape where there are no free atoms on the surface is shown in the inset of the same figure on the left. It is a circular cluster of a predetermined size. During the evolution of the cluster to its equilibrium size and shape, the cluster and surface atoms are subject to the kinetic processes (i), (ii), and (iii) mentioned in Section 2. When the simulation begins, particles are released to the surface from the perimeter of the cluster and some of them are recaptured during the evolution processes and finally the cluster reaches an equilibrium state where the number of atoms in the cluster begins to fluctuate around an average value. The last state of the surface consisting of the cluster and the free atoms on the surface are taken together as the initial starting surface for our simulations and this morphology is simulated at least five times independently using different random number sequences and then the average values of relevant quantities are taken. The simulation of clusters of a predetermined size is carried out in the way just explained for all clusters considered in the present study. A typical equilibrated cluster shape is shown in the inset of Fig.  1 on the right. As one can see the final shape looks more like a square whose corners are rounded. This is an expected shape because of the underlying geometry. Ejection of particles from a straight edge of the cluster is less probable than from the corners where the edge is rougher or more kinky. Therefore most of the attachment/detachment processes take place near the corners.

From Eq. ( 3 ), one can see that the slope of log(〈 r ² 〉)–log(〈 t 〉) plot must be unity. Similarly, the plot of 〈 r ² 〉/〈 t 〉 as a function of time is expected to approach a constant value at a long time. These two behaviours are consistently checked in our simulations. A typical log(〈 r ² 〉)–log(〈 t 〉) graph is shown in Fig.  2 . The corresponding (〈 r ² 〉)/〈 t 〉–〈 t 〉 graph of the same plot is shown in the inset (on the right at the bottom) of Fig.  2 . The slope of the former is close to unity (1.03) and the latter one approaches a constant value as t increases. The constant value in the latter plot is the diffusion constant of the cluster. That same quantity can be found from the slope of (〈 r ² 〉)–〈 t 〉 plot as discussed below, also shown in the inset (upper left corner) of Fig.  2 is the variation of (〈 r ² 〉) as a function of time which is expected to be a straight line as shown in Fig.  2 .

Fig. 2.

Figure Option View Download New Window

Fig. 2. A typical log(〈 r 2 〉) versus log( t ) graph whose slope is expected to be unity. The least squares best fit line (dashed lines) has a slope of 1.03 in the figure. The inset on the upper left corner shows the variation of mean square displacement of the center of the cluster as a function of time whose slope gives the diffusion coefficient D defined in Eq. ( 3 ). The inset on the lower right corner shows the variation of (〈 r 2 〉)/ t as a function of time t which is expected to be a constant at long time as seen in Eq. ( 3 ) which also gives the diffusion coefficient.

We now consider the diffusion coefficient defined in Eq. ( 3 ) as a function of the size of the cluster at various temperature values. The diffusion coefficient is found from the slope of the 〈 r ² 〉– t plot by fitting the data to a straight line. The theoretically predicted dependence of diffusion coefficient on the cluster size N is given in Eq. ( 4 ). The diffusion coefficients of clusters of various sizes at four different temperatures are found from the slopes of the mean square displacement of the center of the cluster as a function of time. The log–log plots of the diffusion coefficients found in this manner are shown as a function of N in Fig.  3 at four different temperatures. The best least squares line fits to the data are also shown in the figure. The slopes of the curves ( α values) are as follows: –0.64 (170 K), –0.67 (180 K), –0.73 (200 K), and –0.65 (220 K). These values are equal to neither the –3/2 value expected for the displacement of the center of mass of the cluster due to the peripheral motion of particles, nor –1/2 value expected for the motion of the center of mass due to the attachment/detachment of the particles to/from the cluster edge. Rather these values are intermediate between the expected exponent values. Therefore one can assume that both kinetic processes operate in the mobility of the cluster. The value of exponent α at 200 K deviates from the others as can be seen in Fig.  3 . There are two processes that compete in the motion of the center of the cluster: the peripheral motion of particles and the attachment/detachment of particles to/from cluster edge. Depending on which of these processes is more dominant compared to the other, the exponent α would be expected to vary in the range of [–3/2, –1/2]. The deviation of α at 200 K from the others should be judged from this perspective.

Fig. 3.

	Figure Option View Download New Window
	Fig. 3. The log( D )–log( N ) plot of the diffusion coefficient as a function of cluster size at various temperatures. The slopes of the curves are as follows: –0.64 (170 K), –0.67 (180 K), –0.73 (200 K), and –0.65 (220 K).

Next, we consider the effect of peripheral motion of particles along the cluster edge on the mobility of the clusters. The energy barrier to diffusion along the cluster edge is chosen this time to be much bigger than both the binding energy E _b and the surface diffusion energy E _d . This choice prohibits the motion of particles along the cluster edge and therefore the motion of the cluster occurs only by attachment/detachment of the particles to/from the cluster edge, the so called evaporation/condensation mechanism. The log–log plots of surface diffusion coefficient D as a function of cluster size N for this case is shown in Fig.  4 for two temperature values. The exponents for α calculated for the present case are as follows: –0.53 (180 K) and –0.50 (200 K). In Fig.  4 , the E _ds is chosen to be approximately 10 times bigger than E _b and E _d . The value found for the exponent α in Eq. ( 3 ) is very close to the predicted value of –1/2.

Fig. 4.

	Figure Option View Download New Window
	Fig. 4. The same as Fig.  3 but here E ds ≫ E d , E b , hence the diffusion of particles along the cluster edge is almost prohibited. The slopes of the curves are as follows: –0.53 (180 K), –0.50 (200 K).

The dynamics of clusters has been studied in Refs. [ 20 ] and [ 21 ] from different perspectives. Kodambaka et al. ^{[ 20 ]} studied the role of step energy anisotropy in the motion and final equilibrium shapes of two-dimensional clusters. In addition to determining the orientation-dependent step-free energies and stiffnesses from the analysis of anisotropic island shapes and time-dependent fluctuations, they also determine the final shape of clusters after coalescence of two-dimensional clusters on surfaces by a combination of analytic and numerical methods. Khare and Einstein ^{[ 21 ]} studied the motion of vacancies of a certain radius and adsorbed atom islands and found that the diffusion constant scales with the radius of vacancy. Our findings on the dependence of cluster diffusion constant on the number of particles in the cluster is similar to their findings of its dependence on the radius of vacancy.

Finally, we show the behavior of the diffusion coefficient D as a function of the cluster size (number of particles in the cluster) in Fig.  5 . We also show a fit of Eq. ( 4 ) to the same data at two temperatures in the same figure. The fits appear to be satisfactory. For small clusters (2–12 atoms in a cluster) the diffusion coefficient is found to be oscillatory. ^{[ 3 , 11 ]} This is an expected behaviour because for small clusters the geometry becomes important. For example a cluster with three atoms may be expected to be more mobile than a cluster with four atoms on a square grid because four atoms may complete more incomplete bonds. We did not choose small clusters in our study so such a behavior does not arise. The mobility of the clusters is simulated in a lattice model in Ref. [ 7 ]. They found that the diffusion coefficients for small clusters ( N < 100) and larger ones depend on the number of clusters differently where they found that the exponent is α ∼ −1.5 for small clusters and it is α ∼ –0.5 for large ones. The clusters we have chosen are large ( N > 100) and we did not observe such a behaviour in our simulations. The values of diffusion coefficient found in our simulations vary an order of magnitude as the cluster size increases. For example at T = 200 K, if one assumes a lattice size of a = 3 Å, the diffusion coefficient for a cluster of size N ∼ 200 is about D = 6 × 10 ⁻¹⁴  cm ² /s, and it is about D = 7 × 10 ⁻¹⁵  cm ² /s for a cluster of size N ∼ 3500 particles.

Fig. 5.

	Figure Option View Download New Window
	Fig. 5. The variation of diffusion coefficient D as a function of the number of particles N in the cluster for T = 200 and 220 K (see Fig.  3 ). A fit of Eq. ( 3 ) to the data points is also shown. The fits for T = 170 and 180 K are not shown for a clear representation of the results.

We have considered the evolution of a single layer cluster located on a square grid with size of 350×350 in units of lattice constant a . The evolution is studied as a function of temperature, cluster size, and the energy barrier to diffusion along the edge of the cluster. When both the particle attachment/detachment to/from the cluster perimeter and the diffusion of particles along the cluster edge are in effect, the exponent in Eq. ( 4 ) is between –0.64 < α < –0.75. The exponent is neither –3/2 as expected for the cluster edge motion only, nor –1/2 as expected for process of particle exchange of the cluster with the surface. When the motion of particles along the cluster edge is prohibited by assigning a large value for the energy barrier along the cluster edge, the exponent approaches the value of –1/2 as expected. The value of diffusion coefficient varies by an order of magnitude as the size of the cluster increases. Typical values found for the diffusion coefficient are of the order of 10 ⁻¹⁵ –10 ⁻¹⁴  cm ² /s.

The initially circularly shaped clusters evolve to a square-like shape as they approach the equilibrium shape (see Fig.  1 ). This is mainly due to the square lattice chosen in the study. If one uses, for example, a triangular lattice instead, one would expect to obtain a more rounded equilibrium shape of cluster with more corners. This situation will affect the probability of motion of particles at the perimeter of the cluster because of the change of geometry but one would not expect the overall behaviour to change considerably.