Computer Physics Communications ELSEVIER
Computer Physics Communications 97 (1996) 101 - 105
Computer simulation of dynamics of surface growth Andrzej Z. Maksymowicz a, Mafia Magdofi b, Jeremy S.S. Whiting c a Faculty of Physics and Nuclear Techniques, University of Mining and Metallurgy AI. Mickiewicza 30, 30-059 Krak6w, Poland b Department of Physics and Computer Science, Pedagogical University Ul. Podchorq~ych 2, 30-084 Krak6w, Poland c Department of Physics, Unioersity of York York Y01 5DD, UK
It is shown that limited surface diffusion may strongly modify properties of the film surface. For linear samples, when some analytical results are known for random deposition, the fractal dimension D of the surface is discussed. Model predictions and computer simulation shows that D < 1 for the static case, when assembly of clusters is scanned. The surface is smoother for positive binding energy between atoms, and more rough otherwise. With zero binding energy, analytical results of random growth is recovered.
A surface may be seen as a boundary between different phases in d-dimensional space. The topology of a surface, its dynamics, dimensionality and many other aspects may be considered. Crystal formation, sedimentation, biological cell growth or division of living cells, wetting phenomena and electrodeposition are a few examples of the growth of the surfaces. This, of course, corresponds to a great variety of models, approaches (analytical calculations, numerical calculations, computer simulation, etc.) and approximations. It covers a vast area of topics and there are hundreds of publications devoted to the subject of surface growth. For a review, see for example . Some exact results are known which can be compared with results obtained from computer simulation, thus testing the simulation technique. Simulations may be applied as a complementary method which is sometimes still efficient while analytical treatment fails. Here we confine ourselves to thin films obtained by evaporation onto a flat substrate. A full description of the surface morphol-
ogy is given by the height function h(r,t), the number of atoms located on site r at time t. However, such information is neither necessary nor available. Instead, we characterize the surface in terms of some parameters such as film thickness h, surface roughness described by its fractal dimension D, height correlation function F(s) between different sites s-apart or correlation length Q, then surface thickness 3. Film thickness h is defined as the average over sites r, h-~ ( h ( ~ ' , t ) ) .
The roughness of the film surface may be characterized by its fractal dimension D from the postulated power law between surface length L and base length B (definitions of B and L are given in the next section) of the clusters formed in the deposition process, L = c - B ° + L o,
where, as it is usually assumed, the residual L o = O. For flat surfaces D = 1.
0010-4655/96/$15.00 Copyright © 1996 Elsevier Science B.V. All rights reserved. PII S 0 0 1 0 - 4 6 5 5 ( 9 6 ) 0 0 0 2 4 - 0
A.Z. Maksymowicz et a L / Computer Physics Communications 97 (1996) 101-105
The surface roughness can also be interpreted in terms of the correlation length Q. It is defined from the standard pair correlation function F(-~) of the local thickness h(r-9, F ( s ' ) = ( h ( 7 ) . h( ~'+ 3) ) - h. h,
where the shift vector ~' is in the film plane. It is seen that the correlation function F(~') is zero for randomly distributed atoms. It is negative for rough surfaces and positive for flat ones. The correlation length Q is then the average of ~" with weights F(-g'), Q = ~] s-F( s~.
The fluctuations 6 in film thickness around its average value h may be interpreted as the thickness of the surface itself, 62 = F ( 0 ) = ( h( 7 ) - h(~') ) - h . h .
Each of the parameters D, Q and 6 may be thought of as a function of the film deposition time t or average film thickness h. The correlation function F(E) and correlation length Q also may depend on the direction ~' for anisotropic surfaces which may be the case for non-perpendicular incidence of atoms. In principle, all parameters are of a statistical nature and one can distinguished two typical approaches which lead to different values. For example, for the fractal dimension D we may assume a static picture of the surface at a given time t, which corresponds to a given film thickness h, then we fit the power law (2) for the assembly of clusters. Alternatively, in the dynamical case we watch one cluster growth, and again fitting the power law for the L(t), B(t) dependence of this cluster. We confine ourselves to d = 1, the one-dimensional case when some analytical results are known for random deposition. We concentrate on fractal properties of the film surface and discuss the effect of surface diffusion on roughness.
is a formal definition• The base line, which is the substrate, is divided into unit cells randomly filled with square particles. The ith cluster itself extents between empty cells. The base length B i = 1,2,3 .... is defined as the distance in the geometrical sense along the base of a cluster• Then the surface length L i > B i is measured along the cluster surface, and so L i is the number of all the sides we pass when travelling over the surface between the two points. For the smallest cluster which is one particle, B = 1 and L = 3. For B = 1, and two atoms sitting one on top of the other, L = 5. Each cluster is characterized by two numbers, B and L, which represent a point in the (B, L) plane. Fitting the exponential law, one gets the dimension D from the computer experiment. As it was mentioned, in the l d case and for random deposition the analytical results for D are known . For a given B, and for a given configuration of m~, m 2, m s . . . . . m s particles at each site in the cluster, we may work out the formulae for L(B, ml, m2, m 3 . . . . . m s ) . There are many possible configurations m t , m 2, m 3 . . . . . m B and we may introduce their distribution function p ( B , m~, m 2, m s . . . . . m e ) . The next step is to replace the m ' s by their average values, (m) = E P(B,m). L(B,m),
which corresponds to the operation of replacing all the points along the L-axis of a given B by a single value point in the ( L - B ) plane. In other words, we reduce the scattered points to a L(B) line. Within assumed average number of particles per cell h, and for random deposition, one gets after some simple yet tedious calculations, L(B)=[1-1(2h).e
• [l(2h)-e -2h-
where 2. Fractal dimension for 1D surfaces
J ( h ) = (1 - e - h ) 2 / 2 h
Let us consider the 1D case and discuss the fractal dimension D which may be extracted from the L(B) power law dependence, L = c- B D, from which O = ( dL/dB)/(L/B)
and I(2h) is expressed in terms of modified Bessel functions I,,
l ( 2 h ) = 10(2h) + l l ( 2 h ) .
A.Z. Maksymowicz et al. / Computer Physics Communications 97 (1996) 101-105
The fractal dimension D is then obtained from the formal definition (6). In fact, we need to specify how the definition is to be applied. D depends e x p l i c i t on film thickness h and the cluster size B. However, these two are related and as the film grows the average cluster size must also get larger. For random deposition we have clusters of the average base length (see the appendix) ( B ) = e h, a monotically increasing function of h, as expected. In other words, for given h most clusters are of size B such that h = l o g ( B ) . T h e idea is to replace h in the L ( B ) dependence in (8) by the corresponding l o g ( B ) . T h e fractal dimension D is then always larger than one and tends to unity for large h. The asymptotic behaviour is
7r-~/4 h ~
This result is very close to the one obtained by Tanaka  who got similar result using the saddle point method. His result becomes identical with form (11) after replacing ~" with e. Computer simulations on growing clusters [3-5] confirms that D is always larger than unity, decreases with h, and tends to one for thicker films. This was a dynamical case when we watch one growing cluster in a process of deposition, thus obtaining the time sequence of L (t) a n d B ( t ) for this cluster, to which we fit D. Alternatively, one may consider the static case when all clusters are scanned at a given time. It means that the thickness h is fixed. From our model calculations (6), and L ( B ) given by (8), we rigorously expect D < 1, contrary to the dynamical case. Some computer simulation [6-8] support the predictions and agreement between numbers obtained from analytical calculations and computer simulation is remarkably good . For small h we get 1 4h D = 3 + 9
54 + . . . .
3. Example: limited surface diffusion model In order to get closer to real systems, it is necessary to account for atomic ordering and diffusion processes. The process is assumed of two consecutive steps. First, each particle falls vertically (no oblique direction) at a random site. If the particle sticks here, we get the simplest but most unrealistic random deposition model. In the next step we allow for surface diffusion. In the surface diffusion model proposed by Family , the particle diffuses until it reaches a local minimum heights. In this paper we assume limited diffusion when the particle may move to the nearest neighbour and come to rest there. The mobility is controlled by a diffusion energy barrier E so that the probability of this move is proportional to the Boltzmann factor e x p ( - E / T ) , where T is temperature. For large E the particle sticks where it has landed and we believe to recover the standard 'hit and sit' random deposition model. For the other limit, E = 0, the particle may move left, stay where it landed or move right with same probability 1 / 3 unless we account for a preferential anisotropic diffusion. Such anisotropy may be related to a surface tension connected with a difference in the number of particles on neighbouring sites. With n h = number of particles at the hit site, before our particle gets fixed, n n -- n t, n r, the number of particles in the nearest, left or right cell, we assume the simplest rule that the particle gains an energy V if the jump takes place to the nearest neighbour site with n, > nh, and so the probability of such a jump is again proportional to the Boltzmann factor e x p ( - V J T ) , where V~ = V • s i g n ( n , - nh). ThUS the binding energy V is another model parameter which describes surface tension leading to the wetting phenomena . Combining the two effects, mobility and surface tension, we assume that the probability p(left,stay, right) that the particle moves left, stays or moves to the right, is given by p(left) = A . e ( - E - V,)l r),
p(stay) = A- e (°),
p(right) = A . e ( - e - v , ) / r ) ,
and in the dense limit D=
l _ e-h~'-~
with the proper normalization factor A. In simulation we applied the E = 0 limit to allow for maximum
A.Z. M a k s y m o w i c z et al. / Computer Physics Communications 97 (1996) 1 0 1 - 1 0 5
mobility and a range of different V to see its effect on the surface roughness. The random deposition model corresponds to large E.
4. Evaluation of the memory required In order to achieve the necessary accuracy for D of the order of say 3 decimal digits, just above the accuracy we may expect from an experiment, it is necessary to account for a minimum of at least 10 6 or so clusters, nCluster. For a given number of cells nCell we estimate the number of clusters (see the Appendix) nCluster = nCell, e -h • ( 1 -
Table 1 Dependence of the fractal dimension D on the average number of atoms h for different values of the binding energy V and cut-off minB = 2
0.01 0.04 0.10 0.50 1.00 1.50 2.00 3.00 5.00 10.00 12.00
V = -2
Error for D
0.54 0.54 0.57 0.63 0.673 0.700 0.724 0.758 0.809 0.884 0.904
0.54 0.56 0.59 0.68 0.755 0.797 0.826 0.872 0.939 0.998 1.001
0.55 0.56 0.59 0.72 0.827 0.881 0.914 0.960 1.000 -
0.2 0.02 0.01 0.01 0.001 0.001 0.002 0.004 0.01 0.15 0.5
for h atoms per cell on the average and for randomly distributed atoms. As expected, in both limits of a few atoms or very thick films, nCluster goes to zero and we have to increase the nCell significantly to reach the assumed nCluster. It is only n e a r e h --= 0.5 (or h = 0.69) that nCluster is comparable with nCell. With that, if we intend to scan clusters over a wide range on h below the percolation limit, say h = 0.01, to a dense limit of h = 10 or 20, it is quite easy to run with a memory of several hundred Mbytes. Perhaps it is reasonable to recall that the average size of a cluster ( B ) = e h, and so it is sensible to declare the nCell and B variables as long for larger h.
a smoother surface with smaller D. In each case D
5. Results and summary First, we recovered the random deposition results for the surface diffusion model with binding energy V = 0, as expected. That is, results of computer simulation are independent on whether we use a large diffusion energy barrier V, when the atoms are truly fixed at the landing position, or a small V when atoms may randomly move to the left or to the right. In Table 1 we collected results of our simulations for the static case of D(h) dependence for different values of V. As we mentioned, the V = 0 case agrees reasonably well with the calculated values. Case V = 2 produces a more flat surface, as expected. On the opposite, negative binding energy V = - 2 gives
Table 2 Dependence of the fractal dimension D on the binding energy V for different values of cut-off, minB = 1, 2, 4 and 10. The film thickness h = 1.5. V
minB = 1
minB = 2
minB = 4
minB = 10
-2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0
0.571 0.610 0.649 0.688 0.722 0.753 0.779 0.801 0.817
0.700 0.725 0.747 0.772 0.797 0.822 0.845 0.865 0.881
0.833 0.838 0.846 0.854 0.863 0.875 0.887 0.900 0.909
0.923 0.899 0.927 0.924 0.926 0.930 0.933 0.942 0.944
A.Z. Maksymowic z et a L / Computer Physics Communications 97 (1996) 101-105
and should be continued in one of the following directions: • simulations of the dynamic behaviour in which we observe one growing cluster; • the case when the film substrate is treated differently from the deposited atoms, a necessary step if aggregates are to be investigated on the monolayer-or-so scale; • going over from the one-dimensional case (when theoretical predictions at least for random deposition exist) to higher dimensions study of surface morphology for thicker films, above the percolation limit when all substrate is covered by a single cluster; • study of anisotropic correlations due, for example, to oblique incident of atoms or epitaxy; • further on we may simulate deposition of binary alloys with atomic order, to obtain concentration profiles or concentration inhomogeneities.
from which the average cluster size is oo
= e h=- 1 / p .
For an average number h of particles per cell, and for a random deposition of the particles among the cells, we apply Poisson distribution for a number k of particles in a cell, p ( k ) = e - h . h k / k ! , k = 0,1,2,3 . . . .
in the continuous limit of large number of cells nCell. Then p = p(0) = e "h is a probability that the
cell is empty, and with probability (1 - p ) the cell is occupied, k > 0. To find a cluster of base length B = 1,2,3 . . . . . we claim the situation when, starting from the first occupied site, we accept a sequence of the remaining ( B - 1) sites also occupied (1 - p is the probability) followed by the terminating empty cell (with probability p). This leads to the geometrical distribution p(B) for a cluster of length B, p ( B ) = e - h • (1 - e - h ) (1~- 1),
Now we may estimate the number of clusters nCluster. We expect a small number of clusters in
the dilute limit, most of them are of size B = 1. In the dense limit, the number of clusters must also vanish since nearly all cells are occupied and individual clusters merge. If we consider all the clusters on the nCell long base line, then the total number of cells the clusters occupy is nCluster. < B > . The ratio of the total number of occupied cells, nCluster. < B > , to the number of cells nCell must be the probability (1 - p) that the cell is occupied. So we have 1 - p = nCluster. ( B ) / n C e l l ,
from which nCluster = nCell, e - h . ( 1 -- e - h ) .
For the diluted case h ~ O, nCluster = n C e l l , h, as expected. In the dense limit h ~ 0% the number of clusters vanishes exponentially, n C l u s t e r = nCell. e "h. The maximum of nCluster is reached for an intermediate region e h = 1/2.
References        
M. Kotrla, Czech. J. Phys. 42 (1992) 449. T. Tanaka, J. Phys. Soc. Japan 59 (1990) 3898. M. Nakamura, Phys. Rev. B 40 (1989) 2549. M. Nakamura, Phys. Rev. B 40 (1989) 3358. M. Nakamura, Phys. Rev. B 41 (1990) 12268. A.Z. Maksymowicz, J. Rafa and K. Zakrzewska, Computers in Physics 6 (1992) 317. A.Z. Maksymowicz, J. Rafa and K. Zakrzewska, Phys. Rev. B 50 (1994) 8901. A.Z. Maksymowicz, M. Magdofi, J.S.S. Whiting and I.L. Maksymowicz, paper presented at First School of Large Scale Computer Calculations, 21-23 February 1994, Krak6w. F. Family, J. Phys. A 19 (1986) L441. P.G. de Gennes, Rev. Modern Phys. 57 (1985) 827.