Graef M. Introduction to conventional transmission electron microscopy

Подождите немного. Документ загружается.

308 Dynamical electron scattering in perfect crystals

When we substitute this back into equation (5.8) we ﬁnd:





dψ

− i2π

− (k

+ g)

2|k

+ g|cos α



2πi(k

+g)·r



iπ







g−g



+ iU



g−g





+ g



|cos α



2πi(k

+g)·r

. (5.10)

Let us brieﬂy summarize what we have done so far: starting from equation (5.5)

we have introduced a plane wave expansion for the wave function, based on Bragg’s

law; we have introduced a Fourier expansion for the complex lattice potential

U + iU



; we have applied the high-energy approximation; and we have assumed

that we have a perfect crystal. Equation (5.10) must be valid for every location

r, therefore the equality must be valid for each of the individual terms in the

summation:

dψ

− i2π

− (k

+ g)

2|k

+ g|cos α

= iπ







g−g



+ iU



g−g





+ g



|cos α



. (5.11)

This set of equations, one for each g, states that the change of the amplitude of

the diffracted wave ψ

with depth z in the crystal (ﬁrst term on the left-hand

side) is determined by: (1) the current value of the amplitude (second term on the

left-hand side) and (2) the amplitudes of all other waves ψ



(summation on the

right-hand side). The strength of the interaction between the waves is determined

by the Fourier coefﬁcients of the complex electrostatic lattice potential. In other

words, the amplitude of any given diffracted wave is determined by the amplitudes

of all other waves, and the electrostatic lattice potential determines how strong these

interactions or couplings are.

We have seen in Chapter 2 (on page 119) that for a ﬁnite crystal, diffraction may

occur even when the reciprocal lattice point is not precisely located on the Ewald

sphere. We have introduced a geometrical parameter s

, the excitation error or devi-

ation parameter, which measures the distance between the reciprocal lattice point

and the Ewald sphere, along the direction of the foil normal (or, equivalently, along

the relrod). Replacing K

by k

to allow for refraction, we recognize the expression

for s

in the second term of the left-hand side of equation (5.11), and we write:

dψ

− 2π is

= iπ







g−g



+ iU



g−g





+ g



|cos α



. (5.12)

From our discussion of the relation between reciprocal space and the objective

lens back focal plane in Chapter 2, we remember that we can move the reciprocal

lattice points through the Ewald sphere (i.e. change their distance to the Ewald

5.3 General derivation of the Darwin–Howie–Whelan equations 309

sphere) by changing the orientation of the crystal, or by tilting the incident beam.

Therefore, the parameters s

are experimental variables, and their values can be

set by the microscope operator. One should note that the various excitation errors

are not independent of each other; in fact, when the excitation errors of two

non-collinear vectors g and g



are selected (by tilting the crystal in a particular

orientation), all other excitation errors are ﬁxed for a given electron wavelength.

Since the ﬁrst term in equation (5.12) has the dimensions of reciprocal length,

the right-hand side must have the same dimensions. Therefore, we introduce two

new parameters with the dimensions of length: the extinction distance ξ

and the

absorption length ξ



. The Fourier coefﬁcients of the real and imaginary part of the

electrostatic lattice potential can be written in terms of their modulus and phase

angle, as follows:

=|U

iθ

;



=|U



iθ



For a centrosymmetric crystal (which must have real Fourier coefﬁcients), the phase

factor θ

can only take on the values 0 and π, corresponding to positive and negative

values for the Fourier coefﬁcients U

. There is no a priori reason why, for a non-

centrosymmetric crystal, the phase factors of U

and U



should be the same. We

will address this issue in some detail in Chapter 6.

We can use the moduli of these Fourier coefﬁcients to deﬁne the extinction

distance and the absorption length:

(extinction distance)

−1

→

≡

+ g|cos α

; (5.13)

(absorption length)

−1

→



≡



+ g|cos α

. (5.14)

Noting that, to a very good approximation, we have |k

+ g|≈|k

+ g



|≈|k

|,so

we can write

g−g



+ iU



g−g



)

|cos α

iθ

g−g



+ i

iθ



g−g



For notational purposes it is convenient, following [GBA66], to introduce the fol-

lowing complex quantity:

≡

+ i

iβ



, (5.15)

with

= θ



− θ

. (5.16)

310 Dynamical electron scattering in perfect crystals

Note that, since U

= 0 (and hence ξ

=∞), and U



is real, we also have

+ i

iβ



. (5.17)

Finally, we can rewrite equation (5.12) as

dψ

− 2π is

= iπ





iθ

g−g



g−g



. (5.18)

Equations (5.18) are known as the multibeam Darwin–Howie–Whelan equations

of dynamical electron diffraction theory [HW61]. They are valid for an arbitrary

