- Email: [email protected]
A lattice kinetic Monte Carlo code for the description of vacancy diffusion and self-organization in Si
Nuclear Instruments and Methods in Physics Research B 148 (1999) 262±267
A lattice kinetic Monte Carlo code for the description of vacancy diusion and self-organization in Si A. La Magna b
a,*
, S. Coa a, L. Colombo
b
a CNR-IMETEM, Stradale Primosole 50, I-95121 Catania, Italy INFM and Dipartimento di Scienza dei Materiali, Universita degli Studi di Milano, via Emanueli 15, I-20126 Milano, Italy
Abstract We have investigated the mechanism driving void evolution in crystalline Si by means of 'on lattice' kinetic Monte Carlo (MC) simulations. The implementation of an ecient algorithm allowed us to compare, on a macroscopic time scale, the predictions of simulation codes using three dierent binding models. We demonstrate that the use of a modi®ed Ising model, taking also into account second neighbour interaction, results in a MC code embodying the main physical features on vacancy (V) interaction, derived by Tight Binding Molecular Dynamics (TBMD) simulations. In particular the dependence of cluster stability and reactivity on size and geometrical con®guration can be eciently taken into account. We have found that, when this model is used, the ripening process is also driven by the migration of small V clusters, and not only by free Vs as predicted by less re®ned models. This feature strongly modi®es the kinetics of growth and the basic mechanism of void ripening. Ó 1999 Published by Elsevier Science B.V. All rights reserved. PACS: 66.30.L; 61.72.Q; 02.50.N; 71.15.D Keywords: Voids; Ripening; Monte Carlo methods; Defects
1. Introduction The nucleation and growth of vacancy (V) aggregates in crystalline Si have recently drawn a large scienti®c and technological interest. For instance, formation of voids is observed in as-grown Czochralski Si ingots [1] as a result of V supersaturation occurring when the crystal is cooled down from the melting point. Voids can also be formed by light ion (H,He) implantation at high
* Corresponding author. Tel.: +39 095 591212; fax: +39 095 7139154; e-mail: [email protected]
¯uences and post-annealing [2], and have several applications in semiconductor processing. However the mechanism ruling the formation of voids from a given supersaturation of point defects is still under debate. Proper combination of Molecular Dynamics (MD) calculations and Monte Carlo (MC) simulations could provide precious insights into this topic. Indeed a MC approach [3,4] can be applied to bridge the gap between the microscopic length and time scales of MD calculation and the domain of experiments. However, this task can be reliably met only if the main physical features, coming from MD, are implemented in the MC code in a feasible way.
0168-583X/98/$ ± see front matter Ó 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 5 8 3 X ( 9 8 ) 0 0 7 9 8 - 8
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
Recently, several features [5,6] of V±V interaction have arisen from MD. First of all, accurate Tight Binding Molecular Dynamics (TBMD) [5] simulations show that the binding energy for the last added particle Eb (n) (i.e. the energy gain obtained by adding one V to an existing cluster with (n)1)V) is non-monotonic and depends on the assumed pattern of growth. Furthermore, analysis of the stability of a V±V complex by MD calculations [6] based on the Stillinger±Weber potential shows that the interaction between Vs extends at least to the second neighbours in c-Si. This last result is also con®rmed by the fact that annealing and reorientation for a Si divacancy [7] can only be understood when this extended interaction is taken into account. In this work we present a kinetic MC (KMC) code in which the main features of V energetics, including the eect of lattice structure, cluster topology and extended V interactions, are taken into account. The predictions of this code are compared to those obtained when using a MC program with a less re®ned V binding model.
263
2. A re®ned model for V interaction We have implemented a novel Ôon latticeÕ KMC code with the following eective atomic description of V diusion and clustering in Si: (i) the particles are allowed to randomly jump to any of their 4 neighbour sites. The jumping rate is tuned to reproduce the estimated diusivity; (ii) a bound state with binding energy Eb1 forms between two Vs sitting at two nearest neighbour (nn) sites (Fig. 1(a)). This V interaction follows from the reduction of the total number of dangling bonds and hence of the total energy; (iii) a bound state with binding energy Eb2 forms between two Vs sitting at two second neighbour (sn) sites (Fig. 1(b)). This interaction, predicted by MD calculations [6], is also necessary to explain experimental data on divacancy annealing and reorientation [7]. The interaction might be caused by a rearrangement of the two dangling bonds
Fig 1. Schematic description of the interaction between Vs sitting: (a) at nearest neighbour position; (b) at second neighbour position; (c) for a system containing three bound vacancies. E1b , E2b and 2E1b are the total binding energies of the three systems.
264
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
associated to the Si atom in the lattice location between the two Vs; (iv) this sn interaction is screened when a third V approaches the lattice site between the two, since the dangling bond-mediated interaction disappears (Fig. 1(c)). Formally, the model maps an extended Ising model (EIM) described by a binding energy: X E1 X E2 Sk b b Si Sj 1 Sl Sm ; Eb 2 2 hi;ji hhl;mii where hi; ji indicates that i and j are nearest neighbour sites while hhl; mii indicates that l and m are second neighbour sites. The variables Si are equal to 1 if the site i is occupied by a vacancy and equal to 0 otherwise. Since sn interaction could be screened, as described in (iv), we have 0 if Sk 1; 2: Eb2 Sk 2 Eb if Sk 0;
The reliability of the three models presented has been tested by comparing the binding energy calculated using the CBM (dashed line in Fig. 2(a)), the IBM (dashed line in Fig. 2(b)) and the EIM (solid line in Fig. 2(b)) with those derived by accurate TBMD calculations, taken by Ref. [5] (+ and solid line in Fig. 2(a)). In the comparison we have considered the case of hexagonal ring cluster pattern of growth (i.e. V clusters grown removing Si atoms from adjacent 6-membered bond rings present in c-Si). The following parameters have been chosen: Eb1 1:4 eV and Eb2 0:8 eV. The monotonic behaviour of CBM interpolates the values obtained from the TBMD simulation, but it fails in describing the oscillations of Eb (n) caused by the insertion of the growing cluster in the Si lattice and by the relaxation of the internal surface. Using the IBM we obtain a non-monotonic be-
where k identi®es the lattice site between l and m. In order to verify if the proper description of V interaction signi®cantly aects V agglomeration, we have also implemented other MC codes using less re®ned V binding models: the Ising binding model (IBM) and the continuum binding model (CBM). The IBM is obtained from the expression (1) simply neglecting sn interaction (i.e. Eb2 0). In the CBM the topology of the cluster and its structural relation with the underlying lattice is neglected [4]. In this case Eb is proportional to the change in the surface energy when a vacancy is added to an existing cluster: Ebc n K1 K2 n2=3 ÿ n ÿ 1
2=3
;
3
where the two constants K1 and K2 are usually derived by a smoothed ®t to MD data for clusters of dierent sizes. The implementation of CBM [4] diers from that of on lattice models since: (i) free Vs move through isotropic random jumps with a jumping distance Rj ; (ii) both Vs and vacancy clusters are point-like object and a capture radius Rc is introduced to describe their interaction. In the following we use the same parameter for the CBM reported in Ref. [4]: K1 3.65 eV, K2 5.15 eV while Rc and Rj are equal to the second neighbour distance in Si.
Fig 2. Binding energy for a cluster of n-vacancies calculated by TBMD (+ and solid line) and by the CBM (dashed line); (b) binding energy for a cluster calculated by the EIM (solid line) and according to the IBM (dashed line).
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
haviour of Eb (n) which is linked to the diamond structure of the underlying lattice. The stability of the clusters with closed rings (n 6, 10, 14, 18) is recovered but the great stability of the clusters with n 5, 8, 12, 16 members in not accounted for. Finally we note that these last features emerge through EIM which gives the best approximation for the result of TBMD simulation. 3. Eects on vacancy evolution and ripening Our MC codes implements proper accelerated algorithms as described in detail in Ref. [8]. The (cubic) simulation box contain 103 ±106 atoms and cyclic boundary conditions have been imposed. For a given system, after the same initialization has been chosen by picking randomly the Vs in the simulation cell, the dynamical evolution has been followed using the codes based on the dierent models.
265
The amount of information that can be obtained from MC simulations is signi®cantly higher when using on lattice models. With the CBM only cluster sizes and positions are accessible while the internal cluster arrangement could only be investigated when EIM or IBM are considered. In order to obtain a reliable graphic representations of the dierences between the kinetic evolution in the IBM and EIM cases, in Fig. 3 we report a comparison between the simulation at high temperature (T 800°C) of two small and dense V systems ruled by the two dierent on lattice binding models. In both cases no free Vs are present at the times shown here. The system ruled by the IBM (snapshots a-b-c) is frozen in the coalescence stage: note that the con®guration remain identical in the time interval explored apart for a movement of a V in two equivalent positions. On the other hand, the system ruled by the EIM (snapshots d-e-f) shows a clear tendency toward a self-organization of Vs in bigger clusters.
Fig 3. Snapshots of V clusters obtained by the IBM (a-b-c) and the EIM (d-e-f) after 106 , 5 ´ 106 and 107 Monte Carlo steps (MCS).
266
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
Two clear features already emerge from the EIM simulation: (i) the second neighbour interaction results in an eective mobility of small clusters, due to the successive detachment and recapture of single Vs lying at the cluster periphery (note the divacancy joining the bigger aggregate in the sequence of the ®gures); (ii) the tendency of clusters with intermediate dimensions to quickly arrange their surface along h1 1 1i planes. In order to investigate the evolution in the cluster size distribution we have performed a long time simulation on larger systems containing 106 sites, with dierent V densities q in the range 1018 ± 1020 and for temperatures up to 1000°C. In Fig. 4 the distributions calculated using the three dierent models after 1012 MC steps are reported for a system with q 1020 cmÿ3 and T 900°C. Signi®cant dierences in the distributions obtained for the three dierent models are evident in this ®gure. In the case of CBM (Fig. 4(a)) the distribution is almost symmetric and shows a peak corresponding to the cluster size dominant at the time considered. Both the spread of the distribution and the size of biggest V cluster are small when they are compared with those obtained using on lattice models. The emerging features of the distribution for the IBM (Fig. 4(b)) are the clear peaks at sizes (6, 10, 14,. . .) corresponding to the close ring aggregates. This behaviour, which produces a severe bottleneck for ripening of voids, appears during the time evolution soon after the dissolution of smaller clusters (till V5 ) when the clusters with size n 7, 8, 9, 11, 12 become the main source of free particles which drive the ripening process. The kinetics appears strongly modi®ed in the EIM case (Fig. 4(c)). The distribution shows an increased spreading in the cluster size and reveals the presence of very large clusters (containing as much as 61 Vs). We stress again that the cause of this features is the mobility of small V aggregates which play an essential role in the cluster self-organization and which boosts considerably the ripening process. The results that we have presented have several implications on our current understanding and description of the ripening of voids in Si. First of all, the ripening does not occur only by the stan-
dard Ostwald ripening [3] mechanism, but rather through a self-organization of microvoids mediated by the eective mobility of small clusters. Furthermore, the presence of particularly stable V clusters, which produces a severe bottleneck for void evolution using IBM, is completely overwhelmed by the motion of small clusters induced
Fig 4. Cluster size distribution after 1012 MCS using: (a) the CBM; (b) the IBM; (c) the EIM.
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
by the sn interaction. This implies that the use of extended interactions not only produces a better description of the static V cluster energetics, in agreement with the TBMD calculations, but has also important consequences on the dynamic evolution. Acknowledgements We would like to acknowledge E. Rimini for several useful discussions. This work has been partially supported by Progetto 5% Microelettronica.
267
References [1] H. Yamagishi, I. Fusegawa, N. Fujimaki, M. Katayama, Semicond. Sci. Technol. A 135 (1992) 1. [2] S.M. Myers, G.A. Petersen, C.H. Seager, J. Appl. Phys. 80 (1996) 3717; V. Raineri, M. Saggio, Appl. Phys. Lett. 71 (1997) 1673. [3] M. Strobel, K.-H. Heinig, W. Moller, Computational Materials Science 494 (1997) 1. [4] M. Jaraiz, G.H. Gilmer, J.M. Poate, T. Diaz de la Rubia, Appl. Phys. Lett. 68 (1996) 409. [5] A. Buongiorno, L. Colombo, T. Diaz de la Rubia, Europhysics Letters 43 (1998) 695. [6] A. Buongiorno, L. Colombo, Phys. Rev. B 57 (1997) 8767. [7] J.A. Van Vechten, Phys. Rev. B 33 (1986) 2674. [8] A. La Magna, S. Coa, Comp. Mat. Sci., submitted.
A lattice kinetic Monte Carlo code for the description of vacancy diusion and self-organization in Si A. La Magna b
a,*
, S. Coa a, L. Colombo
b
a CNR-IMETEM, Stradale Primosole 50, I-95121 Catania, Italy INFM and Dipartimento di Scienza dei Materiali, Universita degli Studi di Milano, via Emanueli 15, I-20126 Milano, Italy
Abstract We have investigated the mechanism driving void evolution in crystalline Si by means of 'on lattice' kinetic Monte Carlo (MC) simulations. The implementation of an ecient algorithm allowed us to compare, on a macroscopic time scale, the predictions of simulation codes using three dierent binding models. We demonstrate that the use of a modi®ed Ising model, taking also into account second neighbour interaction, results in a MC code embodying the main physical features on vacancy (V) interaction, derived by Tight Binding Molecular Dynamics (TBMD) simulations. In particular the dependence of cluster stability and reactivity on size and geometrical con®guration can be eciently taken into account. We have found that, when this model is used, the ripening process is also driven by the migration of small V clusters, and not only by free Vs as predicted by less re®ned models. This feature strongly modi®es the kinetics of growth and the basic mechanism of void ripening. Ó 1999 Published by Elsevier Science B.V. All rights reserved. PACS: 66.30.L; 61.72.Q; 02.50.N; 71.15.D Keywords: Voids; Ripening; Monte Carlo methods; Defects
1. Introduction The nucleation and growth of vacancy (V) aggregates in crystalline Si have recently drawn a large scienti®c and technological interest. For instance, formation of voids is observed in as-grown Czochralski Si ingots [1] as a result of V supersaturation occurring when the crystal is cooled down from the melting point. Voids can also be formed by light ion (H,He) implantation at high
* Corresponding author. Tel.: +39 095 591212; fax: +39 095 7139154; e-mail: [email protected]
¯uences and post-annealing [2], and have several applications in semiconductor processing. However the mechanism ruling the formation of voids from a given supersaturation of point defects is still under debate. Proper combination of Molecular Dynamics (MD) calculations and Monte Carlo (MC) simulations could provide precious insights into this topic. Indeed a MC approach [3,4] can be applied to bridge the gap between the microscopic length and time scales of MD calculation and the domain of experiments. However, this task can be reliably met only if the main physical features, coming from MD, are implemented in the MC code in a feasible way.
0168-583X/98/$ ± see front matter Ó 1999 Published by Elsevier Science B.V. All rights reserved. PII: S 0 1 6 8 - 5 8 3 X ( 9 8 ) 0 0 7 9 8 - 8
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
Recently, several features [5,6] of V±V interaction have arisen from MD. First of all, accurate Tight Binding Molecular Dynamics (TBMD) [5] simulations show that the binding energy for the last added particle Eb (n) (i.e. the energy gain obtained by adding one V to an existing cluster with (n)1)V) is non-monotonic and depends on the assumed pattern of growth. Furthermore, analysis of the stability of a V±V complex by MD calculations [6] based on the Stillinger±Weber potential shows that the interaction between Vs extends at least to the second neighbours in c-Si. This last result is also con®rmed by the fact that annealing and reorientation for a Si divacancy [7] can only be understood when this extended interaction is taken into account. In this work we present a kinetic MC (KMC) code in which the main features of V energetics, including the eect of lattice structure, cluster topology and extended V interactions, are taken into account. The predictions of this code are compared to those obtained when using a MC program with a less re®ned V binding model.
263
2. A re®ned model for V interaction We have implemented a novel Ôon latticeÕ KMC code with the following eective atomic description of V diusion and clustering in Si: (i) the particles are allowed to randomly jump to any of their 4 neighbour sites. The jumping rate is tuned to reproduce the estimated diusivity; (ii) a bound state with binding energy Eb1 forms between two Vs sitting at two nearest neighbour (nn) sites (Fig. 1(a)). This V interaction follows from the reduction of the total number of dangling bonds and hence of the total energy; (iii) a bound state with binding energy Eb2 forms between two Vs sitting at two second neighbour (sn) sites (Fig. 1(b)). This interaction, predicted by MD calculations [6], is also necessary to explain experimental data on divacancy annealing and reorientation [7]. The interaction might be caused by a rearrangement of the two dangling bonds
Fig 1. Schematic description of the interaction between Vs sitting: (a) at nearest neighbour position; (b) at second neighbour position; (c) for a system containing three bound vacancies. E1b , E2b and 2E1b are the total binding energies of the three systems.
264
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
associated to the Si atom in the lattice location between the two Vs; (iv) this sn interaction is screened when a third V approaches the lattice site between the two, since the dangling bond-mediated interaction disappears (Fig. 1(c)). Formally, the model maps an extended Ising model (EIM) described by a binding energy: X E1 X E2 Sk b b Si Sj 1 Sl Sm ; Eb 2 2 hi;ji hhl;mii where hi; ji indicates that i and j are nearest neighbour sites while hhl; mii indicates that l and m are second neighbour sites. The variables Si are equal to 1 if the site i is occupied by a vacancy and equal to 0 otherwise. Since sn interaction could be screened, as described in (iv), we have 0 if Sk 1; 2: Eb2 Sk 2 Eb if Sk 0;
The reliability of the three models presented has been tested by comparing the binding energy calculated using the CBM (dashed line in Fig. 2(a)), the IBM (dashed line in Fig. 2(b)) and the EIM (solid line in Fig. 2(b)) with those derived by accurate TBMD calculations, taken by Ref. [5] (+ and solid line in Fig. 2(a)). In the comparison we have considered the case of hexagonal ring cluster pattern of growth (i.e. V clusters grown removing Si atoms from adjacent 6-membered bond rings present in c-Si). The following parameters have been chosen: Eb1 1:4 eV and Eb2 0:8 eV. The monotonic behaviour of CBM interpolates the values obtained from the TBMD simulation, but it fails in describing the oscillations of Eb (n) caused by the insertion of the growing cluster in the Si lattice and by the relaxation of the internal surface. Using the IBM we obtain a non-monotonic be-
where k identi®es the lattice site between l and m. In order to verify if the proper description of V interaction signi®cantly aects V agglomeration, we have also implemented other MC codes using less re®ned V binding models: the Ising binding model (IBM) and the continuum binding model (CBM). The IBM is obtained from the expression (1) simply neglecting sn interaction (i.e. Eb2 0). In the CBM the topology of the cluster and its structural relation with the underlying lattice is neglected [4]. In this case Eb is proportional to the change in the surface energy when a vacancy is added to an existing cluster: Ebc n K1 K2 n2=3 ÿ n ÿ 1
2=3
;
3
where the two constants K1 and K2 are usually derived by a smoothed ®t to MD data for clusters of dierent sizes. The implementation of CBM [4] diers from that of on lattice models since: (i) free Vs move through isotropic random jumps with a jumping distance Rj ; (ii) both Vs and vacancy clusters are point-like object and a capture radius Rc is introduced to describe their interaction. In the following we use the same parameter for the CBM reported in Ref. [4]: K1 3.65 eV, K2 5.15 eV while Rc and Rj are equal to the second neighbour distance in Si.
Fig 2. Binding energy for a cluster of n-vacancies calculated by TBMD (+ and solid line) and by the CBM (dashed line); (b) binding energy for a cluster calculated by the EIM (solid line) and according to the IBM (dashed line).
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
haviour of Eb (n) which is linked to the diamond structure of the underlying lattice. The stability of the clusters with closed rings (n 6, 10, 14, 18) is recovered but the great stability of the clusters with n 5, 8, 12, 16 members in not accounted for. Finally we note that these last features emerge through EIM which gives the best approximation for the result of TBMD simulation. 3. Eects on vacancy evolution and ripening Our MC codes implements proper accelerated algorithms as described in detail in Ref. [8]. The (cubic) simulation box contain 103 ±106 atoms and cyclic boundary conditions have been imposed. For a given system, after the same initialization has been chosen by picking randomly the Vs in the simulation cell, the dynamical evolution has been followed using the codes based on the dierent models.
265
The amount of information that can be obtained from MC simulations is signi®cantly higher when using on lattice models. With the CBM only cluster sizes and positions are accessible while the internal cluster arrangement could only be investigated when EIM or IBM are considered. In order to obtain a reliable graphic representations of the dierences between the kinetic evolution in the IBM and EIM cases, in Fig. 3 we report a comparison between the simulation at high temperature (T 800°C) of two small and dense V systems ruled by the two dierent on lattice binding models. In both cases no free Vs are present at the times shown here. The system ruled by the IBM (snapshots a-b-c) is frozen in the coalescence stage: note that the con®guration remain identical in the time interval explored apart for a movement of a V in two equivalent positions. On the other hand, the system ruled by the EIM (snapshots d-e-f) shows a clear tendency toward a self-organization of Vs in bigger clusters.
Fig 3. Snapshots of V clusters obtained by the IBM (a-b-c) and the EIM (d-e-f) after 106 , 5 ´ 106 and 107 Monte Carlo steps (MCS).
266
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
Two clear features already emerge from the EIM simulation: (i) the second neighbour interaction results in an eective mobility of small clusters, due to the successive detachment and recapture of single Vs lying at the cluster periphery (note the divacancy joining the bigger aggregate in the sequence of the ®gures); (ii) the tendency of clusters with intermediate dimensions to quickly arrange their surface along h1 1 1i planes. In order to investigate the evolution in the cluster size distribution we have performed a long time simulation on larger systems containing 106 sites, with dierent V densities q in the range 1018 ± 1020 and for temperatures up to 1000°C. In Fig. 4 the distributions calculated using the three dierent models after 1012 MC steps are reported for a system with q 1020 cmÿ3 and T 900°C. Signi®cant dierences in the distributions obtained for the three dierent models are evident in this ®gure. In the case of CBM (Fig. 4(a)) the distribution is almost symmetric and shows a peak corresponding to the cluster size dominant at the time considered. Both the spread of the distribution and the size of biggest V cluster are small when they are compared with those obtained using on lattice models. The emerging features of the distribution for the IBM (Fig. 4(b)) are the clear peaks at sizes (6, 10, 14,. . .) corresponding to the close ring aggregates. This behaviour, which produces a severe bottleneck for ripening of voids, appears during the time evolution soon after the dissolution of smaller clusters (till V5 ) when the clusters with size n 7, 8, 9, 11, 12 become the main source of free particles which drive the ripening process. The kinetics appears strongly modi®ed in the EIM case (Fig. 4(c)). The distribution shows an increased spreading in the cluster size and reveals the presence of very large clusters (containing as much as 61 Vs). We stress again that the cause of this features is the mobility of small V aggregates which play an essential role in the cluster self-organization and which boosts considerably the ripening process. The results that we have presented have several implications on our current understanding and description of the ripening of voids in Si. First of all, the ripening does not occur only by the stan-
dard Ostwald ripening [3] mechanism, but rather through a self-organization of microvoids mediated by the eective mobility of small clusters. Furthermore, the presence of particularly stable V clusters, which produces a severe bottleneck for void evolution using IBM, is completely overwhelmed by the motion of small clusters induced
Fig 4. Cluster size distribution after 1012 MCS using: (a) the CBM; (b) the IBM; (c) the EIM.
A.La Magna et al. / Nucl. Instr. and Meth. in Phys. Res. B 148 (1999) 262±267
by the sn interaction. This implies that the use of extended interactions not only produces a better description of the static V cluster energetics, in agreement with the TBMD calculations, but has also important consequences on the dynamic evolution. Acknowledgements We would like to acknowledge E. Rimini for several useful discussions. This work has been partially supported by Progetto 5% Microelettronica.
267
References [1] H. Yamagishi, I. Fusegawa, N. Fujimaki, M. Katayama, Semicond. Sci. Technol. A 135 (1992) 1. [2] S.M. Myers, G.A. Petersen, C.H. Seager, J. Appl. Phys. 80 (1996) 3717; V. Raineri, M. Saggio, Appl. Phys. Lett. 71 (1997) 1673. [3] M. Strobel, K.-H. Heinig, W. Moller, Computational Materials Science 494 (1997) 1. [4] M. Jaraiz, G.H. Gilmer, J.M. Poate, T. Diaz de la Rubia, Appl. Phys. Lett. 68 (1996) 409. [5] A. Buongiorno, L. Colombo, T. Diaz de la Rubia, Europhysics Letters 43 (1998) 695. [6] A. Buongiorno, L. Colombo, Phys. Rev. B 57 (1997) 8767. [7] J.A. Van Vechten, Phys. Rev. B 33 (1986) 2674. [8] A. La Magna, S. Coa, Comp. Mat. Sci., submitted.