High Performance Computing on Hewlett-Packard Systems


Computer Simulation of Anisotropic Thin Film Growth

A.Maksymowicz1,*, K.Malarz1,**, M.Magdon2,***, S.Thompson3 and J.Whiting3

1 Department of Theoretical and Computational Physcics (ZFTiK),
Faculty of Physics and Nuclear Techniques (WFiTJ),
University of Mining and Metallurgy (AGH),
al. Mickiewicza 30, PL-30059 Krakow, Poland

2 Department of Mathematical Statistics,
Faculty of Agriculture,
Agricultural University,
al. Mickiewicza 21, PL-31120 Krakow, Poland

3 Department of Physics, University of York,
York YO1 5DD, England

e-mail: *amax@agh.edu.pl, **malarz@agh.edu.pl, ***rrmagdon@cyf-kr.edu.pl
*phone: +48 12 6173539, fax: +48 12 6340010


Computer simulations, based on a simplified model of thin film grown by deposition of atoms, and using Monte Carlo method, are carried out to study morphology of the surface. We consider how binding energy of atom pairs J and diffusion energy barrier V influence the anisotropic thin film growth. The morphology of the surface is described by its roughness, dispersion of film height and ellipticity. Preliminary calculations and tests of prepared computer programs were performed on HP715 machine, then for the fully operational version on Convex Exemplar SPP. Typical running time for a sample of 500*500*16 deposited atoms was about a couple of hours during regular sessions, CPU time was 3 min, and the required memory size was about 100 MB. The model, implemented algorithm, and results of numerical calculations are presented and discussed.


1. Introduction

The growth of surfaces is considered in many theoretical models [1] using different analytical and numerical approaches and approximations. Exact solutions are known only for some limiting cases, for which the computer simulations may be tested. Then we may analyse the problem in question by computer experiment with some confidence, and compare results of calculations with experimental data [2].
In this paper, which is a continuation of our previous investigations [3, 4] of anisotropic growth of thin films, we intend to study the influence of limited surface diffusion on morphology of the surface, against the migration barrier energy V and the local binding energy of particles J. To characterise roughness, we define thickness of the surface simply as the dispersion sigma of film heights z=h(x,y), over the film plane (x,y), then h-h correlation function F between nearest neighbours. We define the correlation as configurational averages <...>,

Fs = <(h i - < h >)(h i+s - < h >)>, EQ.1.

where h i = h(x,y) and < h > denotes the average film thickness, then s is the vector from the atom i, localised at co-ordinates (x,y) to its nearest neighbour. Dispersion sigma2 = F corresponds to s=0. The correlations Fx, Fy allows us also to describe the possible anisotropy of film growth in x- or y-direction, as provoked either by anisotropic diffusion barriers Vx, Vy or by the anisotropy in binding energies Jx, Jy. Then Fx != Fy and an ellipticity parameter e, defined as
e = (Fx - Fy)/sigma2, EQ.2.
may be used as a measure of the anisotropy. For random deposition we expect sigma = < h >1/2, where < h > is the film thickness, and F = 0. For flatter surfaces sigma is smaller and correlation function F is positive. In this paper we systematically study the growth of film surface and its morphology in terms of dispersion sigma, correlations functions Fx, Fy, and ellipticity e. The Monte Carlo method and computer simulations for different sets of parameters of the proposed model are used to analyse how binding energy of atom pairs J and diffusion energy barrier V influence the anisotropic film surface growth.

2. Model

Each particle is represented as a unit volume cubic, which is randomly placed at co-ordinates (x,y) on two-dimensional n*m lattice. If two cubics have one common wall we call them nearest neighbours. We assume limited surface diffusion when the particle may stay at randomly chosen place (x,y,z(x,y)) or diffuse to one of its nearest neighbours position in one of the s directions to minimise the local energy in terms of the binding energy. Thus we have two factors which determine the particle behaviour, a tendency to find the energetically most suitable site in vicinity of (x,y) site, then some diffusion barriers Vs that prevent particles from too easy migration.
The mobility of particles is controlled by the diffusion energy barrier Vs, which may be anisotropic when it depends on the direction s. The probability of such a jump to a neighbouring site is proportional to the Boltzmans factor:

p~exp(-Vs/kBT), EQ.3

