Discrete simulation models of surface growth

(1)

Theoretical Physics

Discrete simulation models of surface growth

Martin Bj¨ork

martbj@kth.se Erik Deng

edeng@kth.se

SA104X Degree Project in Engineering Physics, First Level Department of Theoretical Physics

Royal Institute of Technology (KTH) Supervisor: Jack Lidmar

May 22, 2014

(2)

Abstract

In this thesis the time evolution and scaling properties of different discrete models of surface growth using computer simulation is studied. The models are mainly using random deposition at a perpendicular angle to the substrates to model the adsorption process, and both one dimensional and two dimensional surfaces are considered. The Edwards-Wilkinson, Kardar-Parisi-Zhang and Mullins equations are also studied as analytical methods to describe the growth of surfaces. The scaling exponents derived from these equations are used as reference when analysing the exponents calculated from the simulation models studied in this thesis.

We have found that the simulation models do not correspond perfectly with the analytical models for surface growth, suggesting possible flaws in our models or definitions. Despite the possible flaws, the models prove to be powerful tools for analysing the time evolution of surface growth. Furthermore, we have shown that most of the simulation models exhibit the expected scaling properties, which indicates that the surfaces do have the self-affine structure they are presumed to have.

Sammanfattning

I denna avhandling har tidsutvecklingen och skalningsegenskaperna hos olika diskreta modeller för yttillväxt studerats med hjälp av datorsimuleringar. De flesta modellerna använder slumpmässig deposition vinkelrätt mot substraten för att modellera adsorption, och b˚ade en- och tv˚adimensionella modeller har studerats. Edward-Wilkingson-, Kardar-Parisi-Zhang- och Mullinsekvationen har även studerats som analytiska modeller för att beskriva yttillväxt.

Skalningsexponenterna som har erh˚allits fr˚an dessa ekvationer har anv¨ants som referenser vid analysen av exponenterna som har r¨aknats ut fr˚an de simuleringsmodeller som har studerats i denna avhandling.

Vi har kommit fram till att simuleringsmodellerna inte st¨ammer ¨overrens perfekt med de

analytiska modellerna, vilket tyder p˚a m¨ojliga brister i v˚ara modeller eller definitioner. Trots de

m¨ojliga bristerna har modellerna visat sig vara kraftfulla verktyg vid analys av tidsutvecklingen

av yttillväxt. Vidare har vi visat att de flesta av simuleringsmodellerna uppvisar de förväntade

skalningsegenskaperna, vilket ¨ar ett tecken p˚a att ytorna har den sj¨alvaffina struktur de antas ha.

(3)

Introduction

Surface growth is often associated with the accretion of a physical surface, such as growing crystals and metals.

This is also largely the focus of the area, and there are many applications of this approach, both in academics and industries [1], including crystal growth [2], biological growth [3] and growing snow layers [4]. However, there are also other areas where the same general ideas and concepts are applicable, such as fluid flows [5], fire fronts [6] and bacterial growth [7][8].

Generally, surface growth can be considered to be the time evolution of an interface - the interface representing the growing surface. It has been observed that these interfaces share some common properties in a wide range of applications, for example the time evolution of the interface roughness and the self-affine structure of the interface [9].

These common properties makes this a very interesting subject to study because of the wide range of possible applications of the results.

In this thesis, we have chosen to study different discrete simulation models of surface growth. The main focus has been the implementation of the different models, but we have also investigated the time evolution of the roughness of each model, as well as their scaling properties. Furthermore, we have analysed a number of analytical stochastic growth equations as a comparison to the results obtained from the simulation models. The equations we have analysed are the Edward-Wilkinson equation, the Kardar-Parisi-Zhang equation and the Mullins equation - each having a different approach accounting for the processes occurring in surface growth. The scaling properties of each equation are compared to the simulation models to determine the validity of the different models.

We begin by presenting the background theory of the subject, introducing the most important concepts of surface growth as well as definitions of the properties we study in this thesis and analyses of the stochastic growth equations. We continue by describing the simulation models we have considered and our implementations of them, including a short summary of the strengths and weaknesses of each model. Subsequently, we present our results and a discussion, giving a full analysis of each model and their relation with the stochastic growth equations. Lastly, a summary of our findings is presented, concluding this thesis.

Background theory

In this section, we present the background theory of the subject at hand. This theory is presented in its entirety in [9]

and [10], but is reiterated here for convenience. We begin by introducing the concept of fractals and self-affinity, followed by a presentation of important topics, such as roughness and scaling exponents. We conclude this section by presenting different analytical methods for describing surface growth, i.e. the Edwards-Wilkinson equation (EW), the Kardar-Parisi-Zhang equation (KPZ) and the Mullins diffusion equation. The results derived from these equations will be used as references when comparing with the models we have used.

When observed at different scales, the morphology of fractals does not change, meaning that they look the same regardless of the scale at which they are being observed.

For fractals, the scaling factor is the same in all directions, but there are also objects that share a lot of properties with fractals, the sole difference being that they have different scaling factors in each direction. Objects with this property are called self-affine. Surfaces are an example of self-affine objects; when scaled correctly it is virtually impossible to tell at which magnification the surface is observed.

However, they are not self-affine in the classical sense since they do not have identical patterns repeating on different scales, but they are statistically self-affine; the statistical properties of the surface are conserved when scaling - an example of this is shown in 1. This self-affinity can be used to determine the scaling properties of the surface.

Figure 1: An example of a statistically self-affine structure; without a reference it is virtually impossible to tell which image is shown at a higher magnification. Image adapted from [11].

(5)

The surface roughness, also sometimes called the interface width, is defined as the standard deviation of the time dependent height, as

w(L, t)⌘ vu ut 1

L^d

X

i=1

[h(i, t) ¯h(t)]², (1) where w is the surface roughness, L is the system size, d is the dimension, h(i, t) is the height of the surface at lattice site i and time t, and ¯h(t) is the mean height of the surface at time t. For a general surface the roughness increases as a power of time up to a crossover time t_⇥, sometimes called the saturation time,

w(L, t)⇠ t , [t⌧ t⇥] . (2) The exponent is called the growth exponent which describes the time-dependent roughening dynamics. After the crossover time is reached the roughening saturates, giving the saturation value, wsat. The saturation value increases with increased system size L and the dependence also follows a power law,

