- Email: [email protected]
The melting behaviour of small silicon clusters
15 August 1994 PHYSICS ELSEVIER
LETTERS
A
Physics Letters A 191 (1994) 339-345
The melting behaviour of small silicon clusters P. Tchofo Dinda a, G. Vlastou-Tsinganos b, N. Flytzanis b, A.D. Mistriotis b~c a Laboratoire de Physique de I’UniversitPde Bourgogne, 6 Boulevard Gabriel, 21000 Dijon, France b Department of Physics, University of Crete, Heraklion, Crete, Greece 71409 c Research Center of Crete, P.O. Box 1527, Heraklion. Crete, Greece 71110 Received 23 February
1994; revised manuscript received 22 June 1944; accepted for publication Communicated by L.J. Sham
25 June 1994
Abstract We report an analysis of the melting behaviour of small silicon clusters interacting via a nonlinear interatomic potential with four-body terms. The analysis shows, by means of Monte Carlo and molecular dynamics simulations, that the small silicon clusters undergo, in a vacuum, structural changes from a solid rigid state to a liquid-like state. The melting temperature exhibits a strong variation with cluster size.
Small clusters of atoms, namely semiconductor clusters, have received a great deal of attention owing to their interest in condensed matter physics, and to some potential industrial applications. For example, the area of application of semiconductor clusters includes the production of amorphous films by the deposition of clusters on a substrate [ 1,2], or the production of materials with some specific properties by incorporating small clusters in a matrix [ 31. Important information can also be obtained for nucleation and cluster growth, as well as for the structure of the amorphous materials. Two problems are generally considered in the theoretical developments: first, the static and dynamic properties of the cluster, and more particularly the geometrical structure of the atoms in the cluster; and second, the cluster behaviour under some specified external conditions of pressure and temperature. In particular the study of melting properties of small clusters can give us a more complete description of its structural properties (global and local minimum, and dynamical stability of different isomers), but also information on cluster
formation and fragmentation. The nature and character of the melting process in small noble gas clusters modelled by a Lcnnard-Jones potential has been the subject of extensive investigation [4,7]. In Ref. [ 51 an analysis of the melting transition is performed using the method of molecular dynamics coupled with local quenching. Thus a mechanistic description is given of the phase change in terms of the local potential minima accessed in the transition region, the isomerisation pathways and frequencies of interwell transitions. Semiconductor clusters, however, require the introduction of more complicated potentials that take into account the directional character of covalent bonding. A very successful potential to study the dynamical properties of the cluster is the Stillinger-Weber potential [ 8 1, which has been shown to describe satisfactorily the behaviour of the silicon crystal and amorphous materials. Using this potential Blaisten-Barjoras and Levesque [ 91 studied small Si clusters (neutral and positively charged) to determine their minimum energy configuration. The melting was also studied by
0375-9601/94/%07.00 @ 1994 Elsevier Science B.V. All rights reserved SsDIO375-9601(94)00542-7
340
P. TchofoDinda et al. /Physics LettersA 191(1994) 339-345
monitoring the average potential energy as a function of temperature. In most cases however no strong signatures of melting were observed but rather a continuous transition with fluctuations. In order to suitably describe the small clusters of silicon, Mistriotis, Flytzanis and Farantos (MFF) [ lo] have modified the Stillinger-Weber potential. The MFF potential [lo] consists of two-body, three-body, and fourbody terms. The two-body term is constructed to simulate the interaction between the Si (silicon) atoms in the diamond structure, b = A(B/rt
- 1) exp[a/(rij
- R)],
(1)
where rij = ]ri - rj], with rk the coordinate of the kth particle. The physical meaning of R and the values of the parameters A and B will be made more precise later on. The three-body and four-body terms are introduced to describe the directionality of the covalent Si-Si bond. Specifically the form of a three-body term is v, = kjik + kijk + kikja
(2)
where kjik =
lXhjf;:k{
1 -
exP[-Q(cOS8jik
+ f.)‘]}.
(3)
The form of the four-body term is v4 =
gijki
+
gjikl
+ gkijl
+
(4)
glijk,
where
+
the melting point), and also to describe satisfactorily the geometries and energies for the ground states and low metastable structures of Si clusters. The two main improvements over the Stillinger-Weber potential [8] are the possibility to vary independently the amplitude of the three-body term Asand the curvature near the minimum. This is achieved by only an extra parameter, but it allows stronger deviations from the tetrahedral symmetry. The introduced four-body terms (with no extra parameters) become important as the cluster size increases, since the three-body term increases the cluster energy more slowly than expected from ab initio calculations. The purpose of the present Letter is to extend these studies to the melting behaviour of small Si clusters of various sizes, from N = 7 atoms upward, indeed it was shown in Ref. [lo] that this model potential provides the best description of the properties of Si clusters for N > 7. For smaller clusters, linear structures with strong n-bonds are predominant, while the Stillinger-Weber potential [ 8 ] is better suited for bulk structures. In order to facilitate the calculation we use a Monte Carlo simulation algorithm introduced by Creutz [ 111, in which an energy transfer process takes place between the system and an extra degree of freedom, a “demon” with energy ED. The demon is allowed to exchange energy with the system subject to the restriction that ED > 0. One sets up a random walk through configurational space while maintaining a constraint on the total energy E. That is, EC + ED = E = cons&
(COSeji/
+
f)’
+
(WSekil
+
f12]}).
(5)
The functions hj = exp[y/(rij
- R)l
(6)
are cutoff functions, so that the interaction vanishes smoothly for rij b R. It was shown in Ref. [lo] that the parameters A = 16.3 eV, B = 11.581 A4, (Y = 2.095 A, R = 3.771 A, A3 = 4.0 eV, A4 = 47.0 eV, Q = 5.0, y = 2.4 A, allow one to reproduce several geometric and dynamical properties of the Si crystal (namely the energy per atom, the structure of the crystal ground state (diamond structure), the lattice constant, and
(71
where EC is the configurational energy of the system. For each configuration generated, one must check whether the configurational energy EC is reduced or not. If EC is reduced, the new configuration is accepted and the excess energy is given to the demon. However if EC is increased, the new configuration is only accepted if the demon has enough energy. In such a case the energy is taken from the demon. Otherwise no change occurs. In thermal equilibrium, the temperature T is determined from the ensemble average of the demon energy with the following relation, T =
%A
=
(&)/h
(8)
where the subscript CA refers to the Creutz algorithm, kB is the Boltzmann constant. This algorithm has been
P. Tchofo Did
Lo 0.1 0.05
0
M
0
0.1
et al. /Physics Letters A 191(1994) 339-345
Si, 0.2
0.3
0.15 Lo 0.1
0.05 0 ii 0
Si, 1
2
Fig. 1. The bond length fluctuation parameter 6 as a funo tion of the effective temperature for Sir, obtained from the Monte Carlo Creutz algorithm (a), and from the molecular dynamics simulation (b ) . applied successfully by Grimson [ 121 to the study of the melting behaviour of small Lennard-Jones clusters. During the simulation we measure the degree of rigidity of the cluster by using the dimensionless bondlength fluctuation parameter (BLF) defined by ((a) -
2
6 =
N(N-
1) c
i
341
tions. This figure exhibits clearly three temperature regions in which the cluster behaves differently: (i) At low temperatures the BLF is relatively small and increases only slightly and smoothly as the temperature increases, which indicates that the cluster is in a rigid solid state, where the configuration remains near the global minimum (as we will see in the forthcoming discussion); (ii) Next, in a relatively small temperature region, the BLF undergoes a sharp increase, thus indicating a strong increase in the particle mobility. Such a behaviour is usually identified with a cluster melting transition. We define the melting temperature Tm as the temperature at which the melting behaviour begins to manifest itself, that is, Tm (Sir ) = TCA a 1590 K for the case of Fig. la; (iii) Upon heating, the cluster attains a steady state in which the BLF is relatively high, but increases again slightly and smoothly as the temperature increases. This state can therefore be associated with a liquidlike state ’ ’ 2, where the cluster configuration is often away from the global minimum. In order to confirm these results we have also examined the evolution of the BLF 6 versus temperature via the classical MD.5 (molecular dynamics simulations), by putting initial velocities on the atoms in such a way that the initial momentum of the system is zero. Consequently for obtaining the temperature from the MDS we disregarded the contribution of the three translational degrees of freedom to the total kinetic energy of the system. So the temperature was obtained from the MDS in the following simple way, T = TMDSE
2(h) 3kB(N - 1)’
(10)
(rij)')lj2
(rij)
’
(9)
where ( ) denotes an ensemble average over an entire simulation after equilibration. So, 6 serves to indicate structural changes in the cluster. Fig. 1a shows the results obtained via the Creutz algorithm for a cluster of size N = 7, in the vacuum. The starting configuration was the ground state obtained in Ref. [ lo] The data reported are obtained by averaging over 7 x lo6 equilibrium configurations after an initial equilibration period of lo6 configura-
1 In the range of temperature that we consider, the energy introduced in the cluster is not large enough to cause the fragmentation of the cluster; consequently no evaporation of the atoms occurs. 2 We willshowin a forthcoming paper that after the melting transition the distribution of the interatomic distances for the atoms that were initially the nearest neighbours (in the solid state) becomes essentially identical to the distribution for the atoms that were initially situated at the farthest dis tance; which indicates a liquid-state behaviour in the sense that in the liquid the mobility of the atoms is such that any two atoms can lie, from time to time, at the nearest or the farthest distance.
342
P. Tchofo Dinda et al. /Physics Letters A 191(1994)
where the subscript MDS refers to the molecular dynamics simulation. Fig. 1b, which gives the result obtained for 7 x lo6 time steps 3, shows that all general features mentioned for the Creutz algorithm are preserved. The MDS gives us T, (Si, ) = TMDS M 1550 K, which agrees with the results obtained via the Creutz algorithm within 3%. However, an important point to note in doing this comparison is that the Creutz algorithm requires only a small fraction of the amount of calculations required when using the MDS. Consequently for obtaining the melting temperature for larger clusters we have used only the Creutz algorithm. Fig. 2 shows the results obtained for Sis, Si9, Silo, Sill, and Siir, with an increasing number of accepted configurations (that is, respectively 9 x 106, 10 x 106, 11 x 106, 12 x 10” and 13 x 106), so that the results are reproducible. For each of these clusters the melting behaviour manifests itself by an abrupt increase in the BLF in a relatively small temperature region. The melting temperatures for these clusters are shown in Fig. 3. Two outstanding points emerge from this figure: first, we observe that the melting temperature for the clusters is always significantly lower than the melting temperature of the perfect Si crystal (w 2000 K). This is quite conceivable; indeed the large percentage of surface atoms in the cluster (compared to the crystal) leads to significant differences for some interatomic bondings in the cluster, so that one can regard the cluster as a sort of crystal with a maximum of physical imperfections. As it is furthermore well known that the defects in a crystal provide nucleation sites that favour structural phase transitions, it is therefore conceivable that the melting temperature of the cluster is lower than that of the crystal. Second, we observe the surprising behaviour that the melting temperatures for Sis, Sig, Sin,, and Sill are very low compared to the melting temperatures for Sir, and S&z; this result can be understood by invoking a similar behaviour already observed in previous works; for example it is now well recognized that there exist so-called “magic numbers”, which correspond to highly stable clusters with respect to the fragmentation [ 13- 15 1. Similarly our results show that there exist some clusters which are particularly stable with 3 The dynamical simulation time step is chosen to conserve the energy of the system to an accuracy better than 0.1% for the length of the simulation.
339-345
m (a)
Si
0.2 co
0.1
0
0.1
0
0
0.05
0.15
0.1
0
0.05
0.1
0.15
_;f+-jjj;fjy,-j
0
0.05
0.1
0.15
0
0.05
0.2 Lo 0.1
0
0.15
Si
(e) 0.3
0.1
Fl.li 0.06
0.1 0.15 0.2
Fig. 2. The bond length fluctuation parameter S as a function of the effective temperature, obtained from the Monte Carlo Creutz algorithm, for (a) Sis, (b) Sig, (c) Silo, (d) Sill, and (e) Silz.
respect to the melting transition, two of which are Si,, and S&r, as shown in Fig. 3. Such clusters could be called “highly rigid clusters”. In this respect we point out that the analysis of the fragmentation experiments leads to the existence of magic numbers which theoretically appear as having the highest binding energy per atom, where the reference state is that with all atoms separated. In the melting process however the cluster does not break, but visits a part of its configurational space. The rigidity of the cluster therefore depends on the separation of the lowest energy isomers and their
343
P. Tchofo Dinda et al. /Physics Letters A 191(1994) 339-345
1
0.2 0.6 0.4 0.2
0 6
8
10
12
-24.6 -24 -23.6 McAL2NmGYMmMA 1
n Fig. 3. The melting temperature Tmversus cluster size.
0.2
-24.6 -24 -23.6 mcALENE2GYMmMA
T=2639 K cc>
* 0.8
accessibility. The above analysis shows that upon heating the small silicon clusters undergo, in a vacuum, structural changes from a solid rigid state to a liquid-like state (see footnote 2). However an interesting question remains to be considered: what is the physical meaning of the cluster melting phenomenon for a collective entity - such as a cluster - which consists of only a few number of atoms. In order to gain insight into this problem, and more generally into the physical processes which govern the structural changes in the cluster, it is necessary to examine in detail the dynamics of a cluster of size N. A simple way to understand this problem is to regard the cluster as a collective entity embedded within a configurational potential with characteristic troughs (wells) and crests (barriers). Then the question becomes: how does the cluster, as a whole entity, evolve within this configurational potential, as the temperature increases. We, simulate the dynamics by starting with a cluster that is initially relaxed to the ground state. Then we put some initial velocities on the atoms in order to set the cluster in motion. During the dynamics we periodically quench the system, in order to determine the potential wells that the cluster visits and the relative time spent in each well. That is, using the classical pseudo-dynamics relaxation process, we extract energy from the system to relax the cluster by periodically setting the velocities of all atoms to zero after a specified number of time steps. After relaxing, the cluster finds itself in a bound state, which corresponds to a local-energy minimum. Next, the cluster is released in the state where it was found just before the quenching, so that the dynamics continue without being.affected by the quenching process. Fig. 4, which gives the results obtained for
k
0.4
Eli,
0.2 0~ -24.6 -24 -23.6 LmAL2NmGYmmA
Fig. 4. Plot of Fv, the frequency of visitingthe possiblelocal energy minima, for Si7 at (a) T = 576 K, (b) T = 1691 IC,and (c) T = 2639 K. Si7 at different temperatures, shows the local-energy minima corresponding to the wells that the cluster has visited during the dynamics, and gives a fairly good estimation of the the relative time spent in each of these wells. Indeed the total time spent in a well is closely related to the number of times that the cluster finds itself in the bottom of this well, owing to successive relaxation processes. We have shown in Fig. 4, for simplicity, the number of times that the cluster finds itself in the bottom of each well, normalized to the number of times that the cluster finds itself in the ground state. This quantity is denoted in the figures as Fv. Fig. 4a corresponds to the dynamics of the Si7 cluster at T = 576 K, that is, in the solid-like state. We see that throughout the dynamics the cluster visits only the well corresponding to the ground state Ei (Si7 ) = -24.454 eV, which is the global minimum. The cluster is then trapped in the global-minimum well owing to the fact that the initial kinetic energy introduced in the system is not sufficient to permit the cluster to cross the intervening potential barriers. In general, we have found the interesting result that the solid-like state always corresponds to a cluster trapped in the global-minimum well, and that the melting behaviour begins to manifest itself only when the kinetic energy
344
P, Tchofo Dinda et al. /Physics Letters A 191(1994) 339-345
becomes sufficient to permit the cluster to move from the global-minimum well to another well. Fig. 4b, which corresponds to the cluster dynamics during the melting transition (T = 1691 K), shows precisely that the cluster visits not only the globalminimum well but also another well corresponding to the bound state Ez(Sir ) = -24.005 eV. This indicates that the initial kinetic energy introduced in the system is now sufficient to permit the cluster to overcome some potential barriers. The cluster is no longer trapped in the global-minimum well and can therefore visit some other potential wells; the atoms can now execute relatively large-amplitude motions: hence the strong increase in the BLF during the melting, which we described previously. Furthermore we see in Fig. 4b that the cluster spends only a small fraction of its time in the second well; indeed the total time spent in this well represents only 10% of the total time spent in the global-minimum well. Furthermore we have found that the relative time spent in the second well increases as the temperature increases, and that the BLF begins again to increase smoothly (liquid-like state) as soon as the cluster acquires enough kinetic energy to visit a third well. Fig. 4c, which corresponds to a cluster in a liquidlike state at T = 2639 K, shows that during the dynamics the cluster visits four wells. The two first wells correspond to those mentioned above and the two other wells correspond respectively to the bound state energies E3 (Si7) = -23.57966 eV and E4(Si7) = -23.6083 eV. Furthermore we see that the relative time that the cluster spends in the second well has significantly increased compared to the case of Fig. 4b, and this time now represents more than 50% of the total time spent in the global-minimum well. We also note the surprising behaviour that although /Es] < lE41 the time that the cluster spends in the third well represents only a small fraction of the time spent in the fourth well. Such a behaviour is related to the shape of the configuration of the potential energy of the Sir cluster. The most important point that emerges from the above discussion is that only the two or three lowest local energy minima play a prominent role in the melting behaviour. It is therefore clear that the melting temperature is related to the amount of energy that is necessary to permit the cluster to cross the intervening barrier between the global-minimum well (asso-
ciated to the ground state) and the well corresponding to the energy situated just above the ground state. Consequently, it is conceivable that the melting temperature is somehow related to the difference between the two lowest energy minima, that is, the first gap in the spectrum of the local energy minimum. One of course expects small deviations from this criterion depending on the accessibility of the path between two minima. We have found a similar behaviour for the other clusters considered in the present paper; however due to the limited number of clusters examined it would be premature to try to determine a relationship between the first gap of the energy spectrum and the melting temperature; such a relationship would have been very useful for obtaining in a straightforward way the melting temperature for larger clusters without having to go through Monte Carlo simulations. Another interesting problem that remains is to examine how the melting behaviour manifests itself when the clusters are confined within hard spheres, as is the case, for example, when the clusters are embedded inside a solid matrix. Furthermore our results are in qualitative agreement with those of Ref. [ 91 except that the introduction of the bond length fluctuations is much more sensitive in order to monitor structural phase changes. In order to eliminate fluctuations in the 6 versus (ED) plots we had to consider many steps in the Monte Carlo and molecular dynamics calculations. Finally it should be mentioned that in spite of recent advances in the theoretical development of empirical model potentials for small clusters of atoms, such as the MFF potential [ lo] that we have used in the present paper, the results obtained via those model potentials provide in some cases only a crude approximation of the behaviour of real materials. However, as we already mentioned previously, many fundamental properties have been satisfactorily described by those empirical model potentials. In this respect, we think that most of the essential features of the cluster melting phenomenon, that we have described in the present paper, are generic features which should be preserved by some other bona fide empirical model potentials for Si clusters. More specifically, the most important conclusion that emerges from this work is the strong variation of Tm with cluster size. The drop in Fig. 3 is about 900 K (more than 50%) and is clearly outside the possible error in defining Tm (less
P. Tchofo Din& et al. /Physics Letters A 191(1994) 339-345
than 5%). This observation is very important to determine suitably the energy of incident clusters in ionised cluster beam (ICB ) experiments for the fabrication of amorphous films. A suitable choice of cluster masses by differentiation can assure a liquid-like state at the smallest possible temperature, in order to produce uniform films.
345
[41 J. Jellinek, T. L. Beck and R. S. Berry, J. Chem. Phys.
84 (1986) 2783. [51 T. L. Beck and R. S. Berry, J. Chem. Phys. 88 (1988)
3910. [61 F. G. Amar and R. S. Berry, J. Chem. Phys. 85 (1986)
5943. [71 I. L. Garzon, M. Avalos Botja and E. Blaisten-Bajoras,
Phys. Rev. B 40 (1989) 4749. WI F. H. Stillinger and T. A. Weber, Phys. Rev. B 31
(1985) 5262.
One of the authors, P.T.-D., would like to thank the Physics Department of the University of Crete for hospitality. Part of this work is supported by the EEC Science project under grant No. SCl*-CT91-0686.
[lOI A. Mistriotis, N. Flytzanis and S. Farantos, Phys. Rev.
References
1131 L.A. Bloomlield and W. L. Brown Phys. Rev. L&t. 54
[91 E. Blaisten-Bajoras and D. Levesque, Phys. Rev. B 34
(1986) 3910. B 39 (1989) 1212. 1111 M. Creutz, Phys. Rev. L&t. 50 (1983) 1411. 1121 Malcolm J. Grimson, Chem. Phys. L&t. 195 (1992)
92. (1985) 2246.
[ I] T. Takagi, Vat. Sci. Tech. A 2 (1984) 382. [2] I. Yamada, in: Semiconductor and semimetals, Vol. 32A, ed. J. Pankove (1984). [ 31 T. Inoshita and H. Watanabe, in: Microclusters, Proc. First NEC Symposium, Hakone, Japan, eds. S. Sugano, Y. Hishima and S. Ohnishi ( 1986) p. 281.
114 L.A. Bloomfield, M. E. Geusic, T. J. McIlrath, M.
F. Jarrold, R. R. Freeman and W. L. Brown in: Microclusters, Proc. First NEC Symposium, Hakone, Japan eds. S. Sugano, Y. Hishima and S. Ohnishi (1986) p. 238. [15 K. Raghavachari and V. Logovinsky Phys. Rev. Lett. 55 (1985) 28.53.
LETTERS
A
Physics Letters A 191 (1994) 339-345
The melting behaviour of small silicon clusters P. Tchofo Dinda a, G. Vlastou-Tsinganos b, N. Flytzanis b, A.D. Mistriotis b~c a Laboratoire de Physique de I’UniversitPde Bourgogne, 6 Boulevard Gabriel, 21000 Dijon, France b Department of Physics, University of Crete, Heraklion, Crete, Greece 71409 c Research Center of Crete, P.O. Box 1527, Heraklion. Crete, Greece 71110 Received 23 February
1994; revised manuscript received 22 June 1944; accepted for publication Communicated by L.J. Sham
25 June 1994
Abstract We report an analysis of the melting behaviour of small silicon clusters interacting via a nonlinear interatomic potential with four-body terms. The analysis shows, by means of Monte Carlo and molecular dynamics simulations, that the small silicon clusters undergo, in a vacuum, structural changes from a solid rigid state to a liquid-like state. The melting temperature exhibits a strong variation with cluster size.
Small clusters of atoms, namely semiconductor clusters, have received a great deal of attention owing to their interest in condensed matter physics, and to some potential industrial applications. For example, the area of application of semiconductor clusters includes the production of amorphous films by the deposition of clusters on a substrate [ 1,2], or the production of materials with some specific properties by incorporating small clusters in a matrix [ 31. Important information can also be obtained for nucleation and cluster growth, as well as for the structure of the amorphous materials. Two problems are generally considered in the theoretical developments: first, the static and dynamic properties of the cluster, and more particularly the geometrical structure of the atoms in the cluster; and second, the cluster behaviour under some specified external conditions of pressure and temperature. In particular the study of melting properties of small clusters can give us a more complete description of its structural properties (global and local minimum, and dynamical stability of different isomers), but also information on cluster
formation and fragmentation. The nature and character of the melting process in small noble gas clusters modelled by a Lcnnard-Jones potential has been the subject of extensive investigation [4,7]. In Ref. [ 51 an analysis of the melting transition is performed using the method of molecular dynamics coupled with local quenching. Thus a mechanistic description is given of the phase change in terms of the local potential minima accessed in the transition region, the isomerisation pathways and frequencies of interwell transitions. Semiconductor clusters, however, require the introduction of more complicated potentials that take into account the directional character of covalent bonding. A very successful potential to study the dynamical properties of the cluster is the Stillinger-Weber potential [ 8 1, which has been shown to describe satisfactorily the behaviour of the silicon crystal and amorphous materials. Using this potential Blaisten-Barjoras and Levesque [ 91 studied small Si clusters (neutral and positively charged) to determine their minimum energy configuration. The melting was also studied by
0375-9601/94/%07.00 @ 1994 Elsevier Science B.V. All rights reserved SsDIO375-9601(94)00542-7
340
P. TchofoDinda et al. /Physics LettersA 191(1994) 339-345
monitoring the average potential energy as a function of temperature. In most cases however no strong signatures of melting were observed but rather a continuous transition with fluctuations. In order to suitably describe the small clusters of silicon, Mistriotis, Flytzanis and Farantos (MFF) [ lo] have modified the Stillinger-Weber potential. The MFF potential [lo] consists of two-body, three-body, and fourbody terms. The two-body term is constructed to simulate the interaction between the Si (silicon) atoms in the diamond structure, b = A(B/rt
- 1) exp[a/(rij
- R)],
(1)
where rij = ]ri - rj], with rk the coordinate of the kth particle. The physical meaning of R and the values of the parameters A and B will be made more precise later on. The three-body and four-body terms are introduced to describe the directionality of the covalent Si-Si bond. Specifically the form of a three-body term is v, = kjik + kijk + kikja
(2)
where kjik =
lXhjf;:k{
1 -
exP[-Q(cOS8jik
+ f.)‘]}.
(3)
The form of the four-body term is v4 =
gijki
+
gjikl
+ gkijl
+
(4)
glijk,
where
+
the melting point), and also to describe satisfactorily the geometries and energies for the ground states and low metastable structures of Si clusters. The two main improvements over the Stillinger-Weber potential [8] are the possibility to vary independently the amplitude of the three-body term Asand the curvature near the minimum. This is achieved by only an extra parameter, but it allows stronger deviations from the tetrahedral symmetry. The introduced four-body terms (with no extra parameters) become important as the cluster size increases, since the three-body term increases the cluster energy more slowly than expected from ab initio calculations. The purpose of the present Letter is to extend these studies to the melting behaviour of small Si clusters of various sizes, from N = 7 atoms upward, indeed it was shown in Ref. [lo] that this model potential provides the best description of the properties of Si clusters for N > 7. For smaller clusters, linear structures with strong n-bonds are predominant, while the Stillinger-Weber potential [ 8 ] is better suited for bulk structures. In order to facilitate the calculation we use a Monte Carlo simulation algorithm introduced by Creutz [ 111, in which an energy transfer process takes place between the system and an extra degree of freedom, a “demon” with energy ED. The demon is allowed to exchange energy with the system subject to the restriction that ED > 0. One sets up a random walk through configurational space while maintaining a constraint on the total energy E. That is, EC + ED = E = cons&
(COSeji/
+
f)’
+
(WSekil
+
f12]}).
(5)
The functions hj = exp[y/(rij
- R)l
(6)
are cutoff functions, so that the interaction vanishes smoothly for rij b R. It was shown in Ref. [lo] that the parameters A = 16.3 eV, B = 11.581 A4, (Y = 2.095 A, R = 3.771 A, A3 = 4.0 eV, A4 = 47.0 eV, Q = 5.0, y = 2.4 A, allow one to reproduce several geometric and dynamical properties of the Si crystal (namely the energy per atom, the structure of the crystal ground state (diamond structure), the lattice constant, and
(71
where EC is the configurational energy of the system. For each configuration generated, one must check whether the configurational energy EC is reduced or not. If EC is reduced, the new configuration is accepted and the excess energy is given to the demon. However if EC is increased, the new configuration is only accepted if the demon has enough energy. In such a case the energy is taken from the demon. Otherwise no change occurs. In thermal equilibrium, the temperature T is determined from the ensemble average of the demon energy with the following relation, T =
%A
=
(&)/h
(8)
where the subscript CA refers to the Creutz algorithm, kB is the Boltzmann constant. This algorithm has been
P. Tchofo Did
Lo 0.1 0.05
0
M
0
0.1
et al. /Physics Letters A 191(1994) 339-345
Si, 0.2
0.3
0.15 Lo 0.1
0.05 0 ii 0
Si, 1
2
2
6 =
N(N-
1) c
i
341
tions. This figure exhibits clearly three temperature regions in which the cluster behaves differently: (i) At low temperatures the BLF is relatively small and increases only slightly and smoothly as the temperature increases, which indicates that the cluster is in a rigid solid state, where the configuration remains near the global minimum (as we will see in the forthcoming discussion); (ii) Next, in a relatively small temperature region, the BLF undergoes a sharp increase, thus indicating a strong increase in the particle mobility. Such a behaviour is usually identified with a cluster melting transition. We define the melting temperature Tm as the temperature at which the melting behaviour begins to manifest itself, that is, Tm (Sir ) = TCA a 1590 K for the case of Fig. la; (iii) Upon heating, the cluster attains a steady state in which the BLF is relatively high, but increases again slightly and smoothly as the temperature increases. This state can therefore be associated with a liquidlike state ’ ’ 2, where the cluster configuration is often away from the global minimum. In order to confirm these results we have also examined the evolution of the BLF 6 versus temperature via the classical MD.5 (molecular dynamics simulations), by putting initial velocities on the atoms in such a way that the initial momentum of the system is zero. Consequently for obtaining the temperature from the MDS we disregarded the contribution of the three translational degrees of freedom to the total kinetic energy of the system. So the temperature was obtained from the MDS in the following simple way, T = TMDSE
2(h) 3kB(N - 1)’
(10)
(rij)')lj2
(rij)
’
(9)
where ( ) denotes an ensemble average over an entire simulation after equilibration. So, 6 serves to indicate structural changes in the cluster. Fig. 1a shows the results obtained via the Creutz algorithm for a cluster of size N = 7, in the vacuum. The starting configuration was the ground state obtained in Ref. [ lo] The data reported are obtained by averaging over 7 x lo6 equilibrium configurations after an initial equilibration period of lo6 configura-
1 In the range of temperature that we consider, the energy introduced in the cluster is not large enough to cause the fragmentation of the cluster; consequently no evaporation of the atoms occurs. 2 We willshowin a forthcoming paper that after the melting transition the distribution of the interatomic distances for the atoms that were initially the nearest neighbours (in the solid state) becomes essentially identical to the distribution for the atoms that were initially situated at the farthest dis tance; which indicates a liquid-state behaviour in the sense that in the liquid the mobility of the atoms is such that any two atoms can lie, from time to time, at the nearest or the farthest distance.
342
P. Tchofo Dinda et al. /Physics Letters A 191(1994)
where the subscript MDS refers to the molecular dynamics simulation. Fig. 1b, which gives the result obtained for 7 x lo6 time steps 3, shows that all general features mentioned for the Creutz algorithm are preserved. The MDS gives us T, (Si, ) = TMDS M 1550 K, which agrees with the results obtained via the Creutz algorithm within 3%. However, an important point to note in doing this comparison is that the Creutz algorithm requires only a small fraction of the amount of calculations required when using the MDS. Consequently for obtaining the melting temperature for larger clusters we have used only the Creutz algorithm. Fig. 2 shows the results obtained for Sis, Si9, Silo, Sill, and Siir, with an increasing number of accepted configurations (that is, respectively 9 x 106, 10 x 106, 11 x 106, 12 x 10” and 13 x 106), so that the results are reproducible. For each of these clusters the melting behaviour manifests itself by an abrupt increase in the BLF in a relatively small temperature region. The melting temperatures for these clusters are shown in Fig. 3. Two outstanding points emerge from this figure: first, we observe that the melting temperature for the clusters is always significantly lower than the melting temperature of the perfect Si crystal (w 2000 K). This is quite conceivable; indeed the large percentage of surface atoms in the cluster (compared to the crystal) leads to significant differences for some interatomic bondings in the cluster, so that one can regard the cluster as a sort of crystal with a maximum of physical imperfections. As it is furthermore well known that the defects in a crystal provide nucleation sites that favour structural phase transitions, it is therefore conceivable that the melting temperature of the cluster is lower than that of the crystal. Second, we observe the surprising behaviour that the melting temperatures for Sis, Sig, Sin,, and Sill are very low compared to the melting temperatures for Sir, and S&z; this result can be understood by invoking a similar behaviour already observed in previous works; for example it is now well recognized that there exist so-called “magic numbers”, which correspond to highly stable clusters with respect to the fragmentation [ 13- 15 1. Similarly our results show that there exist some clusters which are particularly stable with 3 The dynamical simulation time step is chosen to conserve the energy of the system to an accuracy better than 0.1% for the length of the simulation.
339-345
m (a)
Si
0.2 co
0.1
0
0.1
0
0
0.05
0.15
0.1
0
0.05
0.1
0.15
_;f+-jjj;fjy,-j
0
0.05
0.1
0.15
0
0.05
0.2 Lo 0.1
0
0.15
Si
(e) 0.3
0.1
Fl.li 0.06
0.1 0.15 0.2
Fig. 2. The bond length fluctuation parameter S as a function of the effective temperature, obtained from the Monte Carlo Creutz algorithm, for (a) Sis, (b) Sig, (c) Silo, (d) Sill, and (e) Silz.
respect to the melting transition, two of which are Si,, and S&r, as shown in Fig. 3. Such clusters could be called “highly rigid clusters”. In this respect we point out that the analysis of the fragmentation experiments leads to the existence of magic numbers which theoretically appear as having the highest binding energy per atom, where the reference state is that with all atoms separated. In the melting process however the cluster does not break, but visits a part of its configurational space. The rigidity of the cluster therefore depends on the separation of the lowest energy isomers and their
343
P. Tchofo Dinda et al. /Physics Letters A 191(1994) 339-345
1
0.2 0.6 0.4 0.2
0 6
8
10
12
-24.6 -24 -23.6 McAL2NmGYMmMA 1
n Fig. 3. The melting temperature Tmversus cluster size.
0.2
-24.6 -24 -23.6 mcALENE2GYMmMA
T=2639 K cc>
* 0.8
accessibility. The above analysis shows that upon heating the small silicon clusters undergo, in a vacuum, structural changes from a solid rigid state to a liquid-like state (see footnote 2). However an interesting question remains to be considered: what is the physical meaning of the cluster melting phenomenon for a collective entity - such as a cluster - which consists of only a few number of atoms. In order to gain insight into this problem, and more generally into the physical processes which govern the structural changes in the cluster, it is necessary to examine in detail the dynamics of a cluster of size N. A simple way to understand this problem is to regard the cluster as a collective entity embedded within a configurational potential with characteristic troughs (wells) and crests (barriers). Then the question becomes: how does the cluster, as a whole entity, evolve within this configurational potential, as the temperature increases. We, simulate the dynamics by starting with a cluster that is initially relaxed to the ground state. Then we put some initial velocities on the atoms in order to set the cluster in motion. During the dynamics we periodically quench the system, in order to determine the potential wells that the cluster visits and the relative time spent in each well. That is, using the classical pseudo-dynamics relaxation process, we extract energy from the system to relax the cluster by periodically setting the velocities of all atoms to zero after a specified number of time steps. After relaxing, the cluster finds itself in a bound state, which corresponds to a local-energy minimum. Next, the cluster is released in the state where it was found just before the quenching, so that the dynamics continue without being.affected by the quenching process. Fig. 4, which gives the results obtained for
k
0.4
Eli,
0.2 0~ -24.6 -24 -23.6 LmAL2NmGYmmA
Fig. 4. Plot of Fv, the frequency of visitingthe possiblelocal energy minima, for Si7 at (a) T = 576 K, (b) T = 1691 IC,and (c) T = 2639 K. Si7 at different temperatures, shows the local-energy minima corresponding to the wells that the cluster has visited during the dynamics, and gives a fairly good estimation of the the relative time spent in each of these wells. Indeed the total time spent in a well is closely related to the number of times that the cluster finds itself in the bottom of this well, owing to successive relaxation processes. We have shown in Fig. 4, for simplicity, the number of times that the cluster finds itself in the bottom of each well, normalized to the number of times that the cluster finds itself in the ground state. This quantity is denoted in the figures as Fv. Fig. 4a corresponds to the dynamics of the Si7 cluster at T = 576 K, that is, in the solid-like state. We see that throughout the dynamics the cluster visits only the well corresponding to the ground state Ei (Si7 ) = -24.454 eV, which is the global minimum. The cluster is then trapped in the global-minimum well owing to the fact that the initial kinetic energy introduced in the system is not sufficient to permit the cluster to cross the intervening potential barriers. In general, we have found the interesting result that the solid-like state always corresponds to a cluster trapped in the global-minimum well, and that the melting behaviour begins to manifest itself only when the kinetic energy
344
P, Tchofo Dinda et al. /Physics Letters A 191(1994) 339-345
becomes sufficient to permit the cluster to move from the global-minimum well to another well. Fig. 4b, which corresponds to the cluster dynamics during the melting transition (T = 1691 K), shows precisely that the cluster visits not only the globalminimum well but also another well corresponding to the bound state Ez(Sir ) = -24.005 eV. This indicates that the initial kinetic energy introduced in the system is now sufficient to permit the cluster to overcome some potential barriers. The cluster is no longer trapped in the global-minimum well and can therefore visit some other potential wells; the atoms can now execute relatively large-amplitude motions: hence the strong increase in the BLF during the melting, which we described previously. Furthermore we see in Fig. 4b that the cluster spends only a small fraction of its time in the second well; indeed the total time spent in this well represents only 10% of the total time spent in the global-minimum well. Furthermore we have found that the relative time spent in the second well increases as the temperature increases, and that the BLF begins again to increase smoothly (liquid-like state) as soon as the cluster acquires enough kinetic energy to visit a third well. Fig. 4c, which corresponds to a cluster in a liquidlike state at T = 2639 K, shows that during the dynamics the cluster visits four wells. The two first wells correspond to those mentioned above and the two other wells correspond respectively to the bound state energies E3 (Si7) = -23.57966 eV and E4(Si7) = -23.6083 eV. Furthermore we see that the relative time that the cluster spends in the second well has significantly increased compared to the case of Fig. 4b, and this time now represents more than 50% of the total time spent in the global-minimum well. We also note the surprising behaviour that although /Es] < lE41 the time that the cluster spends in the third well represents only a small fraction of the time spent in the fourth well. Such a behaviour is related to the shape of the configuration of the potential energy of the Sir cluster. The most important point that emerges from the above discussion is that only the two or three lowest local energy minima play a prominent role in the melting behaviour. It is therefore clear that the melting temperature is related to the amount of energy that is necessary to permit the cluster to cross the intervening barrier between the global-minimum well (asso-
ciated to the ground state) and the well corresponding to the energy situated just above the ground state. Consequently, it is conceivable that the melting temperature is somehow related to the difference between the two lowest energy minima, that is, the first gap in the spectrum of the local energy minimum. One of course expects small deviations from this criterion depending on the accessibility of the path between two minima. We have found a similar behaviour for the other clusters considered in the present paper; however due to the limited number of clusters examined it would be premature to try to determine a relationship between the first gap of the energy spectrum and the melting temperature; such a relationship would have been very useful for obtaining in a straightforward way the melting temperature for larger clusters without having to go through Monte Carlo simulations. Another interesting problem that remains is to examine how the melting behaviour manifests itself when the clusters are confined within hard spheres, as is the case, for example, when the clusters are embedded inside a solid matrix. Furthermore our results are in qualitative agreement with those of Ref. [ 91 except that the introduction of the bond length fluctuations is much more sensitive in order to monitor structural phase changes. In order to eliminate fluctuations in the 6 versus (ED) plots we had to consider many steps in the Monte Carlo and molecular dynamics calculations. Finally it should be mentioned that in spite of recent advances in the theoretical development of empirical model potentials for small clusters of atoms, such as the MFF potential [ lo] that we have used in the present paper, the results obtained via those model potentials provide in some cases only a crude approximation of the behaviour of real materials. However, as we already mentioned previously, many fundamental properties have been satisfactorily described by those empirical model potentials. In this respect, we think that most of the essential features of the cluster melting phenomenon, that we have described in the present paper, are generic features which should be preserved by some other bona fide empirical model potentials for Si clusters. More specifically, the most important conclusion that emerges from this work is the strong variation of Tm with cluster size. The drop in Fig. 3 is about 900 K (more than 50%) and is clearly outside the possible error in defining Tm (less
P. Tchofo Din& et al. /Physics Letters A 191(1994) 339-345
than 5%). This observation is very important to determine suitably the energy of incident clusters in ionised cluster beam (ICB ) experiments for the fabrication of amorphous films. A suitable choice of cluster masses by differentiation can assure a liquid-like state at the smallest possible temperature, in order to produce uniform films.
345
[41 J. Jellinek, T. L. Beck and R. S. Berry, J. Chem. Phys.
84 (1986) 2783. [51 T. L. Beck and R. S. Berry, J. Chem. Phys. 88 (1988)
3910. [61 F. G. Amar and R. S. Berry, J. Chem. Phys. 85 (1986)
5943. [71 I. L. Garzon, M. Avalos Botja and E. Blaisten-Bajoras,
Phys. Rev. B 40 (1989) 4749. WI F. H. Stillinger and T. A. Weber, Phys. Rev. B 31
(1985) 5262.
One of the authors, P.T.-D., would like to thank the Physics Department of the University of Crete for hospitality. Part of this work is supported by the EEC Science project under grant No. SCl*-CT91-0686.
[lOI A. Mistriotis, N. Flytzanis and S. Farantos, Phys. Rev.
References
1131 L.A. Bloomlield and W. L. Brown Phys. Rev. L&t. 54
[91 E. Blaisten-Bajoras and D. Levesque, Phys. Rev. B 34
(1986) 3910. B 39 (1989) 1212. 1111 M. Creutz, Phys. Rev. L&t. 50 (1983) 1411. 1121 Malcolm J. Grimson, Chem. Phys. L&t. 195 (1992)
92. (1985) 2246.
[ I] T. Takagi, Vat. Sci. Tech. A 2 (1984) 382. [2] I. Yamada, in: Semiconductor and semimetals, Vol. 32A, ed. J. Pankove (1984). [ 31 T. Inoshita and H. Watanabe, in: Microclusters, Proc. First NEC Symposium, Hakone, Japan, eds. S. Sugano, Y. Hishima and S. Ohnishi ( 1986) p. 281.
114 L.A. Bloomfield, M. E. Geusic, T. J. McIlrath, M.
F. Jarrold, R. R. Freeman and W. L. Brown in: Microclusters, Proc. First NEC Symposium, Hakone, Japan eds. S. Sugano, Y. Hishima and S. Ohnishi (1986) p. 238. [15 K. Raghavachari and V. Logovinsky Phys. Rev. Lett. 55 (1985) 28.53.