An improved kinetic Monte Carlo approach for epitaxial submonolayer growth
1,223, * Robert Deák , Zoltán Néda and Péter B. Barna
1Eötvös Loránd University, Department of Materials Science, Budapest, Hungary 2Babeş-Bolyai University, Department of Theoretical and Computational Physics, Cluj-Napoca, RO-400084, Romania 3Research Institute for Technical Physics and Materials Science, H-1525, Budapest, P.O box 49, Hungary
Two-component submonolayer growth on triangular lattice is qualitatively studied by kinetic Monte Carlo techniques. The hopping barrier governing surface diffusion of the atoms is estimated with an improved formula and using realistic pair interaction potentials. Realistic degrees of freedoms enhancing the surface diffusion of atoms are also introduced. The main advantages of the presented technique are the reduced number of free parameters and the clear diffusion activated mechanism for the segregation of different types of atoms. The potential of this method is exemplified by reproducing (i) vacancy and stacking fault related phase-boundary creation and dynamics; (ii) a special co-deposition and segregation process where the segregated atoms of the second component surrounds the islands formed by the first type of atoms.
Keywords: submonolayer epitaxial growth, kinetic Monte Carlo, pattern formation, two-component systems, segregation
Simultaneous deposition of different types of atoms is widely encountered in experiments and practical applications [1-8]. It can be used for engineering special coating structures (e.g. nanocomposits by inhibitor additive) [9-13], or for improving the quality of the
epitaxially grown films by using one component as a kind of surfactant [14,15]. On the other hand impurities (contaminants), operating generally as inhibitors, are always
present during the deposition process, and even a minute amount of them can drastically modify both the bulk structure and the surface growth morphology of the films [16-20]. Typical topological features related to inhibitor impurity effects are: (i) irregular shapes of monolayer islands, (ii) bunches of growth steps forming hillocks and dents within the surface of crystals, (iii) truncated and rounded crystal shapes, as well as (iv) deep grain boundary grooves decorated by small crystals in polycrystalline films. Appearance of repeated nucleation and islands on the surface indicates directly that the crystal growth is interrupted by a surface covering layer i.e. the crystals became encapsulated by the impurity phase [1,18,19,21,22]. This phenomenon was clearly demonstrated by in situ transmission electron microscopy experiments in carbon contaminated indium films .
* corresponding author: Peter. B. Barna (firstname.lastname@example.org)
The ideal one-component deposition is thus seldom realized and one always encounters the situation where species of several material components participate in the surface atomic processes. This foreign species can control the course of fundamental phenomena and the pathway of structure evolution. To understand the effect of surfactant or inhibitor impurities on the structure evolution and the complex atomic processes taking place on the growth surface a consecrated method is to use kinetic Monte Carlo (MC) simulations [23-26]. This computational approach can help researchers in predicting the pathway of structure evolution and the structures that will form at different experimental conditions. By this way one could understand also the effect of the experimentally controllable parameters and engineer structures with desired practical properties.
It is well known that MC simulations are much weaker approximations to reality than the nowadays fashionable ab-initio Molecular Dynamics simulations . By considering kinetic MC simulations several processes are taken into-account only in a phenomenological manner, without considering a microscopically realistic mechanism for it. The interaction potentials governing the dynamics of the atoms are also heuristic and usually they are not derived from first principles. MC simulations offer however a great advantage (for a review see ): it is fast and one can study thus larger systems and longer time-scales. Due to this advantage it is also more adaptable for moderate computational resources than ab-initio Molecular Dynamics. Quite reasonable number of atoms can be studied on cheep PC type computers. Making MC simulations more realistic is an important task. This could help researchers in elaborating powerful codes for predicting developing structure of thin-films or the topology of auto-organized nano-structures. Here we use a fast and microscopically more realistic kinetic Monte Carlo method for two-component sub-monolayer growths. The presented method can be generalized for several co-deposited components and also for the case of multilayer growth.
First the generally used kinetic Monte Carlo techniques are presented and the specific problem considered in the present work is discussed. Then the improved kinetic Monte Carlo approach is described and applied for the targeted problem.
2. Kinetic Monte Carlo methods for epitaxial monolayer growth
Kinetic (or resident time) Monte Carlo methods are appropriate for simulating those dynamical phenomena where several processes with widely different rates are simultaneously present. In the case of pattern formation during epitaxial growth this is the case: (i) atoms can be deposited on a crystalline surface with a given rate, (ii) atoms can diffuse on the surface (this diffusion being governed again by different rates as a function of the binding energy of the specific atom) and (iii) decohesion of adatoms from the surface are also possible (Fig. 1). When dealing with several material components the different type of neighbouring atoms can exchange by interchanging their positions (Fig. 1) following complicated microscopic mechanisms. This exchange-segregation process is characterized again with a different rate. Sometimes this is the key process for controlling
the structure evolution and, consequently, for engineering specific (e..g. nanocomposite) structures in multi-component systems.
Fig. 1 Processes with widely different rates that govern the dynamics of atoms during epitaxial growth.
In order to save precious computational time in simulating these processes with widely different rates, a Monte Carlo method used for studying equilibrium properties of low temperature systems (the BKL Monte Carlo method ) was adapted and named as kinetic Monte Carlo method (for a review see ). The basic idea is that in each simulation step one process is probabilistically selected (the probability to select one process is proportional with it’s rate) and carried out. The time is than updated non-
uniformly, depending on the rates of all possible processes at that given moment.
Generally the deposition rate is fixed and calculated from the deposition speed
(deposition flux) given as the number of new monolayers deposited in unit time (ML/s).
The diffusion rate (r) of an atom is governed by the thermodynamic temperature (T) X?Y
of the system and the potential barrier (ΔE) that the atom has to overcome between X?Y
the initial (X) and the final (Y) position:
E？X！Y(1) rf exp()？；0X！YkT
In expression (1) k is the Boltzmann factor and f is the attempt rate, which is roughly the 013vibration frequency of atoms in the crystal (f?10 Hz). Since the value of the barrier is 0
not straightforward to estimate (even if the pair-potential between the atoms is known), several simplifying methods are used. The simplest approach is to consider the potential barrier dependent only on the binding energy of the atom in the initial X state [29-32] or by applying the transition state theory . A better, but computationally more costly approach is to consider a realistic pair-potential between the atoms  and estimate the potential in several points between the initial and final state. The embedded-atom method [35,36] offers another possibility for estimating the potential barrier in the hopping process. In such case the difference between the maximum and initial value will yield the
potential barrier. An even more complex approach would be to map the potential in the neighbourhood of the initial and final point not on a line, but on a surface.
The decohesion rate is obtained either by fixing a potential barrier E for this process dec
or by calculating the more realistic potential barrier as the binding energy of the chosen atom at the given site.
Exchange between neighbouring and different type of atoms are microscopically realized through complex vacancy mechanisms. In simulations however an oversimplified geometry is considered where many possible degrees of freedom for diffusion are not considered, so exchange possibilities are mostly blocked. The exchange rate is then
usually postulated in form (1) by assigning a hypothetical E potential barrier for this ex
Simulations are usually performed in a two-dimensional geometry [29-32], the atoms being allowed to occupy only the sites of a pre-defined lattice. By this approach one reproduces an idealized situation where a new layer is growing on a perfect crystalline substrate. The simplest possibility is to consider a square lattice and the sites on the growing layer positioned exactly on the top of the atoms forming the substrate [29-32]. In such manner a non-realistic three-dimensional cubic structure is simulated but approaches on more complex geometries are also possible. One can use lattices with different symmetries and different stacking sequences for positioning the atoms in the growing layer . Simulations can be made more realistic by considering a second layer on the
top of the simulated one so that interchanges between these two layers become possible. This would allow formation of additional defects and vacancies.
Nowadays computationally costly off-lattice kinetic Monte Carlo methods [37,38] are also considered for the case when several types of atoms are simultaneously present and there is lattice constant or symmetry mismatch between the crystalline structures of the components. In such an approach the position of the atoms are computed from an energy minimization procedure and the dynamics of the system is realized with the kinetic Monte Carlo algorithm. The method is an optimal reconciliation between the realistic nature of the Molecular Dynamics simulations and the higher speed of the kinetic Monte Carlo approach. Beside many interesting results obtained with this approach the method is still not usable for reasonably large system sizes in (2+1)D and practically relevant dynamics times.
3. The considered problem and previous results
A problem which was several times studied in the literature by kinetic Monte Carlo
simulations is reconsidered and analyzed in the present work by an improved approach. Atoms of type A (growing materials) and B (impurities) are co-deposited on the planar
surface of a perfect 3D single crystal of A type of atoms [31,32]. There is a deposition
rate for both components (F and F respectively, with F=F for simplicity), and the ABAB
atoms deposited on the substrate are allowed to diffuse there. The conditions necessary to obtain a particular segregation of the A and B atoms where the B impurities will surround
(decorate) the islands formed by the A adatoms are investigated. Earlier kinetic Monte
Carlo simulations [31,32] performed on simple square-lattice geometry and using a first approximation for the ΔE hopping barrier concluded the important result that such X?Y
structures can be obtained only if a direct exchange mechanism between A and B atoms
on neighbouring sites is postulated.
More precisely in the earlier simulations it is considered that the hopping barrier for an
Watom of type W ()depends only on the initial state of the atom. It is the sum of a ？E ！XY
term related to the substrate, Eand a contribution related to each lateral nearest sub ,
neighbour. Contributions depend on the local composition so that for each term we have four possibilities: AA, AB, BA and BB. The hopping barrier is then calculated as
WWWWQWQWQ, (2) ？E？E？(nE;nE)！(XYh0sub1n，QA,B
WWQwhere is 1 if the substrate atom is of type W and 0 otherwise, is the number of nn01
WQnearest-neighbour W-Q pairs, is the corresponding contribution to the barrier En
WQQWWQ(symmetric in W and Q, ) and is the contribution from a free W atom E？EEsubnn
on a substrate atom Q . In order to perform the kinetic Monte Carlo simulation for this
WQWQsystem one has to postulate the and values and the f attempt frequency. The EE0 subn
T temperature of the system and the F=F deposition rate has to be selected. AB
Simulations were performed on a square lattice where the positions of the atoms in the growing layer are on the top of the substrate atoms.
These simulations proved that an energetic bias favouring segregation is not sufficient to obtain configurations with impurities (B) mostly positioned at island edges. To achieve
this peculiar segregation a thermally activated exchange mechanism had to be introduced between A and B atoms on neighbouring sites. The easiest way to realize this was postulating a phenomenological potential barrier E, for this process and to use equation ex
(1) for calculating the exchange rate.
The present work intends to further improve the kinetic Monte Carlo simulation for the two-component system. Our task is to reduce considerably the number of postulated
WQWQparameters (e.g. and ) and to consider a more isotropic geometry which will EEsubn
increase the degrees of freedom for diffusion. As a consequence of this, for a reasonably large parameter set the interchange of A and B atoms will arise automatically and the
phenomenological parameter E is eliminated. In several cases when the energetic bias ex
favours segregation impurity decorated islands will form.
4. The improved kinetic Monte Carlo approach
In the present kinetic MC algorithm less parameters is postulated, diffusion of atoms with increased degrees of freedom and improved potential barriers are considered. At the same time the speed of the algorithm does not drop in considerable manner (the advantage of the MC approach is kept).
The first modification is that the triangular lattice ((111) plane of fcc structure) is used as substrate (filled circles in fig. 2). This leads also to a more compact packing of the atoms. It is assumed that atoms are spheres with the same diameter for both A and B components.
In such manner there are two triangular sub-lattices (empty circles and crosses in fig. 2) on which the adatoms can be deposited forming monolayer lattices of fcc or hcp crystalline phases. Due to geometric restrictions atoms in the growing layer cannot occupy neighbouring sites belonging to different sub-lattices.
Fig. 2 Geometry of the considered lattice. Filled in large circles represent the atoms of the substrate, small empty circles and crosses represents the fcc and hcp lattice sites, respectively, on which the new layer can growth.
Considering a bulk fcc substrate structure stacking fault develops at the interface of the substrate and a growing hcp monolayer island (Fig. 3). By this way phase boundaries will also appear between growing islands of fcc and hcp types. This extra defect mechanism
characteristic for this geometry facilitates the diffusion and interchange of atoms between the islands of the two growing phases. For the case of the two-component system only the growth of the fcc phase will be investigated, however the formation of phase boundaries and their motions in simple homoepitaxy will be also examined.
Diffusion of adatoms on the top of the first fcc type growing monolayer is also considered. These adatoms can also jump down on the substrate.
Apart of geometry a second improvement in the kinetic MC algorithm is in the calculation mode of the hopping barrier for the diffusion. First, realistic pair-potentials between the atoms are used to compute the binding energies of the atoms. The hopping barrier for the diffusion process is then calculated from the binding energies in the initial and final states.
A Lenard–Jones type of pair-potential was considered, although the program permits easily to change it to any other type of accepted form.
Fig. 3.The fcc and hcp phases that can form on a perfect fcc bulk substrate
Assuming the lattice constant as length unit, the interaction potential between atoms of type W and Q separated with distance r can be written with one E parameter: WQ
12UE (3) ？(；)WQWQ126rr
The parameter E will fix the binding energy at the r=1 equilibrium distance. WQ
Interaction between atoms are taken into account up to s=3 lattice site distances. It is
assumed that the hopping barrier from a site X to a site Y should depend not only on the
binding energy in site X, but it should also depend on the change in the binding energy. The following form for calculating the hopping barrier of an atom has been proposed :
XYX (4) ？E？；，E;(1；，)(E；E)！XYnnn
XYIn equation (4) ( ) is the total interaction energy (binding energy) of the atom at EEnn
sites X (Y), respectively. α is a parameter between 0 and 1, whose value will be
determined later. This is the simplest linear form in which the barrier depends both on the binding energy in the initial position and on the difference between the binding energies
Xof the final and initial sites. Moreover, it yields the good barrier for decohesion (-) En
and reasonable values for the self-surface diffusion and edge-diffusion. In order to get positive barriers for each possible process α has to be bounded between 0.3 and 0.6. A
simple exercise using Lenard-Jones type potentials on an fcc structure shows that the ratio of the energy barrier for self-surface diffusion and adsorption energy should be around 0.35. It is immediate to realize that this ratio is exactly the value of α, and a first
estimate is. ，？0.35
The number of postulated parameters in the simulation code is strongly reduced. One has to fix only the values, the F=F deposition rate, the fattempt frequency E,E,EAB0 AAABBB
and the T thermodynamic temperature of the system. Parameters very similar to the one
13used in previous studies [31, 32] were considered. We have fixed f=10 Hz, E 0AA
=0.1eV, E =0.01 eV, so that aggregation of A particles are favoured relative to the BB
aggregation of B particles. E =E is varied in the 0.01 – 0.07eV interval, the ABBA
thermodynamic temperature T was considered in the realistic 300 – 500K interval and the
deposition rate F=F in the 0.0001 – 0.1 ML/s domain. Systems with lattice sizes up to AB
500x500 were easily simulated in a few days on normal PC type computers (Pentium 4, 3.4 GHz). As will be discussed in the next section, the peculiar segregation with the B
impurity atoms decorating the islands formed by the A atoms can be obtained for a quite
reasonable parameter range.
5. Simulation results
In the present paper we will present only some qualitative results for illustrating the applicability of the used method. Specific and practically important problems can be than studied using the described algorithm.
5.1 Evolution and annihilation of stacking-faults and phase-boundaries on an fcc (1,1,1) surface.
The case when A type of atoms are deposited with a low deposition rate (F=0.001ML/s) A
on a planar fcc (1,1,1) surface are simulated. In this case monolayer domains of the two equivalent orientations are nucleated and grown. Formation, motion and annihilation of stacking faults related phase-boundaries appear and can be nicely followed during simulation (fig. 4). Some movies are also given on the home-page dedicated to this study .
Fig.4 Characteristic time evolution and annihilation of stacking faults related phase-boundaries for the case when only A type of atoms are deposited. The pictures from left to right represent some steps in the time-evolution. The F and H f islands correspond to fcc and hcp stacking, respectively. Simulation parameters are: E =0.1eV, T=400K, AA13F=0.001ML/s and f=10 Hz. A0
5.2 Formation of impurity decorated islands in case of two-component deposition.
Simulations with the co-deposition of the A (•) and B (+) type of atoms leads to the 13expected structures. For the fixed parameters (f=10 Hz, E =0.1eV, E =0.01 eV) 0AABB
and as a function of the E parameter two main type of structures are observable: (i) AB
island containing intermixed A and B type of atoms, and islands decorated by B
impurities. As expected, for low E values the impurity decorated islands are stable, AB
while for higher E values the islands containing intermixed A and B type of atoms are AB
observable (fig. 5). Increasing or decreasing the temperature will only shift the boundary between these two type of structures and favour larger or smaller islands of lower or higher number density (nucleation density) respectively, for the same number of deposited atoms.
Fig. 5 Island structures obtained as a function of the E parameter. Simulations with AB13E =0.1eV, E =0.01 eV, T=400K, F=F =0.001ML/s and f=10 Hz. The present AABBAB0
structures are obtained for t=5000 simulation time, and a central part of a much larger
simulation area is presented. For E<0.07 eV impurity decorated islands are formed. AB
The time evolution of the structures for the case of impurity decorated islands is illustrated in Fig. 6. Some movies showing a more complete dynamics are available on the home-page  dedicated to this study.
Fig.6 Snapshots from the time evolution of the system. Simulation parameters are E AA13=0.1eV, E =0.01 eV, E=0.05 eV, T=400K, F=F =0.001ML/s and f=10 Hz. A BBABAB0
central part of a much larger simulation area is presented.
6. Discussion and Conclusions
The aim of the present work was to show that kinetic Monte Carlo methods can easily be improved without loosing it’s main advantage of simulating large systems and reasonable
long dynamics time. The two-component co-deposition process was successfully
simulated by considering two main improvements relative to the simple classical method. First a new hopping barrier formula was used, calculated from the pair interaction potentials of the atoms (assumed to be of Lenard-Jones type in this study). Secondly, more degrees of freedom for the diffusion of particles were allowed by considering a triangular lattice and the growth of a second layer on the top of the simulated one. Interchanges of adatoms between these two layers are also possible. The exchange mechanism between atoms of different types on neighbouring lattice sites appears directly from diffusion without considering an artificially predefined rate for such a process. By this way it is possible to simulate with a reduced number of parameters the
special segregation process in which impurities are decorating the growing islands. Vacancies, stacking faults, phase-boundaries and their dynamics are also successfully reproduced by this improved kinetic Monte Carlo technique. The method can be generalized to computationally study practically important pattern and structure formation problems during the co-deposition of several types of atoms on a given crystalline substrate. The presented method is computationally fast, and systems with tens of thousands of atoms can be simulated in reasonable computational time on PC type computers. The method is limited however to the case when the deposited components have similar crystal structures and similar lattice constants. Whenever there is lattice constant mismatch between the components off-lattice kinetic Monte Carlo methods have to be used [37, 38]. This will reduce however considerably the size of the system and the time length of the dynamics that can be investigated with a given computational resource.
This work was supported by the Hungarian National Science Foundation under contract thNo. OTKA T048699, by INNOVATIAL 6 Framework Program Project IP 515844-21
and by Romanian CNCSIS grant 41/183. The authors acknowledge illuminating discussion with Prof. G. Radnoczi (Research Institute for Technical Physics and Materials Science, Budapest, Hungary).
1. J.F.Pócza, A. Barna, P.B. Barna, I. Pozsgai, G. Radnoczi, Jpn. J. Appl. Phys., Suppl. 2, Part 1, 525, 1974
2 H. A. van der Vegt, J. M. C. Thornton, H. M. van Pinxteren, M. Lohmeier, and E. Vlieg, Phys. Rev. Lett. 68, 3335 ~1992
3. G. Rosenfeld, R. Servaty, C. Teichert, B. Poelsema, G. Comsa, Phys. Rev. Lett. 71,
4. V. Fiorentini, S. Oppo, and M. Scheffler, Appl. Phys. A: Mater. Sci. Process. 60, 399
5. J. A. Meyer, H. A. van der Vegt, J. Vrijmoeth, E. Vlieg, and R. J. Behm, Surf. Sci., 355, L375 ~1996!
6. B. Voigtlander, A. Zinner, T. Weber, H. P. Bonzel, Phys. Rev. B 51, 7583 ~1995!.