Monte Carlo simulations of ion implantation in crystalline targets

Materials Science and Engineering, B2 (1989) 33-40

33

Monte Carlo Simulations of Ion Implantation in Crystalline Targets* A. M. M A Z Z O N E

lstituto LAMEL, ('onsiglio Nazionale delle Recerche, Via Casmgnoli 1, 40126 Bologna (Italy) (Received June 2. 1988)

Abstract In this work, we illustrate different types of Monte Carlo method. The results of such calculations describe the passage of ions with energy in the megaelectronvolt range through a crystalline target. The immediate consequence of such a passage is discussed. 1. Introduction The passage of heavy ions through matter is characterized on a short time (100 ns or less) scale by three main phases. (i) The first phase, which is often called the prompt collisional phase, governs the ion path and is dominated by the energy exchanges by collisions between the incident ion and the target electrons and atoms. In this phase, according to the ion mass and energy, a regime of high electronic excitation is generated and recoils with energy sometimes comparable with the primary energy are set in motion. (ii) The phase successive to the ion penetration is the thermalization of the excited region. In this phase the surplus energy is dissipated via electron-electron, electron-phonon collisions, secondary recoil generation and coupling between low energy recoils and lattice phonons. If the collisional phase is important to the understanding of the relationship between the mechanisms of energy deposition and the formation of simple defects, it has often been pointed out (for a complete list of references see ref. 1) that the more significant evolution of damage occurs during the thermal phase. (iii) Molecular dynamics simulations of low energy damage events show that, at the end of the thermalization phase, even if lattice equilibrium *Paper presented at the Symposium on Deep Implants: I=undamentals and Applications at the E-MRS Spring Meeting, Strasbourg, May 31 -June 2, 1988. 0921-5107/89/$3.50

conditions are regained, interstitials and vacancies are still left in the lattice. The final phase is a process of self-annealing damage due to the thermally activated motions of the defects in the otherwise perfect lattice. Often, there is no neat separation in time between the phases. Rather, they differ in the physical mechanisms which are activated. Our primary aim is to describe, through results of Monte Carlo simulations, how the crystalline nature of the target intervenes in these phases. The collisional phase is by far the best understood of the three and the considerable volume of results available allows us to draw a complete picture which is only flawed by quantitative uncertainties. For the other two phases, we have only fragments of information and the field appears to be open to further developments. However, it will also be seen that one aspect which greatly jeopardizes the understanding of the stages of thermalization and self-annealing is the quantitative uncertainty in the stopping mechanisms active during the collisional phase.

2. The collisional phase 2.1. The ion range Obviously one of the more spectacular effects of the crystalline lattice is the formation of channelled trajectories. It is generally acknowledged that it was the Monte Carlo simulation by Robinson and Oen which in the 1960s initiated the present-day studies of channelling. However, the approach more widely used is the classic theory of Lindhard. According to this theory, the atomic strings forming the crystalline lattice exert two antithetic actions. On one hand, there is a steering effect which governs the ion trajectories through glancing collisions. On the other hand, there is a randomizing action which results in chaotic uncorrelated motions dominated by violent large-angle collisions. The transition from © Elsevier Sequoia/Printed in The Netherlands

34 one type of behaviour to the other is determined by a critical value of the angle formed by the ion trajectory and the direction of the atomic strings. Furthermore the theory of Lindhard identified two fundamental modes of channelling, i.e. axial and planar channelling. The corresponding trajectories consist of reflections at the atomic strings (axial channelling) or of oscillatory motions between planes (planar channelling; axially channelled trajectories confined to the original channel were subsequently defined as hyperchannelled). Early Monte Carlo simulations completed this picture by actually showing the typical enhanced penetration of channelled ion paths for incidence conditions aligned with axes or planes. These simulations were also instrumental in the evaluation of statistically uncorrelated dechannelling caused by processes such as thermal lattice vibrations, electron scattering, native oxide coating and lattice disorder (a list of references can be found in the review paper by Gemmell [2]). More recently, however, great attention is being paid to implants in a nominally random direction. In such a case, according to common expectations, the behaviour of a random and a crystalline target should be identical. On the contrary it has been emphasized [3-6] that in a misaligned crystal the ion distribution shows, with respect to a structureless target, an enhanced indepth penetration which depends in a critical manner on the target offset angles. In silicon the result holds for heavy and light projectiles in an energy range from a few kiloelectronvolts to the near submegaelectronvolt range and leads to straggling parameters up to 60% larger than those of a structureless target [6]. Furthermore, along the same lines, it has been also shown that dechannelling factors (such as native oxide coating or beam misalignment), which tend to suppress channelling for implants into major crystallographic directions, have a reverse effect for nominally random orientations

[7-9]. Two main reasons can be envisaged as accounting for these effects. On one hand, there is the complex structure of a crystalline lattice where many minor channels can be located in the vicinity of a major crystallographic direction (a complete map of the minor channels near the ~001) axis in silicon has been recently reported [10]). On the other hand, the possibility of rechannelling has to be accounted for, i.e. extra scatter-