number of diffracted beams, and they describe how the amplitude of a diffracted

beam changes with depth in the crystal. These changes depend on how well the

Bragg condition is satisﬁed, and on the strengths of the interactions

with the other

diffracted beams.

Using the Dirac delta-function we may rewrite the Darwin–Howie–Whelan (or

DHW) equations as follows:

dψ

= iπ







δ(g − g



) +

iθ

g−g



g−g







. (5.19)

This form of the system of equations will prove to be useful in the following section,

where we convert the dynamical scattering equations to a single matrix equation.

The numbers (q

g−g



)

−1

are related to the probability that an electron will be scattered

from the beam g



into the beam g. In particular, the number (q

g−0

)

−1

= (q

)

−1

refers to scattering from the incident beam into the diffracted beam g. The number

0−g

)

−1

= (q

−g

)

−1

describes the reverse process, i.e. scattering from the diffracted

beam g back into the incident beam direction. In the absence of absorption (ξ



=∞),

and for a centrosymmetric crystal, we have (q

−g

)

−1

= (q

)

−1

In Section 6.5.1, we will take a closer look at the numerical computation of

the extinction distances and absorption lengths for a given crystal structure. For

now, we assume that we know all relevant parameters, and then the only “free”

parameters are the crystal orientation, through the excitation errors s

, the crystal

thickness z

, and the electron wavelength λ.

In principle, we could now use any differential equation solver algorithm to

compute the solutions to these equations, for an arbitrary crystal. At this point, it is

far more useful, however, to rewrite the equations in a number of different forms,

so that we may discover other alternative (and possibly faster) ways to solve the

dynamical scattering equations.

5.4 Formal solution of the DHW multibeam equations 311

5.4 Formal solution of the DHW multibeam equations

We begin with equations (5.19). Instead of labeling every diffracted beam with the

corresponding reciprocal lattice vector g, we will simply number all the beams from

n = 1toN. It is common practice, but not necessary, to label the transmitted beam

with n = 1. The DHW system of equations then becomes

dψ

= iπ









δ(n − n



) +

iθ

n−n



n−n







n = 1,... ,N. (5.20)

Using equation (5.17), these equations can be simpliﬁed substantially by the fol-

lowing substitution:

= S

i(θ

+π z/q

)

= S

iθ

−π z/ξ



. (5.21)

The intensity of the nth scattered beam is given by

=|ψ

= e

−2π z/ξ



i.e. an exponentially damped function. The parameter ξ



is commonly known as

the normal absorption length, and it is common to all beams, including the trans-

mitted beam. For a sample thickness equal to the normal absorption length, the

exponential factor equals e

−2π

= 0.001 8674, which means that only just under

0.2% of the electrons make it through a sample of this thickness. Since this factor is

common to all scattered beams, it is appropriate at this point to remove it from the

equations. Substitution of equation (5.21) into the DHW equations above, results

in the following system of equations:

= iπ









δ(n − n



) +

1 − δ(n − n



)

n−n







. (5.22)

These equations cannot be solved analytically, and we could use, say, the fourth-

order Runge–Kutta method to obtain a numerical solution; we will discuss numeri-

cal solution methods in Chapter 7. For now, we will proceed analytically and rewrite

the system of equations in matrix form. Consider the square matrix A, deﬁned by

= 2πs

;



n−n



n = n





(5.23)

In other words, the matrix A contains along its diagonal the excitation errors of

all the beams contributing to the scattering process. The off-diagonal elements

describe the interactions between the various beams through the complex numbers

n−n



. There is, hence, a clean separation between the diffraction geometry and the

312 Dynamical electron scattering in perfect crystals

coupling between the diffracted beams. We will see later on that we can make good

use of this separation (Section 5.5).

If we write the diffracted amplitudes S

as components of a column vector S, the

DHW equations can be rewritten in matrix form as:

†

= iAS. (5.24)

This equation can be solved formally for a crystal of thickness z

S(z

) = e

iAz

S(0) ≡ SS(0), (5.25)

as can be veriﬁed easily by substitution in the above equation. The matrix A is

the so-called crystal transfer matrix or structure matrix, and S is known as the

scattering matrix [Stu62].

The exponential of a matrix can be deﬁned from the standard Taylor expansion

for the exponential function (e.g. [AS77]):

∞



n=0

= 1 +

x +

+···+

+···. (5.26)

For the scattering matrix S we can write:

S = e

iAz

= 1 + i

−

+···+(i)

+···, (5.27)

where 1 is the identity matrix with the same number of rows and columns as A.We

will discuss an algorithm for the numerical computation of S in Chapter 7.

At this point, it is useful to take a closer look at the two new matrices S and

A. From the theory of differential equations we know that the general solution to