wsat(L)⇠ L^↵, [t t_⇥] , (3) where ↵ is the roughening exponent that describes the roughening after system saturation. The crossover time also depends on a power law

t_⇥⇠ L^z, (4)

where z is called the dynamic exponent. One easy way of estimating t_⇥is shown in figure 2.

Figure 2: Two fitted lines, one for the time dependent growth, the other for the saturated value, are used to estimate the crossover time, t⇥, at their intersection point.

For different system sizes, by plotting w(L, t)/wsat(L) as a function of time, the saturation values are scaled to a

single value, and by plotting the roughness as a function of t/t_⇥, the systems will be scaled to saturate at the same time, as can be seen in figure 3.

Figure 3: Roughness curves for different system sizes, before and after scaling. The roughness is scaled using the saturation value wsat⇠ L^↵and the time is scaled using the crossover time t_⇥⇠ L^z.

This suggests that w(L, t)/wsat(L)is only dependent on t/t_⇥, thus giving the general scaling relation of the roughness, also called the Family-Vicsek scaling relation [12]

w(L, t)⇠ w^sat(L)f

✓ t t_⇥

◆

w(L, t)⇠ L^↵f

✓ t L^z

◆

. (5)

Here the function f(u), with u = t/t⇥, is a scaling function satisfying

f (u)⇠ u , [u⌧ 1]

f (u) = const , [u 1] . (6)

The exponents are also related to each other. Con- sider the fitted straight lines in figure 2; Approaching the crossover point from the left we get w(t_⇥)⇠ t⇥, while approaching from the right gives us w(t_⇥) ⇠ L^↵, according to (2) and (3). Thus t_⇥⇠ L^↵and by using (4) we get that

z = ↵

. (7)

(6)

This relation between the exponents holds for any growth process that obeys the scaling relation (5). These three scaling exponents ↵, and z characterize the growth models and tell us their self-affine structures and morpholo- gies.

Stochastic growth equations

Models based on stochastic continuum growth equations are often able to predict the scaling exponents analytically and have the general form

@h(xxx, t)

@t = (h, xxx, t) + ⌘(xxx, t) . (8) Here h(xxx, t) is the height at position xxx and time t, (h, xxx, t) is a function that reflects the modelled growth process, and ⌘(xxx, t) is the noise term corresponding to the random fluctuations in the deposition of particles. We only consider Gaussian noise, which is the most simple case for which the noise has no correlation between the individual lattice sites and thus averages to zero. The noise therefore satisfies the properties

h⌘(xxx, t)i = 0 (9)

and

h⌘(xxx, t)⌘(xxx⁰, t⁰)i = 2D ^d(xxx xxx⁰) (t t⁰) , (10) where d is the dimension and D is a constant. The form of (h, xxx, t) varies depending on the models used. We will present some of the most commonly used forms in the following sections.

Random deposition

The simplest growth model is when particles are randomly generated at a position and deposited at the top of the column underneath it. This is called the random deposition model (RD). The columns of a generated RD surface has no correlation with each other and thus each column grows independently. Because of this, we cannot expect the surface to reach any saturated state, meaning that the ↵ and z exponents are undefined and irrelevant. The exponent, however, can be determined if we consider a non-equilibrium approach to the model. Since the columns are uncorrelated the growth rate at each site must equal the average number of particles arriving at each site xxx, which implies that the function (h, xxx, t) equals a constant, C. Using this and in- tegrating (8) over time we get

h(xxx, t) = Ct + Z t

0

⌘(xxx, t⁰)dt⁰ (11) ) hh(xxx, t)i = Ct . (12)

Taking the mean of the square of (11) we get

hh²(xxx, t)i = C²t²+ 2Dt . (13) Thus we get that

w²(t) =hh²i hhi²= 2Dt . (14) Since w(t) ⇠ t¹² we have that the roughness exponent for RD is

= 1

2. (15)

As stated previously, the ↵ and z exponents are undefined for this model and can therefore not be calculated.

Edwards-Wilkinson equation

To describe models with correlated lattice sites we need to modify the function (h, xxx, t). By using symmetry arguments for a surface in equilibrium we can deduce general forms of the growth equation.

A surface in equilibrium should be invariant under the transformations

h! h + h xx

x! xxx + xxx (16)

t! t + t .

This means that the surface is independent of the origin of the coordinate system as well as the origin of time, since we should be able to study a surface from any point and any time and it should still behave consistently. For a surface not accounting for empty spaces inside the interface, it should also be symmetric about the origin of the coordinate system and the mean height which suggest invariance under the transformations

h! h

xxx! xxx . (17)

The transformations in (16) rules out any explicit dependence on h, xxx or t, leaving only derivatives of h, except for constants that can safely be ignored. Furthermore, the transformations in (17) rules out any dependence on multiples of odd order derivatives, such as rh and (rh)², leaving only terms on the form r²ⁱhand (r^2jh)(rh)^2k for any combi- nation of positive integers i, j, k. Disregarding higher order terms, the lowest term that satisfies these requirements is r²h, which gives us the Edwards-Wilkinson equation [13]

@h(xxx, t)

@t = ⌫r²h + ⌘(xxx, t), (18) where ⌫ is a diffusion dampening constant. The effect of the Laplacian term is to smooth the profile of the surface,

(7)

while keeping the mean height unchanged. Particles tend to relax from a higher to a lower position, and thus this term is often referred to as the surface relaxation term [9].

To calculate the exponents ↵, , and z we can make use of the symmetry arguments. Assuming that the roughness interface h(xxx, t) is self-affine, it should be statistically identical when rescaling horizontally and vertically as well as rescaling in time,

x x

x! xxx⁰⌘ ✏xxx

h! h⁰⌘ ✏^↵h (19)

t! t⁰⌘ ✏^zt . Using the relation of Dirac’s delta function,

d(✏xxx) = 1

✏^d

d(xxx) , (20)

and (19) in (18) we get

@(✏^↵h)

@(✏^zt) = ⌫r²(✏^↵h) + ⌘(✏xxx, ✏^zt)

✏^{↵ z}@(h)

