A phenomenological electronic stopping power model for molecular dynamics and Monte Carlo simulation of ion implantation into silicon David Cai1, Niels Grønbech-Jensen1, Charles M. Snell2, and Keith M. Beardmore1 1 Theoretical Division, and 2Applied Theoretical and Computational Physics Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545 9 It is crucial to have a good phenomenological model of electronic stopping power for modeling 9 thephysicsofionimplantationintocrystallinesilicon. InthespiritoftheBrandt-Kitagawaeffective 9 chargetheory,wedevelopamodelforelectronicstoppingpowerforanion,which canbefactorized 1 into(i)agloballyaveragedeffectivechargetakingintoaccounteffectsofcloseanddistantcollisions n bytargetelectronswiththeion,and(ii)alocalchargedensitydependentelectronicstoppingpower a foraproton. ThisphenomenologicalmodelisimplementedintobothmoleculardynamicsandMonte J Carlo simulations. There is only one free parameter in the model, namely, the one electron radius 8 rs◦ for unbound electrons. By fine tuning this parameter, it is shown that the model can work 2 successfully for both boron and arsenic implants. We report that the results of the dopant profile simulation for both species are in excellent agreement with the experimental profiles measured by ] secondary-ion mass spectrometry (SIMS) over a wide range of energies and with different incident h directions. Wepointoutthatthemodelhaswideapplicability,foritcapturesthecorrectphysicsof p electronic stoppingin ion implantation. This modelalso providesagood physically-based damping - p mechanismformoleculardynamicssimulationsintheelectronicstoppingpowerregime,asevidenced m bythestriking agreement of dopant profilescalculated in our molecular dynamicssimulations with o the SIMS data. c . s c i s y I. INTRODUCTION h p [ Ion implantation in semiconductors is an important technology in integrated circuit device fabrication1. A reliable description of as-implanted profiles and the resulting damage is needed for technologicaldevelopment, such as device 1 design and modeling, as well as process optimization and control in the fabrication environment. For semiconductor v 6 deviceswhosephysicaldimensionsareoforderofsubmicronsorsmaller,lowimplantenergiesandreductionofthermal 5 processing are necessary, resulting in more prominent channeling effects in the as-implanted profiles and less post- 0 implant diffusion. At these physical dimensions, it is essential to obtain the two- or three-dimensional details of the 1 ever shallower and more compact dopant and damage profiles for post-implant diffusion simulations. 0 Studyoftheenergylossofchanneledparticleshasalonghistory2,forthechannelingfeaturescanbeusedtoelucidate 9 the energy-loss mechanisms. Earlier analytical treatments of the implant profiles based on moment distributions, 9 derived from the Lindhard-Scharff-Schiott theory (LSS)3, preclude channeling because of the amorphous nature of / s the targets assumed in the studies. Later, it was realized that, because of the channeling effect, electronic stopping c i power plays a much more significant role in ion implantation into crystalline solids than otherwise would be deduced s from the application4 of the LSS theory to amorphous materials. It is especially true for heavy ion implants at low y energies, such as arsenic ions in the energy range below 700keV5,6. For implantation into silicon, most Monte Carlo h p (MC) models are only concerned with boron implants, and have not modeled arsenic implants accurately with an : electronicstoppingpowermodelconsistentwiththatusedforboron7,8. Aswillbeshownbelow,thephenomenological v model we developed for electronic stopping power can be implemented into a Monte Carlo simulation program for i X both boron and arsenic implants in different channels with equal success over a wide range of implant energies. r InadditiontoMonteCarlosimulationswiththebinarycollisionapproximation(BCA)9,moleculardynamics(MD) a incorporating multiple-interactions via many-body potentials can also be used to simulate the behavior of energetic ionsinamorphousorcrystallinesilicon. Thismethodisespeciallyapplicableatlowenergies,forwhichmany-body,and multiple interactions are increasingly important10. Although it is well known that the BCA is valid for high incident energies (∼ 0.1keV up to ∼MeV, the upper limit is set by relativistic effects), in a cascade, especially initiated by a relatively low energy ion, the energy of the ions will decrease and eventually reach the lower validity limit of the BCA at which many-body effects become important9,10. For crsytals of high symmetry, the BCA can be modified to account for simultaneous collisions in channels9,11, and MD results can provide good insight into how to successfully modify the BCAinthissituation. Moreover,MDresultscanbecomparedto BCAMonteCarlosimulationsandused to establish the low energy limits of the binary collision approximation. 1 An extremely important issue in deploying molecular dynamics to model collision processes in covalent and ionic solidsishowtoincorporateenergytransfermechanismsbetweenelectronsandions12. Agooddescriptionofdynamical processes in energetic collisions, such as initial displacement damage, relaxation processes, and the cooling phase as the energy dissipates into the ambient medium, requires a theoretical framework that encompasses all interactions between ion-ion, ion-electron, and their interaction with the thermal surroundings. Especially, it should capture the nonequilibrium thermodynamic nature of these physical processes involving a wide range of energy scales, from a low energy electron-phonon interaction regime to a high energy radiation damage regime.13–17. Traditional MD simulations can capture the thermal behavior of an insulator. Since they do not take into account coupling between the phonons and the conduction electron system, obviously, these simulations underestimate the heat-transfer rate for noninsulating materials. In addition to lattice thermal conductivity, the issue ofthe conductivity due to electrons must be addressed. Furthermore, a correct description of the electronic stopping power should be incorporated into MD simulations for high energy implantation. For example, in sputtering processes by particle bombardment, examination of MD simulations with and without inelastic electronic energy loss has established that, independent of the ion’s mass or energy, the inelastic electronic energy losses by target atoms within the collision cascade have greaterinfluence onthe ejected atom yield than the ion’s electronic losses18. This is in contrastto the belief that the electroniclossmechanismisimportantonlyforcascadesinitiatedbylightionsorbyheavyionsathighbombardment energies18. Although a convincing experimental verification of the electronic effects in sputtering is still lacking, the effects should be relevant to defect production rates, defect mobility and annealing, etc.19,20. Also as shown in Ref. 21, traditional MD simulations produce extremely long channeling tails due to the absence of electronic stopping. In order to incorporate the ion-electroninteraction into molecular dynamics simulations, a simple scheme was proposed by adding a phenomenological term, which describes the inelastic electronic stopping in the high energy radiation damage regime, while also capturing the thermal conductivity by coupling low energy ions to a thermal reservoir22. The empirical expression used in Ref. 22 for the strength of the ion-electron coupling is a function of the local electronic density. At the low charge density limit, a density functional result was reproduced23, and at the high chargedensity limit, the linear responseresults were captured. Inthe same spirit, we developa stochastic MD model incorporating the electronic stopping power as a damping mechanism. Our model is based on an effective charge theory24 with the electronic stopping power factorized into two parts. One is the effective charge of the incident ion, which is a globally averaged quantity determined by the average unbound electron density in the medium. The otherfactoristhe electronicstoppingpowerforaproton,forwhichthesamelocaldensityfunctionalresultsareused. Naturally, our damping mechanism incorporates both regimes, i.e., the electronic stopping regime and the electron- phonon interaction regime, into our molecular dynamics simulation, because the inelastic loss for a proton exhibits a similar density dependence as prescribed in Ref. 22, with additional modifications due to the velocity dependence of the effective charge. In the present work, however, we emphasize mainly the electron stopping power in the high energy regime (∼ keV to ∼100keV), i.e., the electrons behave as an energy sink. The validity of the model for the electronic heat conduction regime will be discussed elsewhere. In the following, for boron and arsenic implants into single-crystalsiliconinboththechannelingandoff-axisdirections,wewillshowthataclassicalMDwiththephysically- baseddamping mechanismcangeneratedopantprofilesinexcellentagreementwithexperimentally-measuredprofiles obtained by secondary-ionmass spectrometry (SIMS). As discussed above, the phenomenological model we have developed for electronic stopping power is successfully implemented into both BCA Monte Carlo programs and MD simulations. Wide applicability requires that a model be validfor different implant species overa wide rangeof energies. We emphasize thatthis electronic stopping model is accurate both for boron and arsenic implants, thus providing a crucial test of the generality and validity of the model in capturing the correct physics of electronic stopping. The paper is organizedas follows. In Sec. II, we present the phenomenologicalmodel for electronic stopping power in detail. Atomic units e = h¯ = m = 1 are used throughout the paper unless otherwise specified. In Sec. III, e we briefly discuss different electronic stopping models implemented on the versatile BCA Monte Carlo simulation platform,MARLOWE9,10. Thenthe resultsoftheBCAMonteCarlosimulationsonarare-eventalgorithmenhanced UT-MARLOWE platform25 with our electronic stopping model are summarized. In Sec. IV, the results of the MD with the inelastic electronic energy loss are presented. In Sec. V, we make closing remarks and point out directions for future studies. II. THE MODEL According to the Brandt-Kitagawa (BK) theory24, the electronic stopping power of an ion can be factorized into two components based on an effective charge scaling argument. One is the effective charge of the ion (if not fully ionized), Z∗, which is in general a function of ion velocity v and the charge density of the target ρ, or equivalently, 1 2 the one electronradius r =[3/(4πρ(x))]1/3; the other is the electronic stopping power for a proton, S (v,r ). In the s p s local density approximation, therefore, the total inelastic energy loss ∆E of an ion of constant velocity v is e ∆E = [Z∗(v,r )]2S (v,r )dx, (1) e 1 s p s Z where the integral is along the ion path. Since the effective charge is a continuous function of electronic density, mathematically, it is always possible to find a mean value, r◦ of r , such that Eq. (1) can be rewritten as s s ∆E =[Z∗(v,r◦)]2 S (v,r )dx. (2) e 1 s p s Z Iftheeffectivechargeisaslowlyvaryingfunctionofspace,physically,thismeansthatr◦ describesanaveragenumber s of unbound electrons in the sea and thus can be assumed to determine the Fermi surface. Therefore, we have the relation between the Fermi velocity and r◦: s 1 v = , (3) F αr◦ s where α = [4/(9π)]1/3. We note that this r◦ will be the only tunable parameter in our electronic stopping power s model. Nextweturntoasimplestatisticalmodelforthis partiallyionized,movingprojectile. ForanionwithN =Z −Q 1 bound electrons, where Q is the charge number of the ion of atomic number Z , a radially symmetric charge density 1 N r ρ = exp − (4) e 4πΛ2r Λ (cid:16) (cid:17) isusedintheBKtheory. HereΛistheionsizeparameter,afunctionofthefractionalionization,q =(Z −N)/Z . The 1 1 total energy of the electrons comes from the sum of the kinetic energy estimated by the local density approximation, the electron-electron interaction in the Hartree approximation weighted by a variational parameter λ to account for correlation, and the Coulomb energy of the electrons in the electric field of the nucleus. A variational approach minimizing the total energy leads to the following dependence of the ion size on the ionization fraction q: 2a (1−q)2/3 0 Λ= , (5) Z1/3[1−(1−q)/7] 1 where a =0.24005. In the BK theory, the generalized Lindhard theory of the electronic stopping in a homogeneous 0 electron gas with an electron density n = 3/(4πr¯3) is used. The total electronic stopping is estimated from the sum s of the energy loss in soft, distant collisions, i.e., small momentum transfers with target electrons seeing a charge qZ , and the energy loss to the target electrons experiencing increased nuclear interaction in hard, close collisions 1 corresponding to large momentum transfers. As extensively discussed in the literature (see, e.g., Ref. 24,26,27, and references therein), it is assumed that the charge state of a proton in a solid is unity. Given an ionization fraction q and using the scaling argument for the ratio of ion stopping to the proton stopping at the same velocity, the BK theory produces a simple expression for the fractional effective charge of an ion24,27 2 4Λ γ(r¯ )=q+C(r¯ )(1−q)ln 1+ , (6) s s r¯ " (cid:18) s (cid:19) # where C(r¯ ) is weakly dependent on the target and has a numerical value of about 1/2. We will set C = 0.5 below. s Then, the effective charge is Z∗ =Z γ(r¯ ). (7) 1 1 s For our model, using the procedure (2) outlined above,this dependence of r¯ is identified with the dependence of the s meanvaluer◦. Therefore,theeffectivechargeZ∗ hasanonlocal,i.e.,spatiallyindependent,characteranddependson s 1 the Fermi surface. In the above discussion, as can be seen, q is a parameter which is not fixed by the BK theory. For obtaining this ionization fraction, there are velocity and energy criteria originally proposed by Bohr28 and Lamb29, respectively. Kitagawa also used a statistical argument to justify scaling analyses in terms of the scaling parameter v /(v Z2/3)30. Recently, the issue of which stripping criterion can give rise to a better physical understanding has 1 B 1 beenraised31. However,inlightofthelargeamountofexperimentaldataemployedinRef.27toextractanionization 3 scaling consistent with the Brandt-Kitagawa theory, we will use this empirically verified scaling in our model. As summarized in Ref. 27, a new criterion in the BK approach is proposed24,26, i.e., a relative velocity criterion, which assumes that the electrons of the ion which have an orbital velocity lower than the relative velocity between the ion andthe electronsinthe medium arestripped off. The relativevelocity v is obtainedby averagingoverthe difference r betweentheionvelocityv andtheelectronvelocityv undertheassumptionthattheconductionelectronsareafree 1 e electron gas in the ground state, therefore, whose velocity distribution is isotropic. Performing a further averagingof v over the Fermi sphere leads to26 e v2 v =v 1+ F for v ≥v , (8) r 1 5v2 1 F (cid:18) 1(cid:19) 3v 2v2 v4 v = F 1+ 1 − 1 for v <v . (9) r 4 3v2 15v4 1 F (cid:18) F F(cid:19) For the ionization scaling, a form of the Northcliffe type32 is then assumed for the scaling variable, i.e., the reduced relative velocity: v r y = , (10) r v Z2/3 B 1 where v is the Bohr velocity and v = 1 in our units. The extensive experimental data for ions 3 ≤ Z ≤ 92 are B B 1 used in Ref. 27 to determine q =1−exp[−0.95(y −0.07)]. (11) r In Ref. 27, an ionization scaling fit with even tighter bunching of the experimental data along the fit is presented. However, this approach entails a much more involved computational procedure27. The accuracy level of Eq. (11) is adequate for our present purposes. Inourmodel,theelectronicstoppingpowerforaprotonisderivedfromanonlineardensity-functionalformalism23. In the linear response theory, the energy loss per unit path length of a proton moving at velocity v in the electron gas is obtained by Ritchie33 dE 2v π 1 = ln 1+ − , (12) dx 3π αr 1+αr /π (cid:18) (cid:19)R (cid:20) (cid:18) s(cid:19) s (cid:21) usinganapproximationtothefullrandom-phaseapproximationdielectricfunction,whichamountstotheexponential screeningpotential aroundthe ioninduced by density fluctuations of the electrons. The nonlinear,density-functional calculation based on the formalism of Hohenberg and Kohn, and Kohn and Sham34,35 has been performed23,36,37 to obtainthe chargedensityandscatteringphaseshiftsforthe conductionbandasafunctionofenergyself-consistently. The final stopping power for a proton is obtained via ∞ dE 3v = (l+1)sin2[δ (E )−δ (E )], (13) dx k r3 l F l+1 F F s l=0 X where δ (E ) is the phase shift at the Fermi energy for scattering of an electron of angular momentum l and k is l F F theFermimomentum38. AsshowninRef.27,39,40,comparisonwithexpermentaldatademonstratesthatthedensity functionaltreatmentprovidesanimprovementoverthe linearresponse(dielectric)result41, whichunderestimatesthe stopping powers. In our implementation, the result of the nonlinear, density functional formalism for the electronic stopping power for a proton is used, which can be expressed as dE S (v,r )=− G(r ), (14) p s s dx (cid:18) (cid:19)R where, for computational convenience, the correction factor G(r ) takes the form s G(r )=1.00+0.717r −0.125r2−0.0124r3+0.00212r4 (15) s s s s s for r < 6. We note that a different correction factor was used in Refs. 42,43, which does not have the following s desired behavior for r ≪ 1. Since the density functional result converges to the Ritchie formula as r decreases s s 4 towards values sufficiently small compared to unity23, this requires that the correction factor smoothly tend to unity as r →0. Obviously, the above G(r ) possesses the correct convergence property. s s The last ingredient needed for our model is the charge distribution ρ(x) for silicon atoms in the crystal. We use the solid-state Hartree-Fock atomic charge distribution27, which is spherically symmetric due to the muffin-tin construction. In this approximation, there is about one electron charge unit (0.798 electrons for Si) left outside the muffin-tin. Thissmallamountofchargecanbeeitherdistributedinthevolumebetweenthesphericalatoms,resulting in an interstitialbackgroundchargedensity 0.119e/˚A3, ordistributed betweenthe maximalcollisiondistance used in Monte Carlo simulations and the muffin-tin radius (see details below). III. BCA MONTE CARLO SIMULATION RESULTS First, in comparison with other electronic stopping models used in Monte Carlo simulations based on the MAR- LOWE platform, we stress that in our model the effective charge is a nonlocal quantity, neither explicitly dependent on the impact parameter nor on the charge distribution, and the stopping power for a proton depends on the local charge density of the solid. A purely nonlocal version of the BK theory was implemented into MARLOWE42, in which both the effective charge and the stopping power for a proton depend on a single nonlocal parameter, namely, the averagedone electron radius. Its results demonstrated that energy loss for well-channeled ions in the keV region has high sensitivity to the one-electron radius in the channel. It was pointed out that a correct density distribution is needed to account for the electronic stopping in the channel8,42. Later, a purely local version of the BK theory was developed to take into account the charge distribution of the electrons43,44. Comparison with other electronic stoppingmodels,suchasLindhardandScharff4,Firsov45,andtheabovenonlocalimplementation42,showedamarked improvementinmodelingelectronicstoppinginthechannel7,43,44. Goodagreementbetweensimulateddopantprofiles and the SIMS profiles for boron implants into h100i single crystal silicon was obtained. However, this purely local implementation of the BK theory did not successfully model the electronic stopping for the boron implants into the h110i axial channel and arsenic implants as noted in Ref. 7. In the present work, UT-MARLOWE25 was selected as the platform for our electronic stopping model implemen- tation. UT-MARLOWE is an extension of the MARLOWE code for simulating the behavior of energetic ions in crystalline materials9,10. It has been enhanced with: (i) atomic pair-specific interatomic potentials for B-Si, B-O, As-Si, As-O for nuclear stoppings27, (ii) variance reduction algorithm implemented for rare events, (iii) important implant parameters accounted for, e.g., tilt and rotation angles, the thickness of the native oxide layers, beam diver- gence, and wafer temperature, etc. In our simulations, we have turned off certain options, such as the cumulative damage model in the UT-MARLOWE code, which is a phenomenological model to estimate defect production and recombination rates. Individual ion trajectories were simulated under the BCA and the overlapping of the damage causedbydifferent individualcascadeswasneglected. Inorderto test the electronic stoppingmodel we alsousedlow dose (1013/cm2) implants so that cummulative damage effects do not significantly complicate dopant profiles7. Also, for the simulation results we report below, 16˚A native oxide surface layer, 300K wafer temperature were used. The maximumdistanceforsearchingacollisionpartneris0.35latticeconstant,thedefaultvalueintheUT-MARLOWE5. The excesschargeoutside the muffin-tinsis distributedinthe spacebetweenthis maximumcollisiondistanceandthe muffin-tin radius. In the simulation, the electronic stopping power is evaluated continuously along the path the ion traverses through regions of varying charge density, i.e., the energy loss is given by ∆E = [Z γ(v ,r◦)]2S (v ,r (x))dx. (16) e 1 1 s p 1 s Zion path In the simulations, the free parameter r◦ was adjusted to yield the best results in overall comparison with the s experimental data. The value r◦ = 1.109˚A was used for both boron and arsenic ions for all energies and incident s directions. This value is physically reasonable for silicon. Note that the unbound electronic density in silicon with onlyvalenceelectronstakenintoaccountwillgiverisetoavalueof1.061˚Aforr . Thefactthatourr◦ valueisgreater s s than 1.061˚A indicates that not all valence electrons participate in stopping the ion as unbound electrons. We display the Monte Carlo dopant profile simulation results as follows. We note in passing that the lower and upper limits of energy used in our simulations are determined by the energy range of the SIMS data available to us. InFigs.1,2and3,weshowborondopantprofilesfortheenergies15keV,35keV,and80keValongh100i,h110i,and the off-axis direction with tilt = 7◦ and rotation = 30◦, respectively. It can be seen that the overall agreement with the SIMS data is excellent. In Fig. 1, the simulations show a good fit for the cutoff range. In the high energy regime, the simulated distribution shows a slightly peaked structure. This can be attributed to a strong channeling due to insufficient scatterings in the implants. We have noticed that by increasing,e.g., the native oxide layer thickness, the peak can be reduced. For the h110i channeling case, the distribution indicates a possibility that the total electronic 5 stoppingpoweralongthechannelisalittletoostrongatthehighenergyend. However,itshouldbekeptinmindthat for this channel, the UT-MARLOWE model becomes sensitive to the multiple collision parameter which is employed as an approximate numerical correction for the effect of multiple overlapping nuclear encounters. It is not clear how to separate the contributions from these two different sources. For comparison, in Fig. 4, we also display a low energy (5keV) implant case46. Again the agreement for both the channeling and off-axis directions is striking (thin lines without symbols). In order to illustrate the importance of electronic stopping power at this low energy for boron implants, we have used an artificially reduced electronic stopping power, i.e, multiplying ∆E in Eq. (16) by a factor of 1/10 in the simulation, to genearte dopant profiles in e thechannelingandoff-axisdirections. FromFig.4,evidently,itcanbeconcludedthat,forboronimplantseveninthis lowenergyregime,electronicstoppingpowerhassignificantinfluence onthe channeledtailofthe dopantdistribution and on the cutoff range for both channeling and off-axis directions. Figs.5 and 6 show arsenic dopantprofiles for energiesranging from15keVto 180keValong h100i, andthe off-axis direction with tilt =8◦ and rotation=30◦, respectively. It can readily be concluded that our electronic model works successfully with arsenic as well as boron implants into crystalline silicon. For comparison, a case of arsenic implant intoamorphoussiliconisalsoshowninFig.7. Theimplantenergyis180keV.The effectofelectronicstopping,which shows clearly in the long sloping channeling tail in the crystalline counterpart (see Fig. 5), is less prominent for the amorphous case (Fig. 7). To examine the role that electronic stopping power has on arsenic implants in the low energy regime (∼5keV), we again simulated arsenic dopant profiles with the artificially reduced electronic stopping power. For these low energy implants, the oxide layer thickness 3˚A were used in the BCA simulations on account of the fact that the wafers used fortheseimplantsweretreatedbydiluteHFetchfor30seconds,thenimplantedwithin2hourstopreventnativeoxide regrowth47. In Fig. 8, we show that our electronic stopping power model is successful in both the h100i channeling and off-axis directions (thin lines without symbols). Clearly, the artificial reduction of electronic stopping power leads to incorrectdopantdistributions and cutoff rangesfor both the channeling and off-axisdirections, althoughthe deviationsindicatealesssignificantcontributionfromelectronicstoppingforarsenicimplantsthanforboronimplants at the energy5keV.However,the deviationin the cutoff rangedue to the electronic stopping powerreduction for the channeling case is still significant. Obviously, this reinforces the conclusion that, for channeling implants even at low energies, electronic stopping is not negligible. In summary, the above results demonstrate clearly that our electronic stopping power indeed captures the correct physics of the electronic stopping for ion implants into silicon over a wide range of implant energies. IV. MD SIMULATION RESULTS We have also used classical molecular dynamics simulation to study the electronic stopping power as one of the damping mechanisms in the high energy regime, as discussed above. Here we demonstrate that experimental data, suchasSIMS,canbeusedtotestthevalidityofthisphysically-baseddampingmodel. Theinteractionbetweensilicon atoms are modeled by Tersoff’s empirical potential48: 1 E = f(r )[V (r )−b V (r )], (17) ij R ij ij A ij 2 i6=j X wheref(r )isacutofffunctionthatrestrictsinteractionstonearestneighbors,V (r )andV (r )arepairterms,and ij R ij A ij b is a many-body function that canbe regardedas an effective Pauling bond order. We have modified the repulsive ij partoftheTersoffpotentialbysplinningtotheZBLuniversalpotentialatclose-range.27. TheZBLuniversalpotential isalsousedto modelthe ion-siliconinteractions. InourfullMDsimulationsforthe lowdoseimplantation,the lattice temperature was initialized to 300K and the above electronic stopping model was applied to all the atoms. The only modificationrequiredfor implementationin MD is to takeinto accountthe contributions frommultiple siliconatoms to the local electron density, while ensuring that the background electron density is only counted once. For each individual cascade, all recoils and the accumulation of damage in the ion path are taken into account. Using the parameter value 1.109 for r◦ from the comparison of BCA Monte Carlo simulation results with the SIMS data, we s have simulated the implantation of low energy boron and arsenic ions into the Si {100}(2×1) surface at energies between 0.5keV and 5keV, with both channeling and off-axis directions of incidence. We mention here that, for the h100i channeling case up to 0.16 keV (32%) of 0.5 keV boron implant energy and 0.64 keV (13%) of 5keV arsenic implant energy are lost via electronic stopping in our simulations. Simulations were terminated when the total energy of the ion became less than 5eV, giving typical simulation times of around 0.2ps. Figs. 9, 10 and 11 show the calculated dopant concentration profile for various energies and directions. Each MD profile is generated by a set of between 500 and 1300 individual ion trajectories. Also shown are the profiles obtained using the modified 6 UT-MARLOWE BCA code described in Sec. III. Obviously, the MD calculation results are in very good agreement with the experimental data, and with the BCA results. This demonstrates that our electronic stopping power model provides a good physically-based damping mechanism for MD simulations of ion implantation. V. CONCLUSION We havedeveloped a phenomenologicalelectronic stopping powermodel for the physics of ionimplantation. It has been implemented into MD and BCA Monte Carlo simulations. SIMS data have been used to verify this model in the MD and BCA Monte Carlo platforms. This model has only one free parameter, namely, the one electron radius of unbound electrons in the medium. We have fine tuned this parameter to obtainexcellent results of dopantprofiles compared with SIMS data in both MD and BCA Monte Carlo simulations. We emphasize that this model with a single parameter can equally successfully model both boron and arsenic implants into silicon over a wide range of energiesandindifferentchannelingandoff-axisdirectionsofincidence. Thisversatilityindicateswideapplicabilityof the model in studies of other physical processes involving electronic stopping. As a more stringent test of the model, it should also be applied to implantation of species other than boron and arsenic. Using arsenic implantation as an example, we have also addressed the issue of how significant electronic stopping is for heavy ions in a low energy regime. For instance, to achieve a good quantitative understanding, we still have to take into account the physics of electronic stopping for arsenic implants at 5keV. Asdiscussedabove,itisimportanttoincorporateion-electroncouplingsintoMDsimulationsinboththehighenergy radiation damage regime and the low energy electron-phonon interaction regime. We have demonstrated that this modelprovidesacrucialpieceofphysicsinMDsimulationsformodelingenergeticcollisionsintheelectronicstopping power regime. The agreement of the simulated dopant profiles with the SIMS data shows that the incorporation of this physically-based damping term into MD simulations is a phenomenologically reliable approach in the regime concerned. Under way is an investigation of whether it can be used as a good phenomenological model for electron- phonon coupling in the low energy regime. This agreement also suggests that MD can be used to generate dopant profiles for testing against the low energy BCA results when experimental data is not available. Furthermore, MD simulations incorporating this physically-based damping mechanism can provide valuable insight into how to modify the binary collision approximation. This will enable the validity of the Monte Carlo simulation to be extended further into the lower energy regime, while not destroying computational efficiency required in realistic simulation environments. VI. ACKNOWLEDGMENT We thank Al Tasch for useful discussions and for providing us with SIMS data, which facilitate validation of the model. This work is performed under the auspices of the U.S. Department of Energy. 1For example, Ion implantation: Science and Technology, edited byJ.F. Ziegler (Academic Press, San Diego, 1988). 2M.T. Robinson and O.S. Oen,Phys. Rev. 132, 2385 (1963). 3J. Lindhard, M. Scharff,and H.E. Schiott, Mat. Fys. Medd. Dan. Vid.Selsk. 33, No. 14 (1963). 4J. Lindhard and M. Scharff,Phys. Rev.124, 128 (1961). 5S.H. Yang, Monte Carlo simulation of arsenic ion implantation into single-crystal silicon, (Univ. of Texas at Austin) Ph.D thesis (1995). 6S.M. Sze, Semiconductor Device Physics and Technology (Wiley & Sons, New York,1985). 7S.H.Yang, S.J. Morris, S.Y Tian, K.B. Parab, and A.F. Tasch, IEEE Trans. Semiconductor Manufacturing, 9, 49 (1996). 8C.S. Murthy and G.R. Srinivasan, IEEE Trans. Electron Devices, 39, 264 (1992). 9M.T. Robinson and I.M. Torrens, Phys.Rev. B 9, 5008 (1974). 10W. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer-Verlag, New York,1991). 11M.T. Robinson, Radiation Effects and Defects in Solids, 130-131, 3 (1994). 12A.M. Stoneham, Nucl. Instrum.Methods, B48, 389 (1990) 13P.B. Allen, Phys. Rev.Lett. 59, 1460 (1987). 14C.P. Flynnand R.S.Averback,Phys.Rev. B 38, 7118 (1988). 7 15M.W. Finnis, P. Agnew, and A.J.E. Foreman, Phys. Rev.B 44, 567 (1991). 16W.L. Johnson, Y.T. Cheng, M. van Rossum, and M.A. Nicolet, Nucl.Instrum. Methods B7/8 657 (1985). 17A.Caro, Radiation Effects and Defects in Solids, 130-131, 187 (1994). 18M.M Jakas and D.E. Harrison, Phys. Rev.B 30, 3573 (1984). 19M.T. Robinson, K. Dan.Vidensk.Selek. Mat. Fys. Medd 43, 27 (1993). 20R.M. Nieminen, K. Dan.Vidensk.Selek. Mat. Fys. Medd 43, 81 (1993). 21N. Grønbech-Jensen, T. Germann, P.S. Lomdahl, and D.M. Beazley, IEEE Computational Science & Engineering, 2, 4 (1995). 22A.Caro and M. Victoria, Phys. Rev.A 40, 2287 (1989). 23P.M. Echenique, R.M. Nieminen and R.H.Ritchie, Solid State Commun. 37, 779 (1981). 24W. Brandt and M. Kitagawa, Phys.Rev.B25, 5631 (1982). 25S.H.Yang,S.Morris, S.Tian, K.Parab,M. Morris, B. Obradovich,C.Snell, andA.F.Tasch, UT-MARLOWEVersion 3.0. Microelectronics Research Center, The University of Texas at Austin (1995). 26S. Kreussler, C. Varelas and W. Brandt, Phys. Rev.B23, 82 (1981). 27J.F. Ziegler, J.P. Biersack and U.Littmark, The Stopping and Range of Ions in Solids (Pergamon Press, New York,1985). 28N.Bohr, Phys.Rev. 58, 654 (1940); 59, 270 (1941). 29W.E. Lamb, Phys. Rev. 58, 696 (1940). 30M. Kitagawa, Nucl. Instrum.Methods B2, 123 (1984). 31R.J. Mathar and M. Posselt, Phys.Rev.B 51, 107 (1995). 32L.C. Northcliffe, Phys.Rev. 120, 1744 (1960). 33R.H.Ritchie, Phys. Rev. 114, 664 (1959). 34P. Hohenbergand W. Kohn, Phys.Rev.B136, 964 (1964). 35W. Kohn and L.J. Sham, Phys.Rev.A 140, 1133 (1965). 36P.M. Echenique, R.M. Nieminen, J.C. Ashley and R.H. Ritchie, Phys.Rev.A 33, 897 (1986). 37P.M. Echeniqueand A.Arnau, Phys.Scr. T49, 677 (1993). 38T.L. Ferrell and R.H. Ritchie, Phys.Rev.B 16, 115 (1977). 39A.Mann and W. Brandt, Phys. Rev.B 24, 4999 (1981). 40W. Brandt, Nucl.Instrum. Methods 194, 13 (1982). 41J. Lindhard and W. Winther,K. Dan.Vidensk.Selsk. Mat. Fys. Medd. 34, No. 4 (1964). 42N.Azziz, K.W. Brannon, and G.R. Srinivasan, MRS Symp.Proc. 45, 71 (1985). 43K.M. Klein, C. Park, and A.F. Tasch, IEEE Trans. Electron. Devices, 39, 1614 (1992). 44K.M. Klein, C. Park, and A.F. Tasch, Appl. Phys.Lett. 57, 2701 (1990). 45O.B. Firsov, Sov.Phys. JETP 36, 1076 (1959). 46SIMS data for these cases are from K.Gartner, M. Nitschke,and W.Eckstein, Nucl. Instrum.Methods, B 83, 87 (1993). 47A.F. Tasch, Private Communication. 48J. Tersoff, Phys. Rev.B 39, 5566 (1989). 8 VII. FIGURES 1018 E=15 keV E=80 keV −3m)1017 N (c O TI A R T E=35 keV N E C ON1016 C 1015 0.0 2.5 5.0 7.5 DEP T H ( 1 03 ANGSTROM) FIG. 1. Boron implant profiles for the h100i direction with energies ranging from 15keV to 80keV. Zero tilt and rotation angles. The thick lines are SIMS data. 1018 E=15 keV E=35 keV E=80 keV −3m)1017 N (c O TI A R T N E C ON1016 C 1015 0.0 2.0 4.0 6.0 8.0 10.0 DEP T H ( 1 03 ANGSTROM) FIG. 2. Boron implant profiles for the h110i direction with energies ranging from 15keV to 80keV. Zero tilt and rotation angles. The thick lines are SIMS data. 1019 E=15 keV 1018 E=35 keV −3m) N (c E=80 keV O ATI1017 R T N E C N O C 1016 1015 0.0 2.0 4.0 6.0 DEPTH (103 ANGSTROM) ◦ ◦ FIG.3. Boron implant profiles for the h100i direction with the tilt =7 and rotation =30 . The thick lines are SIMS data. 9 1019 1018 −3m) N (c O ATI1017 <100> channel, angle=0,0 R T N E C N O C 1016 angle=7o,7o 1015 0.0 1.0 2.0 3.0 4.0 DEPTH (103 ANGSTROM) FIG.4. Boron implant profiles for the h100i channeling and for the off-axis direction with the tilt =7◦ and rotation =7◦. The implant energy is 5keV. Thick lines: SIMS data; Thin lines without symbols: BCA Monte Carlo simulation with full electronicstoppingpower;Circles: BCAMonteCarlosimulation withtheartificially reducedelectronicstoppingpowerforthe channeling case; Triangles: thecorresponding case for the off-axis direction (see text). 1019 E=15 keV 1018 −3m) E=50 keV N (c O ATI1017 R T N E=180 keV E C N O C 1016 1015 0.0 5.0 10.0 15.0 DEP T H ( 1 03 ANGSTROM) FIG.5. Arsenic implant profiles for the h100i direction with energies ranging from 15keV to 180keV. Zero tilt and rotation angles. The thick lines are SIMS data. 1019 E=15 keV 1018 E=100 keV −3m) N (c O ATI1017 R T N E C N O C 1016 1015 0.0 1.0 2.0 3.0 DEPTH (103 ANGSTROM) ◦ ◦ FIG.6. Arsenicimplantprofilesfortheh100idirectionwiththetilt=8 androtation=30 . ThethicklinesareSIMSdata. 10