a system of N ﬁrst-order coupled differential equations can be written as a linear

combination of N special solutions. The special solutions can be obtained by stating

N independent boundary conditions at the entrance surface of the crystal. The most

convenient choice for the independent conditions is

S(0) =













;













;













; ... ;













Solving the DHW equations for each of those initial conditions results in N column

vectors S

), n = 1,... ,N , which can be placed in an N × N square matrix S.

†

The order of the entries in S is irrelevant, as long as the same order is used for the matrix A.

5.5 Slice methods 313

If a general incident amplitude is described by the column vector S(0), then the

general solution at the exit plane of the crystal is given by

S(z

) = SS(0).

Hence, the scattering matrix S can be computed by explicitly solving the DHW

equations for N independent initial conditions. However, the exponential expression

for the matrix S in terms of the structure matrix A shows that S depends in a non-

linear way on the crystal thickness z

. This means that the previous equation should

actually read as

S(z

) = S(z

)S(0).

The scattering matrix must, therefore, be recalculated for every crystal thickness.

In Section 5.7 on Bloch waves we will describe how this matrix can be written as

the product of three matrices, with only one of them being dependent on the crystal

thickness.

In the absence of absorption, the sum of the intensities of the beams leaving the

crystal must equal the incident intensity, and the matrix relation above expresses a

unitary transformation, i.e. a transformation that leaves the norm of the vector S

unchanged. Consequently, the dynamical scattering process can be regarded as a

rotation in an N -dimensional space, and S is a hermitian matrix. When absorption

occurs, the transformation is no longer unitary, and the total intensity decreases

with increasing crystal thickness.

5.5 Slice methods

The exponential of a matrix, as deﬁned in equation (5.27), converges after a small

number of terms if the crystal thickness z

is sufﬁciently small. We can make use

of the properties of exponential functions and rewrite relation (5.25) as follows:

S(z

) = S(z

)S(0);

= S(n)S(0);

[

S()

]

S(0).

We have made use of the fact that (x)

= (x

)

. The total crystal thickness equals

= n, and, if we take  sufﬁciently small, then the series expansion for S()

converges rapidly; the scattering matrices for thicknesses that are a multiple of 

are obtained by repeated multiplication of S() with itself, and matrix multiplication

can be carried out efﬁciently on fast processors.

We have subdivided the crystal into a stack of slices of equal thickness ,as

illustrated in Fig. 5.1(a). Now it becomes apparent how we can use the matrix

314 Dynamical electron scattering in perfect crystals

z=n

ε z=n(ε

+ε

)

(a) (b)

Fig. 5.1. (a) Illustration of a crystal divided into parallel slices of equal thickness; (b) more

complex stack of slices, consisting of two different materials. The boundary conditions

at each of the interfaces must be properly taken into account in the computation of the

multibeam scattering process.

formulation to our advantage. Consider a crystal consisting of an alternating se-

quence of layers of different thickness, normal to the electron beam. The layers

could consist of different materials, or identical material with a different orien-

tation. To compute the amplitudes at the exit plane we would then only need to

compute two scattering matrices S

(

) and S

(

), and multiply them together as

many times as needed to obtain the ﬁnal amplitudes:

S(z

) =

[

(

)

]

S(0), (5.28)

where n is the number of bilayers in the stack. The matrix formalism turns this type

of computation into a very simple sequence of matrix multiplications;

†

it would

be very difﬁcult to employ analytical methods to obtain the same results. For a

practical implementation of the matrix method we refer to Chapter 7.

The structure matrix A can be written as the sum of two matrices, one diagonal,

and one with zeros along the diagonal (recall that the term 1/q

was separated out

in (5.21) and therefore does not contribute to V):

= 2πs

;

(5.29)



n−n



n = n



The formal solution to the DHW equations is given by

S(z

) =

i(T +V)

S(0). (5.30)

Matrix multiplication is in general non-commutative and therefore we cannot write:

i(T +V)

= e

iT 

iV

← this is incorrect! (5.31)

It is easy to verify that this statement is incorrect by writing down and comparing

the ﬁrst few terms in the Taylor expansions of left- and right-hand sides. Instead,

†

One must take care to specify the proper boundary conditions at each interface plane.

5.6 The direct space multi-beam equations 315

one has to use the so-called Zassenhaus theorem (from the theory of Lie algebras)

which states the following [Suz77]:

For every pair of square matrices A and B of dimension N there

exists a series of unique matrices C

of degree m in A and B such

that

A+B

and the product series converges. The ﬁrst three elements are

given by C

= A, C

= B and C

=−

[A, B], where [a, b] indicates

the commutator ab − ba.

Hence we can write

i(T +V)

= e

iT 

× e

iV

× e

−

(TV−VT )

×···.

For small values of the slice thickness , higher-order terms can be ignored and the

