Graef M. Introduction to conventional transmission electron microscopy

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

438 Systematic row and zone axis orientations

Table 7.2. Zone axis symmetries for several independent zone axis directions in all

32 crystal point groups [based on Table 4 in [BESR76]].

Point groups 111100110uv0uuw [uvw]

3m 6

4mm1

2mm1

43m 3m 4

432 3m

Point groups 111100uv0 [uvw]

m3 6

2mm1

23 3 2m

Point groups [00.1] 11.01

1.0 [uv.0] [uu.w][u

u.w ][uv.w]

6/mmm 6mm1

2mm1

6m2 3m1

2mm m m

6mm 6mm m1

mm1

622 6m

Point groups [00.1] [uv.0] [uv.w]

6/m61

631

66 m

Point groups [00.1] [11.0] [u

u.w ][uvw]

3m 6

3m 3m 1

32 3m

Point groups [00.1] [uv.w]

33 1

Point groups [001] 100110 [u0w ][uv0] [uuw][uvw]

4/mmm 4mm1

2mm1

42m 4

4mm 4mm m1

422 4m

Point groups [001] [uv0] [uvw]

4/m41

44 m

Point groups [001] 100 [u0w][uv0] [uvw]

mmm 2mm1

2mm1

mm2 2mm m1

222 2m

Point groups [010] [u0w ][uvw]

2/m21

22 m

Point groups [uvw]

7.4 Computation of the exit plane wave function 439

derived from Table 4 in [BESR76], all possible zone axis symmetries for different

zone axis directions for all 32 point groups.

7.4 Computation of the exit plane wave function

In Chapter 10, we will discuss various types of phase contrast microscopy, in

particular high-resolution TEM (HRTEM) and Lorentz microscopy. The main goal

of the HRTEM method is to obtain structural information about the material at the

length scale of the interatomic distances. As for every other observation method that

we have described so far, the resulting experimental high-resolution images must

be compared with theoretical, simulated images. Such a simulation requires two

types of input: the wave function of the electrons at the exit plane of the crystal, and

the transfer function of the microscope. We postpone a discussion of the transfer

function until Chapter 10; in this section, we will describe brieﬂy how the exit plane

wave function can be computed in an efﬁcient way.

As we have seen in Chapter 5, there are three major computational methods which

can be used to obtain this wave function: Bloch waves, the multi-slice method,

and the real-space method. We will assume that the thin foil is oriented so that

a low-index zone axis [uvw] is parallel to the electron beam, and that the foil is

not bent. This means that we are effectively only considering the exact zone axis

orientation, i.e. the point at the center of the bend center. For many crystal structures,

parallel projection drawings along low-index zone axes can be interpreted readily

in terms of the structural elements, such as atom columns, coordination polyhedra,

and so on.

On the

website we provide a single program, BWEW.f90, to compute the exit wave

using the Bloch wave formalism. For a numerical implementation of the multi-

slice formalism we refer the interested reader to the book Advanced Computing

in Electron Microscopy by E.J. Kirkland [Kir98]; this text describes programming

details for the technique and all the underlying theory. The book also provides

extensive C-source code that implements the basic algorithms of the multi-slice

method. The source code for the third method, the real-space method, has been

made available by its authors and can also be downloaded from the

website. For a

critical review of several simulation methods we refer the reader to [SOBS83].

7.4.1 The multi-slice and real-space approaches

We begin with the basic equation (5.35):

∂ψ

∂z



 +



ψ. (7.17)

440 Systematic row and zone axis orientations

If we subdivide the crystal into slices of thickness  normal to the beam direction,

then the formal solution to this equation is given by

ψ(x, y, z = n) =

n−1

...

ψ(x, y, z = 0), (7.18)

with

(x, y) ≡ e



(

+

)

, (7.19)

and

iσ



j

( j−1)

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

This solution is exact in the limit of vanishing slice thickness . In the simplest

case, the slice thickness is equal to the lattice periodicity c along the z-axis and all

slices are identical:

V =

iσ

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

In Chapter 5, we have seen that the exponential of a sum of operators is not simply

equal to the product of the exponentials. If we apply the Zassenhaus theorem to

second order in the slice thickness we ﬁnd for the operator

