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.
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.
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:
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.
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.
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.
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.
< h > | 1.000 | 4.000 | 9.000 | 16.000 | 25.000 | 64.000 |
sqrt(< h >) | 1.000 | 2.000 | 3.000 | 4.000 | 5.000 | 8.000 |
sigma | 0.987 | 1.994 | 2.998 | 3.994 | 5.004 | 7.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.
V | -50 | -10 | -2 | -1 | 0 | 1 | 2 | 10 | 50 |
sigma | 4.0 | 4.0 | 4.0 | 4.0 | 4.0 | 4.0 | 4.0 | 4.0 | 4.0 |
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.
J | -10.0 | -3.0 | -1.5 | -1.0 | -0.5 | -0.25 | 0.0 | 0.25 | 0.5 | 0.75 | 1.0 | 1.5 | 3.0 | 10.0 |
sigma | 1.0 | 1.0 | 1.1 | 1.2 | 1.7 | 2.3 | 4.0 | 7.2 | 10.9 | 14.1 | 16.3 | 18.8 | 21.2 | 21.9 |
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.
V | -50 | -10 | -2 | -1 | 0 | 1 | 2 | 3 | 5 | 10 | 50 |
e | -1.262 | -1.263 | -1.269 | -1.269 | -1.251 | -1.158 | -0.900 | -0.508 | -0.093 | 0.003 | 0.002 |
J | 0 | 0.1 | 0.25 | 1 | 1.5 | 2 | 3 | 5 |
e | -0.002 | -0.283 | -0.718 | -1.251 | -1.552 | -1.595 | -1.597 | -1.597 |
V | 0 | 0.5 | 1 | 1.5 | 2 | 2.5 | 3 | 5 | 10 | 50 |
e | 0.001 | 0.011 | 0.010 | -0.010 | -0.046 | -0.091 | -0.146 | -0.336 | -0.479 | -0.475 |
V | 0 | 0.5 | 1 | 1.5 | 2 | 2.5 | 3 | 5 | 10 | 50 |
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.
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).
range | 1 | 2 | 3 | 4 | 5 |
F | 0.715 | 0.475 | 0.275 | 0.132 | 0.041 |
range | 1 | 2 | 3 | 4 | 5 |
F | -0.472 | 0.061 | 0.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.
range | 1 | 2 | 3 | 4 | 5 |
F(u) | -0.778 | 0.428 | -0.186 | 0.062 | -0.010 |
range | 1 | 2 | 3 | 4 | 5 |
F(v) | 0.814 | 0.587 | 0.380 | 0.222 | 0.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.
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.