The effect of correlated thermal vibrations on crystal blocking patterns

Nuclear Instruments and Methods in Physics Research B28 (1987) lo-16 North-Holland, Amsterdam

10

THE EFFECT OF CORRELATED

THERMAL VIBRATIONS

ON CRYSTAL BLOCKING

PATIXRNS

S. ALLINEY Dipartimento di Matematica dell’liniversitd! and INFN, Sezione di Bologna, Italy

F. MALAGUTI

and E. VERONDINI

Dipartimento di Fisica deN’Universitci and INFN, Sezione di Bologna, Italy

Received 5 March 1987

We discuss a simple way to include correlated thermal vibrations in Monte Carlo simulations of crystal blocking dips. The results are very similar to those given by conventional (i.e. independent) calculations. The reason is that the correlation length is too short compared to the string length needed to generate the blocking pattern.

1. Introduction Since the pioneering work of Ryabov [l], there has been a continuous interest [2] in the role that correlated thermal vibrations might play during propagation of charged particles through crystals. Calculations show that including correlations produces only minor effects in directional phenomena, apart from surface backscattering yields [2]. The reason for such a failure is not completely understood. The aim of the present work is to show that the correlation along a string extends for no more than 2-3 adjacent atoms, while many more atoms are in fact necessary to deflect the ion path in a sizeable way. In other words, the correlation length is too small to produce relevant effects. To simulate correlated ion scattering within a crystal, we will assume a frame of reference with the z-axis along the atomic string. The problem is then to generate a series of atomic thermal displacements orthogonal to the string 5

=

(Xl>

YA

5=(x2,

Ye),...,

m=(x,,

Y,). (1)

Here, we will assume that x and y vibrations are independent, so they can be sampled separately from an equal time [l] distribution that we assume as an multivariate normal, described by a covariance matrix M with elements mik = c*c,k-_i, and correlation

(2) coefficients

that are invariant for translations and reflections along the string. Here u2 is the usual one-dimensional mean-square vibration amplitude, while ci can be calculated in the framework of the Debye theory using e.g. the formalism of Nelson et al. [3] *. The mathematical background needed to describe correlated vibrations is sketched in sect. 2, while the results are discussed in sect. 3. 2. The multivariate normal distribution The n-variate normal distribution N,(O, M) with zero mean and variance-covariance (positive definite) matrix M is defined by the probability density function: f(x)

(4)

to generate

an

2.1. Cholesky factorization As is well known [S], any symmetric n x n matrix M can be factorized as

positive

M = LLT,

definite (5)

where L is a lower triangular matrix, whose diagonal elements are the square roots of the leading principal determinants of M. Since M is positive definite, such determinants are greater than zero. The Cholesky factorization can be obtained using in-place computations [6]. Now, it follows from (5) that M-l=

B.V.

IM-‘/2 1exp[ - QxTMW1x].

Two different methods can be employed n-vector x having this distribution.

*

0168-583X/87/$03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)

= (2?p2

L-TL-1,

The fancier approach by Jackson et al. [4] seems sary for present purposes.

(6) unneces-

11

S. Alliney et al. / Correlated thermal vibrations and crystal blocking L-’ being again a lower triangular matrix, [S]. According to eq. (6), the pdf (4) can be rewritten as f(x)

= (2a)-n’z

]M]-“*

exp[ - $xTL-rL-lx].