= e



(

+

)

≈ e







−

(



−



)



;

= e









. (7.22)

The second expression can be shown to be equal to the ﬁrst by means of a simple

Taylor expansion to second order in . The exit wave can then be written as

ψ(x, y, z = n) = e

−













ψ(x, y, z = 0). (7.23)

The ﬁrst and last terms simply propagate the wave function by half a slice thickness.

The terms between the square brackets are identical to the expression (5.42) on

page 319. This expression is correct to second order in the slice thickness.

Now that we have established a theoretical equation for the exit wave we must

answer the following questions:

how many image points (pixels) do we need to have an accurate representation of the exit

wave in both direct and reciprocal space?

how do we compute the phase grating term?

how do we compute the propagator term?

7.4 Computation of the exit plane wave function 441

0 1 2 3 4 5 6 7

01234-1-2-3-4

-1

-2

-3

-4

(a)

(b)

(c)

01234

-4

-1

-2

-3

-4

-1-2-3

Fig. 7.27. Schematic representation of the relation between real-space (a) and reciprocal

space (b) sampling for a square unit cell. In (c) the origin of the reciprocal lattice is shifted

to the lower left-hand corner.

7.4.1.1 The computational grid

The sampling theorem, introduced in Section 3.11.1, states that a continuous func-

tion can be represented by a discrete sample, provided the sampling frequency is

sufﬁciently high and the function is bandwidth limited. A continuous function f (r)

with an inﬁnite domain has a Fourier transform f (q) that also exists on an inﬁnite

domain. If we select a ﬁnite region of the domain of f (r), and sample it discretely,

then its discrete Fourier transform no longer has an inﬁnite domain but becomes

periodic. Typically, the ﬁnite domain will be the projected unit cell of the structure

of interest. If the sampling grid has a sampling interval δ, and we use N points

(or pixels) to cover the side of the unit cell, i.e. a = N δ (see Fig. 7.27a), then the

Nyquist frequency q

is given by

2δ

≡ g

max

This means that the maximum reciprocal lattice point that must be taken into account

in the computation is

N a

∗

. The same argument can be made for any direction in

the real-space unit cell, so that the relevant portion of reciprocal space is a circle,

as shown in Fig. 7.27(b). There is no point in taking more reciprocal lattice points

because the sampling interval δ is too large to correctly represent the contributions

of these higher spatial frequencies. The reciprocal lattice points indicated in gray

in Fig. 7.27(b) should therefore have zero contribution. In numerical work it is

customary to shift the origin of reciprocal space to the lower left-hand corner of the

array, as shown in Fig. 7.27(c). Each corner of the resulting array is then identical as

a consequence of the periodicity of the Fourier transform. The points with frequency

N/2 and −N/2 are also identical.

It is customary to select the number of grid points N

along each axis to be a power

of 2, i.e. N

= 2

, with n

integers. The main reason for this choice is the fact that the

442 Systematic row and zone axis orientations

max

1/ε

max

1/ε

(a) (b)

Fig. 7.28. Graphical representation of the relation between slice thickness and g

max

. The

slice thickness should be sufﬁciently small, so that the corresponding FOLZ layer intersects

the Ewald sphere far from g

max

HRTEM image simulation makes extensive use of the fast Fourier transform (FFT),

introduced in Section 2.5.4, and the traditional FFT algorithms require a power of

2 for all dimensions. This is no longer a necessary restriction, since modern FFT

algorithms can deal with any number N

, not just powers of 2 or other small integers.

The

fftw-package used for the computations in this textbook can handle transforms

of arbitrary N

; when N

is a large prime number the computational speed will be

somewhat slow, but it is rather easy to avoid such numbers. This means that we can

select the grid spacing δ to be the same along both axes, even when the axes have

different lengths.

The number of sampling points is important because it also determines the

maximum slice thickness that can be considered for multi-slice-type simulations.

Figure 7.28 shows schematically how the slice thickness in real-space relates to

reciprocal space. A slice thickness of  creates a periodic structure along the

beam direction which has a FOLZ layer at a distance 1/ from the ZOLZ layer

(e.g. [CVD84a]). To avoid such situations, the point P in Fig. 7.28 must be kept