product may be truncated after the second term. We ﬁnd for the approximate wave

function components at the exit plane of a crystal of thickness z

= n:

S(z

) = e

iT 

iV

···e

iT 

iV

iT 

iV

S(0), (5.32)

where the product of the pair of exponentials is executed n times. It appears that

we have not really gained anything by splitting the structure matrix into two parts.

However, we know that the matrix T is a diagonal matrix; it is straightforward

to show that the exponential of a diagonal matrix is equal to the diagonal matrix

containing the exponentials of the original diagonal elements. Therefore we have



iT 





= e

2πis





δ(n − n



), (5.33)

and this is rather straightforward to compute numerically. The exponential of the

second matrix, V, is much more demanding computationally, and we will next

derive a method to convert the exponential of V into an expression that is much

easier to compute.

5.6 The direct space multi-beam equations

We start again from the basic Schr¨odinger equation (5.5) on page 305. This time, we

will not make use of Fourier expansions at all, and describe the dynamical scattering

process entirely in direct space, following [VDC84]. The incident electron is repre-

sented by a plane wave with free-space wave vector K

, and after passing through the

316 Dynamical electron scattering in perfect crystals

crystal we assume that the wave function can be written as a modulated plane wave:

(r) = ψ(r)e

2πiK

·r

. (5.34)

Substitution into equation (5.5) results in (dropping the argument r):

ψ + i4π k

·∇ψ =−4π

(U + iU



)ψ,

where we have once again included refraction in the wave vector k

. Since the

z-component of the wave vector k

is much longer than the x and y components,

we can separate the terms in z from the others, and we ﬁnd

∂

∂z

+ i4π k

∂ψ

∂z

+ 

ψ + i4π k

·∇

ψ =−4π

(U + iU



)ψ.

The differential operators 

and ∇

are two-dimensional operators, and contain

only x and y derivatives. Once again we can make use of the high-energy approxima-

tion (see page 306) and ignore the second-order derivative along the beam direction.

That this is again a valid approximation has been shown in [VDC84]. This results in

∂ψ

∂z

4πk

$



+ i4π k

·∇



ψ + 4π

(U + iU



)ψ

Next, we introduce a shorthand notation for the terms in this equation:

 ≡

4πk





+ i4π k

·∇



;

V ≡

iπ

(U + iU



The resulting equation is

∂ψ

∂z



 +



ψ. (5.35)

This equation has the same form as the reciprocal space equation (5.24) after

splitting the crystal transfer matrix A into its two component matrices:

= i(T + V)S.

Again we have separated the geometry of the problem (

) from the electron

interactions (

V ). Let us take a closer look at the two components of this equation.

5.6.1 The phase grating equation

Let us assume ﬁrst that

 = 0, which means that there is no lateral scattering of

the electron; the main equation then reduces to

∂ψ

∂z

V ψ. (5.36)

5.6 The direct space multi-beam equations 317

For a crystal of thickness z

, this equation is integrated readily to give

ψ(x, y, z

) = e

;

V (x,y,z)dz

where we have written the explicit coordinate dependence. The argument of the

exponential can be rewritten in more familiar terms as follows:

V =

iπ

2me

= iσ V

where we have used the deﬁnition of the interaction constant σ (2.35) and k

= 1/λ.

The solution to the ﬁrst component equation is then rewritten as

ψ(x, y, z

) = e

iσ

;

(x ,y,z)dz

≡ e

iσ V

;

is commonly known as the projected potential and is deﬁned by

≡

(x, y, z)dz. (5.37)

Equation (5.36) can be interpreted easily if we consider a thin crystal with z

= ;

in the absence of absorption, the projected potential V

changes only the phase of the

electron wave function, not the amplitude, and hence the thin slice of crystal acts as

a phase shifting layer or a phase grating. If the phase change is small (much smaller

than a radian), then the exponential can be approximated by its Taylor expansion

iσ V

≈ 1 + iσ V

+···

and the crystal is referred to as a weak phase object . We will see in Chapter 10

that ferromagnetic thin foils give rise to large phase shifts, and are hence strong

phase objects.

5.6.2 The propagator equation

Next, we assume that

V = 0 (this is similar to the empty crystal approximation)

and we solve the following equation:

∂ψ

∂z

ψ. (5.38)

This equation is ﬁrst order in z, and second order in x and y, hence it belongs

to the class of diffusion equations, with z taking on the role of time t. The main

difference between this equation and a standard diffusion equation is the fact that



is a complex differential operator. If the incident (refraction-corrected) wave vector

does not have x and y components, i.e. if k

= 0, then the differential operator is

purely imaginary. The formal solution for a slice of thickness  is, as before, given by

ψ(x, y,) = e





ψ(x, y, 0).