@(t) = ✏^{↵ 2}⌫r²(h) + ⌘(✏xxx, ✏^zt) . (21) The second moment of the noise term becomes

h⌘(✏xxx, ✏^zt)⌘(✏xxx⁰, ✏^zt⁰)i = 2D ^d(✏(xxx xxx⁰)) (✏^z(t t⁰))

= 2D✏ ^{(d+z) d}(xxx xxx⁰) (t t⁰) , (22) which implies that

⌘(✏xxx, ✏^zt) = ✏ ^(d+z)/2⌘(xxx, t) . (23) The scaled equation therefore becomes

✏^{↵ z}@(h)

@(t) = ✏^{↵ 2}⌫r²(h) + ✏ ^(d+z)/2⌘(xxx, t)

@(h)

@(t) = ✏^{z 2}⌫r²(h) + ✏^{(z d)/2 ↵}⌘(xxx, t) . (24) Since we require invariance under transformation, this equation must be identical to (18). Therefore, the ✏ factors must equal 1 and the exponents are

↵ = 2 d

2 ; z = 2; =↵

z =2 d

4 . (25)

For d 2the predicted values for ↵,  0, which suggests a non-exponential behaviour. Derivations using power spec- tral density functions as well as simulations has shown that the behaviour is logarithmic for larger time spans [9][10].

Kardar-Parisi-Zhang equation

EW does not account for growth that occurs at local normals of the surface. For a growth rate v along the surface normal, the increase in height is found by using the Pythagorean the- orem

h =p

(v t)²+ (v trh)²

= v tp

1 + (rh)²

= v t

✓

1 +(rh)² 2 + . . .

◆

, (26)

where we have expanded using the assumption |rh| ⌧ 1.

This suggests that a non-linear term (rh)²should be added to account for lateral growth. The new equation becomes

@h(xxx, t)

@t = ⌫r²h +

2(rh)²+ ⌘(xxx, t) (27) and is called the Kardar-Parisi-Zhang equation [14].

Using renormalization group theory, which is outside the scope of this thesis, the exponents can be derived only in 1+1 dimensions.

↵ = 1

2; z = 3

2; =↵

z =1

3. (28)

The values for 2+1 dimensions have only been found through simulation methods, resulting in the values [9]

↵ = 0.38; z = 1.58; = ↵

z = 0.24 . (29)

Mullins equation

To account for surface diffusion, we may consider a macro- scopic current of particles on the surface, jjj(xxx, t). The diffusion process does not change the number of particles and therefore it must satisfy the continuity equation

@h

@t = r · jjj(xxx, t) . (30) The current will flow from higher to lower potentials, which suggests that the current must be related to the local chemical potential,

jjj(xxx, t)/ rµ(xxx, t) . (31) The chemical potential is related to the number of bonds formed between particles. A particle landing in valleys in the surface, corresponding to a positive surface curvature, have more neighbours, making that site a favourable place to stay at with many bonds between particles. Conversely, hills on the surface, corresponding to a negative surface curvature, has few neighbours and thus are unfavourable sites with few bonds. These conditions are satisfied when

µ(xxx, t)/ r²h(xxx, t) . (32)

(8)

Combining the results (30)- (32) we get

@h

@t = r · [ r( r²h(xxx, t))]

= r⁴h(xxx, t) . (33) This term is sometimes called the Mullins diffusion term [15]

and can be added to other stochastic growth equation to account for surface diffusion. This term models a tendency of the surface to smooth out local valleys in favour of growing vertically. To extract the scaling exponents from the Mullins equation,

@h(xxx, t)

@t = r⁴h + ⌘(xxx, t) , (34) we make use of symmetry principles again to rescale and use the same method as with EW (18) to obtain the equation

@(h)

@(t) = ✏^{z 4}r⁴(h) + ✏^{(z d 2↵)/2}⌘(xxx, t) , (35) giving the exponents

↵ = 4 d

2 ; z = 4; =↵

z =4 d

8 . (36)

Models of surface growth

In this section, we present the different models of surface growth that we have considered. Our main focus lies on simulation models in one and two dimensions; we have created the algorithms for the models with inspiration from the sources mentioned below and implemented the algorithms in MATLAB. For the Random deposition (RD), Ballistic deposition (BD) and Simple surface diffusion (SD) algorithms, the inspiration was drawn from [9], while the inspiration for the kinetic Monte Carlo algorithm (KMC) was drawn from [16].

One dimensional models are easy to implement and have theoretical value in fields like fluid flows through porous me- dia [5] and the progress of forest fires [6]. Two dimensional models are slightly more complicated to implement, but are important in the study of physical surfaces, such as in crystal growth [2], and the formation of snow layers [4]. It is con- venient to refer to the surface growth dimension as ”d + 1 ” dimensions where ”d ” denotes the substrate surface dimension and the ” + 1 ” is for the growth occurring in an extra dimension [10].

We begin by describing the simple models, RD and BD, and continue with the slightly more advanced model, SD, before moving on to the most advanced model we have implemented, KMC. In each model description, we start by presenting the general concept of the model, as well as a short explanation of our implementation in 1+1 dimension.

In the cases where we have also implemented a 2+1 dimensional version, we present a short explanation of the 2+1 dimensional implementation as well.

Random deposition

The first model we studied was a simple model where particles collide perpendicular to the surface at random positions and attach irreversibly where they are deposited, ignoring local geometry. This is achieved by repeatedly selecting positions along the surface at random and increment the surface height at that point; that is, select a position i at random, and set

height(i) = height(i) + 1.

Figure 4 shows how particles interact with the surface.

As seen in figure 5, this results in a very rough surface with narrow spikes. This model is not particularly physically plausible since all columns are uncorrelated and grow independently, but it serves as a good starting point for further development of a better model. Particles depositing at random positions each time step is a good model for a homogeneous surface adsorption, and because of its simplicity, it is easy to further develop models that account for more properties of a surface.

(9)

Figure 4: Visualisation of how particles stick to the surface in RD.

Both particle A and B fall as far as they can, without interacting with the other columns on the surface.

Figure 5: Simulation of random depostion in 1+1 dimension.

10 000 particles were deposited, and the system width was 200 lattice units.

The expansion to 2+1 dimensions of this model is quite straight-forward. As in 1+1 dimensions, a position is chosen at random along the surface, and the height at that position is incremented; that is, a position (i,j) is chosen at random, and we set