far enough from the origin so that the excitation error of the g

max

reﬂection will be

signiﬁcantly smaller than 1/:

max





, (7.24)

 

λg

max

8δ

. (7.25)

If the sampling conditions are selected so that this condition is satisﬁed, then HOLZ

effects will be negligible. If the unit cell size along the beam direction is too large

for a given sampling interval δ, then the unit cell must be cut into thinner slices and

the computation time increases, proportional to the number of slices needed.

The value of g

max

to be used for a multi-slice-type computation depends on the

particular atom type(s) present in the unit cell. Heavier atoms scatter more strongly

7.4 Computation of the exit plane wave function 443

and will need more reciprocal lattice points to represent their scattering potentials.

Coene and Van Dyck [CVD84a] have proposed that g

max

should be somewhere

between Z (the atomic number) and Z

1/3

(expressed in nm

−1

), where Z is taken

for the heaviest atom in the unit cell. For pure gold a value of 50 nm

−1

is suggested,

whereas for Cu g

max

≈ 39 nm

−1

. Repeating a computation for several values of g

max

may give an indication as to the minimum value of this parameter (or equivalently,

the maximum sampling interval δ).

7.4.1.2 Computation of the phase grating term

The phase grating exponential can be computed in a number of different ways. We

will only discuss the convolution method introduced by D.F. Lynch [Lyn74]. Other

more advanced methods can be found in [SOBS83] and [dJCVD87]. Consider the

phase grating exponential e



. It is easy to show by means of Taylor expansions

that the following equality holds:



= lim

s→∞



1 +





. (7.26)

Replacing

by iσ V

we have

iσ V

= lim

s→∞

)

1 + i

σ

. (7.27)

This is a slowly converging product and typically we have s ≈ 10

to have a rea-

sonable accuracy of 10

−5

[Lyn74]. If we take the Fourier transform of the s-fold

product, we end up with an s-fold convolution product of the Fourier transforms:

iσ V

δ(h, k) + iσ V



⊗

δ(h, k) + iσ V



⊗···. (7.28)

While the actual product series converges slowly, the convolution product can be

computed in a reasonably short time by taking s to be a power of 2, i.e. s = 2

with typically r = 10–12. If we represent a single term between square brackets

on the right-hand side of the above equation by Q

, then the ﬁrst convolution

produces Q

= Q

⊗ Q

, the second produces Q

= Q

⊗ Q

and so on until

r − 1 convolutions have been carried out.

The coefﬁcients V

are the Fourier coefﬁcients of the lattice potential in the

plane normal to the beam direction. The range of coefﬁcients is again limited

by g

max

, so that Q

is a bandwidth-limited function. Convolutions of bandwidth-

limited functions can be computed by means of the convolution theorem, but require

special care to avoid aliasing artifacts. It is not difﬁcult to show that a doubling

of the number of sampling points from N to 2N in combination with a scaled

Nyquist limit of

will completely eliminate aliasing. This means

that Q

must be bandwidth-limited to two-thirds of the Nyquist frequency for the

444 Systematic row and zone axis orientations

double-sampled function. The resulting convolution product Q

will have amplitude

in all frequencies out to

; application of a mask can then set all frequencies

between

and

of q

to zero before the next convolution is carried out. For more

details we refer to Section 6.7 in [Kir98]. At the end of the r-fold convolution

product, we simply Fourier transform Q

(sampled over N points instead of 2N )

to real-space to obtain the phase grating term. If we represent the mask by M, then

a general step of the algorithm can be summarized as follows:

i+1

≡ F



−1

[MQ

]



. (7.29)

This recipe should be carried out r − 1 times to obtain Q

. The implementation of

this algorithm is left as an exercise for the reader.

The computation of the phase grating term makes the implicit assumption that the

atoms inside a slice are located at the center of the slice. This may not be a realistic

assumption, and more accurate algorithms for the computation of the projected

potential have been proposed in the literature. Van Dyck and Coene deﬁned the so-

called potential eccentricity [VDC84], which explicitly allows the potential inside

a slice to be asymmetric with respect to the center of the slice. This leads to a

slightly modiﬁed solution (7.23), and therefore a slightly slower but more accurate