ing centres may become active which feed ions into the channels. As the governing action of the atomic strings consists of coherent nuclear scattering, we expect the disappearance of channelling at high energies when the stopping regime is dominated by electrons. To elucidate the full meaning of this complex problem, we report some indicative calculations (Figs. 1-4~. For the simulation method used in such calculations and in those in "Fables 1 and 2 see ref. 9. In Fig. l a comparison is made between depth profiles of 0.5 and 5 MeV As + ions for different orientations of the silicon target. (Unless otherwise stated, "random" corresponds to a nominally random direction, '~planar" corresponds to a beam direction aligned with the (110) planes and "axial" corresponds to the (100) axis. The beam has" 0 ° divergency.) At the lower energy the axial peak is located at a depth which is approximately four times the depth of the random peak and the planar profile shows the oscillations induced by the reflections at atomic planes. At 5 MeV the three peaks fall at a shorter spatial separation and the oscillations disappear from the planar profile. Similar dechannelling effects are observed in the lateral distribution. For the cases in Fig. 1 the direction of the incident beam is normal to the target surface. Consequently an ion can penetrate along the lateral direction only after a large-angle collision which consistently deflects the ion path. This requirement is opposite to channelling and a long lateral range is indicative of early dechannelling. Figures 2 and 3 report the depth distribution of As + ions in sections located at a different radial distance R from the ion incidence point. It is seen that, for the axial and planar case at 0.5 MeV, most of the ions are concentrated near the incidence point, whereas a non-negligible fraction with longer lateral range is noticed in the random case. At 5 MeV the differences in the lateral distributions are attenuated. In Table 1 the first-order moments of the ion distribution in the crystalline target are compared with those in a structureless target. It can be seen that in the crystalline target the ion distribution is constantly broader along the in-depth direction and narrower along the radial direction. Both effects, which indicate the formation of channelled paths, become less evident at high energies. In these calculations, following a common line [3, 11, 12], the theory of Lindhard is adopted for electronic losses. This theory does not include

35

Fig. 1. Monte Carlo calculations of depth profiles of As + ions in a crystalline silicon target for the conditions shown.

Fig. 2. Monte Carlo calculations of depth profiles of As + ions in different lateral sections. R is the radial distance from the ion incidence point.

36

Fig. 3. Monte Carlo calculations of depth profiles of As* ions in different lateral sections. R is the rzldial distance, fr~m~ the i~I~ incidence point.

TABLE 1 Average projected range Rp, in-depth straggling parameter dRp and lateral straggling parameter dR r (cylindrical coordinates) of As + ions implanted in a crystalline and in a structureless silicon target Target orientation

Rp (/~m)

(ktmdRp)

dR, (/~m)

As + ions." energy, 0.5 M e V Axial Planar Random

1.16 0.67 0.29

0.33 0.39 0.116

0.028 0.057 0.075

Structureless target

0.28

0.09

0.085

As + ions: energy, 5 M e V Axial Planar Random

4.08 3.67 2.8

0.82 0.62 0.39

0.20 0.27 0.34

Structureless target

2.7

(I.34

0.37

any correlation between the stopping mechanisms and the lattice structure. Consequently the disappearance of channelling effects is hardly surprising. However, reduced electronic stopping for channelled particles is a well-known fact. For intermediate ion velocities the effect has been described in terms of the dependence of the

stopping mechanisms on the impact parameter according to the Firsov formulation [13]. At high energies, this reduction may be heuristically explained as due to a reduction in the close collision contribution associated with the decrease in the electron density at the centre of the channel. A representative, if not strictly correct, idea of such effects is given in Fig. 4. In the figure a comparison is made between a random and axially channelled ((t 10) axis of a silicon target) B + profile. For the axial cases the random stopping has been scaled to different amounts to convey the idea of the different theoretical approaches (in the figure, a is the ratio of channelled to random stopping). The value a = 0.5 corresponds to a naive application of Lindhard's equipartition rule. The value a = 0.8 is obtained from a dielectric formalism which accounts for the non-uniform distribution of the valence electrons within the channel section [14] (m both references quoted in ref. 14, good agreement is shown between the theories and the experimental energy losses of protons in silicon [15 ]). The use of the Firsov formula leads to an effective a value in the range of 0.8 and the corresponding profile has

37

Fig. 4. Monte Carlo calculations of depth profiles of B + ions in crystalline silicon (incidence conditions, random and aligned with the ~110) axis).

been omitted for clarity in the figure. The profiles in Fig. 4 lead to the obvious result that the lattice effects are fully restored if a consistent stopping reduction is operated. However, the uncertain limits of validity of the different theories and the complex evaluation of the effective charge render a realistic estimate of a problematic. The problem is further aggravated by the scarcity of experimental data which are limited to relatively low energies [16] or to very light particles [15]. Consequently, in the megaelectronvolt range, the presence of channelling for either major or minor channels as well as rechannelling effects appear to be important unanswered questions. In the following, within the aim of a worst-case estimate, we use random stopping irrespective of the target orientation. As will be seen, non-negligible channelling effects are still perceivable if the impact region is examined on a fine scale. 2.2. High energy recoils During the collisional phase, recoils are generated whose energy is sometimes comparable with the primary energy. These atoms, travelling to a great distance from the ion track, may often split

the cascade into a more complex structure containing independent subcascades (or branches). Great attention is being paid to the formation of branches because of their influence on the crystalline-to-amorphous transition. In fact, the formation of branches, by giving a broader spatial extension to the cascade, alters the distribution of elastic energy density. Such a density is the fundamental factor in the process of amorphization and we expect different amorphous 'structures according to the branching yield. Obviously a channelled trajectory matches ideally the requirement of a long path and Monte Carlo simulations have dwelt on the subject of recoil channelling. Recoil channelling occurs because the displaced atom is guided by collisions within a crystal channel. Consequently the effect has a complex dependence on the initial energy and direction of the moving particle and is hardly predictable. Although instances of recoil channelling are shown in many Monte Carlo simulations, the calculations point to a more frequent occurrence of quasi-channelling [17, 18]. In a quasi-channelled motion an atom experiences considerable energy loss with the atoms forming

38

the channel walls (and the range is shorter than in proper channelling) but in spite of these losses it retains its original direction. Results representative of the occurrence of recoil channelling and their dependence on the target orientation are illustrated in Table 2. (In the table, "branches" indicate recoils with initial energy greater than 40 keV. The visual inspection of the cascade shows, in fact, that this energy is generally sufficient to guarantee the formation of an independent subcascade.) In the table a comparison is also made between a crystalline and a structureless target. Two points must be emphasized. First, we recall that, when an ion travels under channelling conditions for most of its path, the probability of a collision with a large energy exchange is remote and the formation of one branch is a rare event. As shown in the table, the branching yield for the axial and planar incidence remains appreciably lower with respect to the random case even for energies in the megaelectronvolt range. This result leads to the conclusion that channelling effects, although less conspicuous in the ion path, still intervene by reducing the few collisions which determine the formation of branches and the cascade structure. The second noticeable point is that recoils with a lower energy sometimes have a longer range (see for instance the random and the planar case for 5 MeV). The result is clearly indicative of the formation of channelled trajectories. This suggestion is also supported by the comparison with the structureless target where the recoil range is regularly shorter than in the crystalline target even if the path is initiated with an equal energy. TABLE 2 Structure of the average cascade in a silicon target. Formation of high energy (E = 40 keV) recoils ("branches") Target orientation

N u m b e r of branches

Maximum energy of branches (MeV)

Maximum vector range o f branches (~m)

A s + ions: energy, 0.5 M e V Axial 0.1 Planar 0.63 Random 1.7

0.12 0.15 0.18

0.17 0.2 0.45

Structureless target

0.18

0.35

As + ions:energy, 5 M e V Axial 0.5 Planar 2.2 Random 5.0

0.7 0.8 0.9

0.5 1.8 1.5

Structureless target

0.9

0.9

1.6

6.5

3. The thermalization of the excited region It is by now well-accepted modelling that the primary energy is partitioned between nuclear and electronic losses. Nuclear losses are absorbed by the formation of successive generations of recoils and by the excitation of lattice phonons, which ultimately drain the surplus energy. Very few studies are available on the interaction between phonons and cascade atoms [19, 2(1] so that the processes of de-excitation of the atomic system still have many obscure points. However, in spite of these uncertainties, nuclear losses call be calculated with great accuracy and this removes a fundamental element of ambiguity. For implants in the megaelectronvolt range, electronic losses generally prevail and the problem acquires a far greater complexity. In such a case, in fact, the dissipation of the energy in the medium depends on the competing effect of the electron-electron and the electron-phonon interactions respectively. The rate of such processes is determined by the initial energy of the excited electron, electron-phonon collisions tending to be the more efficacious for a very low electron energy. Furthermore, while electron-electron collisions rapidly disperse the energy in the lattice, possibly contributing to the formation of charged defects and to Coulomb explosions, electronphonon collisions may lead sometimes tt) a noticeable increase in the local temperature [21, 22]. For highly excited and dense electron systems, more complex processes, such as electron-hole collisions [23 l, may also intervene. In addition to the possibility of different relaxation modes, relevant differences are expected to arise from the target orientation. In fact, under channelling conditions, owing to the smaller oscillations of the ion paths, the deposited energy is compactly located in a volume thinner than for random incidence. For instance, the detailed output of the calculations indicates that, if a silicon target bombarded with 2 MeV B + is subdivided to a depth of approximately 0.5 ¢~m in volume elements of 10 A x 10 A x l(t A, the energy deposited in the w~lumes along the ion path is approximately 100 eV for axial channel-. ling. This value decreases by approximately one order of magnitude for random incidence. In such a case, however, the lower level of electronic excitation is balanced by the increase in the number of volumes where the deposited energy is not vanishingly small.

39

To state a quantitative relationship between these figures and the relaxation mechanisms is a doubtful operation as we should accurately know the number of the electrons intervening in each collision and where such electrons are actually located. The discussion in Section 2 on the evaluation of close collisions in channelled stopping clearly indicates an unsatisfactory knowledge of such elements. On this basis, we can only speculate that, owing to the lower level of electronic excitation for random incidence, electron-phonon coupling will be more frequent on random than on channelled paths. In this latter case, especially for paths distant from the empty regions in the channel centre, the possibility of collisions between excited electrons is not to be discarded.

4. Damage self-annealing As shown by molecular dynamics simulations, on completion of the thermalization phase, simple defects or small clusters are left in the lattice. As defects are often mobile at room temperature, a phase of cascade quenching (i.e. processes such as defect recombination and clustering) occurs through normal, thermally activated diffusion. To describe the defect migration within a Monte Carlo approach, two lines have been

TABLE 3

followed. In one case the Monte Carlo simulation describes the annealing of the isolated cascade following, jump by jump, the motions of the defects in the otherwise perfect lattice. The second type of Monte Carlo simulation is a stochastic approach where the exact dynamics of the entire damaged region is constructed by averaging over a statistical selection of representative states from the appropriate thermodynamic ensemble. In Table 3 a comparison is made between the behaviour of the isolated cascade and the behaviour of the entire damaged region for representative values of the fractional disorder, indicated in the table as displacements per atom (the computational procedures are those reported in refs. 24 and 25 respectively, where also a detailed list of references on the simulation methods can be found). For the isolated cascade the table reports the two limiting cases when only interstitials or only vacancies are mobile. The comparison between them shows that the mobility of either of these species is essential for cluster formation, which appears to be enhanced when vacancies are mobile. We recall that a typical cascade structure has a depleted core while the interstitials form a diffuse halo around the vacancy centre. Consequently, when defects are allowed to move, vacancies cluster more easily because of their close proximity to each other.

Damage self-annealing cascades formed by 5 M e V As + implanted in a silicon target

Target orientation

Lsolated cascade quenching (% of same character as isolated displacements)

Maximum cluster size mumber of atoms of the mobile species)

Recombination

Clustering

0.59 0.60 0.63

0.11 0.13 0.13

3.0 4.0 4.0

0.60 0.61 0.61

0.2l) 0.30 0.30

6.0 6.0 7.0

Only interstitials mobile Axial Planar Random

Only vacancies mobile Axial Planar Random

Damaged region Displacements per atom

Maximum cluster size (number oJ'atoms of the mobile species)

0.01 O. 1

9.0 1800.0

40

The defect mobility is a very delicate point for implants in the megaelectronvolt range. In fact, many theoretical and experimental studies indicate that close collisions with the excitation of the outer and inner shells lead to the formation of multiply charged defects [26, 27]. As the defect mobility depends on the defect charge state, we may expect a prevailing interstitial or vacancy clustering according to the defect charge state generated during a close collision or its immediate aftermath. Furthermore, as close collisions are more likely to occur on random than on channelled paths, we expect the target orientation (which according to the table has no effect on the cascade quenching) to be important in damage annealing through its action on the defect mobility. The simulation of the entire damaged region (Table 3) emphasizes the importance of the damage evolution occurring within the isolated cascade. In fact the maximum cluster size appears to increase appreciably above that in the isolated cascade only for a disorder level approaching the range of the crystalline-to-amorphous transition. The conclusions obtained from this are similar to those from the thermalization, i.e. the main source of doubts can be traced back to the occurrence of close collisions, which is a fundamental element of the stopping mechanisms active during the collisional phase. 5. Conclusion In this paper the peculiar features of a crystalline target have been briefly reviewed and discussed. It has been shown that, even if the immediate consequence of the different atomic arrangements corresponding to different orientations can be fully accounted for, the uneven distribution of electrons within a crystal poses many questions about the mode of energy deposition and the mechanisms of energy relaxation. Although a complete formalism is available for problems of this type, for conditions close to equilibrium a great deal of understanding has still to be gained about the perturbed conditions generated by a high energy implant.

References 1 V. I. Protasov and V. G. Chudinov, Radia~ Ef], ,~'~

(1984) 185. 2 D.S. Gemmell, Rev. Mod. l ' h w 16, 19741 t. 3 M. Haulala, in B. R. Appleton. F. H. Eisen and i. W Sigmon (eds.), Materials Research .Soeiet3 ,~vmp. t'ro~. Vol. 45, Material Research SocieLx,. Pittsburgh, PA, t~85, 17. 1()5. 4 0 . S. Oen, Nucl. lnstrum. Methods B, /3 (198,6) 495, 5 M. 1). Giles and J. F. Gibbons, I[,EE Tran~, Elecmm Devices, 32 ( I 0) (1985) 191S. 6 R.G. Wilson, AppI. l'hvs. Leu., ,')d (1980) 2797. 7 A. Desalvo, R. Galloni, R. Rosa and F. Zignani, J, Appl. l'hss., 51 (1980) 1994= 8 M. Hautala, NucL lnstrmn. .~lethods B, 15 (1986) 75. 9 A. M. Mazzone, Philos, Mag. 1,ctl., 55(t 987) 235. 10 J. F. Ziegler and R. F. Lever, ,~ppl. l'hvs, l.eu., 46!1985) ~58. 11 ,1. de Pontcharra, P. Spinelli and M. Bruel, in B. R, Appleton, F. H. Eisen and T. W. Sigmon (eds.',, .~laterial.~ Research Socieo' Syrup, l'roe. '~ol. 45, Material Research Society, Pittsburgh, PA, 1985, p. 115. 12 tl. Wong, E. Deng, N. W. Cheung, P. K. Chu, E. M. Strathman and M. D. Strathman, Nuel. lnstrur~l. Methods B, 21 (1987) 447. 13 A. H. EI-Hoshy and J. F. (;ibbons. l'hs~. Rev.. 17.? 1968) 454. 14 A. l)esalvo and R. Rosa, J. Phys. (, /0(1977) 1595. A. t:. Burenkov, F. lb. Komarov and A. M. Kumakhov. l'hys. Status Solidi B, 20 ( 198(! 417, 15 A. Carnera, G. Della Mea, A, V. l)rigo, S. l,o Russo, P. Mazzoldi and G. G. gentini, Ptn,~ Rev. B, 17 978) ~432. 16 F. Eisen, ('an. J. Phys., 46( 1965; 570. 17 ,I. R. Beeler, Phys. Rev., 150 (1966) 47tl. 18 G. P. Mueller, N. D. Wilsey and M. Rosen, ll;t;I; 7)'arts. :\,ucl. Sci., 29 (1982) 1493. 19 B. W. Dodson and P. A. Taylor, ,I, Mater. Res., 211987) 8/!5. 20 A. M. Mazzone, Nucl. Inslrltm. Methods, in the press. 21 H. Blank, l'hys. Status Solidi A, 10 ( 1975 ) 465. 22 T. A. Tombrello, Nucl, lnstrmn. ~lethods B, 2 1t98,4) 355, 23 M. Combescol, Solid-),lale Eleclron., 3l (1988) 657. 24 A.M. Mazzone, Radiat. t"2.I][2Express, 1 (1987) 129. 25 A. M. Mazzone. I E E E Tran,~. ('omput. Assistal Des.. 4 ii985) 369. 26 I 1. J. Pabst and D. W. Palmer, m S. l)atz, B. R. Appleton and C, l). Moak (eds.), Atomic ('ollisions in Solids, Plenum, New York, 1975, p. 141. 27 ('. H. Wie, T. A. Tombrello and T. Vreeland. Jr,, Phys. Rev. B, 33 (1986) 4083,