height(i,j) = height(i,j) + 1.

The result can be seen in figure 6. As can be seen, the result is a dense surface with a lot of spikes.

Ballistic deposition

As in the previous model we let the particles trajectory be perpendicular to the surface, but this time the particles attach as soon as they come in contact with another particle that is already on the surface. A random lattice site is chosen, and a particle is deposited each time step. If the surface is higher at the points next to this chosen point, the particle will stick irreversibly next to the highest of these points, otherwise, it will stick to the surface as in the previous model; that is, a position i is chosen randomly, and we set

Figure 6: Simulation of random deposition in 2+1 dimensions.

500 000 particles were deposited, and the system size was 200 ⇥ 200 lattice units.

height(i) = max(height(i-1), height(i) + 1, height(i+1)).

Figure 7: Visualisation of how particles stick to the surface in BD.

Both particle A and B stick to the first particle they encounter; particle A sticks to the top of the neighbouring column, whereas particle B sticks to the column directly underneath itself.

In figure 7, the particles interaction with the surface is shown. This yields a slightly more interesting result with treelike patterns, as can be seen in figure 8. The structure obtained from this model resembles some phenomena in nature, for example how snow layers form [4], but the lack of surface diffusion, and the fact that the resulting surfaces are very porous, are unwanted properties for our more advanced models.

As with random deposition, the 2+1 dimensional ballistic deposition model is a direct expansion of the corresponding 1+1 dimensional model. A position is chosen at random, and the deposited particle is attached at the same height as the highest of the neighbouring particles; that is, a position (i,j)is chosen at random, and we set the height at that

(10)

Figure 8: Simulation of ballistic deposition in 1+1 dimensions.

10 000 particles were deposited, and the system width was 200 units.

position, h(i,j), to

h(i,j) = max(h(i-1,j-1), h(i,j-1), h(i+1,j-1), h(i-1,j), h(i,j) + 1, h(i+1,j), h(i-1,j+1), h(i,j+1), h(i+1,j+1)).

The result can be seen in figure 9. The number of particles are the same as in figure 6, but the surface is much higher - apparently the surface is much more porous. This is the same effect that was seen in the 1+1 dimensional case, in figures 5 and 8. Note that in the 2+1 dimensional case we accounted for the nearest neighbours as well as the next- nearest neighbours.

Figure 9: Simulation of ballistic deposition in 2+1 dimensions.

500 000 particles were deposited, the system size is 200 ⇥ 200 lattice units. Note that, for convenience, only the top layer of each column is shown.

Simple surface diffusion

In reality, particles are not stuck irreversibly to the surface;

they can diffuse along the surface or even detach. A very simple model to account for this behaviour is to choose a particle from the top layer of the surface at random each time step and let it try to move one step in an arbitrary direction along the surface. The particle is allowed to move along the surface, but not up or down freely - if it reaches an ascent it must overcome an energy barrier to jump up due to the high binding energy, and if it reaches a descent its movement is restrained due to a reflective barrier existing at the edge. These ascents and descents are governed by probabilities where it is much easier to descend than to ascend, whereas moving along a flat surface is always allowed. The model does not account for the number of neighbour each particle has. A typical result from a simulation is shown in figure 10.

Figure 10: Simulation of SD in 1+1 dimensions. 20 000 particles were deposited, and the system size was 200 lattice units. One particle was deposited and 10 particles were moved each time step.

The probabilities for ascending was set to 5% and for descending 20%.

Kinetic Monte Carlo

The probabilities are in reality dependent on the binding energies from neighbouring particles, as well as surface potential barriers and the system temperature. To account for this, we use a much more advanced model, the kinetic Monte Carlo method (KMC). It takes into account the binding energies and temperature dependence, but also simulates the real time evolution of the growth process.

By using correct values for the different parameters and implementing all possible events that may occur, this model should simulate a physically realistic growth process in a realistic time frame. In our model we only consider the events adsorption, diffusion with different number of neighbours, and desorption with different number of neighbours, all of which are depicted in figure 11. Due to the nature of KMC, there is no point in discussing the number of deposited particles; the interesting quantity is instead time.

(11)

Figure 11: The possible events that are considered in the KMC algorithm. Particle A is performing adsorption, particle B is performing diffusion and particle C is performing desorption. N illustrates the neighbours of particles B and C; B has two neighbours and C has one.

All temperature dependent factors are based on the Ar- rhenius equation,

f / exp E0

kBT , (37)

where E0is the binding energy, kBis Boltzmann’s constant and T is the temperature in Kelvin. Since the exact value of the binding energy is unknown to us, and considering that the probabilities are dependent on the ratio of the temperature and the energy, the numerical value of the temperature in Kelvin is not particularly enlightening. We have therefore chosen to define the temperature as

T = T⁰E0

kB

, (38)

where T is the temperature in Kelvin and T⁰ is the temperature in E0/kB, the latter being used throughout this thesis. The binding energy was chosen to E0 = 1 eV, giving T ⇡ 11 600 T⁰. During simulations, we were mainly using T for convenience, so that value is often also given in figures for clarity.

The KMC algorithm is basically made up of four steps:

1. Calculate the rates of the events.

2. Choose a random event based on the rates and perform it.

3. Increment the time with a random number.

4. Repeat from step 1.

Initially, the rate of each possible event is calculated. The rate of an event is calculated as

revent= feventnevent, (39) where reventis the rate of the event, feventis the expected value of the frequency at which the event occurs for one

single particle, and neventis the number of particles able to carry out the event. In a general KMC algorithm, feventis recalculated each simulation step, but we have chosen them to be constants; the choice of the constants is discussed later in this section. The number of particles, however, changes every simulation step and is therefore recalculated.

The calculated rates are mapped onto a probability scale, as depicted in figure 12, where the total length of the scale is equal to the sum of all the rates, and the length of each scale segment is equal to the corresponding rate. An event is chosen using a uniformly distributed random number in the interval of the probability scale and performed.

Figure 12: The rates for arbitrary events 1 to 4 mapped on a line of length Rtot. The length of Ricorresponds to the total number of particles that may execute event i multiplied with the rate for that event.