multi-slice/real-space computation. We refer the interested reader to [KOK87] for

a detailed comparison of various multi-slice implementations.

7.4.1.3 Computation of the propagator term

The difference between the multi-slice and real-space methods lies principally in the

way both methods deal with the computation of the propagator e





. Both methods

attempt to circumvent the direct computation of this term.

The multi-slice method. In the multi-slice method, the propagator is evaluated in recip-

rocal space, using the following relation (derived in Chapter 5):





ψ = F

−1

−iπλq

F[ψ]

. (7.30)

This expression assumes that the incident beam travels along the zone axis. For a tilted

beam the following expression can be derived:





ψ = F

−1

−iπλ(q

+2q·k

)

F[ψ]

, (7.31)

with k

the normal component of the incident wave vector.

The multi-slice algorithm then proceeds as

follows:

ﬁrst, select a slice thickness  and

compute the phase grating (in real-space) and the propagator (in reciprocal space). We will

denote these two factors by p and q, respectively. Then initialize the incident wave φ(0)

in real-space by assigning a unit complex amplitude to each sampling point. Propagate

7.4 Computation of the exit plane wave function 445

this initial wave by half a slice thickness to obtain φ



(0). The wave amplitude at z =  is

then given by

φ() = F

−1

[qF[ pφ



(0)]]. (7.32)

Repeated use of this single step then results in the wave amplitude at multiples of the

slice thickness . The general step is given by

= F

−1

[

n−1

]]

, (7.33)

where p

is the phase grating term for the nth slice. At the end, the wave function must

be propagated for

half a slice again, as described in equation (7.23). As long as

the slice

thickness is constant, the propagator term q must be computed only twice, once for the

half-slice thickness and once for the whole slice thickness. In principle, each slice could

have a different projected potential V

, which is reﬂected in the use of p

. For a perfect

crystal structure with a small unit cell one would typically use a single phase grating term,

but more complex defect arrangements with

a structure which changes along the beam

direction require the use of possibly a large number of different phase grating terms.

There are many implementations of the multi-slice algorithm. The EMS suite of pro-

grams is particularly well suited for simulations that require many different phase grating

terms [Sta87].

†

The C source code accompanying the book Advanced Computing in Elec-

tron

Microscopy

[Kir98]

provides a good starting point for the reader who wishes to write

specialized multi-slice code. A basic multi-slice implementation (multis.f, in Fortran-77)

can be found in [SZ92];

‡

on the website the reader can ﬁnd an ION-routine ms.pro

which implements the essential components

of the

multis.f program (see examples in

Section 7.4.3).

The real-space method. The real-space method [VDC84, CVD84a, CVD84b] assumes

that the value of the wave function in a given point in a slice is inﬂuenced only by

the values of a small neighborhood of points in the previous slice. In other words, the

propagator is considered in real-space (hence the name of the method) and is written as a

small convolution kernel. The explicit derivation of this kernel is somewhat complicated

and we refer the reader to the original paper for all details [CVD84b]. If we represent the

real-space propagator by q

, then the general step of the real-space method is given by

= q

⊗

[

n−1

]

; (7.34)

there are no Fourier transforms, and the propagator kernel typically is a small circular

array (15–20 pixels across). The convolution product is carried out entirely in real-space,

using periodic boundary conditions. The Fortran-77 source code of the real-space method

can be downloaded from the

website; it is available as a compressed UNIX archive

realspace.tar.gz.

†

An on-line version of EMS is available at http://cimesg1.epfl.ch/CIOL/ems.html.

‡

The multis.f program was written by J.M. Zuo; the source code is listed on pages 311–322 in Electron Micro-

diffraction by Spence and Zuo [SZ92].

The author would like to thank Professor D. Van Dyck for permission to distribute the 1990 version of the

real-space source code, written by W. Coene in standard Fortran-77 and modiﬁed by M. Op de Beeck.

446 Systematic row and zone axis orientations

An important advantage of the real-space method is the fact that the ﬁnite propagator

allows for the patching together of many neighboring images into a larger region. This

would be difﬁcult to achieve with the other methods without the need for very large arrays.

An example of the patching method can be found in [CVDVTVL85].

