- Email: [email protected]
Stopping powers for channeled helium ions in silicon using electron densities from bandstructure calculations
Nuclear Instruments and Methods in Physics Research B 85 (1994) 551-555 North-Holland
RIOMI B
Beam Interactions with Materials &Atoms
Stopping powers for channeled helium ions in silicon using electron densities from bandstructure calculations P.W.L. van Dijk, L.J. van IJzendoorn W. van Haeringen, M.J.A. de Voigt
*, M. de Koning, P. Bobbert,
Faculty of Technical Physics, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands
Stopping powers for channeled He ions have been calculated with a modified version of the Monte Carlo code LAROSE [J.H. Barrett, Phys. Rev. B 3 (1971) 15271.The spatial distribution of the valence electron density in Si was obtained from bandstructure calculations. The stopping power was calculated using Lindhard’s free-electron gas approach within the framework of the local density approximation. Spatial variations of the electron density along individual trajectories produce a significant contribution to the energy loss distribution. The average energy loss of 4 MeV He ions channeled in the axial (lOO), (111) and (110) directions have been calculated and compared with measured values. The contribution of the core electrons to the energy loss is investigated by adding the spherically symmetric Hartree-Fock-Slater electron densities of the closed shells to the valence electron density. Calculations show a high energy loss tail in the spectrum qualitatively in agreement with published experimental results.
1. Introduction
Modelling the electronic energy loss of fast ions in matter has been subject of many theoretical investigations. The framework was set by Bohr [l] who used the impulse approximation for the momentum transfer in a semiclassical treatment. The orbit of the ion was treated classically and the excitations of the electrons at large impact parameters was described by the excitation of a harmonic oscillator. Later, Bloch [2] derived the mean energy loss quantummechanically using the impact parameter formalism. The interaction of the moving ion with the electrons in the atomic potential was calculated in the dipole limit of the Born approximation. Basically for each electron a H-like excitation ladder was assumed with an appropriately scaled Bohr radius. Dettmann and Robinson [3,4] refined this scheme by avoiding the dipole approximation and derived analytical expressions for large and small impact parameters using an average exitation energy for each shell. Sigmund and coworkers [5] persisted the use of the model in which the electronic excitations are described with a set of harmonic oscillators even beyond the dipole approximation and including higher Born terms. An alternative first-principle theory was proposed by Lindhard [6] who discussed the properties of a free electron gas. Lindhard described the dynamic behaviour of the electron gas exclusively in terms of the
* Corresponding 438 060.
author, phone +31 40 474 341, fax +31 40
electromagnetic properties of the system which can be justified by the idea that the behaviour of the particles is revealed only through the field to which they give rise. The motion of a particle is then no longer determined by the individual interactions with numerous other particles but only by the field. By introducing a longitudinal and transversal dielectric constant Lindhard derived an expression for the energy loss in terms of the local electron density. Lindhard noted that his free electron gas treatment was particularly useful for the loosely bound electrons in a solid. The application of these theoretical models is specifically interesting for the description of the anomalously low energy loss of channeled ions [7,8]. Although Dettmann [4] derived the energy loss of hyper-channeled ions, his approach was based on the impact parameter formalism with a lattice comprised of independent atoms. The contribution of the valence electrons was treated separately using Lindhard’s approximation for high velocity ions (i.e. uion >> uFermi). Melvin and Tombrello [9] used separate contributions of the so called local electrons with characteristic high momentum transfer (close to the ion trajectory), valence electrons using Lindhard’s model and core electrons. They compared experiments with calculations in a more phenomenological approach using a distribution of channeled particles according to the continuum approximation. However, in these calculations the valence electron density was assumed to be constant. Pathak and coworkers [lo] derived energy loss along the axial (110) direction by introducing explicitly a position dependent electron density. Using Slater or-
0168-583X/94/$07.00 0 1994 - Elsevier Science B.V. All rights reserved SSDI 0168-583X(93)E0858-E
IX. IBA THEORY
552
P. W.L. van Dijk et al. / Nucl. Instr. and Meth. in Phys. Res. B 85 (1994) 551-555
bitals, the spherically symmetric core electron densities of each atom were averaged in the hexagonal shape of the channel, applying geometrical arguments. The valence electron density was again assumed to be constant. Smulders et al. [ll] introduced (electronic) stopping in Monte Carlo trajectory calculations developed to simulate high energy backscattering experiments. Electronic stopping was introduced along the lines of Melvin and Tombrello with Dettmann’s approach for the core electrons. Simulation of backscattering profiles showed a quite satisfactory agreement with experiments. Murthy and Srinivasan [17] implemented electronic energy loss in the binary-collision-cascade code MARLOWE and used the Fourier components of the Si electron density distribution obtained from both calculated and measured X-ray structure factor data [18]. They calculated the energy loss of high energy alpha particles in both channeled and random directions using the energy loss model of Burenkov, Komarov and Kumakhov [19] which is valid for a non-uniform valence electron density distribution. This paper describes an alternative approach to calculate electronic energy loss which uses the explicit spatial electron density distribution from bandstructure calculations in combination with Lindhard’s free electron gas theory. This allows a direct evaluation of dE/dx using the local electron density along an individual ion trajectory. Calculations in the axial (ill), (100) and (110) directions in Si have been performed with a modified version of Barrett’s Monte Carlo binary-collision-code LAROSE [12,13]. First results for 4 MeV He ions show that stopping power differences for the three major axes ((loo), (110) and (111)) can be reproduced. As expected, variations of the energy loss along individual trajectories prove to contribute significantly to the energy loss distribution. In order to investigate the influence of core electrons to the shape of the calculated energy loss distribution, calculations have been carried out in which the core electron density from Hartree-Fock-Slater orbitals was added to the valence electron densitity.
2. Electronic charge distributions
and energy loss cal-
culations in LAROSE
The spatial charge density distribution of the valence electrons in Si has been calculated with the computer code of Denteneer and van Haeringen [14]. Their program is based on density functional theory which describes a system of many electrons with mutual interactions in an external potential. For the latter a so-called pseudo-potential was applied describing the overall potential that the valence electrons feel from both nuclei and core electrons. The ground state prop-
erties of silicon at T = 0 have been used and a spatial distribution of the valence charge density for the conventional unit cell was obtained. In essence their computer code gives the valence electron density for every point in the silicon unit cell. In order to relate the electron density to energy loss for each channeling direction, a stopping-power grid was defined closely related to the unit cell as defined in LAROSE. Matrix transformations were obtained for conversion of the stopping power grid of each particular channeling direction to the conventional Si unit cell. Each volume element in the stopping-power grid has therefore a unique relation to a position in the Si unit cell. The calculated electron density at the centre of each volume element was taken as the local valence electron density of the corresponding volume element in the stopping-power grid. Therefore the choice of the grid size has two restrictions. The volume elements should be small enough to result in a reasonable approximation of the continuous spatial valence electron density but large enough to save computer time. The lower limit for the gridsize is determined by the restriction that summation of all individual valence charge densities over the whole unit cell should be 32, which is the number of valence electrons contained in the conventional unit cell of silicon. This resulted in a lower limit of 512 volume elements for the conventional unit cell which implies 8 volume elements per edge in the (100) direction. Consequently the upperlimit of the volume element edges of the stopping power grid is 0.0679 nm. To calculate the energy loss in a channeling direction using the local density approximation (LDA) we proceed as follows. At first the stopping power grid proportions should be defined for the channeling direction involved. The valence electron densities were calculated at the centre of each volume element and stored in a one-dimensional array. From these electron densities local stopping powers for each volume element were determined using the Lindhard stopping power formula [6]: dE -= dx
(1)
where m is the electron mass, p the electron density, Z, and u are the atomic number and velocity of the ion. Note that the equation above is in SI units. The stopping number L(p, U) was determined from a sixth order polynomal fit of the digitized graphical representation of L(p, u) at a constant velocity (u = 1000 keV/amu) as published by Iafrate and Ziegler [15]. LAROSE was used to calculate 5000 ion trajectories through the crystal. The intersection length of the ion trajectory with each volume element of the stopping power grid is then multiplied by the stopping power of
P. W.L. uan D#k et al. /Nud.
the corresponding volume element. The sum of the local energy losses resulted in the total electronic energy loss of the ion along the trajectory. Note that the calculated local stopping powers are always average values. Energy straggling caused by statistical fluctuations in the local stopping powers can be introduced separately.
1
*
,
,
,,
,
,
*
3
80
5 $
60
,
,
,
,
,
4 MeV He ions only valence
100
,
thickness
electrons
517
nm
40
Energy-loss spectra have been calculated for 5000 trajectories of 4 MeV He ions traversing through 1100 atomic layers in the (ill), (110) and (100) axial channeling directions in silicon. In all cases exact alignment with the crystal axes was used and angular beam divergence was assumed to be absent. Figs. 1 and 2 show the energy-loss distributions calculated with only valence electrons, for the Si(ll0) and Si(ll1) axial channeling directions. Note that for the (110) and (111) axial directions 1100 layers correspond to a thicknesses of 211.2 nm and 517.3 nm respectively. The size of the edge of the volume elements in the stopping-power grid was approximately 0.03 nm for both directions. Table 1 lists the calculated average energy loss (due to valence electrons only) for the three axial directions and the corresponding channeled-to-random ratio’s. Because of the symmetrical distribution the average energy loss is equivalent to the most probable energyloss defined as the “peak value” in ref. [8]. The random stopping power was taken from Ziegler [16] to be 162 eV/nm for 4 MeV He ions. Note that the ions channeled in the (110) direction have the smallest
4 MeV He ions only valence thickness
I
211
eiectrons nm
n I
-I
40 20
~~
t‘,,,,,*
I
3. Results
c
553
Znstr. and Meth. in Phys. Res. B 85 (1994) 551-555
\
ot,..dl*‘,*“,‘\c*‘,*‘,‘**‘.“,.‘-J 40 10 0 Ezrgy
I!&
50
60
(keV)
Fig. 1. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (110) direction using the valence electron density. The initial energy was 4 MeV and 1100 atomic layers have been traversed corresponding to = 211 nm. The accumulated counts correspond to energy intervals of 0.2 keV.
0
20
40
Energy
60
100
720
LOSS (kzt)
Fig. 2. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (111) direction using the valence electron density. The initial energy was 4 MeV and the 1100 atomic layers have been traversed corresponding to = 517 nm. The accumulated counts correspond to energy intervals of 0.2 keV.
loss. The stopping power for ions in the (100) direction is signi~~tly larger than the stopping power for the (110) and (111) direction. Table 1 shows that the implementation of the valence electron density distribution in the calculations accounts for the trends in the experimentally observed stopping power differences. Insertion of a constant valence electron density of 32 electrons per unit cell in formula (1) would result in a channeled-to-random ratio of 58% for all three directions. Comparing the differences between the measured [8] and calculated channeled-to-random ratio’s, we observe that the difference is significantly larger for the (111) direction (11%) than for the (110) direction (3%). Assuming that these differences are caused by the neglect of core electrons, it seems plausible that the discrepancy for the (110) direction is smallest since from a geometrical point of view, core electrons have the smallest relative energy
Table 1 Calculated average energy loss along the (loo), (110) and (111) axes in Si using only the valence electron density. The measured channeled to random ratio’s originate from from Eisen et al. [S] and correspond to the “peak values” equivalent to the most probable energy loss Axial direction
Calculated average energy loss (only valence electrons)
Calculated channeled-torandom ratio (only valence electrons)
Measured channeled-torandom ratio (peak values)
(110)
73.1 eV/nm 82.2 eV/nm 93.7 eV/nm
46% 51% 58%
49% 62%
(111) (100)
IX. IBA THEORY
554
P. W.L. uan Dijk et al. /Nucl.
Instr. and Meth. in Phys. Rex B 85 (1994) 551-555
contribution to the average stopping power for that direction. The calculated widths of the energy loss distributions (see Figs. 1 and 2) are caused by differences in energy loss for trajectories moving through different regions of the channels and probing different electron densities. The widths we observe in the calculated energy-loss spectra are of the order of the Bohr energy straggling. For instance the straggling as given by the Bohr formula for a He ion travelling 500 nm through silicon is 3.2 keV. If we assume a Gaussian distribution in the energy loss this will lead to a FWHM in the energy-loss spectrum of 7.5 keV. It would be of interest, to unravel the contributions to the energy-loss distributions of energy straggling (i.e. statistical fluctuations in the actual ion-electron interaction) and energy-loss differences due to variations in electron density. In order to achieve this, a local stopping-power distribution has to be defined for each volume element in the stopping-power grid rather than only average stopping powers. If we assume these distributions to be Gaussian, we can specify a local straggling per unit length for each volume element. Although implementation of such a procedure introduces another sampling routine and consequently another drastic increase in the necessary computer time, the ever continuing increase in processing speed will bring this within the realm of possibilities soon. It is important to note that path-length differences between various channeled trajectories do not contribute significantly to the calculated widths in Figs. 1 and 2. Relative path-length differences with respect to the total travelled distance are of the order of 0.1%. This is also illustrated in Fig. 3 which is a plot showing the number of trajectories as a function of the energy loss and pathlength. No correlation is found between the absolute pathlength and the magnitude of the stopping power which confirms the possibility of probing specific areas of the channels keeping the overall path-length constant.
only valence thickness
211
electrons nm
Fig. 3. Number of trajectories as a function of energy loss of the 4 MeV He ions and pathlength for the Si(ll0) direction. The number is calculated for energy intervals of 0.8 keV and pathlength intervals of 0.01 nm.
The spectrum shows an asymmetric energy loss distribution with a high energy loss tail which is in qualitative agreement with experimental results [8,21]. The tail is caused by the ions which traverse oscillating trajectories in the channels thus probing the regions with high electron densities near the lattice rows. Although the average energy loss increases after addition of the core electron density, the most probable energy loss has not changed significantly. This is in disagreement with our expectation that the incorporation of core electrons would increase the channeled to random ratio to the experimental observed peak values (Table
4. Core electrons In order to determine the influence of the interaction of ions with core electrons on the shape of the energy loss distributions, the spherically averaged core electron density calculated by Herman and Skillman [20] was added to the valence electron density of the band-structure calculations. This results in an increased electron density within a radius of 0.7 nm of the target atom positions. The spherically symmetric charge density of the core electrons is thus approximated in two surrounding layers of volume elements. The size of the edges of one volume element was 0.33 A. Fig. 4 shows the result of a calculation for the axial (110) direction without changing the stopping model.
:;t,,r’,,,u 0
10
20
Energy
30
40
50
60
Loss (keV)
Fig. 4. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (110) direction with incorporation of the core electron density. The initial energy was 4 MeV and 1100 atomic layers have been traversed correspondcounts correspond to ing to - 211 nm. The accumulated energy intervals of 0.2 keV.
P. WL. van Dijk et al. /Nucl.
Instr. and Meth. in Phys. Res. B 85 (1994) 551-555
1). This behaviour is most probably due to the limited applicability of the Lindhard free electron gas model for strongly bound electrons. In general, also the question can be raised to what extent Lindhard’s formula can be combined with the local density approximation in these calculations. The stopping number L contains contributions of close collisions and distant collisions and the latter obviously probe the electron density over impact parameters larger than the size of a stoppingpower volume element. It might be interesting to evaluate both terms separately where only the close collisions feel the local electron density. Alternatively, the contribution of the core electrons to the stopping power can be calculated separately with a second stopping model (e.g. the impact-parameter dependent model of Sigmund [5]). These questions are subject of further study.
5. Conclusions With the introduction of the spatial valence electron density distribution in Monte Carlo trajectory calculations we are able to reproduce average stopping power differences in the axial (ill), (110) and (100) directions. However, a quantitative comparison of the average energy loss as well as the detailed profile of the energy loss distribution with experiments can only be expected once the contribution of core electrons to the stopping power has been included properly in the calculations. The natural separation of valence and core electrons introduced by the bandstructure calculations allows a simple addition of the two electron densities. The shape of the energy loss distribution without core electrons is caused by spatially varying electron densities along different ion trajectories, it is nearly symmetrical and has a width slightly larger than expected from Bohr straggling. So far, addition of core electron densities resulted in a high energy loss tail which is in qualitative agreement with experiments.
55.5
Acknowledgements
We would like to acknowledge dr. P.D. Jarvis (University of Tasmania) for stimulating discussions during his visit to Eindhoven. References [l] N. Bohr, Phil. Mag. 20 (1913) 10. [2] H. Bethe, Ann. Physik 5 (1930) 325; and F. Bloch, Ann. Physik 16 (1933) 285. [3] K. Dettmann and M.T. Robinson, Phys. Rev. B 10 (1974) [4] k. Dettmann, 2. Physik A 272 (1975) 227. [5] E.H. Mortensen, H.H. Mikkelsen and P. Sigmund, Nucl. Instr. and Meth. B 61 (1991) 139. [6] J. Lindhard, Kgl. Danske. Videnskab. Selskab., Mat. Fys. Medd. 28(8) (1954) 1. [7] B.R. Appleton, C. Erginsoy and W.M. Gibson, Phys. Rev. 161 (1967) 330. [8] F.H. Eisen, G.J. Clark, J. Bottiger, J.M. Poate, Radiat. Eff. 13 (1972) 93. [9] J.D. Melvin and T.A. Tombrello, Radiat. Eff. 26 (1975) 113. [lo] R. Agnihotri and A.P. Pathak, Nucl. Instr. and Meth. B 67 (1992) 39. [ll] P.J.M. Smulders and D.O. Boerma, Nucl. Instr. and Meth. B 29 (1987) 471. [12] J.H. Barrett, Phys. Rev. B 3 (1971) 1527. [13] J.H. Barrett, Nucl. Instr. and Meth. B 44 (1990) 367. [14] P.J.H. Denteneer and W. van Haeringen, J. Phys. C 18 (1985) 4127. [15] G.J. Iafrate and J.F. Ziegler, J. Appl. Phys. 50 (1979) 5579. 1161 J.F. Ziegler, B.F. Biersack, U. Littmark, The Stopping and Range of Ions in Solids, vol. 4 (Pergamon, 1977). [17] C.S. Murthy and G.R. Srinivasan, Phys. Rev. B. 47 (1993) 1256. [18] C.M. Bertoni, V. Bortolani, C. Calandra and F. Nizzoli, J. Phys. C 6 (1973) 3612. [19] A.F. Burenkov, F.F. Komarov and M.A. Kumakhov, Sov. Phys. JETP 51 (1980) 4. [20] F. Herman, S. Skillman, Atomic Structure Calculations (Prentice-Hall, 1963). [21] M. Cholewa, G. Bench, A. Saint, G.J.F. Legge and L. Wielunski, Nucl. Instr. and Meth. B 56/57 (1991) 795.
IX. IBA THEORY
RIOMI B
Beam Interactions with Materials &Atoms
Stopping powers for channeled helium ions in silicon using electron densities from bandstructure calculations P.W.L. van Dijk, L.J. van IJzendoorn W. van Haeringen, M.J.A. de Voigt
*, M. de Koning, P. Bobbert,
Faculty of Technical Physics, Eindhoven University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands
Stopping powers for channeled He ions have been calculated with a modified version of the Monte Carlo code LAROSE [J.H. Barrett, Phys. Rev. B 3 (1971) 15271.The spatial distribution of the valence electron density in Si was obtained from bandstructure calculations. The stopping power was calculated using Lindhard’s free-electron gas approach within the framework of the local density approximation. Spatial variations of the electron density along individual trajectories produce a significant contribution to the energy loss distribution. The average energy loss of 4 MeV He ions channeled in the axial (lOO), (111) and (110) directions have been calculated and compared with measured values. The contribution of the core electrons to the energy loss is investigated by adding the spherically symmetric Hartree-Fock-Slater electron densities of the closed shells to the valence electron density. Calculations show a high energy loss tail in the spectrum qualitatively in agreement with published experimental results.
1. Introduction
Modelling the electronic energy loss of fast ions in matter has been subject of many theoretical investigations. The framework was set by Bohr [l] who used the impulse approximation for the momentum transfer in a semiclassical treatment. The orbit of the ion was treated classically and the excitations of the electrons at large impact parameters was described by the excitation of a harmonic oscillator. Later, Bloch [2] derived the mean energy loss quantummechanically using the impact parameter formalism. The interaction of the moving ion with the electrons in the atomic potential was calculated in the dipole limit of the Born approximation. Basically for each electron a H-like excitation ladder was assumed with an appropriately scaled Bohr radius. Dettmann and Robinson [3,4] refined this scheme by avoiding the dipole approximation and derived analytical expressions for large and small impact parameters using an average exitation energy for each shell. Sigmund and coworkers [5] persisted the use of the model in which the electronic excitations are described with a set of harmonic oscillators even beyond the dipole approximation and including higher Born terms. An alternative first-principle theory was proposed by Lindhard [6] who discussed the properties of a free electron gas. Lindhard described the dynamic behaviour of the electron gas exclusively in terms of the
* Corresponding 438 060.
author, phone +31 40 474 341, fax +31 40
electromagnetic properties of the system which can be justified by the idea that the behaviour of the particles is revealed only through the field to which they give rise. The motion of a particle is then no longer determined by the individual interactions with numerous other particles but only by the field. By introducing a longitudinal and transversal dielectric constant Lindhard derived an expression for the energy loss in terms of the local electron density. Lindhard noted that his free electron gas treatment was particularly useful for the loosely bound electrons in a solid. The application of these theoretical models is specifically interesting for the description of the anomalously low energy loss of channeled ions [7,8]. Although Dettmann [4] derived the energy loss of hyper-channeled ions, his approach was based on the impact parameter formalism with a lattice comprised of independent atoms. The contribution of the valence electrons was treated separately using Lindhard’s approximation for high velocity ions (i.e. uion >> uFermi). Melvin and Tombrello [9] used separate contributions of the so called local electrons with characteristic high momentum transfer (close to the ion trajectory), valence electrons using Lindhard’s model and core electrons. They compared experiments with calculations in a more phenomenological approach using a distribution of channeled particles according to the continuum approximation. However, in these calculations the valence electron density was assumed to be constant. Pathak and coworkers [lo] derived energy loss along the axial (110) direction by introducing explicitly a position dependent electron density. Using Slater or-
0168-583X/94/$07.00 0 1994 - Elsevier Science B.V. All rights reserved SSDI 0168-583X(93)E0858-E
IX. IBA THEORY
552
P. W.L. van Dijk et al. / Nucl. Instr. and Meth. in Phys. Res. B 85 (1994) 551-555
bitals, the spherically symmetric core electron densities of each atom were averaged in the hexagonal shape of the channel, applying geometrical arguments. The valence electron density was again assumed to be constant. Smulders et al. [ll] introduced (electronic) stopping in Monte Carlo trajectory calculations developed to simulate high energy backscattering experiments. Electronic stopping was introduced along the lines of Melvin and Tombrello with Dettmann’s approach for the core electrons. Simulation of backscattering profiles showed a quite satisfactory agreement with experiments. Murthy and Srinivasan [17] implemented electronic energy loss in the binary-collision-cascade code MARLOWE and used the Fourier components of the Si electron density distribution obtained from both calculated and measured X-ray structure factor data [18]. They calculated the energy loss of high energy alpha particles in both channeled and random directions using the energy loss model of Burenkov, Komarov and Kumakhov [19] which is valid for a non-uniform valence electron density distribution. This paper describes an alternative approach to calculate electronic energy loss which uses the explicit spatial electron density distribution from bandstructure calculations in combination with Lindhard’s free electron gas theory. This allows a direct evaluation of dE/dx using the local electron density along an individual ion trajectory. Calculations in the axial (ill), (100) and (110) directions in Si have been performed with a modified version of Barrett’s Monte Carlo binary-collision-code LAROSE [12,13]. First results for 4 MeV He ions show that stopping power differences for the three major axes ((loo), (110) and (111)) can be reproduced. As expected, variations of the energy loss along individual trajectories prove to contribute significantly to the energy loss distribution. In order to investigate the influence of core electrons to the shape of the calculated energy loss distribution, calculations have been carried out in which the core electron density from Hartree-Fock-Slater orbitals was added to the valence electron densitity.
2. Electronic charge distributions
and energy loss cal-
culations in LAROSE
The spatial charge density distribution of the valence electrons in Si has been calculated with the computer code of Denteneer and van Haeringen [14]. Their program is based on density functional theory which describes a system of many electrons with mutual interactions in an external potential. For the latter a so-called pseudo-potential was applied describing the overall potential that the valence electrons feel from both nuclei and core electrons. The ground state prop-
erties of silicon at T = 0 have been used and a spatial distribution of the valence charge density for the conventional unit cell was obtained. In essence their computer code gives the valence electron density for every point in the silicon unit cell. In order to relate the electron density to energy loss for each channeling direction, a stopping-power grid was defined closely related to the unit cell as defined in LAROSE. Matrix transformations were obtained for conversion of the stopping power grid of each particular channeling direction to the conventional Si unit cell. Each volume element in the stopping-power grid has therefore a unique relation to a position in the Si unit cell. The calculated electron density at the centre of each volume element was taken as the local valence electron density of the corresponding volume element in the stopping-power grid. Therefore the choice of the grid size has two restrictions. The volume elements should be small enough to result in a reasonable approximation of the continuous spatial valence electron density but large enough to save computer time. The lower limit for the gridsize is determined by the restriction that summation of all individual valence charge densities over the whole unit cell should be 32, which is the number of valence electrons contained in the conventional unit cell of silicon. This resulted in a lower limit of 512 volume elements for the conventional unit cell which implies 8 volume elements per edge in the (100) direction. Consequently the upperlimit of the volume element edges of the stopping power grid is 0.0679 nm. To calculate the energy loss in a channeling direction using the local density approximation (LDA) we proceed as follows. At first the stopping power grid proportions should be defined for the channeling direction involved. The valence electron densities were calculated at the centre of each volume element and stored in a one-dimensional array. From these electron densities local stopping powers for each volume element were determined using the Lindhard stopping power formula [6]: dE -= dx
(1)
where m is the electron mass, p the electron density, Z, and u are the atomic number and velocity of the ion. Note that the equation above is in SI units. The stopping number L(p, U) was determined from a sixth order polynomal fit of the digitized graphical representation of L(p, u) at a constant velocity (u = 1000 keV/amu) as published by Iafrate and Ziegler [15]. LAROSE was used to calculate 5000 ion trajectories through the crystal. The intersection length of the ion trajectory with each volume element of the stopping power grid is then multiplied by the stopping power of
P. W.L. uan D#k et al. /Nud.
the corresponding volume element. The sum of the local energy losses resulted in the total electronic energy loss of the ion along the trajectory. Note that the calculated local stopping powers are always average values. Energy straggling caused by statistical fluctuations in the local stopping powers can be introduced separately.
1
*
,
,
,,
,
,
*
3
80
5 $
60
,
,
,
,
,
4 MeV He ions only valence
100
,
thickness
electrons
517
nm
40
Energy-loss spectra have been calculated for 5000 trajectories of 4 MeV He ions traversing through 1100 atomic layers in the (ill), (110) and (100) axial channeling directions in silicon. In all cases exact alignment with the crystal axes was used and angular beam divergence was assumed to be absent. Figs. 1 and 2 show the energy-loss distributions calculated with only valence electrons, for the Si(ll0) and Si(ll1) axial channeling directions. Note that for the (110) and (111) axial directions 1100 layers correspond to a thicknesses of 211.2 nm and 517.3 nm respectively. The size of the edge of the volume elements in the stopping-power grid was approximately 0.03 nm for both directions. Table 1 lists the calculated average energy loss (due to valence electrons only) for the three axial directions and the corresponding channeled-to-random ratio’s. Because of the symmetrical distribution the average energy loss is equivalent to the most probable energyloss defined as the “peak value” in ref. [8]. The random stopping power was taken from Ziegler [16] to be 162 eV/nm for 4 MeV He ions. Note that the ions channeled in the (110) direction have the smallest
4 MeV He ions only valence thickness
I
211
eiectrons nm
n I
-I
40 20
~~
t‘,,,,,*
I
3. Results
c
553
Znstr. and Meth. in Phys. Res. B 85 (1994) 551-555
\
ot,..dl*‘,*“,‘\c*‘,*‘,‘**‘.“,.‘-J 40 10 0 Ezrgy
I!&
50
60
(keV)
Fig. 1. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (110) direction using the valence electron density. The initial energy was 4 MeV and 1100 atomic layers have been traversed corresponding to = 211 nm. The accumulated counts correspond to energy intervals of 0.2 keV.
0
20
40
Energy
60
100
720
LOSS (kzt)
Fig. 2. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (111) direction using the valence electron density. The initial energy was 4 MeV and the 1100 atomic layers have been traversed corresponding to = 517 nm. The accumulated counts correspond to energy intervals of 0.2 keV.
loss. The stopping power for ions in the (100) direction is signi~~tly larger than the stopping power for the (110) and (111) direction. Table 1 shows that the implementation of the valence electron density distribution in the calculations accounts for the trends in the experimentally observed stopping power differences. Insertion of a constant valence electron density of 32 electrons per unit cell in formula (1) would result in a channeled-to-random ratio of 58% for all three directions. Comparing the differences between the measured [8] and calculated channeled-to-random ratio’s, we observe that the difference is significantly larger for the (111) direction (11%) than for the (110) direction (3%). Assuming that these differences are caused by the neglect of core electrons, it seems plausible that the discrepancy for the (110) direction is smallest since from a geometrical point of view, core electrons have the smallest relative energy
Table 1 Calculated average energy loss along the (loo), (110) and (111) axes in Si using only the valence electron density. The measured channeled to random ratio’s originate from from Eisen et al. [S] and correspond to the “peak values” equivalent to the most probable energy loss Axial direction
Calculated average energy loss (only valence electrons)
Calculated channeled-torandom ratio (only valence electrons)
Measured channeled-torandom ratio (peak values)
(110)
73.1 eV/nm 82.2 eV/nm 93.7 eV/nm
46% 51% 58%
49% 62%
(111) (100)
IX. IBA THEORY
554
P. W.L. uan Dijk et al. /Nucl.
Instr. and Meth. in Phys. Rex B 85 (1994) 551-555
contribution to the average stopping power for that direction. The calculated widths of the energy loss distributions (see Figs. 1 and 2) are caused by differences in energy loss for trajectories moving through different regions of the channels and probing different electron densities. The widths we observe in the calculated energy-loss spectra are of the order of the Bohr energy straggling. For instance the straggling as given by the Bohr formula for a He ion travelling 500 nm through silicon is 3.2 keV. If we assume a Gaussian distribution in the energy loss this will lead to a FWHM in the energy-loss spectrum of 7.5 keV. It would be of interest, to unravel the contributions to the energy-loss distributions of energy straggling (i.e. statistical fluctuations in the actual ion-electron interaction) and energy-loss differences due to variations in electron density. In order to achieve this, a local stopping-power distribution has to be defined for each volume element in the stopping-power grid rather than only average stopping powers. If we assume these distributions to be Gaussian, we can specify a local straggling per unit length for each volume element. Although implementation of such a procedure introduces another sampling routine and consequently another drastic increase in the necessary computer time, the ever continuing increase in processing speed will bring this within the realm of possibilities soon. It is important to note that path-length differences between various channeled trajectories do not contribute significantly to the calculated widths in Figs. 1 and 2. Relative path-length differences with respect to the total travelled distance are of the order of 0.1%. This is also illustrated in Fig. 3 which is a plot showing the number of trajectories as a function of the energy loss and pathlength. No correlation is found between the absolute pathlength and the magnitude of the stopping power which confirms the possibility of probing specific areas of the channels keeping the overall path-length constant.
only valence thickness
211
electrons nm
Fig. 3. Number of trajectories as a function of energy loss of the 4 MeV He ions and pathlength for the Si(ll0) direction. The number is calculated for energy intervals of 0.8 keV and pathlength intervals of 0.01 nm.
The spectrum shows an asymmetric energy loss distribution with a high energy loss tail which is in qualitative agreement with experimental results [8,21]. The tail is caused by the ions which traverse oscillating trajectories in the channels thus probing the regions with high electron densities near the lattice rows. Although the average energy loss increases after addition of the core electron density, the most probable energy loss has not changed significantly. This is in disagreement with our expectation that the incorporation of core electrons would increase the channeled to random ratio to the experimental observed peak values (Table
4. Core electrons In order to determine the influence of the interaction of ions with core electrons on the shape of the energy loss distributions, the spherically averaged core electron density calculated by Herman and Skillman [20] was added to the valence electron density of the band-structure calculations. This results in an increased electron density within a radius of 0.7 nm of the target atom positions. The spherically symmetric charge density of the core electrons is thus approximated in two surrounding layers of volume elements. The size of the edges of one volume element was 0.33 A. Fig. 4 shows the result of a calculation for the axial (110) direction without changing the stopping model.
:;t,,r’,,,u 0
10
20
Energy
30
40
50
60
Loss (keV)
Fig. 4. Energy loss distribution calculated for 5000 trajectories of channeled He ions along the (110) direction with incorporation of the core electron density. The initial energy was 4 MeV and 1100 atomic layers have been traversed correspondcounts correspond to ing to - 211 nm. The accumulated energy intervals of 0.2 keV.
P. WL. van Dijk et al. /Nucl.
Instr. and Meth. in Phys. Res. B 85 (1994) 551-555
1). This behaviour is most probably due to the limited applicability of the Lindhard free electron gas model for strongly bound electrons. In general, also the question can be raised to what extent Lindhard’s formula can be combined with the local density approximation in these calculations. The stopping number L contains contributions of close collisions and distant collisions and the latter obviously probe the electron density over impact parameters larger than the size of a stoppingpower volume element. It might be interesting to evaluate both terms separately where only the close collisions feel the local electron density. Alternatively, the contribution of the core electrons to the stopping power can be calculated separately with a second stopping model (e.g. the impact-parameter dependent model of Sigmund [5]). These questions are subject of further study.
5. Conclusions With the introduction of the spatial valence electron density distribution in Monte Carlo trajectory calculations we are able to reproduce average stopping power differences in the axial (ill), (110) and (100) directions. However, a quantitative comparison of the average energy loss as well as the detailed profile of the energy loss distribution with experiments can only be expected once the contribution of core electrons to the stopping power has been included properly in the calculations. The natural separation of valence and core electrons introduced by the bandstructure calculations allows a simple addition of the two electron densities. The shape of the energy loss distribution without core electrons is caused by spatially varying electron densities along different ion trajectories, it is nearly symmetrical and has a width slightly larger than expected from Bohr straggling. So far, addition of core electron densities resulted in a high energy loss tail which is in qualitative agreement with experiments.
55.5
Acknowledgements
We would like to acknowledge dr. P.D. Jarvis (University of Tasmania) for stimulating discussions during his visit to Eindhoven. References [l] N. Bohr, Phil. Mag. 20 (1913) 10. [2] H. Bethe, Ann. Physik 5 (1930) 325; and F. Bloch, Ann. Physik 16 (1933) 285. [3] K. Dettmann and M.T. Robinson, Phys. Rev. B 10 (1974) [4] k. Dettmann, 2. Physik A 272 (1975) 227. [5] E.H. Mortensen, H.H. Mikkelsen and P. Sigmund, Nucl. Instr. and Meth. B 61 (1991) 139. [6] J. Lindhard, Kgl. Danske. Videnskab. Selskab., Mat. Fys. Medd. 28(8) (1954) 1. [7] B.R. Appleton, C. Erginsoy and W.M. Gibson, Phys. Rev. 161 (1967) 330. [8] F.H. Eisen, G.J. Clark, J. Bottiger, J.M. Poate, Radiat. Eff. 13 (1972) 93. [9] J.D. Melvin and T.A. Tombrello, Radiat. Eff. 26 (1975) 113. [lo] R. Agnihotri and A.P. Pathak, Nucl. Instr. and Meth. B 67 (1992) 39. [ll] P.J.M. Smulders and D.O. Boerma, Nucl. Instr. and Meth. B 29 (1987) 471. [12] J.H. Barrett, Phys. Rev. B 3 (1971) 1527. [13] J.H. Barrett, Nucl. Instr. and Meth. B 44 (1990) 367. [14] P.J.H. Denteneer and W. van Haeringen, J. Phys. C 18 (1985) 4127. [15] G.J. Iafrate and J.F. Ziegler, J. Appl. Phys. 50 (1979) 5579. 1161 J.F. Ziegler, B.F. Biersack, U. Littmark, The Stopping and Range of Ions in Solids, vol. 4 (Pergamon, 1977). [17] C.S. Murthy and G.R. Srinivasan, Phys. Rev. B. 47 (1993) 1256. [18] C.M. Bertoni, V. Bortolani, C. Calandra and F. Nizzoli, J. Phys. C 6 (1973) 3612. [19] A.F. Burenkov, F.F. Komarov and M.A. Kumakhov, Sov. Phys. JETP 51 (1980) 4. [20] F. Herman, S. Skillman, Atomic Structure Calculations (Prentice-Hall, 1963). [21] M. Cholewa, G. Bench, A. Saint, G.J.F. Legge and L. Wielunski, Nucl. Instr. and Meth. B 56/57 (1991) 795.
IX. IBA THEORY