- Email: [email protected]
Fast formulae for the time integral in binary collision calculations
Nuclear Instruments North-Holland
NIOMI B
and Methods in Physics Research B 84 ( 1994) 425-433
Beam Interactions with Materials 8 Atoms
Fast formulae for the time integral in binary collision calculations Giles Body and Roger Smith ’ Department of Mathematical Sciences, Loughborough University, Leicestershire LEI I 3TU. UK Received 23 August 1993
Asymptotic and numerical formulae for the time integral T which is required in the simulation of a binary collision encounter are developed for some screened Coulomb potentials. These formulae are much faster to evaluate than those calculated by numerical integration and more accurate than the commonly used hard sphere approximation.
1. Introduction In a binary collision computer program, such as TRIM [ 1 ] or MARLOWE [ 21, for simulating ion impact phenomena, the particle paths are replaced by their asymptotes before and after the (assumed binary) collision. The determination of these asymptotes involves the calculation of the scattering angle and the so-called time integral 7. The formula for these asymptotes (which are not always correctly given in the literature) after a two-body collision of an atom of initial kinetic energy Es and mass Mi with a stationary particle of mass A& are given by z,,,
Af sin & ’ = 1 + Afc0sq5,~
A
sin&,,
+ (1 + A)(1 + Afcos&)
+
(A2p’f + A (P’ + of 1~0s &I + P) (1 -I- Afcos4,) (1 + A) ’ (1)
for the scattered atom, and for the recoil -f sin &x Y = 1 -fCOS&n
p’(f
-
+
cos&)
-
7msin&
+
~(1 -fcos+,)
(1 + A)(1 -fcos$,)
,
(2)
where the co-ordinate system is as given in Fig. 1. In these expressions A = M2/M~, f2 = ( EO- Q) /Eo, where Q is the inelastic energy loss in the collision, C#J~is the scattering angle in the centre of mass system, p is the impact parameter, p’ is the modified impact parameter defined below and 7 and 7’ are the time integrals with zm = fr + 7’. Usually the inelastic loss is small, of the order of 1OeV for ions scattering in the energy range 0.1-10’ keV [ 3 ] and so f is often taken to be unity. To evaluate these expressions it is necessary to know the interatomic potential function V(r) between two particles separated by a scalar distance r. Other quantities required include the impact parameters p and p’ and the distance of closest approach R given by the solution to the non-linear equation 1- V(R)/E-p2/R2
= 0;
E = &EC,.
The modified impact parameter p’ is defined in terms of p by p’ = {Ek’/ (E - Q)} ‘12p, where k = { (E - V (R) - Q) / (E - V(R) )} li2. The scattering angle in the centre of mass system C#J~is given by
*
Corresponding
author, tel. +44 509 223192, fax +44 509 211869, e-mail [email protected].
0168-583X/94/$07.00 0 1994 - El sevier Science B.V. All rights reserved SSDZ 0168-583X(93)E0839-9
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
426
I ’ pmjecuk
Fig. 1. The co-ordinate
&I
system in a binary collision between two particles.
=n-pTrm2g(r) dr-p’Jre2g’(r) dr,
R
(4)
R
and the time integrals T and r’ are z =
dw-
I
m(g(r)
-g(r))
dr,
R
-i?(r)) dr, J k’(r) co
7’=
@qx
(5)
R
where g(r)
= { 1 - V(r)/E
-p2/r2}-1’2,
g’(r) = { 1 - V(r)/(E g(r)
- Q) -p’2/r2}-1’2,
= { 1 - p2/r2}-1’2,
g’(r) = { 1 -p12/r2}-1’2. The integrals for for 7 and 7’ can integrals [ 5 ] for The potential V(r)
&,, can be calculated numerically by Gauss-Legendre quadrature [4]. The integrals also be evaluated by numerical quadrature. Analytic expansions of the scattering small scattering angles are not especially useful because they are slowly convergent. function V(r) is often taken to be of the form
= 2 ai exp(-jljr)/r, i=o
(7)
where cyi > 0, Cy=‘=,CK~= 1, /3i > 0, fii < pi+,, i = 0, .a. n - 1, and is derived from numerical fits to the solution of the Thomas-Fermi equation or experimental data. These are known as screened Coulomb potentials and various choices of rz, oi and pi define different interactions. Usually in BC simulations the Moliere [ 61 or Universal potential [ 71 is used. For these screened Coulomb potentials the scattering angles can be evaluated by using Biersack’s magic formula [8], which is a factor of 7 faster than 10 point Gauss-Legendre numerical quadrature [ 91. In most cases the contribution of the time integrals in the accurate calculation of the asymptotic trajectories is not important, but this is not always true, and the hard sphere approximation 7 = p tan &/2 which it is suggested could be used for speed [ 91, can be considerably inaccurate for the screened Coulomb potentials. The hard sphere approximation always systematically overestimates 7, but despite this inaccuracy the effect of
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Rex B 84 (1994) 425-433
421
r on the positions of the asymptotes is usually small but could be important in certain circumstances, such as when the simulation of ion scattering spectroscopy is used in the determination of the surface atomic positions. Consider, for example, a binary collision (with no inelastic energy loss) involving a 3735 eV Ne atom hitting a stationary Cu atom with an impact parameter of p = 0.11 A. Assume the interaction is modelled by the Moliere potential with no inelastic energy loss in the collision. Let the equation of the asymptote for the scattered particle after collision be y = Ax + B. The correct values of &, and r are 106” and 0.044 A, respectively. This gives A = 32.25 and B = 3.04 A. The hard sphere approximation gives r = 0.15 A and B = 4.53 A. So the asymptotes are displaced through x 0.05 A. For Si-Si interactions at 23 1 eV and an impact parameter of p = 0.12 A, correct values of 4m, r and B are 153”, 0.346 A and 1.577 A, respectively. The hard sphere approximation gives T = 0.5 13 A and B = 2.27 A and the asymptotes are displaced through M 0.16 1 A. This may not be so much of a problem in computer simulation of amorphous materials where the positions of the target particles are determined from a statistical distribution, but it could be a problem in the simulation of atomic collisions in crystalline materials where the atoms are in fixed positions. Thus although it is much less important to have an accurate evaluation for T than it is for 4m, there may be circumstances where inaccuracies can lead to errors. In such cases, the direct evaluation of T by numerical integration, where thousands of trajectories may be needed to be run to obtain good statistics, can be time consuming. The purpose of this note is therefore to give some analytic and asymptotic fits to r which are faster than a direct evaluation by numerical quadrature yet more accurate than the hard sphere approximation.
2. Asymptotic formulae for z We first consider four possible asymptotic formulae for r. In this analysis it is assumed that the screened Coulomb potentials are defined by Eq. (7 ) . (i) Elixed,p+oo.FromEq. (3),R =p+p~((~)/(2E)+O(pl/~(p))andr+O.Thescattering angle I$,,, -, 0 and r > 0. (ii) E + co and p # 0 fixed. From Eq. (3) R = p + pV(p)/(2E) + U( l/E’). The scattering angle $,,, + 0, r + 0, and if p is sufficiently small it can be shown that r < 0. (iii) E + 0, p fixed. In this case r + 00, R -+ cm and 13-+ K. (iv) E + cc and &, fixed. In this case p must be of the form p = Z/E, R = c/E + 0( 1/E2 ), where c2 - F2 = c, r < 0, r + 0. Let
J
Z = m(g(r)
-g(r))
(8)
dr.
R
Substituting
r = R( 1 + s) into Eq. (8) gives Co
Z= R
J
(l-(1
+S)-2)-i~cj(l-(1
0
whereto
= 1, ci = -$,..*
q(S) = (1 -P2/R2)[(l
+s)-2)-jq(S)‘-(l-p2R-2(l
+,,~)-~)-f
&,
(9)
0
,cj = (-1)‘(2j)!/(j!222j). + SIm2 -gaiexp(-PiR(l i=O
The function
q(s)
is given by
+ S))/(kaiexp(-/3iR)(1
+ s))].
(10)
i=O
In cases (i) and (ii) 1q (s) 1< 1 and the binomial series in Eq. (9) are convergent. Z is just dm and the second term gives
The j = 0 term for
428
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
7 M
-qR
Jy1-
(1 + s)-W&)
d&r.
(11)
0
In case (ii) for the screened Coulomb potentials, co
J
r=;(l-R2/R2)
CLWexPt-PiR(1 CY=oaiexP[-PiRl(l
(l-(1+&s)_2)-f[(l+,s)-2-
0
+O(R(l
-p2/R2)2).
The integral can be evaluated analytically imation, is r M Y& kaj(RB,Ki(PiR)
and the approximate
+ S)lI ds + S)
(12) form for r, to the order of this approx-
-KO(BiR)).
(13)
i=O
However, in case (i) the limiting process yields r=
t(l
-p2/~2)exp(PoR)(PoRKl(P&)-K~(PoR)),
(14)
since we need only consider the exponential in I/ with the smallest exponent. To the order of magnitude of the approximation, this is just a simplified form of Eq. (12) and simplifies further using the asymptotic form of the Bessel functions for large argument to give 7=
$1
-P~/R~)(/~OR-
;)dm
(15)
An asymptotic formula for &, valid to the same order as Eqs. ( 12) and ( 14) is well known and is called the momentum approximation [ lo]. The asymptotic series for 4m converge only slowly and the second term in the expansion contributes little to improved accuracy for $m [ 5 1. To be consistent with previous work we therefore term Eq. ( 13) as the momentum approximation for 7. . . ... For limit (111) V M 00 e-80’ /r and Eq. (8) yields a0
IX
(1
-p2/X2)-t-jbj(
e-BOX
Ex
(16)
)‘&
where the bj = (- 1 )‘Cj are the binomial coefficients in the expansion of ( 1 -x )-‘i2 - 1, bl = i, b2 = ;, . . .. Some algebraic manipulation of this integral leads to 7 =
(2ln2//3e)
dm-
+ 0(l/R2).
(17) Now we consider limit (iv). This is the most difficult of the limits to evaluate. If p = Z/E then from Eq. (3) 1 - Z2/(E2R2) - CCX~ e-psR /(ER)
= 0,
so R M c/E where 1 - t2/c2 = l/c. The scattering angle &,, is given by &,-,(E,p
= Z/E) = K - 2cE-’
J
O”l Fz(l -C”‘/(E2r2)
- V(r)/E)-1/2
dr.
(18)
clE
Substituting &(E,p
w = rE/c gives = C/E) = K - F
J$ (
1-
1
which has asymptotic
form
c?~/(c~w~)
-
Cai
e-plwclE
/(cw))-‘I2
dw,
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433 Table 1 A comparison between some asymptotic and exact values for a) limit (i) E 1 0.1 10 0.1 0.1 0.1 b) limit (iv) E 20 50 100 10 20 30 100
P
&I
[c-W
P
4m
[deal
10 15 5 5 4 3
0.24 0.42 0.17 14.7 22.4 35.8
0 0 0 0.01 0.05 0.033333 0.01
wz(l
- Z*/(c*w*)
7
FL0149 0.0442 0.0372 0.350 0.385 0.409
7approx0%. (15)) 0.0149 0.0442 0.035 0.344 0.364 0.350
7 -0.069 -0.037 -0.0224 -0.068 -0.0509 -0.0408 -0.018
7approx -0.077 -0.039 -0.0230 -0.088 -0.0583 -0.0447 -0.019
180 180 180 49.07 51.13 51.84 52.77
429
- l/(wc))-‘/*
dw x 2sin-’
7appmxWt. (13)) 0.0149 0.0442 0.0370 0.353 0.390 0.415
(19)
d_&,
1
as E + CC (this just means that V approximates to Coulomb potential fromV(R) = EwehaveR = l/E-~~~‘=,~~~i/E2+O(E-3).NoteI(E,0) - 1 dr. So I(E,O)
=
m(l - Cai
J
e-~sE-‘s-1)-1/2
- 1 ds
i
R
7 gbj(Tctiem”i)j (Em’s_‘)
=
at small distances). If p = 0 = J,“(l-V(r)/E)-‘I*
l/E+...
d.~.
(20)
j=l
The integral Jp” emx/x dx = -y - lnlel + O(lel), where y = Euler’s constant IntegrationbypartsgivesJc’e-X/xj dx=~l-j/(j-l)(l+O(~ln~))forj>l.Thej= is, after substituting Spi = w,
Jm E-‘ai
blk
e-”
dw = blk
/W
i=Op,/E+...
E-‘ai(-y
- ln(Pi/E)
+ O(E-‘))
i=O
=bl
+ E-‘(-~-~~ih(/3i))
E-‘In(E)
i=O
The j > 1 term is by integrating E-ibj
(
[,I-j(l
+
by parts equal to
e-sPi)j]EE
_j)-l(xai
Tsl-j(zai
e-BiS)j-1
1/E
(z 0.57721 ..a). lterm
Caipi
e-BiS
&
)
1 .
430
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
= E-‘bj
(j - 1)-l
+ higher order terms.
Let A = CE”=, bi/(i - 1) = # (1 - x)-‘/~x-~ - C2 - 0.5x-’ dx = In2 + l/2 . Hence Z(E,O) 0.5E-‘In(E) + E-'(A- 0.5~ - 0.5CyCOailnpi) + 0(Em2)so
=
R-Z(E,O) = ; + 0(E-2)-Z(E,0)
r(E, O)=
= -&lnE
+ A
i + Jjy + kCailnpi-ln2
Suppose now p # 0 in limit (iv) so
1 +
0(Em2).
(21)
03cc
ZW,P) =
JC
bj(l - ~2/(E2x2))-j-1’2V(~)jE-j
dx,
(22)
cj~ i=l
where as found above p = xE/c)
c/E,R = c/E,1 - Z2/c2 = l/c. The jth terms equals (substituting w =
00
bjE_’
J
(1 - ;2/ (c2w2) )-j-W
(Tai
~-h~~/~)j,-j,-j+l
&,,.
1
As for case p = 0 we consider j = 1 and j > 1 separately. First the j = 1 term equals Ccl ;E-1 (1 - c”2/(C2W2))-3’2 - 11 Ccti e-BicwlE /W dW + ;E-’ ICcti e-bfcW/E /W dw.
J[
i
1
The second integral gives gives
1
i
SE-'(In E - In(c) - C ai In pi - y ) + higher order terms. The first integral
J[
1
1 M (1 - 2;2/(c2w2))-3/2 2E
- 1 ;
1
dw.
(23)
The j > 1 term to leading order is
J
E-tbj
m(l _ ;2/(c2w2))-i-Vw-jc-j+l
dw.
(24)
1
So the sum of all terms j > 1 is 1 -c
E
J
ca (1 - P/(c2w2)
- l/(wc))-‘/2
- (1 - ,-2/(c22112))-1’2
-;--&
(1 - z,2/(c2w2))-3’2
dw.
1
(25)
Adding Eq. (25) to Eq. (23) a term cancels and we obtain M
J
c( 1 - Z2/(c2w2)
- ~/(wc))-‘/~
I
1 +fi+Zln2-Zln Therefore
1
-c(1-C2/(c2w2))-1/2 -& dw
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. 3 84 (1994) 425-433
431
Table 2 The coefficients IY@and ‘i_nfor the Molihe and Universal polentiah i 1
0”
2
3
a) Molibre potential: yo = 1.269, ye = 0.255, y2 = -0.188, y3 = -0.052,y4 = 0.016, ys = 0.011, y6 = -0.0002. 0 7.292 -0.093 -0.179 -0.083 -0,057 -0.251 1 6.689 -0.i90 6.047 -0.025 2 -0.324 -0.150 3 5.297 -0.323 -0.002 -0.096 4 4.387 -0.248 -0.054 0.005 -0.144 3.384 -0.023 0.003 5 2.451 6 0.078 0.002 0.000 7 1.684 -0.049 -0.001 0.014 1.095 -0.026 -0.002 8 0.013 0.666 9 -0.009 -#.003 0.010 10 0.376 -0.002 0.000 0.006 b) Universal potential: yo = 1.162, yl = 0.275, yz = -0.137, ~3 = -0.064,~4 = 0.001, ~5 = 0.016, ~6 = 0.01. 7,049 -0.073 0 -0.025 -0.168 -0.046 6.115 -0.120 -0.154 1 2 5.239 -0.168 -0.126 -0.024 4.414 3 -0.155 -0.086 -0.009 4 3.607 -0.134 -0.050 -0,001 5 2.844 -0.119 -0.020 0.001 6 2.157 -0.099 -0.000 0.001 7 1.554 -0.065 0.009 0.000 1.046 -0.031 8 0.011 -0.002 9 0.648 -0.009 0.009 -0.002 0.370 10 0.002 -0.002 0.005
I(E,;/E)
= OSE-‘1nE i-0.5y+
+E-‘(-O&niIn~ t’=O
$6+1n2-4.%-0.51n(2c-
1)) +*--,
and hence r(E,?/E)
4 = E
- I
=--OX?‘InE
+ E-r
O.&&la_ [
~+O*~~-~n2+0.5+~*~1~~2~-~)
idl
+ . ..* Some
I
(27)
of the asymptotic values for z are compared to the exact values in Table 1.
3. Chebyshev fits for z The asymptotic formulae whilst useful have a Iimited range of validity. A simulation code requires a fast evaluation of r which is valid in all possible cases in a similar way to which the “magic- formula can be used for +=. In order to obtain such a formula for r we first write T in terms of non-d~mens~ona~~sed units. Let 5 = a?, where a is the screening length. Screened Coulomb potentials can be written in the form V(r) = Ka/rf (r/a). Writing@ = p/a, R = R/a, E = E/K, and v(u) = f (u)/u then
432
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
00
Fig. 2. The time integral Z as a function of non-dimensionalised
(1 - V/E -@2/~2)-t
’ -(/J
energy i and scaled input parameter 5~ l/*/2.3.
- (1 -p2/~‘)-t)
du.
R
Now ? is dependent on the normalised energy E and @. We define a possible scheme for evaluating ? using Chebyshev polynomials, valid in the range 10e3 < E < 10, & (p, E) > &in. Here &in is taken as 10”. To do this define E’ = 0.25(log,,,e + 3.0) and i = 2; - 1. The range of values for < and i are O
j = 2.3EI-+[7&)y0/2
+ f:
T&y,].
(28)
n=l
Now let 6 = 2p/p” - 1 so that the range of values for J? is - 1 5 d 5 1. Define ii by ii = i/5 - 1 where i = Integer Part [ 10; 1. Define 1 = (i - ii) / (ii+ 1 - ii). The formula for Z is given by Z z
,/m
- Iapproxwhere
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
Zapprox = !j~O(i)(Wtii)(l
-A)
+ %(ii+l)~)
+
5
~n~d~~~n~~i~~l
-1)
+ ~n(ii+l)~).
433
(29)
n=l
The coefficients are given in Table 2 for the Moliere and Universal potentials. A comparison of these formulae with an n point Gauss quadrature for r has shown that there is an increase in the speed of evaluation by a factor of about 0.88~ Since a ten point quadrature formula is about the minimum that would normally be required, the Chebyshev fits show a considerable increase in the speed of evaluation of r. The values of r as a function of non-dimensionalised energy and input parameter are plotted in Fig. 2. The maximum error for r derived from these tables is 0.034 where a is the screening length. For the range of values 10V3 < E < 10, & (p, E ) < 10”) the contribution that the term involving 7 makes to the trajectory asymptotes is small. The lower limit for E of 0.001 is chosen because below this energy, the binary collision approximation is not valid. For the Moliere potential /iI < 0.05. Suppose 1 < E < 10, the error in the position of the asymptotes is given by where $r is the true scattering angle. This is negligible so we may assume 7 M 0. For the range of values 10m3 < E < 1, C#I~ (p, E) < 1O”, it is possible to show that I?approx- ?I < 0.05 and the contribution to the asymptote intercepts is also small and may be ignored, although Eq. (13) may be used if higher accuracy is required.
References [ 1 ] J.P. Biersack and W.D. Eckstein, Appl. Phys. A34 (1984) 73. [2] M.T. Robinson, in: Sputtering by Particle Bombardment 1, ed. R. Behrisch (Springer, Berlin, 1981) chap. 3, pp. 73-144. [3] M. Aono and R. Souda, Nucl. Instr. and Meth. B 27 (1988) 55. [4] F.J. Smith, Physica 30 (1964) 497. [S] C. Lehmann and G. Liebfried, Z. Phys. 172 (1963) 465. [6] G. Moliere, Z. Naturforsch. 2a (1947) 133. [7] J.F. Ziegler, J.P. Biersack and U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985). [8] J.P. Biersack and L.G. Haggmark, Nucl. Instr. and Meth. B 174 (1980) 257. [ 91 W.D. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer, Berlin, 1991). [lo] I.M. Torrens, Interatomic Potentials (Academic Press, New York, 1972) p. 154.
NIOMI B
and Methods in Physics Research B 84 ( 1994) 425-433
Beam Interactions with Materials 8 Atoms
Fast formulae for the time integral in binary collision calculations Giles Body and Roger Smith ’ Department of Mathematical Sciences, Loughborough University, Leicestershire LEI I 3TU. UK Received 23 August 1993
Asymptotic and numerical formulae for the time integral T which is required in the simulation of a binary collision encounter are developed for some screened Coulomb potentials. These formulae are much faster to evaluate than those calculated by numerical integration and more accurate than the commonly used hard sphere approximation.
1. Introduction In a binary collision computer program, such as TRIM [ 1 ] or MARLOWE [ 21, for simulating ion impact phenomena, the particle paths are replaced by their asymptotes before and after the (assumed binary) collision. The determination of these asymptotes involves the calculation of the scattering angle and the so-called time integral 7. The formula for these asymptotes (which are not always correctly given in the literature) after a two-body collision of an atom of initial kinetic energy Es and mass Mi with a stationary particle of mass A& are given by z,,,
Af sin & ’ = 1 + Afc0sq5,~
A
sin&,,
+ (1 + A)(1 + Afcos&)
+
(A2p’f + A (P’ + of 1~0s &I + P) (1 -I- Afcos4,) (1 + A) ’ (1)
for the scattered atom, and for the recoil -f sin &x Y = 1 -fCOS&n
p’(f
-
+
cos&)
-
7msin&
+
~(1 -fcos+,)
(1 + A)(1 -fcos$,)
,
(2)
where the co-ordinate system is as given in Fig. 1. In these expressions A = M2/M~, f2 = ( EO- Q) /Eo, where Q is the inelastic energy loss in the collision, C#J~is the scattering angle in the centre of mass system, p is the impact parameter, p’ is the modified impact parameter defined below and 7 and 7’ are the time integrals with zm = fr + 7’. Usually the inelastic loss is small, of the order of 1OeV for ions scattering in the energy range 0.1-10’ keV [ 3 ] and so f is often taken to be unity. To evaluate these expressions it is necessary to know the interatomic potential function V(r) between two particles separated by a scalar distance r. Other quantities required include the impact parameters p and p’ and the distance of closest approach R given by the solution to the non-linear equation 1- V(R)/E-p2/R2
= 0;
E = &EC,.
The modified impact parameter p’ is defined in terms of p by p’ = {Ek’/ (E - Q)} ‘12p, where k = { (E - V (R) - Q) / (E - V(R) )} li2. The scattering angle in the centre of mass system C#J~is given by
*
Corresponding
author, tel. +44 509 223192, fax +44 509 211869, e-mail [email protected].
0168-583X/94/$07.00 0 1994 - El sevier Science B.V. All rights reserved SSDZ 0168-583X(93)E0839-9
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
426
I ’ pmjecuk
Fig. 1. The co-ordinate
&I
system in a binary collision between two particles.
=n-pTrm2g(r) dr-p’Jre2g’(r) dr,
R
(4)
R
and the time integrals T and r’ are z =
dw-
I
m(g(r)
-g(r))
dr,
R
-i?(r)) dr, J k’(r) co
7’=
@qx
(5)
R
where g(r)
= { 1 - V(r)/E
-p2/r2}-1’2,
g’(r) = { 1 - V(r)/(E g(r)
- Q) -p’2/r2}-1’2,
= { 1 - p2/r2}-1’2,
g’(r) = { 1 -p12/r2}-1’2. The integrals for for 7 and 7’ can integrals [ 5 ] for The potential V(r)
&,, can be calculated numerically by Gauss-Legendre quadrature [4]. The integrals also be evaluated by numerical quadrature. Analytic expansions of the scattering small scattering angles are not especially useful because they are slowly convergent. function V(r) is often taken to be of the form
= 2 ai exp(-jljr)/r, i=o
(7)
where cyi > 0, Cy=‘=,CK~= 1, /3i > 0, fii < pi+,, i = 0, .a. n - 1, and is derived from numerical fits to the solution of the Thomas-Fermi equation or experimental data. These are known as screened Coulomb potentials and various choices of rz, oi and pi define different interactions. Usually in BC simulations the Moliere [ 61 or Universal potential [ 71 is used. For these screened Coulomb potentials the scattering angles can be evaluated by using Biersack’s magic formula [8], which is a factor of 7 faster than 10 point Gauss-Legendre numerical quadrature [ 91. In most cases the contribution of the time integrals in the accurate calculation of the asymptotic trajectories is not important, but this is not always true, and the hard sphere approximation 7 = p tan &/2 which it is suggested could be used for speed [ 91, can be considerably inaccurate for the screened Coulomb potentials. The hard sphere approximation always systematically overestimates 7, but despite this inaccuracy the effect of
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Rex B 84 (1994) 425-433
421
r on the positions of the asymptotes is usually small but could be important in certain circumstances, such as when the simulation of ion scattering spectroscopy is used in the determination of the surface atomic positions. Consider, for example, a binary collision (with no inelastic energy loss) involving a 3735 eV Ne atom hitting a stationary Cu atom with an impact parameter of p = 0.11 A. Assume the interaction is modelled by the Moliere potential with no inelastic energy loss in the collision. Let the equation of the asymptote for the scattered particle after collision be y = Ax + B. The correct values of &, and r are 106” and 0.044 A, respectively. This gives A = 32.25 and B = 3.04 A. The hard sphere approximation gives r = 0.15 A and B = 4.53 A. So the asymptotes are displaced through x 0.05 A. For Si-Si interactions at 23 1 eV and an impact parameter of p = 0.12 A, correct values of 4m, r and B are 153”, 0.346 A and 1.577 A, respectively. The hard sphere approximation gives T = 0.5 13 A and B = 2.27 A and the asymptotes are displaced through M 0.16 1 A. This may not be so much of a problem in computer simulation of amorphous materials where the positions of the target particles are determined from a statistical distribution, but it could be a problem in the simulation of atomic collisions in crystalline materials where the atoms are in fixed positions. Thus although it is much less important to have an accurate evaluation for T than it is for 4m, there may be circumstances where inaccuracies can lead to errors. In such cases, the direct evaluation of T by numerical integration, where thousands of trajectories may be needed to be run to obtain good statistics, can be time consuming. The purpose of this note is therefore to give some analytic and asymptotic fits to r which are faster than a direct evaluation by numerical quadrature yet more accurate than the hard sphere approximation.
2. Asymptotic formulae for z We first consider four possible asymptotic formulae for r. In this analysis it is assumed that the screened Coulomb potentials are defined by Eq. (7 ) . (i) Elixed,p+oo.FromEq. (3),R =p+p~((~)/(2E)+O(pl/~(p))andr+O.Thescattering angle I$,,, -, 0 and r > 0. (ii) E + co and p # 0 fixed. From Eq. (3) R = p + pV(p)/(2E) + U( l/E’). The scattering angle $,,, + 0, r + 0, and if p is sufficiently small it can be shown that r < 0. (iii) E + 0, p fixed. In this case r + 00, R -+ cm and 13-+ K. (iv) E + cc and &, fixed. In this case p must be of the form p = Z/E, R = c/E + 0( 1/E2 ), where c2 - F2 = c, r < 0, r + 0. Let
J
Z = m(g(r)
-g(r))
(8)
dr.
R
Substituting
r = R( 1 + s) into Eq. (8) gives Co
Z= R
J
(l-(1
+S)-2)-i~cj(l-(1
0
whereto
= 1, ci = -$,..*
q(S) = (1 -P2/R2)[(l
+s)-2)-jq(S)‘-(l-p2R-2(l
+,,~)-~)-f
&,
(9)
0
,cj = (-1)‘(2j)!/(j!222j). + SIm2 -gaiexp(-PiR(l i=O
The function
q(s)
is given by
+ S))/(kaiexp(-/3iR)(1
+ s))].
(10)
i=O
In cases (i) and (ii) 1q (s) 1< 1 and the binomial series in Eq. (9) are convergent. Z is just dm and the second term gives
The j = 0 term for
428
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
7 M
-qR
Jy1-
(1 + s)-W&)
d&r.
(11)
0
In case (ii) for the screened Coulomb potentials, co
J
r=;(l-R2/R2)
CLWexPt-PiR(1 CY=oaiexP[-PiRl(l
(l-(1+&s)_2)-f[(l+,s)-2-
0
+O(R(l
-p2/R2)2).
The integral can be evaluated analytically imation, is r M Y& kaj(RB,Ki(PiR)
and the approximate
+ S)lI ds + S)
(12) form for r, to the order of this approx-
-KO(BiR)).
(13)
i=O
However, in case (i) the limiting process yields r=
t(l
-p2/~2)exp(PoR)(PoRKl(P&)-K~(PoR)),
(14)
since we need only consider the exponential in I/ with the smallest exponent. To the order of magnitude of the approximation, this is just a simplified form of Eq. (12) and simplifies further using the asymptotic form of the Bessel functions for large argument to give 7=
$1
-P~/R~)(/~OR-
;)dm
(15)
An asymptotic formula for &, valid to the same order as Eqs. ( 12) and ( 14) is well known and is called the momentum approximation [ lo]. The asymptotic series for 4m converge only slowly and the second term in the expansion contributes little to improved accuracy for $m [ 5 1. To be consistent with previous work we therefore term Eq. ( 13) as the momentum approximation for 7. . . ... For limit (111) V M 00 e-80’ /r and Eq. (8) yields a0
IX
(1
-p2/X2)-t-jbj(
e-BOX
Ex
(16)
)‘&
where the bj = (- 1 )‘Cj are the binomial coefficients in the expansion of ( 1 -x )-‘i2 - 1, bl = i, b2 = ;, . . .. Some algebraic manipulation of this integral leads to 7 =
(2ln2//3e)
dm-
+ 0(l/R2).
(17) Now we consider limit (iv). This is the most difficult of the limits to evaluate. If p = Z/E then from Eq. (3) 1 - Z2/(E2R2) - CCX~ e-psR /(ER)
= 0,
so R M c/E where 1 - t2/c2 = l/c. The scattering angle &,, is given by &,-,(E,p
= Z/E) = K - 2cE-’
J
O”l Fz(l -C”‘/(E2r2)
- V(r)/E)-1/2
dr.
(18)
clE
Substituting &(E,p
w = rE/c gives = C/E) = K - F
J$ (
1-
1
which has asymptotic
form
c?~/(c~w~)
-
Cai
e-plwclE
/(cw))-‘I2
dw,
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433 Table 1 A comparison between some asymptotic and exact values for a) limit (i) E 1 0.1 10 0.1 0.1 0.1 b) limit (iv) E 20 50 100 10 20 30 100
P
&I
[c-W
P
4m
[deal
10 15 5 5 4 3
0.24 0.42 0.17 14.7 22.4 35.8
0 0 0 0.01 0.05 0.033333 0.01
wz(l
- Z*/(c*w*)
7
FL0149 0.0442 0.0372 0.350 0.385 0.409
7approx0%. (15)) 0.0149 0.0442 0.035 0.344 0.364 0.350
7 -0.069 -0.037 -0.0224 -0.068 -0.0509 -0.0408 -0.018
7approx -0.077 -0.039 -0.0230 -0.088 -0.0583 -0.0447 -0.019
180 180 180 49.07 51.13 51.84 52.77
429
- l/(wc))-‘/*
dw x 2sin-’
7appmxWt. (13)) 0.0149 0.0442 0.0370 0.353 0.390 0.415
(19)
d_&,
1
as E + CC (this just means that V approximates to Coulomb potential fromV(R) = EwehaveR = l/E-~~~‘=,~~~i/E2+O(E-3).NoteI(E,0) - 1 dr. So I(E,O)
=
m(l - Cai
J
e-~sE-‘s-1)-1/2
- 1 ds
i
R
7 gbj(Tctiem”i)j (Em’s_‘)
=
at small distances). If p = 0 = J,“(l-V(r)/E)-‘I*
l/E+...
d.~.
(20)
j=l
The integral Jp” emx/x dx = -y - lnlel + O(lel), where y = Euler’s constant IntegrationbypartsgivesJc’e-X/xj dx=~l-j/(j-l)(l+O(~ln~))forj>l.Thej= is, after substituting Spi = w,
Jm E-‘ai
blk
e-”
dw = blk
/W
i=Op,/E+...
E-‘ai(-y
- ln(Pi/E)
+ O(E-‘))
i=O
=bl
+ E-‘(-~-~~ih(/3i))
E-‘In(E)
i=O
The j > 1 term is by integrating E-ibj
(
[,I-j(l
+
by parts equal to
e-sPi)j]EE
_j)-l(xai
Tsl-j(zai
e-BiS)j-1
1/E
(z 0.57721 ..a). lterm
Caipi
e-BiS
&
)
1 .
430
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
= E-‘bj
(j - 1)-l
+ higher order terms.
Let A = CE”=, bi/(i - 1) = # (1 - x)-‘/~x-~ - C2 - 0.5x-’ dx = In2 + l/2 . Hence Z(E,O) 0.5E-‘In(E) + E-'(A- 0.5~ - 0.5CyCOailnpi) + 0(Em2)so
=
R-Z(E,O) = ; + 0(E-2)-Z(E,0)
r(E, O)=
= -&lnE
+ A
i + Jjy + kCailnpi-ln2
Suppose now p # 0 in limit (iv) so
1 +
0(Em2).
(21)
03cc
ZW,P) =
JC
bj(l - ~2/(E2x2))-j-1’2V(~)jE-j
dx,
(22)
cj~ i=l
where as found above p = xE/c)
c/E,R = c/E,1 - Z2/c2 = l/c. The jth terms equals (substituting w =
00
bjE_’
J
(1 - ;2/ (c2w2) )-j-W
(Tai
~-h~~/~)j,-j,-j+l
&,,.
1
As for case p = 0 we consider j = 1 and j > 1 separately. First the j = 1 term equals Ccl ;E-1 (1 - c”2/(C2W2))-3’2 - 11 Ccti e-BicwlE /W dW + ;E-’ ICcti e-bfcW/E /W dw.
J[
i
1
The second integral gives gives
1
i
SE-'(In E - In(c) - C ai In pi - y ) + higher order terms. The first integral
J[
1
1 M (1 - 2;2/(c2w2))-3/2 2E
- 1 ;
1
dw.
(23)
The j > 1 term to leading order is
J
E-tbj
m(l _ ;2/(c2w2))-i-Vw-jc-j+l
dw.
(24)
1
So the sum of all terms j > 1 is 1 -c
E
J
ca (1 - P/(c2w2)
- l/(wc))-‘/2
- (1 - ,-2/(c22112))-1’2
-;--&
(1 - z,2/(c2w2))-3’2
dw.
1
(25)
Adding Eq. (25) to Eq. (23) a term cancels and we obtain M
J
c( 1 - Z2/(c2w2)
- ~/(wc))-‘/~
I
1 +fi+Zln2-Zln Therefore
1
-c(1-C2/(c2w2))-1/2 -& dw
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. 3 84 (1994) 425-433
431
Table 2 The coefficients IY@and ‘i_nfor the Molihe and Universal polentiah i 1
0”
2
3
a) Molibre potential: yo = 1.269, ye = 0.255, y2 = -0.188, y3 = -0.052,y4 = 0.016, ys = 0.011, y6 = -0.0002. 0 7.292 -0.093 -0.179 -0.083 -0,057 -0.251 1 6.689 -0.i90 6.047 -0.025 2 -0.324 -0.150 3 5.297 -0.323 -0.002 -0.096 4 4.387 -0.248 -0.054 0.005 -0.144 3.384 -0.023 0.003 5 2.451 6 0.078 0.002 0.000 7 1.684 -0.049 -0.001 0.014 1.095 -0.026 -0.002 8 0.013 0.666 9 -0.009 -#.003 0.010 10 0.376 -0.002 0.000 0.006 b) Universal potential: yo = 1.162, yl = 0.275, yz = -0.137, ~3 = -0.064,~4 = 0.001, ~5 = 0.016, ~6 = 0.01. 7,049 -0.073 0 -0.025 -0.168 -0.046 6.115 -0.120 -0.154 1 2 5.239 -0.168 -0.126 -0.024 4.414 3 -0.155 -0.086 -0.009 4 3.607 -0.134 -0.050 -0,001 5 2.844 -0.119 -0.020 0.001 6 2.157 -0.099 -0.000 0.001 7 1.554 -0.065 0.009 0.000 1.046 -0.031 8 0.011 -0.002 9 0.648 -0.009 0.009 -0.002 0.370 10 0.002 -0.002 0.005
I(E,;/E)
= OSE-‘1nE i-0.5y+
+E-‘(-O&niIn~ t’=O
$6+1n2-4.%-0.51n(2c-
1)) +*--,
and hence r(E,?/E)
4 = E
- I
=--OX?‘InE
+ E-r
O.&&la_ [
~+O*~~-~n2+0.5+~*~1~~2~-~)
idl
+ . ..* Some
I
(27)
of the asymptotic values for z are compared to the exact values in Table 1.
3. Chebyshev fits for z The asymptotic formulae whilst useful have a Iimited range of validity. A simulation code requires a fast evaluation of r which is valid in all possible cases in a similar way to which the “magic- formula can be used for +=. In order to obtain such a formula for r we first write T in terms of non-d~mens~ona~~sed units. Let 5 = a?, where a is the screening length. Screened Coulomb potentials can be written in the form V(r) = Ka/rf (r/a). Writing@ = p/a, R = R/a, E = E/K, and v(u) = f (u)/u then
432
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
00
Fig. 2. The time integral Z as a function of non-dimensionalised
(1 - V/E -@2/~2)-t
’ -(/J
energy i and scaled input parameter 5~ l/*/2.3.
- (1 -p2/~‘)-t)
du.
R
Now ? is dependent on the normalised energy E and @. We define a possible scheme for evaluating ? using Chebyshev polynomials, valid in the range 10e3 < E < 10, & (p, E) > &in. Here &in is taken as 10”. To do this define E’ = 0.25(log,,,e + 3.0) and i = 2; - 1. The range of values for < and i are O
j = 2.3EI-+[7&)y0/2
+ f:
T&y,].
(28)
n=l
Now let 6 = 2p/p” - 1 so that the range of values for J? is - 1 5 d 5 1. Define ii by ii = i/5 - 1 where i = Integer Part [ 10; 1. Define 1 = (i - ii) / (ii+ 1 - ii). The formula for Z is given by Z z
,/m
- Iapproxwhere
G. Body, R. Smith / Nucl. Instr. and Meth. in Phys. Res. B 84 (1994) 425-433
Zapprox = !j~O(i)(Wtii)(l
-A)
+ %(ii+l)~)
+
5
~n~d~~~n~~i~~l
-1)
+ ~n(ii+l)~).
433
(29)
n=l
The coefficients are given in Table 2 for the Moliere and Universal potentials. A comparison of these formulae with an n point Gauss quadrature for r has shown that there is an increase in the speed of evaluation by a factor of about 0.88~ Since a ten point quadrature formula is about the minimum that would normally be required, the Chebyshev fits show a considerable increase in the speed of evaluation of r. The values of r as a function of non-dimensionalised energy and input parameter are plotted in Fig. 2. The maximum error for r derived from these tables is 0.034 where a is the screening length. For the range of values 10V3 < E < 10, & (p, E ) < 10”) the contribution that the term involving 7 makes to the trajectory asymptotes is small. The lower limit for E of 0.001 is chosen because below this energy, the binary collision approximation is not valid. For the Moliere potential /iI < 0.05. Suppose 1 < E < 10, the error in the position of the asymptotes is given by where $r is the true scattering angle. This is negligible so we may assume 7 M 0. For the range of values 10m3 < E < 1, C#I~ (p, E) < 1O”, it is possible to show that I?approx- ?I < 0.05 and the contribution to the asymptote intercepts is also small and may be ignored, although Eq. (13) may be used if higher accuracy is required.
References [ 1 ] J.P. Biersack and W.D. Eckstein, Appl. Phys. A34 (1984) 73. [2] M.T. Robinson, in: Sputtering by Particle Bombardment 1, ed. R. Behrisch (Springer, Berlin, 1981) chap. 3, pp. 73-144. [3] M. Aono and R. Souda, Nucl. Instr. and Meth. B 27 (1988) 55. [4] F.J. Smith, Physica 30 (1964) 497. [S] C. Lehmann and G. Liebfried, Z. Phys. 172 (1963) 465. [6] G. Moliere, Z. Naturforsch. 2a (1947) 133. [7] J.F. Ziegler, J.P. Biersack and U. Littmark, The Stopping and Range of Ions in Solids (Pergamon, New York, 1985). [8] J.P. Biersack and L.G. Haggmark, Nucl. Instr. and Meth. B 174 (1980) 257. [ 91 W.D. Eckstein, Computer Simulation of Ion-Solid Interactions (Springer, Berlin, 1991). [lo] I.M. Torrens, Interatomic Potentials (Academic Press, New York, 1972) p. 154.