High Performance Computing on Hewlett-Packard Systems


Structural Correlations in Growing Thin Film Surfaces

A.Maksymowicz1,*, K.Malarz1,**, G.Nagel1 and M.Magdon2,***

1 Department of Theoretical and Computational Physics (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

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


In this paper we study the influence of binding energy and energy diffusion barrier on surface morphology of thin film growth. The random deposition model with limited surface diffusion, and its computer simulation is considered. Anisotropic growth, due to either binding energy or diffusion barrier anisotropy, is also analysed. Main results in terms of the height-height correlation function, which discloses the anisotropy, are discussed. The numerical calculations were carried out on Convex Exemplar SPP/1600XA-32. Timing, evaluation of memory required and computational complexity are also discussed.


1. Introduction

Growth processes are very common in nature and they also take place in almost all of modern technologies. This implicates hundreds of publications devoted to the subject of surface growth (see [1], [2] for a review).

In this paper we present a simple model for anisotropic thin film growth. Surface roughness, controlled by the particle-particle interaction and by the diffusion, is measured by the standard deviation sigma of film thickness. The value sigma=0 is for perfect wetting and sigma>0 for rough surface. However, this single-site characteristics is insufficient for anisotropy consideration. Then a two-site properties such as correlation function F, which will be defined in the next section, are necessary. Such a function calculated in different directions may be used as a measure of anisotropy.

Results of simulations for sigma and F will also be compared with analytical results for some limiting cases when the exact solution are known.

2. Model

In the simplest yet rather unrealistic case of random deposition model we can assume that particles, represented by a unit volume cubics simply fall vertically until they reach the top of the column where they were dropped or they reach the substrate. At this site they freeze and become part of the aggregate. The substrate is presented by two dimensions square lattice with lattice. This model is oversimplified in the sense that it ignores correlations between the column heights as it the case for truly random processes that follow Poisson distribution. However, this limiting case serves as a test for the model and proposed computer algorithms.

In a more realistic approach, local surface migration must be included. The particle after landing may diffuse on the surface within a finite distance to find the column with the minimal height [3], or for one offering the strongest binding (i.e. the most occupied neighbouring sites [4]). Such a model should be adequate for films grown by molecular beam epitaxy (MBE). In this work a model of random deposition with limited diffusion is applied. Particles may move one step to the nearest neighbour site and then they freeze. Other models devoted to MBE techniques can be found in [1].

Full description of surface morphology is given by the height function h(r,t), the number of atoms located on site r=[x,y] at time t. However, such information is neither necessary nor available. Instead, we define some parameters which describes surface morphology, such as:

where the shift vector s=[Delta x,Delta y] lies in the film plane [5].

The correlation function F(s) is zero for s not equal [0,0] for randomly distributed atoms. F>0 for smooth surfaces and F<0 for rough ones. Then the anisotropy coefficient
e:=F([0,1])-F([1,0]) EQ.4.
may be used to indicate anisotropic growth.

The single site correlation, s=[0,0], is dispersion
sigma2=Phi(0) EQ.5.
and also characterizes the surface roughness.

3. The Monte Carlo algorithm and rules for film growth

After picking-up, at random, initial position of the particle, r=[x,y] on (n times m) large square lattice, we consider the particle total energy in five possible final positions mu0..4={[x,y,h(x,y)], [x-1,y,h(x-1,y)], [x+1,y,h(x+1,y)], [x,y-1,h(x,y-1)], [x,y+1,h(x,y+1)]} (the volume diffusion [x,y,h(x,y)-1] and evaporation [x,y,h(x,y)+1] is not allowed in our model).

There are two components of the mentioned above particle total energy: diffusion energy barrier V which control the particle mobility (movement resistance), and particle-particle binding energy J, positive for attractive forces, and negative when particles repel each other. We assume that the probability of the particle to stay where it landed is given by the Boltzmann factor,
p(stay)= exp(-Jmu/kT) / summu exp(-(Jmu+Vmu)/kT), EQ.6.
and the probability of diffusion to one of the neighbouring site is:
p(left, right, forward, back)= exp(-(Jmu+Vmu)/kT) / summu exp(-(Jmu+Vmu)/kT), EQ.7.
where k is the Boltzmann constant and T denotes temperature.

The simulation of this system is performed using the standard Monte Carlo method. We used a (500 times 500) square lattice with periodic boundary conditions. The size of the lattice offers sufficient statistics to get accurate results, and also to minimize the boundary effects. The total number (N=500 times 500 times 16) particles were placed during each simulation.

4. Evaluation of memory required and computational complexity

The process of preparing and collecting data was divided into two independent stages. First we grow the thin film and the data are stored in memory (or on disk). Then these data are analysed to extract the relevant numbers that characterize the film and its surface. If we store information about the film in a binary file, for a film of average thickness h, we need (n times m times 4 times h) long array for sputtered (n times m times h) particles represented as integer. So for architectures where size of an integer variable is 4 bytes, we have for memory allocation and disc space about (500 times 500 times 64 times 4) bytes = 64 MB required.

Computational complexity is proportional to problem size N=(n times m times h). Main part of time of simulation is spent on evaluation of the probability of possible future movements. It means calculation the particle total energy in possible virtual five position. It requires N times (25 multiples + 5 adds + 5 powers). The alternative approach takes advantage of a pre-calculated table of particle energies for all possible local configurations, if there are not too many of them.

At this stage of simulation, one can easily adopt the parallelisation technique of computations, which may bring considerable time saving, as predicated the Amdahl law.

5. Main results

For given number of deposited atoms, we consider the influence of the control parameters J and V on surface morphology. As we mentioned earlier, this information, which also includes possible anisotropy of the morphology, must be represented by a 2-site characteristics such as the function F(s) of correlation of the local thickness h(r).

5.1. Roughness

At first we recovered known results for random deposition, among them surface roughness sigma2=Phi(0), which is in this case equal to h= < h(r) >,
sigma=sqrt(h), EQ.8.
as predicated by random deposition model, see tab.1.

TAB.1. Dependence of surface thickness sigma against films thicknes s h for random deposition, see eq.8.
< h > 1.0004.0009.00016.00025.00064.000
sqrt(< h >) 1.0002.0003.0004.0005.0008.000
sigma 0.9871.9942.9983.9945.0047.991

This case is adequate to the situation when binding energy J=0, then results are independent on diffusion barrier as it has been shown in tab.2, or when diffusion barrier V is large enough to prevent migration of atoms.

TAB.2. Test for independence of the roughness parameter sigma on the barrier diffusion V for random deposition, J=0.
V -50-10-2-10121050
sigma 4.04.04.0 4.0 4.0 4.04.04.04.0
FIG.1. Film heights profile for large negative values of J.
FIG.2. Film heights profile for random deposition model, J=0.
FIG.3. Film heights profile for big positive values of J.

We also found that for larger J>0 the surface is rougher (fig.3). In the other limiting case, when J -> - infinity, a perfect wetting take place (fig.1), for which sigma -> 0, as expected.

A characteristic saturation for sigma in the limit of large values of J, is observed in tab.3.

TAB.3. Roughness parameter sigma as a function of the binding energy J for V=0.
J -10.0-3.0-1.5-1.0-0.5-0.250.00.250.50.751.01.53.010.0
sigma 1.01.01.11.21.72.34.07.210.914.116.318.821.221.9

5.2. Anisotropy

The differences in function F(s) for s=[0,1] and [1,0], which corresponds to correlations along x- and y-axis, describes the anisotropy, e=F([0,1])-F([1,0]). Ellipticity e=0 for isotropic surface, then anisotropy manifests itself as e not equal 0.

The anisotropic growth of the film surface could be provoked either by anisotropic binding energies Jx not equal Jy, or by anisotropy in the diffusion energy barrier Vx not equal Vy for nonzero J's which also generates non-randomness. (Of course, we always may prevent any non-randomness for large V -> infinity when no migration takes place as it has been shown on tab.4.) As an example of the first possibility for anisotropy, we used Jx=-Jy and V=0 (tab.5), the latter case of Vx not equal Vy and J not equal 0 is presented in tab.6 and tab.7.

TAB.4. Ellipticity parameter e against diffusion barrier V, for Jx=-Jy=0.5.
V -50-10-2-1012351050
e -1.262-1.263-1.269-1.269-1.251-1.158-0.900-0.508-0.0930.0030.002
TAB.5. Ellipticity parameter e against binding energy Jx=-Jy, for V=0.
J 00.10.2511.5235
e -0.002-0.283-0.718-1.251-1.552-1.595-1.597 -1.597
TAB.6. Ellipticity parameter e against diffusion barrier Vx=-Vy=V, for J=+2.
V 00.511.522.5351050
e 0.0010.0110.010-0.010-0.046-0.091-0.146-0.336-0.479-0.475
TAB.7. Ellipticity parameter e against diffusion barrier Vx=-Vy=V, for J=-2.
V 00.511.522.5351050
e 0.000-0.034-0.069-0.114-0.167-0.227-0.289-0.395-0.037-0.004

We can see the effect of saturation of ellipticity parameter e as function of J and/or V. The limit J=0 corresponds to isotropic growth while the other limit J -> infinity to situation when all deposited particles may move only along one (easy) direction.

5.3. Spatial correlation

The sign of correlation function F(s) offers some information on surface morphology [5]. F<0 corresponds to rough spiky surfaces (see fig.3 and tab.9), F>0 produces terrace-like, smooth surface (see fig.1 and tab.8).

TAB.8. Structural correlation for J=-10.
range 12345
F 0.7150.4750.2750.1320.041
TAB.9. Structural correlat ion for J=+10.
range 12345
F -0.4720.0610.023-0.008-0.001

We check correlation F(s) for vectors u=[i,0] and v=[0,i] when i=1..5. Additionally we consider F(w) in direction w=[1,1]. We find F(v)=F(u)=F(w)=0 for random deposition, as expected. When diffusion is accounted for, we expect decreasing correlations with increasing distances i. For anisotropic surfaces, we get different distance dependence of F along x-axis, F(u) in tab. 10, and y-axis, F(v) in fig. 11.

TAB.10. Structural correlation along x-axis for Jx=-Jy=5 and V=-2.
range 12345
F(u) -0.7780.428-0.1860.062-0.010
TAB.11. Structural correlation along y-axis for Jx=-Jy=5 and V=-2.
range 12345
F(v) 0.8140.5870.3800.2220.117

The three typical situations of smooth, rough and anisotropic surfaces (which is smooth in x-direction and spiky in y-direction) are illustrated in figs. 4, 5 and 6. It should be noted that the diagonal F(w) correlation cannot distinguish the smooth and rough isotropic cases, in figs. 4 and 5, respectively, as in both cases F(w) is positive. On the contrary, anisotropic case in fig. 6 yields F(w)<0.

FIG.4. A sketch of a smooth surface given in fig.8.
FIG.5. A sketch of a rough surface given in fig.9.
FIG.6. A sketch of an anisotropic surface given in tab.10 and tab.11.

6. Conclusion

In this paper we recovered the random deposition results for the surface diffusion model with binding energy J=0, as expected. The results of computer simulation are then independent on the diffusion energy barrier V. This applies both in the limit of large V, when the atoms are truly fixed at the landing position, or small V, when atoms may freely move.

We also found expected asymptotic behaviour of sigma(J,V) and e(J,V) for the limit of perfect wetting, J -> -infinity, or particles diffusion allowed along only one direction (and blocked in the other direction) for large differences between control parameters in different directions.

The influence of J and V for function of height-height correlation was presented.

Acknowledgement

We thanks staff of Academic Computer Centre CYFRONET in Krakow for help with prepare results of calculations. The main calculations was carried on Convex Exemplar SPP/1600XA-32 in ACC-CYFRONET-KRAKOW.

Bibliography

  1. Kotrla,M., Review "Growth of Rough Surfaces", Czechoslovak Journal of Physics (1992) 42 449-544.
  2. Kotrla,M. and Levi,C., Review "Theory and simulation of crystal growth", J. Phys.: Condens. Matter (1997) 9 299-344.
  3. Family,F., J.Phys. A (1986) 19 L441.
  4. Wolf,D.E. and Villain,J., Europhys. Lett. (1990) 13 389.
  5. Maksymowicz,A.Z., Magdon,M.S. and Whiting,J.S.S., Computer Physics Communications (1996) 97 101-105.
  6. Maksymowicz,A.Z., Suliga,J., Kulakowski,K. and Magdon,M.S., Proceedings of the 8th Joint EPS - APS International Conference on Physics Computing (1996) 121-124.

AGH Copyright K.Malarz, 1997