- Email: [email protected]
Computer simulation of nucleation on (001) Kossel crystal surface with and without diffusion
surface s c i e n c e EI.SEVIER
Applied Surface Science 92 (1996) 257-260
Computer simulation of nucleation on (001) Kossel crystal surface with and without diffusion Tong B. Tang * Physics Department, H.K. Baptist University, Kowloon Tong, Kowloon, Hong Kong
Received 14 December 1994; acceptedfor publication2 March 1995
Abstract
The Monte Carlo technique was applied to nucleation-controlled growth of a (001) Kossel crystal surface, either with or without surface diffusion. The sizes of critical nucleus, no, and of supercritical cluster, n+, have been obtained at different temperatures and supersaturations. With the inclusion of diffusion, both sizes are reduced, and this explains why the nucleation dip of nucleation growth is shortened. However, n+/n o -- 2, as predicted by the theory of Burton, Carbrera and Frank, only in the absence of diffusion. Also, when diffusion operates the kink site density along steps is no longer independent of supersaturation. It is concluded that surface diffusion brings the system further away from thermodynamic equilibrium.
1. Introduction Crystal growth involves the transportation of masses and energy from one phase to another and cannot be adequately treated by equilibrium thermodynamics. However, kinetic theories suffer from a problem with spatial inhomogeneity, namely the instantaneous distribution of clusters with different sizes. It is natural to resort therefore to computer simulation, for which the Monte Carlo method proved most effective, since we may with sufficient generality model the elementary steps of absorption, evaporation and surface migration of the growth units by a sequence of stochastic events. Indeed, since the pioneering works by Chernov and Lewis [1] and others,
' Tel.: + 852 2 339 7036; fax: + 852 2 304 6558; e-mail: [email protected].
Monte Carlo simulations of crystal growth have reached increasing levels of sophistication in dealing with greater degrees of freedom and more complicated boundary conditions [2,3]. Nevertheless, early theoretical considerations based on reversible thermodynamics have succeeded in determining a number of equilibrium properties of the crystal surface, that serve as parameters in its growth process, but a few basic questions remain concerning their validity under various conditions. The growth of a low-index, defect-free surface at low temperature proceeds from layer to layer and is preceded by the continuous formation and disappearance of two-dimensional nuclei (the Zeldovich fluctuation regime). The critical nucleus size n o is defined as the size at which the probabilities to grow and to dissolve are equal. Supercritical clusters, however, are those that sustain steady average growth.
0169-4332/96/$15.00 © 1996 Elsevier Science B.V. All fights reserved SSDI 0169-433 2(95)00238-3
258
T.B. Tang~Applied Surface Science 92 (1996) 257-260
Burton, Cabrera and Frank (BCF) have [4] made predictions about their shape, and size: n+ = 2n 0. In this paper, this and other deductions will be examined in the light of results from a Monte Carlo study on crystal growth with and without surface diffusion. Conversely, some observations will also be offered concerning the general implication of diffusion, first introduced into simulation studies by Gilmer and Benema [5].
2. Simulation parameters Temperature T and supersaturation A p, were expressed by the dimensionless relations:
to = 4,ss/(2kr), with ~bs~ representing the bond energy between two nearest neighbors on the surface, and
/3 = Alz/kT= ln(Cf/Cf0 ), where Cf denotes the concentration of the fluid and Cf0 its equilibrium concentration. Simulation runs were carried out at to--1.0 and 2.0, below the roughening temperature toR = 0.77. Nucleation dips are reduced by surface diffusion, to /3--0.05 for to = 1.0 and /3 = 1.0 for to--2.0 [8]. On the other hand, to avoid kinetic roughening we need fl < to so that A/z ~: ~bss. Thus fl values have been chosen between 0.05 and 0.25 for to = 1.0 and in the range 1.0-1.8 for to-- 2.0. We used the kinetic Ising model with a stress-free (100) simple cubic Kossel crystal face under the solid-on-solid condition [6,7]. This is the simplest possible model which stillcontains the essential configurational degrees of freedom for building up new crystal layers. In runs with diffusion, the jump length X~ was set to one lattice unit for consistency with BCF theory [8] and had 8 possible directions, namely, to any of the 4 lateral or 4 diagonal neighbors; going up or down a step was also permitted.
3. Simulation procedures It was assumed that the rate of deposition corresponds to the impingement frequency of the fluid
particles on the surface and is therefore proportional to
p + = Ce/Cfo = exp/3. For evaporation and diffusion, the same activation energy was assumed but they depend on local configuration at the site, such that p( N)d = p ( N ) - = exp(4 - 2 N ) to, where N denotes the number of lateral nearest neighbors [5]. The procedure due to Liu, van der Eerden and Bennema [9], that simulates growth at low temperature with high efficiency, is adopted. The type and the site where an event happened were selected according to the waiting times of different events at all sites, calculated as t += - I n n+/p +,
t-= -In n-/p(N)-,
t d = - I n nd/p(N) d, where n +, n- and n d lie within the interval (0,1] and are drawn from a single pseudo-random number chain (Ref. [10] pp. 16,17). The procedure was implemented in MS FORTRAN 5.0 to double-precision format on a microcomputer. The event with the minimum waiting time was accepted and, after the appropriate changes of site configurations, the program calculated new waiting times for those sites and sorted them again with the others, to maintain the event with the minimum time at the top of the to-do list. After a preset number of Monte Carlo cycles, an output record containing the surface configuration, the numbers of deposition, evaporation and diffusion events, the numbers of kink and step sites, as well as cumulative waiting time would be written to a file. The cycles of waiting time generation and to-do list were updated and then re-commenced. A 30 × 30 surface imposed with the periodic boundary condition (Ref. [10] pp. 22-24) was used. The mean diffusion steps at /3 = 0 have previously been determined via extrapolation [I1] as do--8.1 for to = 1.0 and d0 = 8.7 for to = 2.0. Within the ranges of /3, therefore, the surface used suffices in size to prevent any severe coupling of a lone nucleus with its image across a periodic boundary. However, if more clusters had formed, the dynamic behavior of the system would be distorted [12]. Hence the pro-
T.B. Tang/Applied Surface Science 92 (1996) 257-260
gram always halted at 50% coverage, below which no second cluster would appear [11]. In total 180 simulation runs were executed. Ninety proceeded with surface diffusion (X~ = 1), ninety without (X~ = 0). At least three runs were conducted for either case at each value of /3, and their outputs discarded unless the data deduced from them agree with one another to within 10%.
0.6
•=
259
(a)
0.5
;
0.4
~
:
1.2
0.8
:
:
1.6
beta
,~- 0.5 i(b)
4. R e s u l t s a n d d i s c u s s i o n
Qualitatively, I examined the time development of coverage in each run by viewing its output records of surface configuration on the computer monitor as consecutive pages of a Word document, in quick successions by pressing the 'Page Down' key. The coverage was then plotted against cumulative waiting time for representative runs. The size n+ of the growth cluster was taken from the smallest cluster at the onset of rapid growth, and n o of the critical nucleus from the largest that had dissolved before
8O
(a)
600
.~ 4 0 " 20
:
0
:
:
:
:
1ib80 0.8
1.2
1.6
beta
.N
60
ca 40
20
0 0
0.2
0.1
0.3
beta Fig. 1. Sizes of growth cluster (<>. ×) and critical nucleus ( n , zx) for (a) to = 2.0 and (b) to = 1.0 (<>, [] : X~ = 0 and ×, ,x : X s = 1).
0.2
• 0
:
t
0,1
I 0.2
t 0.3
beta
Fig. 2. Time average of kinks per step site for (a) to = 2.0 and (b) to= 1.0([]: Xs=O and n: Xs= 1).
this onset. Data for to = 1.0 and 2.0 are summarized in Figs. l a and lb, which show that, in the presence of diffusion, both n+ and n o were reduced, the latter to a greater degree. This result matches the previous observation of Van Leeuwen and Van der Eerden [12]. The Wulff theorem and the G i b b s - T h o m s o n formula [4] state respectively that the product n o/32 should be a constant and that n+ = 2n o. A plot of n o versus /3-2 for to = 2.0 is indeed fairly linear and indicates t h a t n 0 / 3 2 = 26 ( X s -- 0) and 10 ( X s = 1). For t o = 1.0, the slopes are n 0 / 3 2 = 0 . 4 2 ( X s = 0 ) and 0.25 ( X s = 1). However, n + / n o = n + / 3 2 / n o / 3 2 = 3.6 and 29 for X s - 1, to = 1.0 and 2.0 respectively. Only the cases of X~ -- 0 satisfy the condition that n + / n o = 2.0. BCF [4] argue that an infinitely long step acts as a continuous-line sink for single particles and its density of kinks during growth, maintained by thermal fluctuation, should be independent of supersaturation. However, our results (Fig. 2) support this deduction only for X~ = 0. In the diffusional cases, the kink densities are lower, especially so at low supersaturation. This result is in qualitative agreement with the prediction from an analytical model [13].
260
T.B. Tang~Applied Surface Science 92 (1996) 257-260
5. Conclusion The size of growth cluster is reduced by the introduction of surface diffusion, and that for critical nucleus, to a more marked extent. This explains the shortening of the nucleation dip found earlier [8]. However, by the same token it is no longer true that one is always half of the other, as derived by BCF, who also predict an invariant density of kinks along surface steps. Again, the presence of diffusion invalidates this prediction, the more so at a higher temperature. At high temperature and supersaturation, kinetic roughening occurs as indicated by a change in slope of the jump frequency [11], and equilibrium structures such as critical nuclei lose meaning. The conclusion is therefore that, contrary to previous views [12], diffusion, by opening up more energy decay paths in the configurational space, brings the surface further away from thermodynamic equilibrium. Lastly, the work reported here may be the first that has been undertaken in its entirety on a personal computer. A 486DX-33 was used, and a single job took overnight to a full day. Since the current price for this or some faster model has fallen under US$1200 (in Hong Kong as of mid 1995), more than one can crunch the numbers independently to cut down the aggregate computational time. On a more advanced level, with the help of communication utilities such as the PVM (Parallel Virtual Machine) v3.5 from Oak Ridge Nat. Lab., we can run the program on a network of PowerPCs or PowerMacs that support Unix. Parallelization methods are of
course applicable to a number of Monte Carlo and molecular dynamic problems [14]. The trend towards 'parallelomicrocomputing' may perhaps extend to other areas of computational physics, that likewise place light demand on math software library.
References [1] A.A. Chernov and J. Lewis, J. Phys. Chem. Solids 28 (1967) 2185. [2] H. Miiiler-Krumhhaar, in: Monte Carlo Methods in Statistical Physics, 2nd ed., Ed. K. Binder (Springer, Berlin, 1986) pp. 282, 290. [3] U. Landman, Dynamics of lnterfaciai Processes: Simulations of Growth and Stability of Films, Junctions and Nanostructurcs, presented at ICSFS-7. [4] W.K. Burton, N. Cabrera and F.C. Frank, Phil. Trans. R. Soc. London A 243 (1951) 299. [5] G.H. Gilmer and P. Bennema, J. Appl. Phys. 43 (1972) 1347. [6] C. van Leeuwen and P. Bennema, Surf. Sci. 51 (1975) 109. [7] J.D. Weeks and G.H. Gilmer, Advan. Chem. Phys. 40 (1979) 157. [8] V.K.W. Cheng, E.C.M. Tang and T.B. Tang, J. Cryst. Growth 96 (1989) 293. [9] G.Z. Liu, J.P. van der Eerden and P. Bennema, J. Cryst. Growth 58 (1982) 152. [10] K. Binder, in: Monte Carlo Methods in Statistical Physics, 2nd ed., Ed. K. Binder (Springer, Berlin, 1986) pp. 1-45. [11] E.C.M. Tang and T.B. Tang, Ferroelectrics 142 (1993) 9. [12] C. van Leeuwen and J.P. van der Eerden, Surf. Sci. 64 (1977) 237. [13] J. Zhang and G.H. Nancollas, J. Cryst. Growth 106 (1990) 181. [14] D.W. Heermann and A.N. Burkitt, Parallel Algorithms in Computational Science (Springer, Berlin, 1991).
Applied Surface Science 92 (1996) 257-260
Computer simulation of nucleation on (001) Kossel crystal surface with and without diffusion Tong B. Tang * Physics Department, H.K. Baptist University, Kowloon Tong, Kowloon, Hong Kong
Received 14 December 1994; acceptedfor publication2 March 1995
Abstract
The Monte Carlo technique was applied to nucleation-controlled growth of a (001) Kossel crystal surface, either with or without surface diffusion. The sizes of critical nucleus, no, and of supercritical cluster, n+, have been obtained at different temperatures and supersaturations. With the inclusion of diffusion, both sizes are reduced, and this explains why the nucleation dip of nucleation growth is shortened. However, n+/n o -- 2, as predicted by the theory of Burton, Carbrera and Frank, only in the absence of diffusion. Also, when diffusion operates the kink site density along steps is no longer independent of supersaturation. It is concluded that surface diffusion brings the system further away from thermodynamic equilibrium.
1. Introduction Crystal growth involves the transportation of masses and energy from one phase to another and cannot be adequately treated by equilibrium thermodynamics. However, kinetic theories suffer from a problem with spatial inhomogeneity, namely the instantaneous distribution of clusters with different sizes. It is natural to resort therefore to computer simulation, for which the Monte Carlo method proved most effective, since we may with sufficient generality model the elementary steps of absorption, evaporation and surface migration of the growth units by a sequence of stochastic events. Indeed, since the pioneering works by Chernov and Lewis [1] and others,
' Tel.: + 852 2 339 7036; fax: + 852 2 304 6558; e-mail: [email protected].
Monte Carlo simulations of crystal growth have reached increasing levels of sophistication in dealing with greater degrees of freedom and more complicated boundary conditions [2,3]. Nevertheless, early theoretical considerations based on reversible thermodynamics have succeeded in determining a number of equilibrium properties of the crystal surface, that serve as parameters in its growth process, but a few basic questions remain concerning their validity under various conditions. The growth of a low-index, defect-free surface at low temperature proceeds from layer to layer and is preceded by the continuous formation and disappearance of two-dimensional nuclei (the Zeldovich fluctuation regime). The critical nucleus size n o is defined as the size at which the probabilities to grow and to dissolve are equal. Supercritical clusters, however, are those that sustain steady average growth.
0169-4332/96/$15.00 © 1996 Elsevier Science B.V. All fights reserved SSDI 0169-433 2(95)00238-3
258
T.B. Tang~Applied Surface Science 92 (1996) 257-260
Burton, Cabrera and Frank (BCF) have [4] made predictions about their shape, and size: n+ = 2n 0. In this paper, this and other deductions will be examined in the light of results from a Monte Carlo study on crystal growth with and without surface diffusion. Conversely, some observations will also be offered concerning the general implication of diffusion, first introduced into simulation studies by Gilmer and Benema [5].
2. Simulation parameters Temperature T and supersaturation A p, were expressed by the dimensionless relations:
to = 4,ss/(2kr), with ~bs~ representing the bond energy between two nearest neighbors on the surface, and
/3 = Alz/kT= ln(Cf/Cf0 ), where Cf denotes the concentration of the fluid and Cf0 its equilibrium concentration. Simulation runs were carried out at to--1.0 and 2.0, below the roughening temperature toR = 0.77. Nucleation dips are reduced by surface diffusion, to /3--0.05 for to = 1.0 and /3 = 1.0 for to--2.0 [8]. On the other hand, to avoid kinetic roughening we need fl < to so that A/z ~: ~bss. Thus fl values have been chosen between 0.05 and 0.25 for to = 1.0 and in the range 1.0-1.8 for to-- 2.0. We used the kinetic Ising model with a stress-free (100) simple cubic Kossel crystal face under the solid-on-solid condition [6,7]. This is the simplest possible model which stillcontains the essential configurational degrees of freedom for building up new crystal layers. In runs with diffusion, the jump length X~ was set to one lattice unit for consistency with BCF theory [8] and had 8 possible directions, namely, to any of the 4 lateral or 4 diagonal neighbors; going up or down a step was also permitted.
3. Simulation procedures It was assumed that the rate of deposition corresponds to the impingement frequency of the fluid
particles on the surface and is therefore proportional to
p + = Ce/Cfo = exp/3. For evaporation and diffusion, the same activation energy was assumed but they depend on local configuration at the site, such that p( N)d = p ( N ) - = exp(4 - 2 N ) to, where N denotes the number of lateral nearest neighbors [5]. The procedure due to Liu, van der Eerden and Bennema [9], that simulates growth at low temperature with high efficiency, is adopted. The type and the site where an event happened were selected according to the waiting times of different events at all sites, calculated as t += - I n n+/p +,
t-= -In n-/p(N)-,
t d = - I n nd/p(N) d, where n +, n- and n d lie within the interval (0,1] and are drawn from a single pseudo-random number chain (Ref. [10] pp. 16,17). The procedure was implemented in MS FORTRAN 5.0 to double-precision format on a microcomputer. The event with the minimum waiting time was accepted and, after the appropriate changes of site configurations, the program calculated new waiting times for those sites and sorted them again with the others, to maintain the event with the minimum time at the top of the to-do list. After a preset number of Monte Carlo cycles, an output record containing the surface configuration, the numbers of deposition, evaporation and diffusion events, the numbers of kink and step sites, as well as cumulative waiting time would be written to a file. The cycles of waiting time generation and to-do list were updated and then re-commenced. A 30 × 30 surface imposed with the periodic boundary condition (Ref. [10] pp. 22-24) was used. The mean diffusion steps at /3 = 0 have previously been determined via extrapolation [I1] as do--8.1 for to = 1.0 and d0 = 8.7 for to = 2.0. Within the ranges of /3, therefore, the surface used suffices in size to prevent any severe coupling of a lone nucleus with its image across a periodic boundary. However, if more clusters had formed, the dynamic behavior of the system would be distorted [12]. Hence the pro-
T.B. Tang/Applied Surface Science 92 (1996) 257-260
gram always halted at 50% coverage, below which no second cluster would appear [11]. In total 180 simulation runs were executed. Ninety proceeded with surface diffusion (X~ = 1), ninety without (X~ = 0). At least three runs were conducted for either case at each value of /3, and their outputs discarded unless the data deduced from them agree with one another to within 10%.
0.6
•=
259
(a)
0.5
;
0.4
~
:
1.2
0.8
:
:
1.6
beta
,~- 0.5 i(b)
4. R e s u l t s a n d d i s c u s s i o n
Qualitatively, I examined the time development of coverage in each run by viewing its output records of surface configuration on the computer monitor as consecutive pages of a Word document, in quick successions by pressing the 'Page Down' key. The coverage was then plotted against cumulative waiting time for representative runs. The size n+ of the growth cluster was taken from the smallest cluster at the onset of rapid growth, and n o of the critical nucleus from the largest that had dissolved before
8O
(a)
600
.~ 4 0 " 20
:
0
:
:
:
:
1ib80 0.8
1.2
1.6
beta
.N
60
ca 40
20
0 0
0.2
0.1
0.3
beta Fig. 1. Sizes of growth cluster (<>. ×) and critical nucleus ( n , zx) for (a) to = 2.0 and (b) to = 1.0 (<>, [] : X~ = 0 and ×, ,x : X s = 1).
0.2
• 0
:
t
0,1
I 0.2
t 0.3
beta
Fig. 2. Time average of kinks per step site for (a) to = 2.0 and (b) to= 1.0([]: Xs=O and n: Xs= 1).
this onset. Data for to = 1.0 and 2.0 are summarized in Figs. l a and lb, which show that, in the presence of diffusion, both n+ and n o were reduced, the latter to a greater degree. This result matches the previous observation of Van Leeuwen and Van der Eerden [12]. The Wulff theorem and the G i b b s - T h o m s o n formula [4] state respectively that the product n o/32 should be a constant and that n+ = 2n o. A plot of n o versus /3-2 for to = 2.0 is indeed fairly linear and indicates t h a t n 0 / 3 2 = 26 ( X s -- 0) and 10 ( X s = 1). For t o = 1.0, the slopes are n 0 / 3 2 = 0 . 4 2 ( X s = 0 ) and 0.25 ( X s = 1). However, n + / n o = n + / 3 2 / n o / 3 2 = 3.6 and 29 for X s - 1, to = 1.0 and 2.0 respectively. Only the cases of X~ -- 0 satisfy the condition that n + / n o = 2.0. BCF [4] argue that an infinitely long step acts as a continuous-line sink for single particles and its density of kinks during growth, maintained by thermal fluctuation, should be independent of supersaturation. However, our results (Fig. 2) support this deduction only for X~ = 0. In the diffusional cases, the kink densities are lower, especially so at low supersaturation. This result is in qualitative agreement with the prediction from an analytical model [13].
260
T.B. Tang~Applied Surface Science 92 (1996) 257-260
5. Conclusion The size of growth cluster is reduced by the introduction of surface diffusion, and that for critical nucleus, to a more marked extent. This explains the shortening of the nucleation dip found earlier [8]. However, by the same token it is no longer true that one is always half of the other, as derived by BCF, who also predict an invariant density of kinks along surface steps. Again, the presence of diffusion invalidates this prediction, the more so at a higher temperature. At high temperature and supersaturation, kinetic roughening occurs as indicated by a change in slope of the jump frequency [11], and equilibrium structures such as critical nuclei lose meaning. The conclusion is therefore that, contrary to previous views [12], diffusion, by opening up more energy decay paths in the configurational space, brings the surface further away from thermodynamic equilibrium. Lastly, the work reported here may be the first that has been undertaken in its entirety on a personal computer. A 486DX-33 was used, and a single job took overnight to a full day. Since the current price for this or some faster model has fallen under US$1200 (in Hong Kong as of mid 1995), more than one can crunch the numbers independently to cut down the aggregate computational time. On a more advanced level, with the help of communication utilities such as the PVM (Parallel Virtual Machine) v3.5 from Oak Ridge Nat. Lab., we can run the program on a network of PowerPCs or PowerMacs that support Unix. Parallelization methods are of
course applicable to a number of Monte Carlo and molecular dynamic problems [14]. The trend towards 'parallelomicrocomputing' may perhaps extend to other areas of computational physics, that likewise place light demand on math software library.
References [1] A.A. Chernov and J. Lewis, J. Phys. Chem. Solids 28 (1967) 2185. [2] H. Miiiler-Krumhhaar, in: Monte Carlo Methods in Statistical Physics, 2nd ed., Ed. K. Binder (Springer, Berlin, 1986) pp. 282, 290. [3] U. Landman, Dynamics of lnterfaciai Processes: Simulations of Growth and Stability of Films, Junctions and Nanostructurcs, presented at ICSFS-7. [4] W.K. Burton, N. Cabrera and F.C. Frank, Phil. Trans. R. Soc. London A 243 (1951) 299. [5] G.H. Gilmer and P. Bennema, J. Appl. Phys. 43 (1972) 1347. [6] C. van Leeuwen and P. Bennema, Surf. Sci. 51 (1975) 109. [7] J.D. Weeks and G.H. Gilmer, Advan. Chem. Phys. 40 (1979) 157. [8] V.K.W. Cheng, E.C.M. Tang and T.B. Tang, J. Cryst. Growth 96 (1989) 293. [9] G.Z. Liu, J.P. van der Eerden and P. Bennema, J. Cryst. Growth 58 (1982) 152. [10] K. Binder, in: Monte Carlo Methods in Statistical Physics, 2nd ed., Ed. K. Binder (Springer, Berlin, 1986) pp. 1-45. [11] E.C.M. Tang and T.B. Tang, Ferroelectrics 142 (1993) 9. [12] C. van Leeuwen and J.P. van der Eerden, Surf. Sci. 64 (1977) 237. [13] J. Zhang and G.H. Nancollas, J. Cryst. Growth 106 (1990) 181. [14] D.W. Heermann and A.N. Burkitt, Parallel Algorithms in Computational Science (Springer, Berlin, 1991).