The time advances by a time step chosen from a Poisson distribution, dt = log(u)/Rtot, where u is a uniform random number and Rtot is the total sum of all rates.

Notice that since Rtot changes at every simulation step, the distribution changes every time step. In other words, if there are many events that may occur at the given moment the time step will be small, and conversely, if there are few events that may occur at the given moment the time step will be big in comparison. In reality, many events can occur simultaneously in contrast to the sequential nature of computer simulations; the dynamic time steps takes this into account and therefore gives a realistic time evolution of the system.

(12)

The following pseudo-code illustrates our implementation of the KMC algorithm.

// Initialise all variables,

// among others the event frequencies calculate_event_frequencies()

t = 0

while t < tend

// Get the current rates of all // possible events

ri = get_current_rates()

// Calculate the cumulative rates // and the total rate

Ri = Xi k=1

rk

Rtot = RN

// Choose an event using the // cumulative rates

u1 = uniform random number 2 (0, 1]

i = i : Ri 1 < u1Rtot  Rⁱ perform_event(i)

// Increment the time with a random // number from a Poisson distribution u2 = uniform random number 2 (0, 1]

dt = log u2

Rtot

t = t + dt end

The event frequencies calculated during the initialisation are based on the Arrhenius equation, with the exception of the adsorption which occurs at constant rate for a consistent and controlled behaviour of the surface growth. The adsorption frequency, fads, is calculated as

fads = a0· (no of lattice sites) , (40) where a0is a constant that we set to 1. The desorption frequency, fdes, is calculated as

fdes=

