PHYSICAL REVIEW B VOLUME 56,NUMBER 5 1 AUGUST 1997-I Point defect movement and annealing in collision cascades K.Nordlund"and R.S.Averback Materials Research Laboratory.University of Illinois at Urbana-Champaign.Urbana,Illinois 61801 (Received 20 December 1996) The effect of collision cascades on preexisting point defects in crystalline materials was studied by simu- lating 5 keV collision cascades in gold,copper,aluminum,platinum,and silicon.The results indicate that collision cascades do not significantly affect interstitials or vacancies outside the liquid core of the cascade, although in the fcc metals the heating of the crystal due to the cascade causes some thermal migration of the interstitials.Within the liquid cascade core,both interstitials and vacancies move towards the center of the molten region when it resolidifies and recombine or cluster there.At elevated temperatures,random jumps of interstitials during the thermal-spike phase can cause significant additional trapping of interstitials in the liquid. In contrast to the annealing effects of preexisting damage in the fcc metals,in silicon the amount of new damage created by a cascade is roughly independent of the number of initial point defects.The difference is attributed to the nature of the bonding in the materials.[S0163-1829(97)05729-9] L.INTRODUCTION We create the point defects in well-defined distributions and follow their motion during the simulations.Although our Ion irradiation methods are of considerable interest in damage distributions are artificial,they allow us to system- many research and practical applications of materials atically probe the effect of cascades on defects at different processing.2 Since most damage produced in materials dur- locations.We focus on two crystalline materials of contrast- ing ion irradiation derives from a complex process occurring ing properties:gold,a dense material with metallic bonding in collision cascades,much research has been devoted to and close-packed fcc crystal structure,and silicon,a co- studying these events.Experimental work has shed some valently bonded material with the relatively loose-packed light on this problem:however,owing to the difficulty of diamond structure.We also report on simulations of cascades resolving defect structures inside materials,these studies in copper and aluminum to determine whether the differ- have had only limited success.Molecular-dynamics (MD) ences are due to the atom mass or crystal structure,and in computer simulation offers an alternative approach,34 which platinum to determine the effect of the melting point.We has proved successful in providing both a qualitative and believe this choice of materials sets the bounds for possible quantitative description of damage production in solids(see, behavior. e.g.,Refs.5-9). The paper is organized as follows.In Sec.II,we present In most of the simulation studies performed so far,the our simulation procedures in detail.In Sec.III,we first initial state of the lattice in which a collision cascade is ini- present and discuss the results of our simulations for the five tiated has been defect free.However,in practice most ion elements separately,and then compare their common fea- irradiation-induced cascades are produced in regions which tures and differences. have been previously damaged by implanted ions.Therefore, to understand the effect of prolonged implantation on ini- II.SIMULATION AND ANALYSIS METHODS tially crystalline samples,it is important to know how a cas- cade affects preexisting defects. A.General principles Relatively few studies to date have focused on predam- In this study we were primarily interested in elucidating aged sample structures.Sayed e have studied a few defect reactions in collision cascades.Therefore,we set up overlapping cascades in Si at energies of the order of a few our simulation conditions in a manner which can be expected hundred eV.Gao et al.12 very recently worked on multiple to provide an as clear view of the processes involved as overlapping 400 ev to 5 keV cascades in a-iron (a bcc possible.The electronic stopping power is neglected,since it metal),and Foreman et al.3 simulated a few I keV cascades is expected to contribute only little to the slowing down of 5 in Cu with preexisting features like vacancy and interstitial keV self-ions in the materials treated here.18 loops.Kapinos and Bacon have studied the effect of high The initial temperature of the simulation cell was 0 K in vacancy concentrations on vacancy loop formation in Cu,Ni, most runs.Since the interstitial in gold migrates very easily and Fe Both Gao and Foreman report that the preexist- even at low temperatures, 1920 using a 0K ambient cell tem- ing features may be partly annealed when overlapped by a perature is advantageous for clearly distinguishing what part new cascade,a fact well known from experiment.' of the defect motion is due to the collision cascade.Experi- There has been much speculation about the motion of ence in the field shows that the initial state of damage pro- point defects in the heat spike and pressure wave associated duction in collision cascades does generally not differ much with cascades.7 However,no systematic investigation of at temperatures roughly below 100 K,except possibly for these effects has been carried out by molecular dynamics, lengths of replacement collision sequences,21 which are few which is the focus of this work. in number and not of interest here. 0163-1829/97/56(5)/2421(11)/S10.00 56 2421 1997 The American Physical Society
Point defect movement and annealing in collision cascades K. Nordlund* and R. S. Averback Materials Research Laboratory, University of Illinois at Urbana-Champaign, Urbana, Illinois 61801 ~Received 20 December 1996! The effect of collision cascades on preexisting point defects in crystalline materials was studied by simulating 5 keV collision cascades in gold, copper, aluminum, platinum, and silicon. The results indicate that collision cascades do not significantly affect interstitials or vacancies outside the liquid core of the cascade, although in the fcc metals the heating of the crystal due to the cascade causes some thermal migration of the interstitials. Within the liquid cascade core, both interstitials and vacancies move towards the center of the molten region when it resolidifies and recombine or cluster there. At elevated temperatures, random jumps of interstitials during the thermal-spike phase can cause significant additional trapping of interstitials in the liquid. In contrast to the annealing effects of preexisting damage in the fcc metals, in silicon the amount of new damage created by a cascade is roughly independent of the number of initial point defects. The difference is attributed to the nature of the bonding in the materials. @S0163-1829~97!05729-9# I. INTRODUCTION Ion irradiation methods are of considerable interest in many research and practical applications of materials processing.1,2 Since most damage produced in materials during ion irradiation derives from a complex process occurring in collision cascades, much research has been devoted to studying these events.3,4 Experimental work has shed some light on this problem; however, owing to the difficulty of resolving defect structures inside materials, these studies have had only limited success. Molecular-dynamics ~MD! computer simulation offers an alternative approach,3,4 which has proved successful in providing both a qualitative and quantitative description of damage production in solids ~see, e.g., Refs. 5–9!. In most of the simulation studies performed so far, the initial state of the lattice in which a collision cascade is initiated has been defect free. However, in practice most ion irradiation-induced cascades are produced in regions which have been previously damaged by implanted ions. Therefore, to understand the effect of prolonged implantation on initially crystalline samples, it is important to know how a cascade affects preexisting defects. Relatively few studies to date have focused on predamaged sample structures. Sayed et al.10,11 have studied a few overlapping cascades in Si at energies of the order of a few hundred eV. Gao et al.12 very recently worked on multiple overlapping 400 eV to 5 keV cascades in a-iron ~a bcc metal!, and Foreman et al.13 simulated a few 1 keV cascades in Cu with preexisting features like vacancy and interstitial loops. Kapinos and Bacon have studied the effect of high vacancy concentrations on vacancy loop formation in Cu, Ni, and Fe.14,15 Both Gao and Foreman report that the preexisting features may be partly annealed when overlapped by a new cascade, a fact well known from experiment.16 There has been much speculation about the motion of point defects in the heat spike and pressure wave associated with cascades.17 However, no systematic investigation of these effects has been carried out by molecular dynamics, which is the focus of this work. We create the point defects in well-defined distributions and follow their motion during the simulations. Although our damage distributions are artificial, they allow us to systematically probe the effect of cascades on defects at different locations. We focus on two crystalline materials of contrasting properties: gold, a dense material with metallic bonding and close-packed fcc crystal structure, and silicon, a covalently bonded material with the relatively loose-packed diamond structure. We also report on simulations of cascades in copper and aluminum to determine whether the differences are due to the atom mass or crystal structure, and in platinum to determine the effect of the melting point. We believe this choice of materials sets the bounds for possible behavior. The paper is organized as follows. In Sec. II, we present our simulation procedures in detail. In Sec. III, we first present and discuss the results of our simulations for the five elements separately, and then compare their common features and differences. II. SIMULATION AND ANALYSIS METHODS A. General principles In this study we were primarily interested in elucidating defect reactions in collision cascades. Therefore, we set up our simulation conditions in a manner which can be expected to provide an as clear view of the processes involved as possible. The electronic stopping power is neglected, since it is expected to contribute only little to the slowing down of 5 keV self-ions in the materials treated here.18 The initial temperature of the simulation cell was 0 K in most runs. Since the interstitial in gold migrates very easily even at low temperatures,19,20 usinga0K ambient cell temperature is advantageous for clearly distinguishing what part of the defect motion is due to the collision cascade. Experience in the field shows that the initial state of damage production in collision cascades does generally not differ much at temperatures roughly below 100 K, except possibly for lengths of replacement collision sequences,21 which are few in number and not of interest here. PHYSICAL REVIEW B VOLUME 56, NUMBER 5 1 AUGUST 1997-I 0163-1829/97/56~5!/2421~11!/$10.00 2421 © 1997 The American Physical Society 56
2422 K.NORDLUND AND R.S.AVERBACK 56 For a straightforward comparison between the different TABLE I.Some results of the cascade events.The root-mean- materials,we used similar simulation conditions for all ele- square (rms)distance of interstitials from the cell center (Au,Cu, ments.The computational cells had a size of roughly Al,and Pt)or the center of the interstitials (Si)is given for the 120x120X 120 A3,and contained about 90 000 atoms, initial and final distributions of interstitials,based on the Wigner- which were initially arranged in the proper crystal structure Seitz cell analysis.Unless otherwise noted,the defects were ini- of the material.Periodic boundary conditions were used,and tially arranged spherically around the center of the cell,as ex- plained in the text. the temperature of the three outermost atom layers was scaled toward 0 K to dissipate excess energy.The down- Element Event Defects (int.,vac.) rms distance (A) wards scaling factor of the temperature was restricted to be number Initial Final Initial Final less than 0.1 during a single time step.A linkcell calculation22 and variable time step23 were also employed to Au 1 0,0 4,4 46.9 speed up the simulations. 50.0 51,1 44.7 41.3 The primary set of runs examined interstitial motion.For 3 50,50 41,41 48.7 49.8 these,either 50 interstitials and 50 vacancies or only 50 in- Y 50,50 45,45 48.7 50.3 terstitials were introduced in the cell.The number of defects 5a 50,50 39,39 48.7 54.1 chosen corresponds roughly to the number of Frenkel pairs which can be expected to form from a previous 20-30 keV Au° 6 20.0 21,1 15.1 20.3 cascade in the same region.The defects were introduced > 50,0 51,1 23.9 21.0 randomly using a Gaussian probability distribution centered 8 200,0 201,1 35.2 30.5 at the middle of the simulation cell for vacancies and cen- 9 20,20 5,5 25.5 45.8 tered on a spherical cell 35 A from the middle for intersti- tials.The width of the distribution o was 18 A.24 The runs Au 10 200,0 201,1 58.0 57.0 11 0.200 6.206 43.9 with these defects distributions are labeled Au 2-5,Si 2-5, Cu 2-3,Al 2-3,and Pt 2-3 in Table I.In other sets of runs, 12 200.200 168.168 58.2 60.6 we packed between 20 and 200 defects as close to the cell 13 100,100 93,93 57.8 58.7 center as possible (Au 6-9),placed 50-200 defects uni- 14 50,50 60,60 56.4 51.8 formly in the simulation cell (runs Au 10-15)or used el- 15 50,50 52,52 56.4 55.9 evated temperatures(Au 16-20). Aud 16 50,0 52,2 45.2 41.2 For all the defect distributions,a minimum distance of 10 17 50,50 27,27 48.9 53.6 A between defects was used to ensure that defects do not 18 50,50 25,25 48.9 52.1 annihilate or form clusters independent of the cascade.Test 19a 50,50 37,37 48.9 52.3 runs showed this to be important.In all the events,the de- 20a 39,39 48.9 51.8 fects were relaxed for I ps before the collision cascade was 50,50 initiated.In the cases where only interstitials were introduced Si 0,0 84,94 54.1 into the cell,the cell size was relaxed to zero pressure using 2 50,0 133.87 52.7 45.1 a pressure control algorithm25 before the initiation of the 3 50,50 140,143 50.0 44.0 cascade. 50,50 154.160 50.0 41.7 The collision cascades were initiated by giving one of the 5 50.50 151,161 50.0 40.3 atoms in the lattice a recoil energy of 5 keV.We chose the initial recoil atom and its recoil velocity direction so that the Cu 0.0 13.13 38.6 center of the resulting collision cascade roughly overlapped 50,50 50,50 43.1 44.6 the preexisting defect structure in the simulation cell.The 3 50,50 48,48 43.1 41.9 evolution of the cascades was followed for 50 ps. 1 0,0 10,10 21.6 50,50 41,41 50.8 50.5 B.Recognition of defects 2 50,50 43,43 50.8 50.3 During some of the cascade simulations for Au and Si,a Pt 1 0,0 6,6 37.6 defect analysis routine was run every 20th time step to rec- 2 50,50 48.48 46.7 47.8 ognize and follow the motion of interstitials.An approach 3 45.45 48.8 49.6 inspired by the structure-factor analysis in Ref.26 was used 50,50 to recognize interstitials and the liquid region.The structure Event with no cascade at a fixed temperature factor P is calculated for each atom i. bDefects packed around the cell center. Defects distributed uniformly in cell. P()p Σ[,)-U2 High substrate temperature. nnb is determined from the ideal crystal structure,and is 4 for P.0-∑[U)-UP, the diamond structure(Si)and 12 for the fcc structure (Au). p()is the distribution of angles in a perfect lattice and ()=jm/nnb(nnb-1)/2 the uniform angular distribution. where ()is a list of the nb(n-1)/2 angles formed be- Before doing the sum over the angles the (j)lists are tween atom i and its nnb nearest neighbors.The number sorted by magnitude.26
For a straightforward comparison between the different materials, we used similar simulation conditions for all elements. The computational cells had a size of roughly 12031203120 Å3, and contained about 90 000 atoms, which were initially arranged in the proper crystal structure of the material. Periodic boundary conditions were used, and the temperature of the three outermost atom layers was scaled toward 0 K to dissipate excess energy. The downwards scaling factor of the temperature was restricted to be less than 0.1 during a single time step. A linkcell calculation22 and variable time step23 were also employed to speed up the simulations. The primary set of runs examined interstitial motion. For these, either 50 interstitials and 50 vacancies or only 50 interstitials were introduced in the cell. The number of defects chosen corresponds roughly to the number of Frenkel pairs which can be expected to form from a previous 20–30 keV cascade in the same region.4 The defects were introduced randomly using a Gaussian probability distribution centered at the middle of the simulation cell for vacancies and centered on a spherical cell 35 Å from the middle for interstitials. The width of the distribution s was 18 Å.24 The runs with these defects distributions are labeled Au 2–5, Si 2–5, Cu 2–3, Al 2–3, and Pt 2–3 in Table I. In other sets of runs, we packed between 20 and 200 defects as close to the cell center as possible ~Au 6–9!, placed 50–200 defects uniformly in the simulation cell ~runs Au 10–15! or used elevated temperatures ~Au 16–20!. For all the defect distributions, a minimum distance of 10 Å between defects was used to ensure that defects do not annihilate or form clusters independent of the cascade. Test runs showed this to be important. In all the events, the defects were relaxed for 1 ps before the collision cascade was initiated. In the cases where only interstitials were introduced into the cell, the cell size was relaxed to zero pressure using a pressure control algorithm25 before the initiation of the cascade. The collision cascades were initiated by giving one of the atoms in the lattice a recoil energy of 5 keV. We chose the initial recoil atom and its recoil velocity direction so that the center of the resulting collision cascade roughly overlapped the preexisting defect structure in the simulation cell. The evolution of the cascades was followed for 50 ps. B. Recognition of defects During some of the cascade simulations for Au and Si, a defect analysis routine was run every 20th time step to recognize and follow the motion of interstitials. An approach inspired by the structure-factor analysis in Ref. 26 was used to recognize interstitials and the liquid region. The structure factor Pst is calculated for each atom i, Pst~i!5 1 pu~i!S ( j @u i~j!2u i p ~j!# 2 D 1/2 , pu~i!5S ( j @u i u ~j!2u i p ~j!# 2 D 1/2 , where u i(j) is a list of the nnb(nnb21)/2 angles formed between atom i and its nnb nearest neighbors. The number nnb is determined from the ideal crystal structure, and is 4 for the diamond structure ~Si! and 12 for the fcc structure ~Au!. u i p (j) is the distribution of angles in a perfect lattice and u i u (j)5 jp/nnb(nnb21)/2 the uniform angular distribution. Before doing the sum over the angles the u i(j) lists are sorted by magnitude.26 TABLE I. Some results of the cascade events. The root-meansquare ~rms! distance of interstitials from the cell center ~Au, Cu, Al, and Pt! or the center of the interstitials ~Si! is given for the initial and final distributions of interstitials, based on the WignerSeitz cell analysis. Unless otherwise noted, the defects were initially arranged spherically around the center of the cell, as explained in the text. Element Event Defects ~int., vac.! rms distance ~Å! number Initial Final Initial Final Au 1 0, 0 4, 4 - 46.9 2 50, 0 51, 1 44.7 41.3 3 50, 50 41, 41 48.7 49.8 4 50, 50 45, 45 48.7 50.3 5a 50, 50 39, 39 48.7 54.1 Aub 6 20, 0 21, 1 15.1 20.3 7 50, 0 51, 1 23.9 21.0 8 200, 0 201, 1 35.2 30.5 9 20, 20 5, 5 25.5 45.8 Auc 10 200, 0 201, 1 58.0 57.0 11 0, 200 6, 206 - 43.9 12 200, 200 168, 168 58.2 60.6 13 100, 100 93, 93 57.8 58.7 14 50, 50 60, 60 56.4 51.8 15 50, 50 52, 52 56.4 55.9 Aud 16 50, 0 52, 2 45.2 41.2 17 50, 50 27, 27 48.9 53.6 18 50, 50 25, 25 48.9 52.1 19 a 50, 50 37, 37 48.9 52.3 20 a 50, 50 39, 39 48.9 51.8 Si 1 0, 0 84, 94 - 54.1 2 50, 0 133, 87 52.7 45.1 3 50, 50 140, 143 50.0 44.0 4 50, 50 154, 160 50.0 41.7 5 50, 50 151, 161 50.0 40.3 Cu 1 0, 0 13, 13 - 38.6 2 50, 50 50, 50 43.1 44.6 3 50, 50 48, 48 43.1 41.9 Al 1 0, 0 10, 10 - 21.6 2 50, 50 41, 41 50.8 50.5 3 50, 50 43, 43 50.8 50.3 Pt 1 0, 0 6, 6 - 37.6 2 50, 50 48, 48 46.7 47.8 3 50, 50 45, 45 48.8 49.6 a Event with no cascade at a fixed temperature. b Defects packed around the cell center. c Defects distributed uniformly in cell. d High substrate temperature. 2422 K. NORDLUND AND R. S. AVERBACK 56
56 POINT DEFECT MOVEMENT AND ANNEALING IN 2423 By analyzing the structure factor of each atom and its neighbors,combined with a kinetic energy criterion,we were Initial int. able to recognize interstitials and the liquid region during the 35 ◆ Final int. simulation. c 30 Initial vac. The motion of interstitials throughout a simulation was Final vac. followed using an interstitial database which contains histo- 25 ries of all interstitials.When an interstitial structure is recog- 6 nized during one time step,its position is compared to the 20 position of all interstitial structures recognized during previ- 15 ous time steps.In the new time step,if the new interstitial atom is closer than one lattice constant from one of the pre- wnN 10 viously recognized interstitial structures,it is interpreted to 5 be part of the same interstitial.Otherwise,a new interstitial database entry is created. 20 40 60 80 100 The position of each interstitial structure during any given time step is calculated as the average position of all atoms Distance from center(A) which are part of it.If an interstitial structure becomes part of the liquid zone or has not been seen for at least three FIG.1.Distribution of interstitials (int.)and vacancies (vac.)as interstitial analysis steps,it is discarded from the list of cur- a function of the distance from the center of the simulation cell rent interstitial structures.Thus,after the simulation the in- before and after the cascade event.The numbers are the sum over terstitial database will contain information on the positions, the three cascade events discussed in the text.The low number of lifetime,and movement of all interstitial recognized during defects is the reason for the "roughness''of the distributions. the simulation. Independently of the structure-factor analysis,we per- vacancies and interstitials as a function of the distance from formed an analysis of the population of the Wigner-Seitz cell the cell center.summed over the three cascade events. of each perfect lattice atom position for a few time steps in each simulation.Wigner-Seitz cells with more than one atom were interpreted as interstitials and empty cells as vacancies. D.Silicon The numbers of defects obtained for any given time in the For realistic MD simulation of covalently bonded materi- Wigner-Seitz and structure-factor methods were in good als like silicon,it is important to use an interatomic potential agreement with each other. which adequately describes the bonding interactions in the C.Gold material.We used the Tersoff(C)three-body interatomic potential which gives a good description of several proper- For the MD simulations of gold we employed the gold ties of Si,including the different bonding types and elastic embedded-atom method(EAM)potential used previously at moduli.27.28 The potential also gives a reasonable description this laboratory in studies of collision cascades.720 The uni- of point defect energies,although the order of the point de- versal repulsive potentials has been fitted to the potential to fect energies is not correct.2829 Since we are primarily inter- realistically describe strong collisions. ested in the possibility of defect motion,the migration ener- The interstitials introduced into the cell were given the gies of point defects are more significant in any case than (100)fcc dumbbell interstitial structure,which is the lowest- their equilibrium structure.Experimental and density- energy interstitial of fcc metals.This interstitial is known functional results suggest the vacancy migration energy is to migrate very easily.Experimentally no determination of -0.3 eV.30,31 Recent tight-binding molecular-dynamics re- the migration energy in gold has been achieved.19 The EAM sults give a migration energy of 1.4 eV for the dumbbell potential predicts a migration energy of 0.06 eV 20 which is interstitial2 (however,in real Si the interstitial migration realistic for high-temperature migration where quantum ef- may be enhanced by defect charge state processes30).For the fects are no longer important. Tersoff potential,many different values for the interstitial In the automatic interstitial recognition procedure,we and vacancy migration energies have been reported in the evaluated the structure factor P for the 12 nearest neighbors literature,probably due to the difficulty of recognizing the of each atom.We found that different values of P could be lowest-energy migration path.28,34 For the vacancy,the used to distinguish between crystalline atoms,interstitials, lowest migration activation energy calculated for an actual and liquid atoms.By comparing the Pst analysis to visual migration path is 1.6 ev,33 and for the tetrahedral interstitial inspection of moving defects,we ensured that the analysis 1.1 ev.34 Thus,it appears that the Tersoff potential clearly did not lose track of moving interstitials and fine-tuned the overestimates the vacancy migration energy,but gives a recognition criteria to recognize interstitials even in interme- fairly reasonable value for interstitials. diate configurations during migration. Frequently,the high melting point predicted by the Ter- We performed four simulations of 5 keV recoil events in soff potential is used as an argument against its suitability for which our interstitial analysis method was used:a reference collision cascade studies.However,using a new reliable ap- run in an undamaged gold crystal,two simulations with 50 proach to calculate the melting point3s we obtained a value interstitials and 50 vacancies,and one with only 50 intersti- of 2300+100 K,which is much less than the conventionally tials.Figure 1 shows the initial and final distributions of cited value of 3000 K (Refs.27 and 28)and in better agree-
By analyzing the structure factor of each atom and its neighbors, combined with a kinetic energy criterion, we were able to recognize interstitials and the liquid region during the simulation. The motion of interstitials throughout a simulation was followed using an interstitial database which contains histories of all interstitials. When an interstitial structure is recognized during one time step, its position is compared to the position of all interstitial structures recognized during previous time steps. In the new time step, if the new interstitial atom is closer than one lattice constant from one of the previously recognized interstitial structures, it is interpreted to be part of the same interstitial. Otherwise, a new interstitial database entry is created. The position of each interstitial structure during any given time step is calculated as the average position of all atoms which are part of it. If an interstitial structure becomes part of the liquid zone or has not been seen for at least three interstitial analysis steps, it is discarded from the list of current interstitial structures. Thus, after the simulation the interstitial database will contain information on the positions, lifetime, and movement of all interstitial recognized during the simulation. Independently of the structure-factor analysis, we performed an analysis of the population of the Wigner-Seitz cell of each perfect lattice atom position for a few time steps in each simulation. Wigner-Seitz cells with more than one atom were interpreted as interstitials and empty cells as vacancies. The numbers of defects obtained for any given time in the Wigner-Seitz and structure-factor methods were in good agreement with each other. C. Gold For the MD simulations of gold we employed the gold embedded-atom method ~EAM! potential used previously at this laboratory in studies of collision cascades.7,20 The universal repulsive potential18 has been fitted to the potential to realistically describe strong collisions. The interstitials introduced into the cell were given the ~100! fcc dumbbell interstitial structure, which is the lowestenergy interstitial of fcc metals.19 This interstitial is known to migrate very easily. Experimentally no determination of the migration energy in gold has been achieved.19 The EAM potential predicts a migration energy of 0.06 eV,20 which is realistic for high-temperature migration where quantum effects are no longer important. In the automatic interstitial recognition procedure, we evaluated the structure factor Pst for the 12 nearest neighbors of each atom. We found that different values of Pst could be used to distinguish between crystalline atoms, interstitials, and liquid atoms. By comparing the Pst analysis to visual inspection of moving defects, we ensured that the analysis did not lose track of moving interstitials and fine-tuned the recognition criteria to recognize interstitials even in intermediate configurations during migration. We performed four simulations of 5 keV recoil events in which our interstitial analysis method was used: a reference run in an undamaged gold crystal, two simulations with 50 interstitials and 50 vacancies, and one with only 50 interstitials. Figure 1 shows the initial and final distributions of vacancies and interstitials as a function of the distance from the cell center, summed over the three cascade events. D. Silicon For realistic MD simulation of covalently bonded materials like silicon, it is important to use an interatomic potential which adequately describes the bonding interactions in the material. We used the Tersoff (C) three-body interatomic potential which gives a good description of several properties of Si, including the different bonding types and elastic moduli.27,28 The potential also gives a reasonable description of point defect energies, although the order of the point defect energies is not correct.28,29 Since we are primarily interested in the possibility of defect motion, the migration energies of point defects are more significant in any case than their equilibrium structure. Experimental and densityfunctional results suggest the vacancy migration energy is ;0.3 eV.30,31 Recent tight-binding molecular-dynamics results give a migration energy of 1.4 eV for the dumbbell interstitial32 ~however, in real Si the interstitial migration may be enhanced by defect charge state processes30!. For the Tersoff potential, many different values for the interstitial and vacancy migration energies have been reported in the literature, probably due to the difficulty of recognizing the lowest-energy migration path.28,33,34 For the vacancy, the lowest migration activation energy calculated for an actual migration path is 1.6 eV,33 and for the tetrahedral interstitial 1.1 eV.34 Thus, it appears that the Tersoff potential clearly overestimates the vacancy migration energy, but gives a fairly reasonable value for interstitials. Frequently, the high melting point predicted by the Tersoff potential is used as an argument against its suitability for collision cascade studies. However, using a new reliable approach to calculate the melting point35 we obtained a value of 23006100 K, which is much less than the conventionally cited value of 3000 K ~Refs. 27 and 28! and in better agreeFIG. 1. Distribution of interstitials ~int.! and vacancies ~vac.! as a function of the distance from the center of the simulation cell before and after the cascade event. The numbers are the sum over the three cascade events discussed in the text. The low number of defects is the reason for the ‘‘roughness’’ of the distributions. 56 POINT DEFECT MOVEMENT AND ANNEALING IN . . . 2423
2424 K.NORDLUND AND R.S.AVERBACK 56 ment with the experimental value 1685 K.36 Details of this 10 ps (this includes all of the final interstitials)is shown. calculation will be given elsewhere.37 Most of the interstitials seen in the figure actually existed A short-range repulsive part of the potential was deter- throughout the simulation.The open markers indicate the mined from density-functional theory calculations.38 It was initial positions and the solid ones the final positions.The smoothly fitted to the Tersoff potential using a Fermi func- defects move very little during the simulation,even though tion F(r)=(1+eb)-1 with the values 12 A-1 and they clearly have been exposed to the pressure and heat 1.6 A for b and r,respectively.40 These values provide a waves from the cascade.Noteworthy is that the interstitials smooth fit between the two potentials both between two at- move slightly inwards on average. oms in a dimer and two nearby atoms in bulk silicon. These conclusions were verified by a calculation of the The point defect recognition procedure in silicon was total movement of interstitials with respect to the center of similar to the one used in gold.The Pst structure factor was mass of the liquid zone throughout the simulation.The inter- evaluated for the four nearest neighbors of each atom,and stitials which remained after the cascade had on average used to recognize both interstitials and vacancies.Liquid at- moved inwards 1.0,1.5,and 1.6 A in the three cascade oms were recognized using a combined P and kinetic en- events containing preexisting defects.There was,however,a ergy criterion. considerable spread in the distribution of the movement:The We carried out 5 simulations for silicon:one reference maximum inwards movement seen in the three events was 10 run with no initial defects,one run with 50 initial intersti- A and the maximum outward movement 14 A.On average, tials.and three with 50 initial interstitials and 50 initial va- the interstitials outside the liquid region performed about cancies. three lattice site jumps.From Fig.I we see that the overall effect of the interstitial movement on the defect distribution E.Copper,aluminum,and platinum is quite small. For copper,aluminum,and platinum we employed EAM In the events with 50 initial interstitials and vacancies. potentials4142 onto which the universal repulsive potentiall about 40-45 interstitial structures remained after the cascade had been fitted to realistically describe strong collisions2 event.A few of the interstitials produced in our simulations We carried out 3 simulations for these metals.I without formed small clusters,but the majority remained as single initial defects and 2 with 50 initial interstitials and vacancies. dumbbell interstitials.Some of the initially existing intersti- The initial defect distributions were the same as for the gold tials have recombined with vacancies during the event.Be- simulations.The analysis of defects was performed with the cause this is more likely to occur close to the cell center,the Wigner-Seitz method for about 20 selected time steps. root-mean-square distance of interstitials from the cell center may actually be larger for the final than the initial defect distribution (cf.Table D). IIL.RESULTS AND DISCUSSION Vacancies were identified geometrically using Wigner- A.Gold Seitz cells;otherwise,the vacancy analysis employed the same principles as those described above for interstitials.The The gold results will be presented in four subsections.The vacancy movement results are illustrated in Fig.4.We see first deals with defects outside the liquid core of a cascade, that the vacancies outside the liquid region do not move at all the second with defects in the core,the third with uniform defect distributions,and the last one with high-temperature or move very little.The average number of lattice site jumps per vacancy was 0.2.This is readily understood in terms of runs. their large migration energy,0. Inspection of animations of the collision cascades showed 1.Defects outside the liquid core that most of the interstitial movement occurred when the Results of one simulation with 50 initial vacancies and liquid zone was shrinking in size or had already vanished.To interstitials are illustrated in Figs.2 and 3.Figure 2 shows understand the reason for the migration,we calculated tem- snapshots of the positions of all atoms in the system pro- perature and pressure distributions in the cell,which are jected onto the xy plane at 0,1.5,and 50 ps.At 0 ps,the shown in Figs.5 and 6.High temperatures are present in the initial interstitials and lattice relaxation caused by them are cell for a relatively long period of time after the collision visible among the otherwise undisturbed fcc lattice atoms.At cascade,suggesting this might cause thermal migration of 1.5 ps,we see the liquid cascade center and the pressure the interstitials.The pressure decreases towards the center of wave emanating from it.Some of the interstitials can be the cell before and after the cascade due to the presence of discerned even through the pressure wave,showing that they vacancies near the cell center and the presence of interstitials are not destroyed by it.At 50 ps,spike effects are no longer in the region near the borders(see Fig.6).We believe this visible,and the lattice has relaxed back close to its initial contributes to the inward motion of the interstitials. state.Comparison of the 0 and 50 ps figures shows that most In the event which had 50 initial interstitials and no va- interstitials are positioned at roughly the same sites before cancies,a similar pressure gradient was seen,which indicates and after the cascade.In examining the figures,one should the interstitials are the predominant cause of it.In the cas- keep in mind that when a dumbbell interstitial is oriented cade event with no initial defects,in which only four inter- along the =axis,it will not be visible when projected on the stitials existed after the cascade,no reduction of the pressure xy plane. towards the center occurred. Our analysis of interstitial movement enabled us to relate To test whether the defect movement could be explained initial and final interstitials to each other.In Fig.3,the as a thermal migration process,we simulated the same initial movement of all interstitials which have existed for at least defect distribution as in the cascade events but with no cas-
ment with the experimental value 1685 K.36 Details of this calculation will be given elsewhere.37 A short-range repulsive part of the potential was determined from density-functional theory calculations.38,39 It was smoothly fitted to the Tersoff potential using a Fermi function F(r)5(11e2b f(r2rf) )21 with the values 12 Å 21 and 1.6 Å for bf and rf , respectively.40 These values provide a smooth fit between the two potentials both between two atoms in a dimer and two nearby atoms in bulk silicon. The point defect recognition procedure in silicon was similar to the one used in gold. The Pst structure factor was evaluated for the four nearest neighbors of each atom, and used to recognize both interstitials and vacancies. Liquid atoms were recognized using a combined Pst and kinetic energy criterion. We carried out 5 simulations for silicon; one reference run with no initial defects, one run with 50 initial interstitials, and three with 50 initial interstitials and 50 initial vacancies. E. Copper, aluminum, and platinum For copper, aluminum, and platinum we employed EAM potentials41,42 onto which the universal repulsive potential18 had been fitted to realistically describe strong collisions.26 We carried out 3 simulations for these metals, 1 without initial defects and 2 with 50 initial interstitials and vacancies. The initial defect distributions were the same as for the gold simulations. The analysis of defects was performed with the Wigner-Seitz method for about 20 selected time steps. III. RESULTS AND DISCUSSION A. Gold The gold results will be presented in four subsections. The first deals with defects outside the liquid core of a cascade, the second with defects in the core, the third with uniform defect distributions, and the last one with high-temperature runs. 1. Defects outside the liquid core Results of one simulation with 50 initial vacancies and interstitials are illustrated in Figs. 2 and 3. Figure 2 shows snapshots of the positions of all atoms in the system projected onto the xy plane at 0, 1.5, and 50 ps. At 0 ps, the initial interstitials and lattice relaxation caused by them are visible among the otherwise undisturbed fcc lattice atoms. At 1.5 ps, we see the liquid cascade center and the pressure wave emanating from it. Some of the interstitials can be discerned even through the pressure wave, showing that they are not destroyed by it. At 50 ps, spike effects are no longer visible, and the lattice has relaxed back close to its initial state. Comparison of the 0 and 50 ps figures shows that most interstitials are positioned at roughly the same sites before and after the cascade. In examining the figures, one should keep in mind that when a dumbbell interstitial is oriented along the z axis, it will not be visible when projected on the xy plane. Our analysis of interstitial movement enabled us to relate initial and final interstitials to each other. In Fig. 3, the movement of all interstitials which have existed for at least 10 ps ~this includes all of the final interstitials! is shown. Most of the interstitials seen in the figure actually existed throughout the simulation. The open markers indicate the initial positions and the solid ones the final positions. The defects move very little during the simulation, even though they clearly have been exposed to the pressure and heat waves from the cascade. Noteworthy is that the interstitials move slightly inwards on average. These conclusions were verified by a calculation of the total movement of interstitials with respect to the center of mass of the liquid zone throughout the simulation. The interstitials which remained after the cascade had on average moved inwards 1.0, 1.5, and 1.6 Å in the three cascade events containing preexisting defects. There was, however, a considerable spread in the distribution of the movement: The maximum inwards movement seen in the three events was 10 Å and the maximum outward movement 14 Å. On average, the interstitials outside the liquid region performed about three lattice site jumps. From Fig. 1 we see that the overall effect of the interstitial movement on the defect distribution is quite small. In the events with 50 initial interstitials and vacancies, about 40–45 interstitial structures remained after the cascade event. A few of the interstitials produced in our simulations formed small clusters, but the majority remained as single dumbbell interstitials. Some of the initially existing interstitials have recombined with vacancies during the event. Because this is more likely to occur close to the cell center, the root-mean-square distance of interstitials from the cell center may actually be larger for the final than the initial defect distribution ~cf. Table I!. Vacancies were identified geometrically using WignerSeitz cells; otherwise, the vacancy analysis employed the same principles as those described above for interstitials. The vacancy movement results are illustrated in Fig. 4. We see that the vacancies outside the liquid region do not move at all or move very little. The average number of lattice site jumps per vacancy was 0.2. This is readily understood in terms of their large migration energy, 0.8 eV.19 Inspection of animations of the collision cascades showed that most of the interstitial movement occurred when the liquid zone was shrinking in size or had already vanished. To understand the reason for the migration, we calculated temperature and pressure distributions in the cell, which are shown in Figs. 5 and 6. High temperatures are present in the cell for a relatively long period of time after the collision cascade, suggesting this might cause thermal migration of the interstitials. The pressure decreases towards the center of the cell before and after the cascade due to the presence of vacancies near the cell center and the presence of interstitials in the region near the borders ~see Fig. 6!. We believe this contributes to the inward motion of the interstitials. In the event which had 50 initial interstitials and no vacancies, a similar pressure gradient was seen, which indicates the interstitials are the predominant cause of it. In the cascade event with no initial defects, in which only four interstitials existed after the cascade, no reduction of the pressure towards the center occurred. To test whether the defect movement could be explained as a thermal migration process, we simulated the same initial defect distribution as in the cascade events but with no cas- 2424 K. NORDLUND AND R. S. AVERBACK 56
56 POINT DEFECT MOVEMENT AND ANNEALING IN .. 2425 0ps 1.5ps (a) (b) 50 ps (c) FIG.2.Snapshots of a 5 keV Au cascade.All atom position in the (115 A)3 simulation cell are included in the figure,projected on the xy plane.The interstitials and lattice relaxation around them breaks the otherwise regular fec structure in the 0 and 50 ps figures.The liquid zone in the center of the cell and the pressure wave surrounding it can be clearly seen in the 1.5 ps figure.Comparison of the nonregular features in the 0 and 50 ps figures shows that most interstitials have not moved much from their initial positions cade initiated in the cell.In case the movement of the inter- Tcase(t)e(-006eV)kgT(dt stitials is determined solely by the pressure decreasing in- Tm= Te(-0.06eVIkBT(dt (1)） wards and the temperature of the cell,interstitial movement similar to the cascade events should be seen From the average of T in the three interstitial cascade runs A 50 ps simulation was carried out with a cell temperature we obtained T =205 K. T equivalent to the temperature in the cascade simulations. We carried out a 50 ps run for a cell which had 50 inter- We assumed the thermal defect migration probability is pro- stitials and 50 vacancies distributed as in the cascade runs. portional to exp(-E/kgT),where E=0.06 eV.20 Thus the but now heating it to 205 K from the borders.We obtained average migration temperature Tm of interstitials is obtained similar results to the cascade runs.Due to the distribution of from the temperature of the entire simulation cell as a func- defects in the cell,the pressure decreased inwards.The in- tion of time Tease(1)in the cascade runs,weighted by the terstitials existing after the simulation had,on average, migration exponential. moved 0.6 A inwards.As in the case of the interstitials in the
cade initiated in the cell. In case the movement of the interstitials is determined solely by the pressure decreasing inwards and the temperature of the cell, interstitial movement similar to the cascade events should be seen. A 50 ps simulation was carried out with a cell temperature Tm equivalent to the temperature in the cascade simulations. We assumed the thermal defect migration probability is proportional to exp(2EA /kBT), where EA50.06 eV.20 Thus the average migration temperature Tm of interstitials is obtained from the temperature of the entire simulation cell as a function of time Tcasc(t) in the cascade runs, weighted by the migration exponential, Tm5*Tcasc~t!e~20.06 eV!/kBT~t! dt *e~20.06 eV/kBT~t! dt . ~1! From the average of Tm in the three interstitial cascade runs we obtained Tm5205 K. We carried out a 50 ps run for a cell which had 50 interstitials and 50 vacancies distributed as in the cascade runs, but now heating it to 205 K from the borders. We obtained similar results to the cascade runs. Due to the distribution of defects in the cell, the pressure decreased inwards. The interstitials existing after the simulation had, on average, moved 0.6 Å inwards. As in the case of the interstitials in the FIG. 2. Snapshots of a 5 keV Au cascade. All atom position in the ~115 Å! 3 simulation cell are included in the figure, projected on the xy plane. The interstitials and lattice relaxation around them breaks the otherwise regular fcc structure in the 0 and 50 ps figures. The liquid zone in the center of the cell and the pressure wave surrounding it can be clearly seen in the 1.5 ps figure. Comparison of the nonregular features in the 0 and 50 ps figures shows that most interstitials have not moved much from their initial positions. 56 POINT DEFECT MOVEMENT AND ANNEALING IN . . . 2425