exp{ -t[GrWll&

+ x,&E1

+

6TWl2X”

+

5522

I}

(7)

This suggests considering a new vector random variable z defined as z = L-ix

or equivalently:

(8)

(13) The expression in brackets can be put in a more illuminating form by some simple manipulations. First of all, such expression can be rewritten as

and the corresponding pdf j(z)

= (2~)~“”

)M l-l/* exp[ -

tz’z] .

(9)

Here, the variance-covariance matrix reduces to the identity matrix I. Thus, the components zi, z2,. . ., z, are independent zero-mean normal variates. Hence, an obvious method to generate an n-vector x with the required N,,(O, M) distribution is to generate zi, z2,. . ., z, as independent N,(O, 1) variates, forming the n-vector z with N,(O, I), distribution and then deduce x as x= Lt.

(IO)

It should be noted that the given variance-covariance matrix M is to be factorized once and for all at the beginning of the procedure.

+

X” +

1[

~(th2)

w22

[

Now, we can eliminate from of W and give an equivalent original variance-covariance relations follow immediately

x,

+

$(

4&*>

I (14) .

eq. (14) all the submatrices expression in terms of the matrix M. The following from the definitions (11):

MiiW,i +

ml2wl2 T-l -

bw12

+

m12w22

=

0

Wb)

T m12w12

+

m22w22

=

1.

(I5c)

(1%

From (15a) and (15b) we obtain: 2.2. Recursive computation Wi, A second method [7] builds up the n-vector x from univariate conditional distributions. The first result to be recalled is that the marginal distribution of every m-dimensional subset of x1,. . . , x, is m-variate normal [8]. The second is that the conditional distribution of x,, gtven xl, x2,. . . , x,_~ is again normal [9]. According to the previous remarks, the recursive method can be established as follows. We start by partitioning M and M-i as shown below MC

!I!;_-[ ml2

’

m12

,

*22

--

1

}n-1;

=&I’.

(16)

Using again (15b) together with (15~) ~22

=

(m2,

-

mT2MG’m12)-l

(17)

and finally from (15b) w12/w22= -MCi’mi2.

(18)

Thus we can rewrite the original pdf as the product * Ai exp{ -

KMG%, }

x-42 exp( - t( x,, - &rMll’m,2)

1’

X ( rnaz - m~2M~‘m,2)-1

1

II-1

3~124)

X(x,-mT,Mii’I,)).

(11)

(18)

Hence, the distribution N,(O; M) of x can be expressed as the conditional distribution

-n-1

Ni(mT2MG’&;

1

(19)

for x, and the marginal

corresponding to x= [x1, x2,...,

mz2 - mT2M;‘m12)

N,-,(O;

Xn-l/XnlT = [Wnl T.

The pdf (4) then becomes (apart from the normalization coefficient)

(12)

(20)

Ml,)

for 4 = [xi . . . x,_~]~.

The well-known

case n = 2 may

* Al and A, are the standard normalization constants; a classical result states that det(M) = det(Mil) det(m22 GM,’ m12).

12

S. Alliney et al. / Correlated thermal vibrations and crystal blocking

serve as an illustration.

Let us consider

M=a2

(21) From eq. (14) or eq. (18) we have @(X2/% u;

=

> =

Cl%

dvar(x,/x,)

(22) =

a/g.

(23)

2.3. Computational remarks

n-l X,-l}

=

C

c2 = 0.368,

cs = 0.182,

%,lS

a16,14= 0.040,

a,jxj,

(24)

j=l

The results of the preceding section show that the “second neighbour approximation”, eqs. (22) and (23), which has sometimes been used [lo] to describe correlation effects in channeling, appears to be very good when compared with the exact values, eqs. (24) and (26). We therefore introduced such an approximation within our multistring Monte Carlo code Ill], neglecting correlations between adjacent strings. This is justified because during the gentle collision we are concerned with here, the time one ion needs to travel from the proximity of one string to the other is very large compared to that necessary to travel the correlation length, which extends 2-3 interatomic spacings only, as determined by the actual coefficients in eq. (24). The effect of correlations is expected to be higher in the “low energy region” [12] (27)

-w,//w,,.

(25)

b) Its standard deviation is u;=

Var(x,/x,,...,x,_,)

0.3%

3. Discussion

where the drift coefficients anj are ani=

=

cq = 0.109 a16,13= 0.016.

Further discussion on this point is given in Appendix A.

It is not obvious to give a comparative evaluation of the numerical performance of the two methods reported before (for a review of the published data, see ref. [7]). In general, however, the Cholesky factorization seems more efficient: once the factorization of the variancecovariance matrix has been obtained, we do not need any further matrix operation. Using the recursive method, on the other hand, implies the numerical inversion of every leading submatrix IUI’,;’of M. Clearly, this can be done once and for all, but the computational burden may become excessive when the dimension n of the problem increases. In our application, however, the particular structure of the matrix M suggests the choice of the recursive method. In order to discuss this relevant point, we will use eq. (14) instead of eq. (18). We recall that: a) The expected value of x,, given xi, . . . , x, _ 1, is ~{X,/$,...,

of eqs. (25) and (16) starting from correlation coefficients given by the Debye theory [3]. Detailed calculations for different crystals and temperatures shows that for every n: (i) u,,! [eq. (26)] is very near to u; [eq. (23)]. Typical values for Al (110) axis at 300 K are u;/o = 0.930, CI;/U = 0, 928,. . . , o;,/o = 0.927,. . . . (3 a, n-l = c2, a,,n-2 = O.lc,, while the remaining drift coefficients are negligible. For Al (110) axis at 300 K we have e.g.

=w;;‘2.

(26)

Here W denotes the inverse of n X n variance-covariante matrix M. From these results a very simple recursive method follows to sample the displacements [eq. (l)], very convenient for a Monte Carlo simulation: first, choose xi (yi) with zero mean and u as standard deviation; then x2 (y2) using eqs. (22) and (23) together with the chosen value for xi (yi); then x3 (y3) using n = 3 in eqs. (24) and (26) together with the chosen values for xi (yl) and x2 (~~1, and SO on. As n becomes larger and larger, the drift coefficients

(25) stabilize, so that the expectation (24) contains only a limited number of terms. This means that every atom is influenced only by the nearby ones, and not by the distant ones. This is confirmed by numerical evaluation

and d being Lindhard’s critical angle and the (VQL interatomic spacing along the string) than in the “high energy region” Gt,d
(28)

where the continuum model holds [12]. In fact, condition (27) simulates in some way the two-atom model originally proposed by Oen [13], where correlations are of extreme importance. Fig. 1 displays a comparison between correlated and independent Monte Carlo simulations of the blocking dip due to 0.1 MeV a-particles along the Al (110) axis. We are here in the low-energy region, eq. (27), but the difference, though clearly visible, is much smaller than expected (similar calculations in the high-energy region, eq. (28), showed even weaker effects of the correlations). This somewhat surprising result can be understood if we compare the correlation length, which extends 2-3 interatomic spacings, with the number of atomic scatterings needed to produce a sizeable blocking dip. This is

S. Alliney et al. / Correlated thermal vibrations and crystal blocking

13

cal basis by means of a simple model. Let us make a parabolic approximation string potential in the transverse plane F/(r)

to the ion-

= V, - ikr2.

(29)

If we use the Moliere thermally averaged potential, then a straightforward calculations sketched in Appendix B gives k

=

JI2,TF

(30)

7.

d

where T is the ion kinetic energy. The function F, whose explicit expression is given in Appendix B, depends only on the ion and crystal characteristics. One can easily integrate the equations of motion for potential (29) and (30) assuming as initial conditions that the ion starts parallel to the string and at a distance u from it. Then the time needed to get a deflection 2#r (i.e. the full width of the blocking dip) can be evaluated and transformed to a number N of interatomic spacings d. The final result is

101 n (mrad)

Fig. 1. Simulated blocking dips due to 0.1 MeV a-particles on the Al(110) axis (100 planes). The points are the correlated calculation, while the crosses (displayed only where the difference can be seen) are the standard (independent) calculation. GL is Lindhard’s angle, and n is the ordering number of 1 mrad wide rings centered on the dip.

N=&C, shown in fig. 2, where the minimum yield xmin is given for a blocking dip due to a string containing n atoms (the simulation refers to the same case as fig. 1). One sees that some 5 to 10 atoms are necessary, hence definitely more than the 2-3 that are significantly correlated. This number becomes larger and larger as the ion energy increases and is, for example, n = 20-40 for 2-12 MeV o-particles on Al (110) axis. In any case the correlation length is too short to play a significant r6le in the build-up of the blocking dip. These results can receive a semiquantitative analyti-

I

I

I

I

(31)

where the correction

c=

factor

+l-‘(*g).

(32)

The function F (see Appendix B), hence the factor C, depends only on the ratio a/a between the Thomas-Fermi screening length and the one-dimensional thermal vibration amplitude. The behaviour of C vs a/a is shown in fig. 3. Values of the ratio a/a for several ion-crystal com-

I

I

I

E, I

2 12

I

I

k&V

E,.O.l

100 -e

’

I

I

MeV

-- 0 0

n *

ao.I

G

0

.

6040-

0

0 0

20I 1

I 5

I 10

I 20

0 40

zo

n

I 1

I 23

I

,

I 5

Fig. 2. Minimum yield for a blocking dip due to a crystal containing n planes. a-particles on Al(110)

l’

axis.

n

S. Alliney et al. / Correlated thermal vibrations and crystal blocking

14 Table 1 Valuesof the ratio a/o C

Z,

Z,=6

1 2 8 14 32 56 82

5.64 5.29 4.33 3.88 3.20 2.16 2.49

Al

Si

Fe

Ge

MO

Ag

Ta

W

Pt

Au

Pb

13

14

26

32

42

41

13

14

78

19

82

1.75 2.67 1.44 1.32 1.13 0.99 0.90

2.39 2.30 1.99 1.83 1.57 1.38 1.26

2.20 2.14 1.93 1.80 1.59 1.42 1.31

1.66 1.61 1.47 1.38 1.23 1.11 1.02

2.21 2.22 2.05 1.94 1.75 1.59 1.48

1.34 1.32 1.22 1.16 1.05 0.96 0.89

1.70 1.68 1.58 1.52 1.39 1.29 1.21

2.17 2.14 2.01 1.93 1.78 1.65 1.55

1.62 1.59 1.50 1.45 1.33 1.24 1.16

1.22 1.20 1.14 1.09 1.01 0.94 0.88

0.64 0.63 0.60 0.57 0.53 0.49 0.46

binations are collected in table 1. One can seen that 0.5 2 a/u

the kinetic energy dhard angle, as

5 5.5

throughout the table, and, if one excludes the carbon and lead columns 1
(33)

From fig. 3 one sees that this last limitation corresponds to 5 < c < 6.5.

(34)

From the preceding discussion it follows that correlations are effective if N in eq. (31) is less than of say 2, i.e. aC/$,d

< 2

(C/2)0.

(35)

Eq. (35) represents a simple extension of eq. (27), the effect of the correlations being contained in the correction factor C/2. Eq. (35) can be rewritten in terms of

I

0.1

the Lin-

115.2Z,Z,d Tc

T

c=o=

Z, and Z, being the charges of the ion and the lattice atom respytively. d is the atomic distance along the string (in A if the limitation on T is in ev). Assuming C2 = 34 [see eq. (34)], one obtains eventually T < 3.5Z,Z,d a2 For a-particles ture, this gives

(36)

. along Al (llO)-axis

at room tempera-

T-c 23 keV.

which can be written #Ld’

of the ion, by expressing

1

I

I

I

1111111

lo a/0

Fig. 3. Behaviour of the correction factor c (see eq. (29)) versus a /a.

Hence the effect of correlations in the case of cr-partitles along Al (IlO)-axis should appear for energies as low as some ten keV. This prediction is very difficult to check by a computer simulation, because at such low energies #2 = 110 mrad and the whole blocking pattern would extend up to some 200 mrad. To simulate realistically such a situation the directions of the initial velocities of the a-particles should be randomly chosen within a cone of = 250 mrad, so, to obtain meaningful half-aperture statistics, prohibitively long computer times would be needed. Besides, the small-angle approximation would no longer be valid and this would increase the complexity of the mathematical description of the ion-atom interaction. By looking at table 1, one could choose a more favourable situation. Heavy-ion scattering would appear to be attractive in view of the high value of Z,Z,, but it must be remembered that heavy-ion beams are only partially ionized and this in a statistical way. Summarizing, calculations show that correlations in single-crystals are limited to few neighbouring atoms, while many more are needed to build up the whole blocking pattern as we observe both in experiments and in computer simulations. Eq. (36), though being nothing

S. Alliney et al. / Correlated thermal vibrations and crystal blocking

more than an order-of-magnitude estimate, indicates that correlations become important at such low energies that Monte Carlo simulations are in practice not feasible and experiments are very difficult. To end with, the present work lead us to the following main conclusions: a) correlations seem to be ineffective in the majority of blocking phenomena; b) where they are important, as for example in surface back scattering yields, the simple “second neighbour approximation”, eqs. (22)-(23) may be used in place of the exact solution, eqs. (24)-(26). The authors are deeply indebted to Prof. J.U. Andersen for valuable criticism and useful suggestions.

sion: (A.41

Obviously, we set [ .ln = 0 for i = 0 or k = 0. The proof of this can be found in p. 22 of ref. [14]. Eq. (A.4), apart from its practical usefulness, provides a good explanation of the observed effect. Indeed, the last row of M;’ is but a shifted version of the last row of M ;A i, with a correction given by the last term in (A.4). As n increases xi strongly dominates the absolute values Ix,], j=2,3 ,..., n, because of the decay of correlation between very distant objects in the chain.

Appendix B

The Moliere ion-atom Appendix A

4(R)

The stabilization effect observed on the drift coefficients defined in section 2 [see eq. (25)] can be discussed with more precision. We recall that M, is the n X n variance-covariance matrix of the Gaussian process and therefore its elements are of “displacement-type”:

a

= -=ize2

M,x=e,

exp(-&R/a),

being the Thomas-Fermi

@.I)

length, and

{ c$} = {0.1,0.55,0.35}; {pi} = (6, 1.2, 0.3). The corresponding plane

(B.2) string potential

in the transverse

(B.3) (d being the interatomic distance along the string) can be convoluted with a 2-dimensional Gaussian describing atomic vibrations. The result is the thermally averaged string potential V(r)=~Jdb~(p)~exp(-:021r-p12). (B.4)

(‘4.2)

being the first unit vector of the canonical basis of

Introducing x=r/o;,

R”. Remark:

The

first component of X, say xi, is strictly positive. This follows immediately if we apply the Cramer’s rule to the system (A.2). Indeed we have

deW,-d X1 =

cai 1

potential in 3-space is

(Al)

[M.l,j=f,(Ii-jl).

Matrices having this special structure are known in the mathematical literature (see e.g. refs. [14-161) as symmetrical Toeplitz matrices. Among many other properties, such matrices can be inverted with only 0(n2) multiplications, which can be a substantial saving when n is large. Furthermore, any variance-covariance matrix is symmetrical and positive-definite. Such additional property allows an efficient recursive evaluation of the inverses M,‘; n = 1, 2,. . . . The exposition which follows is based mainly on pp. 19-23, ref. [14]. Let the n-vector x denote the first column of M;‘:

e,

15

det(M,)

the dimensionless z=p/a;

quantities

7, = Pi o/a

(B.5)

it can be cast in the form V(OX)=~exp(-$x2)Ca, I

’

(A.3)

where M,_ 1 is the (n - 1) x (n - 1) matrix obtained by removing the last row and column of M,. Notice that through the property (A.l), M,_ 1 is identical to the cofactor of the (1,l) entry of M,. Now, because of the positive-definiteness of M,, both det(M._,) and det(M.) are positive. Theorem: The entries of M,’ can be determined from M;?, once the n-vector x is known by the recur-

X

/

omK,(yiz)I,,(xz) exp(-$z’)z

dz. (B.6)

The second derivative of V at the origin can be found expanding the exponential and I,, functions up to quadratic terms in x: exp( - ix’) &(xz)

= 1 - +x2 + 0(x2),

= 1+ :x2z2+

o(x’z2).

S. Alliney et al. / Correlated thermal vibrations and uystal blocking

16 Then a simple calculation

d2V dr2 ,_O=

yields

-~~~i[gfP(lli)-IfQ(gi)]. i

(B-7) Here P(q)

= Lmd&(2fi)

exp(-q2y),

Q(q)= ~md_vKo(2fi) eXP(-q*Y)

(B-8) (B.9

can be evaluated numerically, while the dimensionless quantities

depend only on the ion-crystal combination through the constant a/o. Therefore F in eq. (31) has the expression F(a/o)

= Coi[$P(vi)

-qfQ(qi)]-

(B.ll)

References [l] V.A. Ryabov, Phys. Status SoIidi 38 (1970) 63; Fir. Tverd. Tela 12 (1970) 2747 @gIish translation: Soviet Phys. Solid State 12 (1971) 22161; Phys. Status SoIidi Sohdi B49 (1972) 467; See also Fir. Tverd. Tela 26 (1984) 941 [English translation: Sov. Phys. Solid State 26 (1984) 575)].

[2] D.P. Jackson and J.H. Barrett, Comput. Phys. Commun. 13 (1977) 157; Phys. Lett. 71A (1979) 359. [3] R.S. Nelson, M.W. Thompson and H. Montgomery, PhiIos. Mag. 1 (1962) 1385, Appendix I. [4] D.P. Jackson, B.M. PowelI and G. Dolling, Phys. Lett. 51A (1975) 87. [5] G.W. Stewart, Introduction to Matrix Computations (Academic Press, New York, 1973). [6] J.H. WiIkinson and C. Reinsch, Linear Algebra (Springer, Berlin, 1971). [7] W.J. Kennedy, Jr. and J.E. Gentle, Statistical Computing (Dekker, New York, 1980). [8] B.W. Lmdgren, Statistical Theory (Collier MacMillan, New York-London 3d ed. (1976) Chap. 10. [9] W. Feller, An Introduction to Probability Theory and its Applications, Vol II (Wiley, New York, 2nd ed., 1971) Chap. III. [lo] I. Stensgaard, L.C. Feldman, and P.J. Silverman, Surf. Sci. 102 (1981) 1. [ll] E. Fuschini, F. MaIaguti, A. Uguzzoni and E. Verondini, Radiat. Eff. 81 (1984) 37. [12] D.V. Morgan (ed.), Channeling; Theory, Observation and Applications (Wiley, New York-London 1973) Section 2.7; J.U. Andersen, private communication. [13] O.S. Oen, Phys. Lett. 19 (1965) 358. [14] G. Heinig and K. Rost, Algebraic Methods for ToephtzIike Matrices and Operators (Birkh%user Verlag, Basel, 1984). [15] I.C. Gohberg and LA. Feldman, Convolution Equations and Projection Methods for their Solutions, in Translations of Mathematical Monographs, Vol. 41 (Amer-Math. Sot., 1974). [16] T. KaiIath, A. Vieira and M. Morf, SIAM Rev. 20 (1978) 106.