(A0exp⇣

nEb+2E0

kBT

⌘ if n < 2d + 1

0 else , (41)

where A0is called the thermal frequency or molecular vibra- tion frequency, Ebis the bonding energy between the particle and its neighbours, E0 is a potential barrier needed to overcome to leave the surface entirely, and d is the dimension. The diffusion rate, fdif, is calculated as

rdif =

(D0exp⇣ _nE

b+E⁰₀ kBT

⌘ if n < 2d + 1

0 else , (42)

where D0 is the diffusion coefficient and E₀⁰ is an energy barrier corresponding to the lattice potential existing between two neighbouring particles [9]. The thermal frequency A0and diffusion coefficient D0were both set to the same value for simplicity, taking A0 = D0 = 10¹³Hz as a common order of magnitude of thermal frequencies [17].

The values for the energies were set to E0 = 1 eV, E₀⁰ = 2eV and Eb = 0.3eV. Even though they are not directly related to known physical values, as we could not find any satisfactory values for all the needed parameters, they are chosen within a reasonable order of magnitude [9][18].

Figure 13: Simulation of temperature dependent surface diffusion in 1+1 dimension using a kinetic Monte Carlo simulation method at T⁰= 0.0517 E0/kB, corresponding to 600 K.

The expansion of KMC to 2+1 dimensions is straight- forward: The same calculations are performed, but particles are allowed to move in 2+1 dimensions instead of one. For simplicity as well as efficiency we have disregarded effects of next-nearest neighbours.

Typical high temperature simulations can be seen in figures 13 and 14 showing relatively smooth surfaces.

(13)

Results and discussion

Results showing the roughness behaviour is presented in this section. The roughness and scaling exponents of each model are discussed in the same order as the models were presented and we conclude this section with a general discussion of our results.

The errors of the results are not calculated due to two reasons. The first and foremost reason is that the curve fits require visual estimation to determine the slopes; the point where the slopes start are not always apparent. The system- atic errors outweigh the calculated errors by far as a minus- cule shift in the estimation can result in a significant change in the calculated values. A second reason is that our definition of roughness may not be entirely correct, as we will show later, which can affect the results and produce erro- neous values.

Random deposition

For RD we predicted in the background theory that the roughening exponent = ¹₂. As seen in figure 15 the fitted curve from a single simulation has the slope = 0.48, which corresponds well with theory. By performing more simulations and averaging over them, a better value could probably have been produced. The values for ↵ and z are undefined since the roughness never saturates, as argued in the background theory. Although the RD is not a very interesting model, it does set a reference for the value of ; that is, for any model that accounts for smoothing effects, the value of should be less than 0.5, since the surface will be smoother than for RD and thus have a smaller slope.

Figure 15: A single simulation of the surface roughness of RD in 1+1 dimensions, depositing one particle at each time step, using a total of 10 000 time steps on a system size of L = 100. The slope of the power law fit is = 0.48 ⇡ 0.5, as predicted by theory, see equation (15).

Ballistic deposition

The time evolution of the roughness for BD in 1+1 dimensions is seen in figure 16. The behaviour of the system was not expected, as the results from our references follows the same curve up to each crossover time, as in figure 3, while our results are clearly not aligned [9][10].

Figure 16: The roughness over time of BD simulations using different system sizes, averaged over 100 simulations.

The roughness was defined as the standard deviation of the mean height in equation (1), which is the general definition commonly used [9][10][19]. Looking at the special case with only one particle deposited at time t = 1 we get

w(L, 1) = vu ut 1

L^d

"

1·

✓

1 1

L^d

◆2

+ (L^d 1)

✓

0 1

L^d

◆2#

= s

1 L^d



1 2

L^d + 1

L^2d +L^d 1 L^2d

= s 1

L^d

L^d 1 L^d

⇡ 1

pL^d , [L^d 1] . (43)

This shows that the initial roughness, which we assumed to grow uniformly at different system sizes, is actually dependent on the system size. The slope of the curve in a loglog- plot is still the same, but the values of the roughness are all shifted with a factor 1/p

L^d. This suggests that our definition of the roughness may not be reliable, as we do not want a size dependent time growth during t ⌧ t⇥. A better definition might be

w(L, t)⌘ vu ut^L

Xd

i=1

[h(i, t) ¯h(t)]², (44) which is the same definition as before but without the normalising system size factor.

(14)

By rescaling the data seen in figure 16 with a factor pL^d we get the expected initial form of the curves from which we fit our scaling exponents, as seen in figures 17 and 18. The final curves, scaled with the scaling exponents, do not match entirely in the crossover time region, due to the fitting of the slope being done from the initial transient growth, where the initial conditions of the surface still affect the growth exponent .

Figure 17: The roughness rescaled with a factor L¹² for similar initial growth. The dashed lines and the rings show the fitted lines and points used to find the scaling exponents.

Figure 18: The scaled curves of all system sizes. The curves do not match entirely due to the non-linear form of the transient growth and perhaps a poor fit of the curves.

The listed values in table 1 show that the values of BD 1+1 do not follow the relation between the scaling exponents (7), which suggests that the original roughness definition indeed is incorrect. The scaling exponents obtained from BD 1+1 (rescaled) do follow relation (7), which indicates that the new definition is a much better choice.

The curves in BD 2+1 were, just as in BD 1+1, not aligned during the initial growth. By rescaling the curves by the factor p

L^d the initial curves are aligned, which is

Table 1: Scaling exponents for BD

Model ↵ z

BD 1+1 0.39 0.53 2.01

BD 1+1 (rescaled) 0.89 0.50 1.77

BD 2+1 2.90 1.51 2.55

BD 2+1 (rescaled) 3.90 1.51 2.55

shown in figure 19. The estimation for the slopes were done entirely by hand due to the non-linear growth. As can be seen in figure 19, all curves start out with a common slope, after which they diverge. The common slope is in the region t ⌧ L², that is, the region where the system has less than one layer of particles, and therefore can be expected to behave like RD; the slope is indeed 0.5 as expected. In figure 20, it seems like the larger systems converge towards a common slope right before the saturation; this slope is most probably the value of we are searching for, and the calculated value matches well with the values of ↵ and z. The aberrant behaviour of the small systems is most probably an artefact of finite size effects. To support this, it would have been preferable to have simulated larger system sizes where the slope would have been more prominent, but due to the long execution times of the simulations and the limited time available, this was unfortunately impossible.

Figure 19: The roughness of BD 2+1 rescaled with a factor L for similar initial growth. The simulations were averaged over 100 times.

The values of the scaling exponents of BD 2+1 shows, just like in the 1+1 dimensional case, that the relationship between them does not follow (7). Once again BD 2+1 (rescaled) shows a correct relationship between the exponents, suggesting that the new definition of the roughness is more reliable. Note that there is no difference in and z between the rescaled and unscaled version in 2+1 dimensions, since these values were measured by hand, as opposed to

(15)

Figure 20: The roughness of BD 2+1 scaled with the scaling exponents. All curves have the same saturation value and crossover time, and the large systems can be seen to converge towards a common slope near the crossover time. The simulations were averaged over 100 times.

the values from 1+1 dimensions that were estimated using simple computer algorithms.

Simple surface diffusion

For equal system sizes, we should expect that the saturation levels only depend on the number of particles that move per time step. Figure 21 shows that this indeed is the case.

Figure 21: SD in 1+1 dimensions using system size L = 2⁶ and different number of particle moves per time step. The probability for ascending was set to 5 % and for descending was set to 20 %.

These are averages of 12 simulations for each different number of particle moves.

As with BD in 1+1 dimensions, the initial roughness for different system sizes is shifted by the same factor in SD, pL^d. The rescaled data of simulations using 100 particle moves per time step is shown in figure 22. These curves are well aligned for growth occurring before the crossover times and the exponents are easily extracted from these curves. By scaling with these exponents we get a single curve, as seen in figure 23.

Figure 22: Simulation of SD in 1+1. The number of movements at each time step was set to 100. The probability for ascending was set to 5% and for descending 20%. These are averages of 12 simulations for each system size.

Figure 23: The scaled curves of SD in 1+1 dimensions with 100 moves at each time step. The probability for ascending was set to 5% and for descending 20%. These are averages of 12 simulations for each system size.

Table 2: Scaling exponents for SD

Model ↵ z

SD 1+1 0.50 0.52 0.95

Kinetic Monte Carlo

The different events that occur in our KMC models are highly dependent on the temperature. Figure 26 shows the fraction of events in 1+1 dimensions. For low temperatures, where T⁰ ⌧ 0.045 E⁰/kB, the events are entirely dominated by adsorption and the growth process is essentially the same as in RD. For temperatures above T⁰= 0.050 E0/kBdiffusion plays a larger role and smooth out the surface. Higher temperatures should decrease the value of the growth exponent and the saturation values as a result of the smoothing effect. The difference in the surface roughness using different temperatures can be seen in figures 24 and 25.

(16)

For temperatures above T⁰ = 0.07 E0/kB, desorption processes are active. From equations (41) and (42), it is easy to predict that the frequency for desorption will be higher than the frequency for adsorption at high temperatures.

From figure 26, a simple extrapolation suggests that this point will be somewhere between T⁰ = 0.085 E0/kB and T⁰ = 0.09 E0/kB. When this happens the particles desorb faster than they adsorb, making surface growth impossible.

Figure 26: The fraction of events of KMC simulations in 1+1 dimensions at different temperatures in E0/kB, corresponding to temperatures 500 1000K. The fraction axis is shown using a logarithmic scale to better visualise the small fraction values of adsorption and desorption at high temperatures.

Figure 27 shows the rescaled roughness of a KMC simulation at different system sizes at T⁰ = 0.0603 E0/kB. When compared to corresponding graphs from the other simulation models, such as figure 22, it is apparent that this model yields completely different results; there is no obvi- ous initial power law slope and the transition to the saturation value is not as swift as in the other models. The smaller system sizes reach the saturation value relatively fast after the crossover time, whereas the roughness of the large systems keep growing slowly for a long time after the apparent crossover time, before at last reaching the saturation value.

It is also noteworthy that the crossover time seems to occur at the same time for all system sizes in both 1+1 and 2+1 dimensions, in contrast to the other models. This all suggests that the scaling relation that we have considered is too simple to apply for KMC. However, the similarity of the curves in figures 27 and 28 suggests that there is some kind of scaling relation that applies to KMC too.

Figure 27: KMC simulation of different system sizes at T⁰ = 0.0603 E0/kB, corresponding to 700 K. Rescaled with a factor pL^d. These are averages of 12 simulations for each system size.

Figure 28: KMC simulation of different system sizes at T⁰ = 0.0603 E0/kB, corresponding to 700 K. Scaled using the calculated scaling exponents. These are averages of 12 simulations for each system size.

(17)

Figures 29 and 30 show the roughness of KMC simulations in 2+1 dimensions at T⁰ = 0.0517 E0/kB. Figure 29 is unscaled, whereas figure 30 is scaled with the calculated scaling exponents. The curves are almost undistinguishable, and saturates at almost the same values. Since low temperatures leads to a behaviour similar to RD, the roughness is expected to approach a straight line with slope 0.5 at low temperatures.

Figure 29: Unscaled curves of KMC in 2+1 dimension at T⁰ = 0.0517 E0/kB, corresponding to 600 K. These are averages of 12 simulations for each system size.

Figure 30: Scaled curves of KMC in 2+1 dimension at T⁰ = 0.0517 E0/kB, corresponding to 600 K. The simulations were averaged over 12 times.

Oscillatory behaviour of the roughness occurs at high temperatures for KMC in 2+1 dimensions, as can be seen in figure 31. The oscillation begins when almost one layer is complete, when t . 1, and oscillates with decreasing amplitude until about ten layers are complete, when t ⇡ 10.

The growth of the islands on top of the current layer is slowed down since the particles tend to fill the valleys first; this serves as a retarding force that suppresses the roughening until the valleys are filled. These two opposing forces, the roughening and the smoothing of the surface, counteract each other and create an oscillating system

with constant period until equilibrium is reached [20][21].

Note that this oscillatory behaviour is not present at lower temperatures, as seen in figure 29. At lower temperatures, the diffusion is not prominent enough to compete with the adsorption, meaning that the oscillatory behaviour never arises.

Figure 31: Simulations showing the oscillatory behaviour of KMC in 2+1 dimensions at temperature T = 0.0603 E0/kB, corresponding to 700 K. These are averages of 12 simulations for each system size.

The values of all scaling exponents obtained, as well as the values calculated from the stochastic growth equations, are presented in table 3.

Table 3: Table of all scaling exponents

Model ↵ z Temp [E0/kB]

BD 1+1 0.39 0.53 2.01 -

BD 1+1 (rescaled) 0.89 0.50 1.77 -

BD 2+1 2.90 1.51 2.55 -

BD 2+1 (rescaled) 3.90 1.51 2.55 -

SD 1+1 0.50 0.52 0.95 -

KMC 1+1 0.46 0.37 1.18 0.0603

KMC 2+1 0.20 0.53 0.04 0.0603

KMC 2+1 0.27 0.38 0.80 0.0517

EW 1+1 0.5 0.25 2 -

KPZ 1+1 0.5 0.33 1.5 -

KPZ 2+1 0.38 0.24 1.58 -

Mullins 1+1 1.5 0.375 4 -

Mullins 2+1 1 0.25 4 -

(18)

The KPZ equation predicts surfaces with non-linear growth and the scaling exponents obtained from it is shown to correspond well with BD in both 1+1 and 2+1 dimensions [9][14][22]. As we can see from table 3, our values do not correspond well with KPZ theory. The BD 1+1 (rescaled) values seem to have a close relation to KPZ 1+1, but the values are all about 0.2 greater than the KPZ values. The exponents from the rescaled BD in both cases reflects the relation (7), but the values, especially in 2+1 dimensions, deviate from the values obtained from KPZ so much that there undoubtedly must have been mistakes or errors done on our part.

EW and Mullins equation account for relaxation and diffusion of particles, which should be able to describe KMC for certain temperatures, as well as SD 1+1. By comparing the values of SD 1+1 with EW 1+1 and Mullins 1+1, we see that the roughness exponent ↵ corresponds well with EW but is far too low for Mullins, even though it models diffusion. Perhaps using different probabilities in SD for ascents and descents would have yielded better results, as the probabilities were set to 5% and 20% respectively. The probabilities affect the smoothing and saturation of the surface and are therefore affecting the dynamic exponent z. The value of the growth exponent is much higher than expected, which might be due to the unlimited vertical jumps that particles are allowed to do since we did not set a limit, as well as the values of the probabilities.

The values from KMC 1+1 is obtained from a high temperature simulation, where the dominant event occurring on the surface is diffusion. The saturation levels predicted by Mullins equation should be sparsely spaced as ↵ is very high. We see that at high temperature KMC does have the same value as Mullins, suggesting that Mullins describe the growth process well before saturation occurs. The other scaling exponents do not correspond well with Mullins equation, but since the temperature is high and diffusion is very dominant, the smoothing effects makes ↵ values lower, which should also affect z.

The estimation of the exponents for KMC 2+1 was done by hand, and the results are not very reliable since the growth does not seem to follow a simple power law.

As in the 1+1 dimensional case, the ↵ and values are low, which is expected of . The value of ↵ shows that the saturation values are closely separated and is not at all predicted by Mullins equation, and the z value is also small, which means that the crossover times occurs close to each other.

The values obtained through simulations show that in a lot of cases, they do not correspond well with the values obtained by the stochastic growth equations. There are sim- ilarities that can be seen, but these might be coincidences since there are many uncertain factors involved, such as the roughness definition, the estimation for calculating the exponents, and maybe bugs in the code. Bigger systems and longer time spans as well as more simulations to average values from would have also helped to yield better and more accurate results, which due to limited time and computing power unfortunately was not possible.

(19)

Summary and conclusions

In our thesis we have studied stochastic growth equations and calculated their scaling exponents as references for our simulation models. We have implemented simple models such as RD and BD in 1+1 and 2+1 dimensions, as well as SD in 1+1 dimensions. A more advanced model using a KMC algorithm was implemented in 1+1 and 2+1 dimensions to account for adsorption, desorption and diffusion, as well as the temperature and the number of neighbouring particles.

The model using KMC used binding energies, which were chosen to be in the same order of magnitude as other references [9][10]. Other physical constants were chosen in the order of magnitudes found for many different substrates [17]. It is possible that they might simulate a physically plausible time evolution of surface growth, although small energy differences give significant changes in temperature and time span.

The exponents obtained from simulations were not always consistent with the values predicted by theory, especially simulations using BD in 1+1 and 2+1 dimensions.

The calculations for the scaling exponents were dependent on visual estimations, where slight adjustments can lead to significant changes in values. The main concern, however, was the unaligned initial slopes of the curves, leading us to question the definition of the surface roughness itself.

Our analysis showed that there indeed might be a need to correct the definition, but this has not been supported by any references found so far.

Another possibility was that there were flaws in our models or algorithms. The well documented BD was simple to implement and plots showed no sign of any unwanted effects. The roughness was calculated using a predefined standard deviation function in MATLAB, but its validity was checked as explicit calculations of standard deviations were done confirming that the definition of roughness was identical. To find where the problem was, it would help to rerun longer simulations as well as trying out different vari- ants of similar models and other definitions of roughness.

In general, we have found that surfaces exhibit statistically self-affine properties as they were scalable, enabling the roughness for different system sizes to be mapped onto a single curve. Even though the calculated exponents were different from theory, there existed values that made scaling possible, which confirms the self-affine structure of the surfaces. For higher temperatures of KMC simulations we found behaviour in the roughness that suggested that the dependence on system sizes and time does not follow a simple power law, and thus is not scalable in the same way

as the other models.

One particularly interesting phenomena was observed for KMC in 2+1 dimensions at high temperatures. The tendency for particles to fill valleys before adding to upper layers created two opposing forces that generated an oscillating system with constant period. The smoothing and roughening effects counteracted each other periodically until equilibrium was reached. This behaviour has also been observed by others in thin film growth simulations [20].

The process for growing thin film of semiconducting substrates and metals using molecular beam epitaxy (MBE) is of importance in material science and electronics industries, and is still an active research field today [23][24]. Our KMC models MBE well, although it only accounts for simple phenomena such as diffusion and desorption. It can easily be improved by taking more types of events into account, giving more detailed and physically correct simulations. Ex- amples of interesting types of events to study include diffusion along any axis, different angles of deposition, ballistic aggregation and sputtering. Accounting for these phenomena could be the subject of further work.

(20)

References

[1] B. Sapoval, J. Gouyet, M. Rosso, A. eds. Bunde, and S. Havlin, Fractals and Disordered Systems (Springer-Verlag Berlin Heidelberg, 1991).

[2] A. C. Levi and M. Kotrla, “Theory and simulation of crystal growth,”

Journal of Physics: Condensed Matter9, 299 (1997).

[3] J. Ganghoffer and J. E. Piraj´an, Thermodynamics of Surface Growth with Application to Bone Remodeling, Thermodynamics - Interaction Studies - Solids, Liquids and Gases (InTech, 2011).

[4] C. Manes, M. Guala, H. L¨owe, S. Bartlett, L. Egli, and M. Lehning,

“Statistical properties of fresh snow roughness,” Water Resources Re- search44 (2008), ISSN 1944-7973.

[5] T. H. Kwon, A. E. Hopkins, and S. E. O’Donnell, “Dynamic scaling behavior of a growing self-affine fractal interface in a paper-towel- wetting experiment,” Phys. Rev. E54, 685 (1996).

[6] J. Zhang, Y.-C. Zhang, P. Alstrøm, and M. Levinsen, “Modeling forest fire by a paper-burning experiment, a realization of the interface growth mechanism,” Physica A: Statistical Mechanics and its Appli- cations189, 383 (1992), ISSN 0378-4371.

[7] T. Vicsek, M. Cserz˝o, and V. K. Horv´ath, “Self-affine growth of bacterial colonies,” Physica A: Statistical Mechanics and its Applications 167, 315 (1990), ISSN 0378-4371.

[8] A. Epstein, A. Hochbaum, P. Kim, and J. Aizenberg, “Control of bacterial biofilm growth on surfaces by nanostructural mechanics and geometry,” IOPScience (2011).

[9] A. Barab´asi and H. Stanley, Fractal Concepts in Surface Growth (Cambridge University Press, 1995).

[10] M. Pellicione and T. Lu, Evolution of Thin Film Morphology - Mod- eling and Simulations (Springer series in material science 108, 2008).

[11] V. A. Yastrebov, “Contact between elastic rough surfaces,” (2011), URL http://www.yastrebov.fr/research_contact.

html.

[12] F. Family and T. Vicsek, “Scaling of the active zone in the eden process on percolation networks and the ballistic deposition model,”

Journal of Physics A: Mathematical and General18, L75 (1985).

[13] S. F. Edwards and D. R. Wilkinson, “The surface statistics of a granu- lar aggregate,” Proceedings of the Royal Society of London. A. Math- ematical and Physical Sciences381, 17 (1982).

[14] M. Kardar, G. Parisi, and Y.-C. Zhang, “Dynamic scaling of growing interfaces,” Phys. Rev. Lett.56, 889 (1986).

[15] W. W. Mullins, “Theory of thermal grooving,” Journal of Applied Physics28, 333 (1957).

[16] D. T. Gillespie, “Exact stochastic simulation of coupled chemical re- actions,” The Journal of Physical Chemistry81, 2340 (1977).

[17] T. Shimanouchi, “Tables of molecular vibrational frequencies consol- idated. volume i,” NIST Chemistry WebBook, NIST Standard Refer- ence Database Number 69 (1972).

[18] S. Das Sarma and P. Tamborenea, “A new universality class for kinetic growth: One-dimensional molecular-beam epitaxy,” Phys. Rev. Lett.

66, 325 (1991).

[19] H. Gould, J. Tobochnik, W. Christian, and E. Ayars, An introduction to computer simulation methods: applications to physical systems, vol. 1 (Addison-Wesley Reading, 1988).

[20] Z. J. Liu and Y. G. Shen, “Oscillating growth of surface roughness in multilayer films,” Applied Physics Letters84, 5121 (2004), ISSN 0003-6951.

[21] D. E. Savage, N. Schimke, Y. Phang, and M. G. Lagally, “Interfacial roughness correlation in multilayer films: Influence of total film and individual layer thicknesses,” Journal of Applied Physics71, 3283 (1992).

[22] P. Meakin, P. Ramanlal, L. M. Sander, and R. Ball, “Ballistic deposition on surfaces,” Physical Review A34, 5091 (1986).

[23] A. Cho and J. Arthur, “Molecular beam epitaxy,” Progress in Solid State Chemistry10, 157 (1975).

[24] J. R. Arthur, “Molecular beam epitaxy,” Surface science500, 189 (2002).

Discrete simulation models of surface growth