where kB is Boltzmans constant, T is temperature. For large Vs the particle stays where it landed and we expect the standard random deposition model results. For the other limit, when Vs, the particle may move to one of four places: to the "left", "right", along x-axis or "up", "down", along y-axis or it may stay at this site which was randomly chosen. Thus we have 5 possibilities corresponding to 5 sites where the particle may finally be fixed on the lattice: (x,y,z), (x-1,y,z), (x+1,y,z), (x,y-1,z), (x,y+1,z), where (x,y) describes the randomly chosen place and z is taken at the actual site where the particle is frozen. Each of the possibilities takes place with the same probability p=1/5 for random deposition, when Vs = 0, while for the opposite limit, Vs -> infinity, probability of any jump is zero and so particles are fixed at the hit location.
The binding energy Jsij for atom pairs (i,j), where s is direction of atom j, is the diffusion causative force which makes system tends to the minimum of energy. Therefore, for nonzero binding energies, Jsij != 0, the final distribution of particles is not random because of the possible migration to the nearest site of lower energy. We assume that the energy gain Ei for the particle settling at a given i site is proportional to the number n(i,j) of particle-particle nearest neighbour pairs (i,j) linked to the given site in question,
Ei = sumj(n(i,j)Jsij), EQ.4
where s is the direction of atom j. So the energy gain Ej depends on the local environment of the central site i and also its next-nearest-neighbours. Again, the probability of picking up one to neighbouring site j is given by the Boltzmans factor
p~exp(-Ej/kBT), EQ.5
where Ej is the energy of the proposed neighbouring site j.
In this simplified picture we take into account both factors, mobility and surface tension. We assume that the probability p(u,v) of a particle falling on site u and settling finally on one of the sites v (to the "left", "right", along x-axis or "up", "down", along y-axis) is given by
p(u,v) = C*exp(-Ej/T) if v = u and p(u,v) = C*exp(-Ej-Vs/T) if v != u, EQ.6
where the normalisation constant C must satisfy the equation: sumv p(u,v) = 1. For the limited surface diffusion model, a random number generator is employed three times per deposited each particle. First two drawings indicate the hit location (x,y)=u and the third one serves the decision making on the diffusion process. We simulate the thin film growth by the random distribution of particles on two-dimensional n*m lattice. To discuss results of computer simulations we use the standard deviation sigma of the height average < h >, auto-correlation function F and s-correlation functions Fs for both two directions s=x,y, described in introduction. The ellipticity parameter e, defined by eq.2, is used as a measure of the anisotropy in surface morphology.

3. Numerical results

Numerical simulations of thin film growth were performed for two dimensional n*m = 500*500 lattice where n*m*h atoms were placed. Here we consider only nearest neighbours interaction Jsij = Js which still, however, may depend on s direction. The morphology of the surface is studied for three characteristic cases:

We used HP715 machine to test in details the prepared algorithm and computer programs in series of preliminary calculations.

The basic flowchart of the algorithm is described below. The final version of the program was prepared and implemented on Convex Exemplar SPP. For numerical simulation of thin film growth with typical values of model parameters and sample of 500*500*16 deposited particles the required memory size was about 100MB, CPU time was about 3 minutes. The main results of our calculations are discussed in paragraphs 4.2 - 4.4.

3.1 Basic flowchart of the algorithm

  1. Declare all variables.
  2. Import initial data from command line and/or from file.
  3. Allocate memory, initialize tables and variables.
  4. Prepare substrate of fixed n*m bottom layer atoms.
  5. Pick randomly the hit (x,y) site on the n*m lattice.
  6. Calculate local energy at the site (x,y,z), where z=h(x,y), and its vicinity.
  7. Draw the coin and decide where the atom freezes.
  8. Update information and repeat steps 5, 6 and 7 until all atoms were deposited, then go to step 9.
  9. Export results to screen and/or output files.

3.2 Testing chaotic surface growth

First, we recovered the random deposition results for the surface diffusion model with zero binding energy, Js = J = 0, as expected. In particular we verified the known result that roughness sigma is a linear function of square root of average film height < h >, sigma = < h >1/2.
Fig.1 shows results of our computer simulations which are represented by full squares, whereas the solid line is given by function sigma = < h >1/2. The agreement with theoretical model is excellent. It should be added that this result is obviously independent on the diffusion energy barrier Vs. The same results of chaotic surface growth were obtained for the limiting case of large value of V, V -> infinity, and finite values of binding energy J.
FIG.1. Dependence of roughness sigma on square root of average film height < h >1/2 for random deposition.

3.3 Surface roughness for isotropic growth

Here we assumed isotropic diffusion energy barrier, Vx = Vy = V, and recorded the influence of isotropic binding energies Jx = Jy = J on the film roughness. Numerical calculations were carried out for n*m = 500*500 lattice and n*m*h deposited atoms with h=16. Results for sigma(J) dependence obtained for V=0 are presented on fig 2. Fig.3 shows the correlation Fx(J), also for V=0. Values of binding energy J and barrier diffusion V are in units of kBT.
As it was expected we found that positive values J, J > 0, generate rougher surfaces. The roughness of the film surface increases rapidly from value sigma = < h >1/2 (for J=0) reaching some saturation for large J. This effect is also clearly seen in the behaviour of the height correlation function Fx (or Fy) which absolute values increase with increasing J. Negative binding energy, J < 0, produce smoother surfaces; in this case roughness of the film decreases as it is shown on fig 2. Again, absolute values of height correlation function Fx increases for smaller J, see fig 3. In each case, the roughness sigma is insensitive to binding energy J for the limits of large positive or negative values.
FIG.2. Dependence of roughness sigma on binding energy J (in units of kBT) for diffusion barrier V=0.
FIG.3. Dependence of the height correlation function Fx on binding energy J (in units of kBT) for V=0.
For nonzero values of barrier energy, Vx = Vy = V > 0, the influence of binding energy J on film roughness sigma diminishes. This is so until large values of V prevent any diffusion, i. e. up to the limit of random deposition.