7.4.2 The Bloch wave approach

We know from Section 5.7 that the exit wave function can be expressed as a sum

of Bloch waves, each with its own wave vector k

( j)

and excitation amplitude α

( j)

The individual Bloch waves for both the systematic row case and the zone axis case

have already been analyzed earlier in this chapter. All that remains to be done is

to combine the individual Bloch waves C

( j)

(r) with the proper relative weights to

obtain the total wave function for each grid point. The grid of points at which the exit

wave function is to be evaluated is chosen in the same way as for the multi-slice and

real-space methods. Bloch waves are typically used for smaller unit cells for which

the number of beams is relatively small (of the order of a few hundred beams).

The zone axis Bloch wave program

BWEW.f90 is available from the website; the

program computes the exit wave for a series of foil thicknesses.

7.4.3 Example exit wave simulations

Figure 7.29 shows the exit wave amplitude and phase for the [100] zone axis

orientation of silicon at 200 kV, computed with the

BWEW.f90 program. The com-

putation includes 373 beams (all beams for which |g| < 40 nm

−1

; note that this

is signiﬁcantly smaller than g

max

≈ 235 nm

−1

) and was carried out for 12 thick-

nesses (multiples of 10 nm). The advantage of the Bloch wave method is that the

thickness can be substituted at the end of the computation; this means that for large

numbers of beams the results of the eigenvalue computation can simply be stored

on disk and used again should other thicknesses be required.

It is apparent from Fig. 7.29 that the wave function remains sharply peaked at the

atom positions for the entire thickness range. We have seen earlier in this chapter

that, for the exact zone axis orientation, typically only one Bloch wave is excited.

For the copper example in Figs 7.21 and 7.22, the ﬁrst Bloch wave is centered

on the atom positions and is strongly excited. The same happens for most crystal

structures in which the projected structure consists of columns of identical atoms.

The lowest-order Bloch wave(s) have the highest excitation coefﬁcients in exact

zone axis orientation and the beam electrons “channel” along the atomic column.

From γ

( j)

dispersion surfaces similar to those shown in Fig. 7.20, we ﬁnd that

one (or a few) Bloch waves correspond to bound states, i.e. the electrons in the

beam “fall” into the potential well created by the atom column and remain there

throughout the thickness of the foil.

7.4 Computation of the exit plane wave function 447

102030405060

70 80 90 100 110 120

ampl.

phase

ampl.

phase

Fig. 7.29. Bloch wave exit wave simulation for the

[100] zone axis orientation of silicon

(373 beams, 200 kV, 256 × 256 pixels). Phase and amplitude are shown for thicknesses

from 10 to 120 nm. The white spots in the amplitude at 10 nm coincide with the projected

atom positions for a single unit cell.

For an isolated column of atoms the bound states are similar to the radially sym-

metrical states of the hydrogen atom, i.e. similar to s states. We can use the

ms.pro

ION-program on the website to study this channeling behavior. The ms.pro pro-

gram takes a limited number of input parameters. The crystal structure is tetragonal,

with a = 0.4 nm, and c is variable. Two identical atoms can be placed anywhere

in the unit cell. The cell is sampled at 64 × 64 pixels. The user can select the slice

thickness , the number of slices, the acceleration voltage E, the value of g

max

and the c lattice parameter. By varying these parameters, the reader can obtain a

good understanding of how channeling works. An example simulation is shown in

Fig. 7.30. The unit cell is shown in Fig. 7.30(a) and consists of two W atoms at

0.142 nm spacing. The unit cell thickness is c = 0.2 nm and a slice thickness of

 = 0.05 nm is used (200 kV). The projected potential and the propagator (origin

shifted to lower left-hand corner) are shown in Figs 7.30(b) and (c), respectively.

Figure 7.30(d) shows the amplitude (top row) and sin(φ) (bottom row) for

10 thicknesses. There are three different intensity features, indicated by small ar-

rows in the 40 nm amplitude image, which occur for all images, though not always

all at the same time. The smallest feature is a bright spot at the atom location. The

second feature is ring-shaped and surrounds the atom position. The third feature

is more diffuse and surrounds each atom spot at a larger distance; this feature is