3.4 Anisotropic growth

Anisotropic thin film growth takes place, when at least one of the two considered energy contributions, i.e. diffusion barriers or/and binding energies had different values for x- and y- directions. It is obvious from eqs. (4) - (6) that either anisotropy in the diffusion energy barrier, Vx != Vy, or anisotropy in the binding energy, Jx != Jy, directly influence the distribution of particles. As result of computer simulations, we get different morphology of the surface observed in different directions. In our investigations we systematically studied each of these two possible sources of anisotropic growth. Anisotropy is measured by ellipticity parameter e, defined by eq.2. Numerical results are presented in fig.4 and tab.1.
First, we put zero diffusion barrier, V=0, and analysed how anisotropy of binding energies, Jx != Jy, influences the film growth. As it is shown on fig 4, the dependence is evident. Ellipticity parameter e strongly depends on the differences between binding energies in x- and y- directions, deltaJ = Jx - Jy.
The similar dependence was obtained if we took into account nonzero and isotropic values of diffusion barrier energies, Vx = Vy = V > 0. As it is expected, the influence of anisotropy of binding energies on ellipticity e is then less pronounced, and in the limit of large barrier energy V we recover isotropy e=0.
FIG.4. Dependence of ellipticity e on binding energy difference deltaJ for diffusion barrier energy V=0.
The influence of anisotropy in diffusion barrier energies, Vx != Vy, on ellipticity is also seen. In order to verify our expectations we assumed isotropic binding energies, Jx = Jy = J, zero diffusion barrier energy along x-direction, Vx = 0, and we scanned the diffusion barrier energy along y-direction, Vy. Again, we observe anisotropy from the e(deltaV) = e(Vy) dependence. Results of numerical simulations, carried out for n*m = 500*500 lattice and h=16 (the same dimension of lattice and its height as in previous calculations) are collected in Tab 1. Values of binding energy J and diffusion barrier energy along y-direction Vy are in units of kBT .

TAB.1. Dependence of ellipticity e on diffusion barrier difference deltaV = Vy - Vx for vx = 0, and different values of isotropic binding energies Jx = Jy = J.
Vx = 0J=-2J=-0.5J=0.5J=2
deltaV=Vye(deltaV)e(deltaV)e(deltaV)e(deltaV)
0-0.0010.002-0.0020.001
0.250.0090.018-0.005-0.006
0.50.0160.041-0.017-0.009
10.0340.094-0.031-0.021
50.3320.443-0.082-0.080

It is seen that the anisotropy of diffusion energy barrier produces anisotropic film growth measured by ellipticity parameter e. The dependence of ellipticity on diffusion barrier difference, e(deltaV), is different for positive and negative values of binding energy J and its influence on surface morphology is pronounced. For smaller values J < 0 we get greater absolute values of e. For large and positive values of binding energy J the role of anisotropic barrier diffusion is greatly reduced.

4. Summary

The main aim of this work was to study the influence of limited surface diffusion on anisotropic thin film growth, as seen in computer experiment. Numerical calculations and computer simulations were based on a simple model of surface growth. In the model discussed the diffusion is controlled by the barrier energy V and by the binding energy J between nearest neighbours. Morphology and anisotropy in the film growth are given by roughness sigma, the spatial height-height correlation functions Fs and ellipticity parameter e. We systematically studied chaotic surface growth, isotropic and anisotropic growth. Main results are presented in figs. 1 - 4 and in tab 1.

Bibliography

  1. A.C. Levi and M. Kotrla, J. Phys.: Condens. Matter 9(1997) 299
  2. N.-H. Cho, K.M. Krishnan, C.A. Lucas, R.F. Farrow, J. Appl. Phys. 72(1992) 5799
  3. A.Z. Maksymowicz, M.S. Magdon and J.S.S. Whiting, Computer Physics Communications 97(1996) 101
  4. A.Z. Maksymowicz, J. Suliga, K. Kulakowski and M.S. Magdon, Proceedings of the 8th Joint EPS - APS International Conference on Physics Computing (1996) 121

AGH Copyright K.Malarz, 1997