Image Calculations in High-Resolution Electron Microscopy: Problems, Progress, and Prospects

ADVANCES IN ELECTRONICS AND ELECTRON PHYSICS. VOL . 65

Image Calculations in High-Resolution Electron Microscopy: Problems. Progress. and Prospects D. VAN DYCK University of Antwerp Rijksuniversitair Centrum Antwerp. Belgium

I. Introduction . . . . . . . . . . . . . . . . . . . . . . I1. Principles of Image Formation . . . . . . . . . . . . . . . 111. Present Status of Many-Beam Electron-Diffraction Calculations . . . A . History . . . . . . . . . . . . . . . . . . . . . . . B. The Multislice Method . . . . . . . . . . . . . . . . . C. Computations . . . . . . . . . . . . . . . . . . . . D . Use of the Fast Fourier Transform . . . . . . . . . . . . E. Special-Purpose Hardware . . . . . . . . . . . . . . . F . Discussion . . . . . . . . . . . . . . . . . . . . . . IV. Systematic Approach to the Electron-Diffraction Problem . . . . . A . General Theory . . . . . . . . . . . . . . . . . . . . B. High-Energy Approximation . . . . . . . . . . . . . . C. Validity of the High-Energy Approximation . . . . . . . . . D. Discussion . . . . . . . . . . . . . . . . . . . . . . V. Study of Existing Many-Beam Diffraction Formulations . . . . . . A . Differential Equations . . . . . . . . . . . . . . . . . . B. MatrixMethods . . . . . . . . . . . . . . . . . . . . C. The Iterative Method . . . . . . . . . . . . . . . . . . D. Direct Integrating Methods . . . . . . . . . . . . . . . E. SliceMethods . . . . . . . . . . . . . . . . . . . . F . Discussion . . . . . . . . . . . . . . . . . . . . . . VI. A New Formulation: The Real-Space Method . . . . . . . . . . A . Principles . . . . . . . . . . . . . . . . . . . . . . B . Comparison with Other Slice Methods . . . . . . . . . . . C. Numerical Procedure . . . . . . . . . . . . . . . . . . D . Analysis of Input Parameters . . . . . . . . . . . . . . E. Discussion . . . . . . . . . . . . . . . . . . . . . . VII . Recent Developments, Unsolved Problems, and Prospects for the Future A . The Real-Space Patching Technique . . . . . . . . . . . . B. Considerations Concerning the Periodic-Continuation Method . . C. Atomic Column Approximation . . . . . . . . . . . . . D. Absorption. Inelastic Scattering, and Thermal Diffuse Scattering . . E . The Reverse Problem: Direct Structure Retrieval . . . . . . . F . Upper Layer Lines and Beam Tilt . . . . . . . . . . . . . G . Application to General Quantum-Mechanical Problems . . . . . H. Conclusion and Prospects for the Future . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . .

. . . 296 . . . . 297 . . . . 299 299

. . .

. . .

. . . .

. . . .

.

. . . . .

. . . .

. . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . .

. . . . .

.

.

. .

. . . . .

. . . . . . . . .

. . . .

301 303 304 305 306 308 308 309 310 313 316 316 319 320 321 322 324 324 324 327 332 334 336 337 337 341 343 344 346 350 351 352 353

295 Copyright Q 1985 by Academic Press. Inc. All rights of reproduction in any form reserved. ISBN 0-12-014665-7

296

D. VAN DYCK

I. INTRODUCTION

High-resolution electron microscopy is definitely breaking through in both materials science and biology. It is now quite well established that the interpretation of the high-resolution electron micrographs in terms of the microstructure of the specimen is more reliable if the experimental images are compared with computer simulations in which the dynamic electron diffraction in the specimen as well as the image formation in the electron microscope are simulated for all plausible trial structures. Today, most of the specialized electron microscopic laboratories in the world have their own image calculation programs. There is now a growing tendency to perform image simulations interactively and closer to the electron microscope, e.g., on a small laboratory computer with display facilities. Interactive programs have been written and complete, dedicated, simulation hard- and software is available. However, the early expectations regarding the simulation technique, which ran rather high after its first successes in the early 1970s, have yet to be fulfilled. Indeed, since its introduction, the simulation technique has hardly undergone any further development. Simulation is sometimes used more as a gadget than as a tool. For instance, the comparison between experimental and simulated images is often done visually, i.e., qualitatively, without even taking account of the characteristics of the photographic material. Furthermore, a number of experimental parameters, such as specimen thickness, absolute value of defocus, and angle of beam convergence, which are needed as input parameters in the calculations, are often not determined, determined with insufficient accuracy, or even treated as adjustable parameters (Howie 1983). On the other hand, the number of experimental and simulated images of the same specimen, such as through-focus or thickness series, is often insufficient. Furthermore, the normalization of the total diffraction intensity is still often used as a criterion for the selection of input parameters, such as slice thickness and number of beams. However, while this criterion is questionable, other simple criteria are still lacking. It is apparent that some of the image simulations reported in the literature are performed using input parameters that obviously lead to erroneous simulations. As a result, the reliability of the high-resolution image simulation technique for the determination of the microstructure is still relatively poor and can only be used with some confidence in cases where the structure is already sufficiently known and slight modifications of the structure are tested by trial and error. Hence the use of the technique is seldom convincing. On the other hand, the simulation technique would be much more powerful if the contrast of the experimental images were to be measured quantitatively, e.g., using microdensitometry or on-line image capturing and digitalization. In that

CALCULATIONS IN

HIGH-RESOLUTION ELECTRON MICROSCOPY 297

case, the experimental and simulated images could be compared by using least-squares techniques, and iterative image calculations could lead to structural refinements. However, at this stage the computational procedures are still too time-consuming and expensive. For example, a typical image calculation performed on a large mainframe computer can cost somewhere between $100 and $1000 of computer time. Hence, least-squares matching or applications to more complicated systems, such as extended defects in crystals, or to convergent beam electron diffraction are often restricted. Although most of the current many-beam diffraction computations are based on the multislice theory, which dates from 1957, this does not mean that the theory itself is completely satisfactory, but rather that it has not been tested to its extremes. Rigorous approaches however are scarce. Furthermore, there are still a number of problems that are not quite solved or understood, or even not recognized, such as the influence of higherorder Laue zones, the effect of absorption, the influence of inelastically scattered electrons, the presence of anomalous fringes in some high-resolution images [e.g., (001) rutile], and the validity of the approximations, such as the high-energy approximation. There is also some hope that in the future, the reverse problem can be solved, i.e., how to retrieve the specimen structure directly from the electron micrographs without trial-and-error procedures. For this purpose the dynamic diffraction must be “reversed,” which seems to be a tedious problem in which questions of uniqueness are expected to occur. Nowadays, the situation of high-resolution electron microscopy with image simulation is comparable to the situation in the early days of x-ray diffraction, when the power of the new technique was recognized but the applications were limited by theoretical and computational problems. In this work, we sketch the principles of image simulation, compare different computational methods, and discuss recent developments and possible future prospects. 11. PRINCIPLES OF IMAGE FORMATION

Since the famous paper of Scherzer (1949), the principles of image formation in the electron microscope have.been well established and will be briefly sketched below (for a more general treatment see, e.g., Van Dyck, 1978a or Spence, 1981). Calling q(x, y ) the wave function at the exit plane of the object, the diffraction pattern, i.e., the wave function in the focal plane of the objective lens is, then, according to Abbe’s theory, given by its Fourier transform

F u, v

CdX,

Y)l

298

D. VAN DYCK

where u, v are the reciprocal coordinates that characterize the diffracted beams. In the second part of the image formation, before undergoing a second Fourier transform, the diffracted beams are selected for imaging by an aperture function A(u, v) and suffer a phase shift ~ ( uu,) due to the spherical aberration of the objective lens and the defocusing of the microscope. The aperture function A(u, v) can eventually be modulated by a damping factor in order to account, in a first approximation, for the effects of beam divergence and chromatic aberration (see, e.g., Spence, 1981). Finally, the wave function in the image plane is given by

In the hypothetical case of an ideal electron microscope without aberrationphase shift, and aperture, the image intensity I+(x, y)I2 of a pure phase object q(x, y ) = exp[i+(x, y)] would show no contrast at all. For the study of phase objects, such as thin biological or material foils, the inevitable occurrence of a phase shift due to spherical aberration, which is inherent in the use of magnetic lenses, is in fact a fortunate coincidence. Indeed, as shown by Scherzer (1949), by proper defocusing of the objective lens (called the “Scherzer” focus), the phase shift of the spherical aberration can be compensated in such a way that the central beam is approximately retarded over n/2 in phase with respect to the major part of the other diffracted beams. Hence, the electron microscope can be considered as a kind of phase contrast microscope with the image intensity roughly proportional to the projected electrostatic potential of the specimen. However, for thicker specimens and other focus values, the interpretation of the images is often less straightforward and must rely on computer simulations. By using fast Fourier transform routines, the computer simulation of the imaging process using Eq. (1) does not require stringent computer demands. Typical calculation times are on the order of seconds. The image intensity I $(x, y)I2 can be displayed as a halftone picture by using several techniques such as overprinting, photowriting, electrostatic printing, storage CRT, and television screen, depending whether the simulations are performed in batch or in real time. The most laborious part of the simulation process is the computation of the electron wave function q(x, y) at the exit plane of the specimen. Due to the complexity of the electron-object interaction, which is a fully dynamic diffraction process in which many beams can be involved, the computer demands and costs can be so high that they put serious restrictions on the use of the simulation technique. This is, for instance, the case when complicated objects such as extended defects in crystals are involved, in which sometimes 10,000 or more beams take part in the diffraction process, or when the convergence of the incident beam is treated rigorously by computing

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

299

separately the images corresponding with all different directions of incidence and adding the image intensities afterwards. In this work we focus our attention primarily on the problems emerging from the computation of q(x, y) and discuss the most recent developments in the theory and the computer algorithms for the simulation of the many-beam electron diffraction process. 111. PRESENT STATUS OF MANY-BEAM ELECTRON-DIFFRACTION CALCULATIONS A . History

The theory of dynamic elastic electron scattering in a perfect crystal has a long, involved, and even devious, history (Moodie, 1981). The first, and still one of the most complete treatments, was already provided by Bethe as early as 1928. Inspired by Ewald’s work on the dynamic scattering of electromagnetic waves and Born’s theory of quantum scattering, he treated the electronscattering problem starting from the Schrodinger equation by using a Fourier expansion for the periodic crystal potential and a plane wave expansion for the electron wave function equivalent to the concept of Bloch waves. In this way he obtained a set of coupled linear dispersion equations for the plane wave expansion coefficients that can be put into matrix form and, in the approximation of forward scattering, can be reformulated as an eigenvalue problem. However, since the dimensions of the dynamic matrix are proportional to the number of plane waves under consideration, this formulation is only appropriate when the number of beams involved is relatively small. Since Bethe’s original approach, (Bethe, 1928) several many-beam electron-diffraction theories have been developed, some of which were adaptations of existing methods of band theory. Before the famous Kyoto Conference (1961), many of the interconnections between the various formulations were not known or obscure (Moodie, 1981). However, they all turned out to be reformulations of the original Bethe theory. In the early 1950s, experiments showed that many-beam effects could be significant and, with increasing accelerating voltage, hundreds and even thousands of simultaneous reflections could occur. However until 1955, the most reliable many-beam diffraction theory available was still the kinematic theory (Sanders and Goodmann, 1981). Hence, the need for alternative formulations became urgent. In 1950, Sturkey presented his scattering matrix theory (Sturkey, 1962).In 1955, Fujimoto, using the work of Niehrs and Wagner (1959), showed that the scattering matrix theory could be derived directly from the Bethe dispersion equations. Howie and Whelan (1961) used a totally different approach.

300

D. VAN DYCK

They started from Darwin’s treatment of x-ray scattering and obtained a system of coupled first-order differential equations known as the system of Howie and Whelan, which is also equivalent to the eigenvalue method derived from Bethe’s dispersion equations. This approach turned out to be extremely valuable for the calculation of the contrast of electron microscopic images in two-beam or three-beam situations. In combination with the column approximation, it has rendered immense service to electron microscopic studies of crystal defects. At the same time, Tournarie (1960, 1961, 1962) started by Fourier transforming the Schrodinger equation using a semireciprocal notation, and, in a short and clear treatment, he derived a similar set of coupled first-order differential equations. However, none of these treatments was particularly suitable for the description and calculation of many-beam diffraction situations. The first (and for a long time only) method capable of treating many-beam diffraction was the well-known multislice method of Cowley and Moodie (1957), the starting point of which differed completely from the other existing theories. Originally, the concept of the multislice theory was inspired by previous work of Cowley and Moodie (1951) on Fourier images in optics. Here a crystal is considered as a succession of two-dimensional slices in the direction of the incident beam. Each slice acts as a pure phase grating. According to the laws of physical optics, the propagation of electrons from one slice to the next is governed by the Fresnel formula. The whole diffraction process can thus be described as a succession of phase-grating transmissions and propagations (see Section 111). The multislice theory gained considerable interest in the early 1970s when the Australian group (see Alpress et al., 1972) succeeded in implementing the multislice theory for computer simulation (Lynch and O’Keefe, 1972; OKeefe, 1973) of high-resolution images. The agreement between the famous experimental images of Iijima (complex oxides) and the first computer simulations was so striking that it could convince the most sceptical theoretician. However, because the multislice theory was in principle derived from a physico-optical formalism, its relationship to the other existing diffraction theories has been rather obscure for some time. Goodman and Moodie (1974) in an excellent paper elucidated the interrelationship between the various existing theories and the multislice theory, starting from Schrodinger’s equation. Independently and using the totally different but equivalent approach of the Feynman path integral formalism of quantum mechanics, Van Dyck (1975) showed the equivalence between the multislice formulas and the system of Howie and Whelan. Furthermore, the path integral concept provides much complementary insight into the diffraction process. The relationship between the multislice and the Darwin methods was established by Spence and Whelan (1975).

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

301

It seems surprising that, since that time, only relatively little progress has been made in the field of many-beam diffraction calculations. The major improvement in multislice calculations was introduced by Ishizuka and Uyeda (1977) by making use of the fast Fourier transforms for the calculation of the convolution products that describe the electron propagation between the slices. Most of the current image simulation programs are still based on this method; other attempts at speeding up the calculations were proposed by Van Dyck (1979). A possible promising alternative, the real-space (RS) method, was first introduced by Van Dyck in 1980. In the beginning, this method was subject to criticism (Self, 1982; Self et al., 1983; Kilaas and Gronsky, 1983) and some computational difficulties needed to be resolved. Recently (Van Dyck and Coene, 1984; Coene and Van Dyck, 1984a,b), the last difficulties were resolved and the method seems to be competitive, especially when very complicated systems are involved. In the following we discuss the different many-beam diffraction theories with respect to their numerical performance.

B. The Multislice Method Although the multislice formula can be derived from quantum-mechanical principles (Goodman and Moodie, 1974; Van Dyck, 1975), we follow a simplified version of the more intuitive original optical approach (Cowley and Moodie, 1957). Part of the discussion has already been presented elsewhere (Van Dyck, 1980); a more rigorous treatment is given in Section V. Consider a plane wave, incident on a thin specimen foil and nearly perpendicular to the incident beam direction z. If the specimen is sufficiently thin, we can assume the electron to move approximately parallel to z so that the specimen acts as a pure phase object with transmission function

with

where E is accelerating voltage, A = l/k is the electron wavelength, and V,(x, y ) = f V(x, y , z) dz is the specimen potential projected along z. A thick specimen can now be subdivided into thin slices, perpendicular to the incident beam direction. The potential of each slice is projected into a plane which acts as a two-dimensional phase object. Each point (x, y ) of the exit plane of the first slice can be considered as a Huyghens source for a secondary spherical wave with amplitude q(x, y ) (Fig. 1).

302

D. VAN DYCK

9

---E

E

E

E

FIG. 1. Schematic representation of the propagation effect of electrons between successive slices of thickness E. [From Van Dyck (1983).]

Now the amplitude 4(x’,y’) at the point x’, y’ of the next slice can be found by the superposition of all spherical waves of the first slice, i.e., by integration over x and y, yielding exp(2nikr) dx dy When Ix - x’ I 6 E and ly - y’( e E, with approximation can be used, i.e.,

E

the slice thickness, the Fresnel

so that

[(x

-

x‘)’

+ (y

-

y’)’]

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

303

which, apart from constant factors, can be written as a convolution product: 4(x, y ) = exp[ioI/,(x, y ) ] * exp[ink(x2

+ y2)/&]

(3) The propagation through the vacuum gap from one slice to the next is thus described by a convolution product in which each point source of the previous slice contributes to the wave function in each point of the next slice. The motion of an electron through the whole specimen can now be described by an alteration of phase object transmissions (multiplications) and vacuum propagations (convolutions). In the limit of the slice thickness E tending to zero, this multislice expression converges to the exact solution of the nonrelativistic Schrodinger equation in the forward-scattering approximation. C. Computations For numerical computations, the wave function and the slices are sampled by a two-dimensional network of N closely spaced points. In order to represent the local fluctuations in the wave function and the phase object, the distances between adjacent points are typically on the order of 10 pm for heavy atoms to 50 pm for light atoms. The required computer memory is proportional to N . Typical slice thicknesses are on the order of 0.1-1 nm. For small N the computations can be performed at 32-bit precision. For large N (e.g., 4096), 48-bit precision is recommended so as to avoid the accumulation of rounding-off errors (Lynch, 1982).Hence, with a computer memory on the order of 64 kBytes, common in most mini- or microcomputers, a multislice calculation can be performed for up to about 1000 points, which corresponds approximately to a real-space area (or projected unit cell in the case of crystals) of about 0.4 x 0.4 nm for heavy materials to 2 x 2 nm for light materials. If a fast disk memory is used to store intermediate data, the number of points can eventually be increased by a factor of two or more. In order to handle larger areas within the same computer memory such as for the simulation of extended defects in crystals, the specimen can be subdivided into partly overlapping patches in which the calculations are done separately by using the method of periodic continuation after which the images are joined together (Olsen and Spence, 1981). Using this technique, the memory requirements are not stringent even for image simulation in larger areas. For most purposes, a 64-kByte computer proves satisfactory. However, the overlap areas between adjacent patches should be large enough to avoid scattering from adjacent patches into the image of the nonoverlap region, although this can depend on whether we want to know all the coordinates of all the atoms in the defect. Since the width of these overlap areas is proportional to the crystal thickness, this method cannot be extended to thicker specimens without increasing the computer memory. A much more stringent requirement is

304

D. VAN DYCK

the computing time necessary for the calculation of Eq. (3). The major part is required for the calculation of the convolution products, with time proportional to N 2 , whereas the time for the direct products with the phase objects is proportional to N and, thus, is not dominant. In practice, most of the computations are carried out in reciprocal space by using the Fourier transform of Eq. (3). Here the sampling points are transformed into diffracted beams and the direct and convolution products are interchanged with the total calculation time remaining on the same order of magnitude. When the specimen is not periodic in the x, y-plane, such as for a crystal containing one or more defects, the diffraction pattern is continuous. By sampling the reciprocal space with discrete beams, an artificial periodicity is created in real space, which is termed the periodic continuation. This approximation is valid when the periodicity is large enough, i.e., when the defect is surrounded by enough perfect crystal area to avoid scattering from adjacent periods (unit cells). The area plays the same role as the overlap area used in the method of Olson and Spence, so that in this case also the validity of the approximation decreases with increasing crystal thickness. For most of the image simulations thus far, the reciprocal version of the multislice method has been used. Various programs, written by different authors (e.g., P. Fejes, K. Ishizuka, M. A. O’Keefe, and A. J. Skarnulis) are available. Typical calculation times on mainframe computers are of the order of 10 sec/slice for 1000 beams (17 sec on IBM 370,7 sec on Amdahl, 90 sec on CDC 3600, 20 sec on CDC 6600, and 3 sec on CDC 7600). These values should be multiplied by a factor between 10 and 100 for minicomputers (e.g., 460 sec on Prime 750). D. Use of the Fast Fourier Transform A very elegant way to increase the speed of the computations was introduced by Ishizuka and Uyeda (1977). Here each propagation in Eq. (3) is computed in reciprocal space and each phase-object transmission is computed in direct space. The links are made by a fast Fourier transform (FFT) with a calculation time proportional to N log, N . The number of beams needs to be a power of two in each direction and must be increased by a factor of two, because of aliasing effects (Ishizuka and Uyeda, 1977). However, the total calculation time can be reduced considerably especially for large N . Calculations have been performed for up to 256 x 256 beams. Typical calculation times for the FTT multislice are on the order of 1 sec/slice for 1000 beams on a large mainframe (2.8 sec on IBM 370,0.6 sec on Amdahl, 0.3 sec on CDC 7600) in comparison to 45 sec on H P 21 MX minicomputer.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

305

E . Special-Purpose Hardware The computation of the FFT multislice can be carried out almost entirely using vector arithmetic. Hence, it is very suitable for a minicomputer equipped with a dedicated hardware F F T processor or an array processor with sufficient dynamic range, by which the computation can be speeded up by another order of magnitude (even faster than a large mainframe). Typical calculation times are on the order of 0.1 sec/slice for 1000 beams C0.3 sec on Prime 750 with Floating point systems A P 120 (Skarnulis et al., 1981); 0.025 sec on P D P 1134 with Analogic A P 400; 0.1 sec on P D P 1134 with SKY MNK (estimated)]. As follows from Fig. 2 (Rez, 1980), the computing time for the A P 400 (as compared with the CDC 7600) is no longer dominated by the Fourier transforms but rather by the setup of the phase object and propagator, which are still calculated by the host computer. The calculation time for the phase object is proportional to the number of atoms times the number of beams and is on the order of 1-100 sec on a mainframe.

SLICES FIG. 2. Calculation time for a 1000-point multislice: solid line, PDP 1103 with AP 400; dotted line CDC 6600 mainframe; dashed line, CDC 7600 mainframe. [From Rez (1980).]

306

D. VAN DYCK

F . Discussion

At this point, a number of problems are not yet elucidated. (i) The image calculation technique is still laborious and expensive when complicated systems are involved. For example, the simulations computed for the single-split interstitial in gold, using periodic continuation, have taken 8 h on a UNIVAC 1100 (Fields & Cowley, 1978). Often the simulation of convergent beam electron diffraction (CBED) is unfeasible because of unreasonably large computer times. Although the image simulation technique can be compared to curve fitting in other scientific work, a least-squares refinement between computed images of trial structures and the corresponding experimental micrographs still belongs to science fiction. The future certainly lies in minicomputers (e.g., 32-bit) equipped with array processing facilities, which are becoming continually cheaper. Probably, dedicated image processing systems can also be used efficiently for image simulation. However, any speeding up of the computational algorithms, including the setting-up time for the phase object, would still be very welcome. (ii) In the propagator convolution, an electron is allowed to travel from any point x on one slice to any point x’ on the next slice (Fig. 1). The number of paths N 2 dominates the computing time. However, for most of these paths the Fresnel approximation ( x - x’( 4 E, ( y - y’( 4 E is violated, so that the form of the propagator in Eq. (3) is, in principle, not valid. According to Berry (1971), the effect of these paths is canceled by their rapidly oscillating phases and should be, in principle, unimportant. However, this statement is not necessarily true for numerical computations. For instance, suppose the object is sampled with a square mesh of size 6 x 6 nm. If the slice thickness E is chosen so that, by accident, it obeys the relation k6’/& = 2n with n integer, then exp[ink(x2 y ’ ) / ~ ]= 1 for all paths in the propagator equation [Eq. (3)], i.e., every point (x, y) contributes to every point (x’,y’) on the next slice with the same weight (unity). Hence, the wave function in the next slice becomes constant and all information of the scattering in previous slices is lost. This is an artificial resonance effect, comparable to what is known in multislice theory as pseudo-upper-layer line reinforcement. The choice of the mesh size and slice thickness thus requires special attention. In practical multislice calculations, we use the requirement that ks2/E < 1, which can be derived in different ways. It is related to Heisenberg’s inequality since l/k6’ is related to the excitation error of the outermost Bragg beam, which in fact is the uncertainty on the wave vector, (see Fig. 5), and E is the uncertainty in the position of the scattering. However, the effect of the Fresnel approximation on the propagator as a function of mesh size and slice thickness should be further elucidated. Furthermore, if the propagator paths are restricted to a small

+

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

307

cone (only in adjacent points), which is a reasonable assumption for fast electrons, then the artificial resonance effects are avoided and at the same time the computing time might be reduced to a value proportional to N . (iii) By projecting the slice potential onto a two-dimensional plane, the information about the atom coordinates along z is lost. When, in the case of a crystal, the slice thickness is a multiple of the unit cell parameter c along z, the potential is in fact considered as independent of z throughout the whole crystal. In reciprocal space, here only the zeroth-order Laue zone is excited. In order to obtain information about the upper-layer lines (ULL) or higher-order Laue zones (HOLZ), the slices must be made much smaller and the computation times consequently become larger, e.g., if the nth layer lines are important, the potential should be sampled at intervals Az = c/2n. However, in practice it is sometimes necessary to take the slice thickness E much less than Az in order to obtain a convergent multislice calculation. A recent comparison between multislice and Bloch wave calculations by Shannon (1979) shows that in the case of heavy atoms, such as gold, slice thicknesses on the order of 0.03 nm, which is less than the atomic diameter, must be used. If atom sectioning proves necessary, the multislice method is not very practical. In that case, it should be worth searching for a slice method that accounts for ULL interactions to such an extent that atom slicing can be avoided. (iv) Another problem that has not been fully elucidated is how to know that the slice thickness E and the interval 6 are adequately chosen. In principle we could reduce E and 6 until satisfactory convergence is obtained. Frequently, the normalization of the total diffracted intensity is used as a test for normalization; however, c and 6 may not be considered independently. Indeed, normalization can be obtained for any 6 (i.e., any number of beams) in the limit for E tending to zero (Anstis, 1977) and for any E in the limit for 6 tending to zero (see Section VI). Shannon (1 979) has proposed a rather complicated convergence criterion in which E and 6 are coupled. Another possibility is that E and 6 be coupled by the relation derived in point (ii). However, further investigation of the problem is required. (v) In the derivation of the multislice formulas, forward scattering of the high-energy electrons is assumed. However, thus far no convincing proof has been given for the validity of this approximation, especially when thicker foils are involved. (vi) With the current many-beam diffraction theories, it is still impossible to include inelastically scattered electrons in the calculations. However, since the mean-free-path for inelastic scattering is on the order of tens of nanometers, this omission might put serious restrictions on the maximum foil thickness in the calculations.

308

D. VAN DYCK

(vii) The effect of the thermal motion of the atoms is usually included by multiplying the atomic scattering factors with appropriate Debye-Waller factors. However, this technique stems from x-ray diffraction in which the kinematically scattered x-ray “sees” a spatially averaged atom. However, at present it is not clear whether this approximation still holds for dynamic electron scattering. In the following section we treat the diffraction problem in a more systematic way in order to search for an optimal computational algorithm and to elucidate some of these problems. Iv. SYSTEMATIC APPROACH TO THE ELECTRON-DIFFRACTION PROBLEM A . General 7heory’

For a general treatment of electron diffraction, we start from the timeindependent Schrodinger equation. Contrary to the practice of describing diffraction in reciprocal space, we treat the problem in real space. This approach will prove to have certain advantages (Van Dyck, 1980). Consider an electron with wave vector k and energy on the order of 100 keV-1 MeV, which is incident on a thin specimen foil (Fig. 3) with a local potential V(r). The wave function $(r) of the electron in the specimen is then governed by At,h(r) 4n2k2$(r) V(r)$(r) = 0 (4) with V(r) = (2rne/ti2)V(r) (5)

+

+

and in which rn and k are the relativistic values for the electron mass and wave vector, respectively (Fujiwara, 1961). The incident electron can be represented by a plane wave $(r) = exp(2nik * r) The problem now consists of computing the wave function of the electron from Eq. (4) at the exit plane of the foil (Fig. 3). For high-energy electrons ( 2 100 keV), the influence of the foil can be considered as a perturbation, i.e., V(r) 4 4n2k2,so that it is convenient to rewrite the wave function as a modulated plane wave $(r) = exp(2nik * r)4(r) (6) Van Dyck (1979,1980,1983); Van Dyck and Coene (1984).

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

0

309

z

FIG.3. Geometry of the problem.

The modulation function &r), which is in fact the real-space representation of a depth-dependent Bloch wave, will vary relatively smoothly as a function of the crystal depth. Substitution of Eq. ( 6 ) into Eq. ( 4 ) yields A+(r)

+ 471ik - V+(r) + V(r)$(r)

=0

(7)

When the incident beam is nearly parallel to the surface normal e , , as is the case in most of the experimental situations, the normal component k, of the wave vector is large compared with the parallel component kx,. Hence, it is appropriate to separate Eq. (1) into

where V,, and Ax), are, respectively, the gradient and Laplacian operators in the coordinates, x, y , which are chosen in the foil plane. This equivalent form of Schrodinger's equation requires two boundary conditions for the wave function both at the entrance ( z = 0) and exit ( z = z ) planes of the foil and appears difficult to solve when the specimen is not a perfect and simple crystal. However, by using the high-energy approximation outlined in Subsection IV.B, this equation can be transformed into a firstorder differential equation that requires only one boundary condition at the entrance face.

B. High-Energy Approximation In high-energy electron diffraction where the specimen potential is small compared with the energy of the incident electrons, it is generally assumed that the z component of wave vector k, is large and that the wave function

310

D. VAN DYCK

4(r) varies smoothly with z so that a 2 4 / d z 2 can be neglected with respect to the term 4nik, d 4 / a z resulting in a first-order differential equation

or in shorthand notation W r ) / d z = ACA

+ V(r)l4(r)

(10)

with il= i/4nkZ A = Axy 4nikXy* V,,

+

This approximation is sometimes called the forward-scattering approximation. To the knowledge of the author, Schiske (1950) in his quantum-mechanical treatment of electron motion in magnetic lenses was the first to approximate the time-dependent Schrodinger equation by a parabolic equation including only the first derivative with regard to z, which he termed “paraxial Schrodinger equation.” Equation (10) can now be integrated readily from the entrance face to the exit face. It is the basic equation for the dynamic diffraction calculations. As will be shown later, nearly all existing diffraction formulations are implicitly or explicitly based on this equation. In case of normal incidence Eq. (10) is equivalent with the time-dependent Schrodinger equation in two dimensions (Van Dyck, 1980; Gratias and Portier, 1983). Hence, methods for solving Eq. (10) can probably be used for the solution of genera1 quantum-mechanical problems (see Section V1I.G). C. Validity ofthe High-Energy Approximation

Although the high-energy approximation is intuitively sound, we give, for completeness, a more rigorous proof of its validity, which is not essential for the remainder of this work and may be skipped by the reader. In the case of normal incidence, we represent the wave function as the sum of two particular solutions of the general differential equation [Eq. (411, one corresponding to forward-scattered electrons and one corresponding to backscattered electrons = e2nikz 4(,.) + ,-Znikz W) (1 1) 4(r) obeys the differential equation [Eq. (8)] in shorthand notation

4’ = A(A + V ) 4 + $6”

(12)

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

31 1

and similarly for 8,

e = +A + v)e

-

AW

with

2 = il4nk Since the high-energy approximation is based on the idea that the wave vector k is large (small wavelength), it is appropriate to expand Eqs. (12) and (13) in powers of the wavelength. If we successively substitute 4‘ [left-hand side of Eq. (12)] into [right-hand side of Eq. (12)], we obtain

4’ = [A(A + v) + 22v

+ 23(v” + (A + v)*)+ o(n4)14

(15)

and similarly for Eq. (1 3),

If the potential varies slowly enough over one wavelength distance, i.e., 23V”

22V’ < 2v

then it may be assumed that the series in powers of A can be truncated after a few terms, and the second-order differential equations [Eqs. (1 5) and (16)] can be approximated by first-order differential equations. If Eq. (15) is truncated after the first-order term, the usual high-energy equation is obtained. In the vacuum space at the entrance side of the crystal (Fig. 2), the wave function can be written as tj

= exp(2xikz)

+ exp( -2nikz)

R(r)

(17)

where exp( -2nikz)R(r) represents the reflected waves. Here R(r) obeys Eq. (1 6) in which V(r) = 0, i.e., R‘ = -AAR

+ o(23)

(18)

This approximation is equivalent to the parabolic approximation for the Ewald sphere. In the vacuum space at the exit side of the crystal, the wave function takes the form $(r)

= exp(2nikz)

T(r)

corresponding to the transmitted electrons. By using the same approximation, we similarly obtain for T(r):

T’ = AAT

+ 0(A3)

(19)

312

D. VAN DYCK

From the continuity requirements of the wave function and its derivative at the boundary planes z = 0 and z = d, we obtain at z = 0 1 2nik[1

-

R(O)]

+ R(O)= +(o) + e(o)

+ R'(0) = 2nik[4(O) - O(O)] + @(O) + W(0)

(20) (21)

atz=d

+

exp(2nikd) 4(d) exp( -2nikd) O(d) = exp(2nikd) T(d) x 2nik[exp(2nikd) 4(d) - exp( - 2nikd) O(d)] exp(2nikd) @(d) + exp( - 27cikd)W(d) = 2nik exp(2nikd) T(d) exp(2nikd) T'(d)

+

+

(22)

(23)

By eliminating R(0) from Eqs. (20) and (21) and using Eqs. (12), (13), (14), and (18), we obtain

+

1 - 4(0) = - 2A2A4(0) - 'A V 4 ( 0 ) A2 VO(0)

+ O(A3)

(24)

and similarly from Eqs. (22) and (23), using Eqs. (12), (13), (14), and (19), we obtain exp(-2nikd) O(d) = -A2Vexp(2nikd) 4(d) + A2(2A + V ) x exp( - 2nikd) O(d) O(A3)

+

(25)

so that finally,

4(0) = 1 + A ~ V ( O+) o(n3) and exp( - 2nikd) e(d) = - A2V(d) exp(2nikd) 4(d) + 0(i3)

(26)

i.e., the amplitude of the backscattered electron wave function is small and on the order of A2V. The wave function at the exit plane of the crystal is then [Eqs. (22) and (26)], $(d) = [I - A2V(d)

+ o(n3)] exp(2nikd) 4(d)

(27)

Similar to the results obtained in reciprocal space (Van Dyck, 1976), the realspace expression [Eq. (27)] shows two differences when compared with the high-energy wave function as obtained from Eq. (10). a. Amplitude Correction - A2 V ( d )

This correction stems from the backscattered electrons. It is dependent on the value of the potential at the boundary plane and its effect is negligible.

CALCULATIONS I N HIGH-RESOLUTION ELECTRON MICROSCOPY

313

b. Corrected Differential Equation The first-order differential equation [Eq. (15 ) ] , which determines $(d), is slightly different from the high-energy differential equation [Eq. (lo)] in which the potential V is replaced by V

+ AV’ + A2V” + A2(A + V)’ + ...

(29)

As follows from Eq. (29), the first-order correction on Vis of the form AV‘ and is, thus, of the same magnitude as if the whole foil were shifted in the z direction over the very small distance I A 1, which can be assumed to be negligible. The same argument holds for the term 1’V”. O n the other hand, the secondorder term A2(A + V)’ accumulates on integration and might become appreciable for thicker foils. In reciprocal space this correction term was estimated to become important for foils on the order of hundreds of nanometers thick. Although a careful analysis of this correction term in real space requires a separate investigation, it can be concluded that the high-energy equation [Eq. (lo)] will be sufficiently accurate for the usual dynamic calculations. This equation will be used in the remainder of this chapter. D. Discussion Equation (10) is the basic equation from which most of the current dynamic diffraction formulations can be deduced. The operator A

=

a2

ax

+

a2 .-

ay2

+ h i k , a + 47ciky a ~

ax

-

aY

(in Cartesian coordinates) describes the propagation of the electron in free space and V(r) describes the scattering of the electron at the potential. Notice that AV(r> # V(r)A, i.e., the operators do not commute. The physical meaning of Eq. (10) becomes clear if we study both parts of the right-hand side of the equation separately, which correspond to two simple but physically different processes. If the propagation term is neglected, the basic equation is reduced to

which can be solved readily

3 14

D. VAN DYCK

and from Eq. (6), we have for the wave function $(r)

= exp(2nik * r)

exp( A j I V ( r ) d r )

This is the well-known phase-grating expression in which the projected crystal potential acts as a phase grating for the incident plane wave and which has also been introduced in the physico-optical approach [Eq. (2)]. If the scattering term is neglected, the basic equation becomes a4(r)/az = AA4(r) (33) This is a diffusion type of differential equation in which time t is replaced by foil depth z , thus expressing the diffusion of the electron on propagating through the foil. This physical process can be illustrated clearly in the one-dimensional case, assuming normal incidence (k, = 0, k, = k), Eq. (33) becomes

This equation describes a kind of one-dimensional Brownian motion (apart from the imaginary factor i). Imagine the foil subdivided into thin slices of equal thickness E parallel to the entrance plane (Fig. 4). The wave function in a point x of the slice at depth z E can be calculated from

+

4(x, z

+

E) =

4(x, z )

+ 84 aZ (x, + O2(E) + ... E -

2)

The second-order derivative term in x can be calculated from the wave function in the adjacent points x + 6, x - 6 by using 2) ax2

- +(x + 6, Z )

+ 4(x - 8, Z ) - ~ h2

( xZ ),

+ OZ(6) +

* * '

FIG.4. Schematic representation of the diffusion of the electrons between successive slices. [From Van Dyck (1980).]

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

315

so that finally,

ie _______

-t 47ck6’

[qqx

+ 6, z) + 4(x - 6, z)] + OZ(E,6) + ... +

The wave function at a point x in slice z E can thus be calculated by using a linear combination of the wave function at three adjacent points x 6, x, and x - 6 of the previous slice. The propagation of the electron from slice to slice can thus be described as a random walk process. Hence, if an electron starts from the origin, in the entrance plane, and diffuses through a large number of thin slices, its wave function at the exit plane is given by the Gaussian distribution function

+

as indicated schematically in Fig. 4. The total wave function at the exit plane is then obtained by adding up the probability amplitudes @(x‘,0) for starting at points x’,each multiplied by its properly shifted Gaussian distribution, yielding

or in convolution notation

$(x, z) =

exp(i7ckx2/z)* $(x, 0)

It is clear that the solution of the general diffusion equation Eq. (33), which can be written formally as

can be constructed in the same way. Hence, the phase grating as well as the Fresnel expression for the propagation used in the physico-optical approach of Subsection 1II.B is established on a more rigorous quantum-mechanical basis. Moreover, it is now clear that phase grating and propagator are intimately coupled in one differential equation [Eq. (lo)] for which the multislice expression [Eq. (3)] is only an approximate solution.

316

D. VAN DYCK

v. STUDY OF

EXISTINGMANY-BEAM FORMULATIONS DIFFRACTION A. Diferential Equations

1. Real-Space Representation The starting point from which most of the current dynamic diffraction formulations can be derived is the high-energy equation (10)

and

A

=

a2 ~

ax2

+

a2 -

+ 4nik, a + 47cik, -

ay2

ax

d

-

aY

(37)

We now search systematically for the most efficient formulations and computational algorithms for the solution of Eq. (35). It is interesting to note that Eq. (35) is also equivalent to the time-dependent Schrodinger equation in which the time t is replaced by the depth z in the specimen (Van Dyck, 1980; Gratias and Portier, 1983). Effective computer algorithms might probably also find applications for the solution of more general time-dependent quantum-mechanical problems. This aspect will be discussed in Subsection V1I.G.

2. Reciprocal-Space Representation If the potential V(r) is periodic in three dimensions, such as is the case in perfect crystals, it is appropriate to expand V(r) in its Fourier series V(r)

=

1 5 exp[hig - r]

(38)

8

where g represents the vectors of the reciprocal lattice. By writing

m)= 1 4,w exp(2nig

*

r>

8

and substituting in Eq. (35), we obtain, after identification of the terms,

(39)

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

3 17

3

I

crystal foil

I

I

sphere

FIG.5. Ewald sphere and extinction distance.

where 4,(z) is the amplitude of the diffracted beam g at depth z and s, is the excitation error, approximately equal to the distance between the reciprocal point g and the Ewald sphere, measured along e, (Fig. 5). Equation (40) is the well-known system of coupled first-order differential equations of Howie and Whelan (1961), originally derived for a distorted crystal. Denoting by @(z) the column vector of diffraction amplitudes $&z), T, the diagonal matrix with elements and V, the matrix with elements Vg,g,= inVg-g’

Eq. (40) can be written in matrix notation as

dW)ldz

=

CT +

~ldm

(43)

or

d+(z)ldz

= Sd4Z)

(44)

with S=T+V

(45)

where S is called the dynamic matrix or sometimes the scattering matrix; T can be considered as the Fourier transform of A, which represents the propagation of the electron in free space, i.e., the electron remains in the beam g but receives a phase shift. V can be considered as the Fourier transform of V(r),

318

D. VAN DYCK

and Vg,grrepresents the matrix element (probability amplitude) for scattering from beam g’ into g. It is interesting to note that, henceforth, each step in the description in real space has its analog in reciprocal space and vice versa, which can be obtained simply by Fourier transforming; direct products and convolution products are then interchanged. This duality can provide much complementary insight into the diffraction process. 3. Mixed Representation In some cases, e.g., when the potential of the foil is not periodic along the z axis, it can be more convenient to expand V(r) and #(g) in a two-dimensional instead of three-dimensional Fourier series V(r)

=

1 VG(z)exp(2niG

*

R)

G

and #(r)

=

c &(z) exp(2niG - R ) G

(47)

+

where r = R ze, and g = G + gzez, i.e., the reciprocal vector G belongs to the [OOl] zone. Note that the Fourier coefficients VG(z)of the potential are now dependent on z so as to include the influence of the higher-order Laue zones (upper layer lines). Hence, no generality is lost as compared to the other representation. Substitution of Eqs. (46) and (47) in the high-energy equation now results in

where sG is the excitation error as indicated in Fig. 5. This is the system of coupled differential equations in mixed representation as originally derived by Tournarie (1960, 1961, 1962). It can also be put into matrix form, similar to the foregoing,

d#(z)ldz = CT + V(Z)l#(Z) (49) in which the matrix V ( z )is now dependent on z. In case only the [OOl] Laue zone is excited, i.e., only the projected potential is considered, the matrix Vis independent of z and Eq. (49) becomes d+(z)ldz = (7-+ ~ > + ( z )

(50)

which is similar to Eq. (43). Another kind of mixed representation has been introduced by Howie and Basinski (1968) for the treatment of electron dif-

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

3 19

fraction in a simple crystal containing a lattice defect. If the crystal is not drastically deformed by the defect, it is appropriate to expand the crystal potential and the wave function in a Fourier series as for the perfect crystal, hereby allowing the Fourier coefficient to vary with position (x, y ) so as to account for the influence of the defect. Hence, a system similar to Eq. (40) is obtained, in which partial derivatives with respect to x and y occur. Spence (1978) treated this problem by using a kinematic approximation for the scattering at the defect, thereby obtaining a drastic gain in computing time. Neither of these methods will be discussed here.

B. Matrix Methods In principle, the solution of Eq. (43) or (44) can formally be written in matrix form as

with '1

1

OI 4(0) = 0

I]

0

and where the exponential of a matrix is defined in the usual sense by its power series. This solution was first introduced by Sturkey in the early 1950s (Sturkey, 1962). Sturkey formulated the diffraction process as a unitary transformation that rotates the initial vector +(O) in the n-dimensional diffraction space to the final vector @(z).If absorption is neglected, the scattering matrix S is imaginary and, from simple matrix theory, it can be shown that

l4(Z>l2

= &)4*(4

=

6(0)4*(0) = l4(0)l2 = 1

i.e., the total diffracted intensity remains constant. Since no intensity is lost, this formulation can be termed self consistent. It is very important to note that, since the total intensity is maintained, irrespective of the number of beams g considered (even in a two-beam case), the normalization test cannot be used as a criterion for determining the number of diffracted beams that contribute to the diffraction process. Another matrix method consists of diagonalizing the scattering matrix S by means of a unitary transformation matrix U , i.e.,

U-'SU=D

or

S = UDU-'

(52)

320

D. VAN DYCK

where D is a diagonal matrix with elements yn. The columns of the matrix U are the eigenvectors On of S , where y, are the corresponding eigenvalues, i.e.,

and by using Eq. (52), we obtain from Eq. (44) +(z)

= exp(Sz)+(O) =

U exp(Dz)U-'+(O)

(54)

with eDz a diagonal matrix with diagonal elements eYnz.Equation (53) is the eigenvalue problem that can be derived directly from Bethe's dispersion equation in the forward-scattering limit (Howie, 1978). It has also been used by Sturkey (1957) and Niehrs and Wagner (1955), whereas Eq. (54) has been introduced by Fujimoto (1959). The eigenvalue method is also self consistent in the sense that the total diffracted intensity remains unity (no absorption) irrespective of the number of beams. The eigenvalue method has the advantage that once the eigenvalues y, have been calculated, b(z) can be obtained directly for all thicknesses z and for all +(O). Furthermore, in situations where only a few beams are excited, the method with concepts such as Bloch waves and dispersion surface offers much physical insight into phenomena such as the critical voltage effect (e.g., Gevers et al., 1974 and many others). With modern diagonalization programs, the eigenvalue problem is still tractable when the number N of diffracted beams does not exceed a few hundred. Eigenvalue methods can be used for image simulation in simple perfect crystals, possibly matched with slice calculations for the crystal parts containing defects. Since the memory requirements for matrix methods is on the order of N 2 and the calculation time at least on the order of N 3 , the matrix methods become unmanageable when the number of beams is very large, as is usually the case for image simulations in more complex systems.

C . The Iterative Method When the number of diffracted beams becomes so large that full matrix expressions cannot be handled, we can expand Eq. (51) by using only one vector at the time, i.e., starting with +(O) and ending with +(z). An example of such a method is the iterative method (Van Dyck, 1978a,b, 1979). It is constructed from the definition of the exponential matrix exp[(T

+ l/)z] = 1 ( T +n!v)nzn OU

n=O

(55)

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

321

By a suitable iterative procedure, it can be computed in vector form b y expanding

where Onis defined from Eq. (55) and (56) as

with From Eq. (57), we have the recurrence relation which can be computed by means of three memory vectors [remember that multiplication with V, which is a convolution, can always be speeded up using fast Fourier transforms (FFTs)]. Since the order in powers of z is, in principle, unlimited, this method is by far the best when very accurate results are required, since it takes only two FFT per order. Furthermore, it is normalized, so that the total intensity of the diffracted beams equals unity for n -+ co.This can serve as a criterion for the accuracy but is worthless as an estimation of the number of beams required. This method has been used successfully for some years in o u r laboratory. Disadvantages of the method are It is only applicable to perfect crystal parts. (2) The method can be applied for one crystal thickness z at a time. (3) Especially when thick crystals and large excitation errors are involved, the values of the elements of @ can rise seriously before converging to their limiting values, thus requiring a larger number of significant digits in the calculation. (1)

D. Direct Integrating Methods

In principle, the differential equations [Eq. (35), (43), (49), or (SO)] can be integrated directly from the entrance plane to the exit plane of the specimen. For this, standard routines such as Runge-Kutta or predictor-corrector methods can be used. Specifically for this purpose, a special Runge-Kutta method has been constructed with minimal computer memory requirements (Van Dyck, 1978a; Van Dyck et al., 1980).

322

D. VAN DYCK

If the excitation errors in the diagonal matrix Tare large, so as to include diffracted beams with a large Bragg angle, it is appropriate to use a kind of perturbation technique by replacing 4(z)

(59)

= exp(W

and substituting in Eq. (43), yielding dO(z)/dz

=

Q(z) O(z)

with Q(z) = exp( - Tz) V exp( Tz) The new differential equation allows larger integration steps to be used and the integration proceeds faster. However, it can be shown (Van Dyck, 1979) that Runge-Kutta methods are somewhat slower and more memory-consuming than the slice methods of the same order. Hence, slice methods are still preferred. E . Slice methods

The best candidates for image simulation are the slice methods that, starting from the boundary condition at the entrance plane, integrate Eq. (35), (43), or (49) or expand Eq. (51) step by step throughout the crystal, thereby keeping the memory requirement proportional to N (the number of beams or sampling points). The accuracy and/or speed of slice methods can best be compared to the expansion of the exact solution in one slice in powers of A and V(or Tand V ) , which are proportional to the wavelength (A = l/k) and converge rapidly for fast electrons (Van Dyck, 1980). Several slice methods are discussed in the literature. a. Multislice Method Here the solution within one slice is separated into scattering and propagation

exp(Vz) exp(Tz) = exp[(T

+ V)z] + &VT - TV)’ + ...

if the projection plane is the front plane of the slice or

exp(Tz) exp(Vz) = exp[(T

+ V)z] + &TI/

-

(62)

+

~ ) z ’

if the projection plane is the back plane of the slice. Both slice methods are, thus, expansions of first order in z. Since they do not contain ULL effects within one slice, it is not clear whether the slice must be made sufficiently thin to take into account these effects.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

6. Improved Multislice Method metrical center of the slice, i.e.,

323

By placing the projection plane in the geo-

exp($Tz) exp( Vz) exp(+Tz) = exp[(T

+ V)z] + 03(Tz, Vz)

exp()Vz) exp(Tz) exp()Vz)

+ V)z] + 03(Tz, Vz)

or = exp[(T

(63)

the multislice formulas are transformed into expansions of second order, just as the first-order rectangular integration rule can be transformed into the second-order trapezium integration rule by starting and ending with one-half of an integration step (Van Dyck, 1978b). This improvement can be easily introduced into the existing multislice programs. However, ULL effects within the slice are not included. c. Third-Order Slice Method Van Dyck (1978b, 1979) introduced a third-order slice method that consists of an alternation of functions of T and V, suitable for the use of fast Fourier transforms, i.e.,

exp($aTz) exp(aVz) exp($Tz) exp(bVz) exp(3bTz) = exp(T + V)z 04(Tz, Vz)

+

(64)

with

a

= (3

f i,,/3)/6

and

b

=

1

-

a

The iterative method [Eq. (57)] and the direct integrating method [Eq. (60)] can also be considered as higher-order slice methods. Higher-order methods provide a higher precision for the same slice thickness; they increase the accuracy for the same slice thickness, but they do not increase the speed when the accuracy of the second-order method is sufficient. However, since the differences between the first-order and the second-order multislice is simply equivalent to the shift of the whole crystal by half a slice, which can hardly be observed experimentally, it can be expected that the accuracy of the second-order method will be sufficient for most experimental situations. Higher-order methods also allow a larger slice thickness to be used for obtaining the same precision. However, this benefit is only hypothetical when the ULLs within the slice are important. It will be better to search for a second-order slice method that accounts for ULL interactions within the slices.

324

D. VAN DYCK

F. Discussion The diffraction process can be described adequately by one of the four equivalent equations

In the projection approximation

a&)/az

= (T

+ V)&)

(50)

In principle all the methods discussed earlier are mathematically equivalent in the sense that in the proper limits (e.g., infinite number of beams, zero slice thickness, etc.), they all yield the same “exact” solution. The main difference between these methods lies in the case of treating many beams with respect to computing time and memory. At present the optimal speed/accuracy compromise can be obtained by the slice methods. With slice methods, the wave function can be obtained for a number of foil thicknesses in the course of one calculation, which is advantageous for image simulations. When ULL effects are negligible, the multislice method is sufficiently accurate. However, for a number of possible applications, the present computational algorithms are still too slow and any speeding up would be very welcome. Furthermore, a number of problems and questions raised in Section 111 are still open. In the next section we follow a rigorous but more systematic approach to finding a slice method and a computational algorithm that is optimized with respect to speed/accuracy. We also try to answer some of the open questions. VI. A NEW FORMULATION: THE REAL-SPACE METHOD A . Principles’

The power of a particular method for solving a quantum-mechanical problem is strongly dependent on the choice of a suitable basis. For instance, when electron diffraction is described in a simple perfect crystal, the crystal

* Van Dyck (1980, 19083); Van Dyck and Coene (1984); Coene and Van Dyck (1984a).

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

325

potential, as well as the wave function of the electron, can be described by using a small number of Fourier terms (plane waves), in which case it is advantageous to describe the diffraction process in reciprocal space. However, this benefit disappears when large and complicated (artificial) unit cells are involved or when the periodicity is lost, as is the case for diffraction at extended crystal defects. By using the fast Fourier transform (see Section 111) in the multislice programs, the phase-grating transmission is described in real space, while the propagation is still described in reciprocal space. However, the use of the fast Fourier technique has some disadvantages, e.g., the aliasing effect, the requirement that the number of sampling points in each direction needs to be a power of 2, the calculation time ( N log, N ) and the memory requirements. Hence, if complicated structures are involved, it can be more attractive to describe the diffraction process entirely in real space, where the forward-scattering character of the electrons during propagation can be exploited to speed up the calculations and reduce the memory requirements. For this purpose we try to construct a slice method as a solution of Eq. (35) in real space. Since, for high-energy electrons, A is a small quantity, it seems logical to expand the solution of Eq. (35) in one slice of thickness E in powers of i, which converges faster than the Taylor expansion in powers of E. For this purpose, we start from the integral equation form of Eq. (35):

4(x, Y, E ) = 4 ( x , Y , 0) + 2

s:

[A

+ V(x, Y , z > l $ ( x , Y , z ) d z

(65)

It is interesting to note that Eq. (65) is, in principle, equivalent to the Green’s function approach in the forward-scattering limit. This approach has been used by Fujiwara (1959) and Ishizuka and Uyeda (1977). By using the stationary-phase approximation, the latter authors have rederived the multislice expression. Expanding Eq. (65) in powers of A by iteration yields

[A

+ V ( x , Y , 211 dz

+ ,I2 j i [ A + V ( X ,y , z ) ] dz By calling

+ V ( x , y , z’)] dz‘

326

D. VAN DYCK

the projected potential of the slice at the point (x,y) and

‘zV(x, y , z ) dz

=qx,y) =

z

the center of potential of the slice at the point ( x , y), and using the identities

1:

j I V ( x , y , z ) dz

V ( x ,y , z’) dz‘

1). j;qx,

y, z’) dz’

and writing in shorthand 4(x, y, E ) +(E)

=

{ + LAE + 1

+ A’A(E

-

=

1 Vi 2

=-

= ( E - T)V,

4(~),we have explicitly from Eq. (66)

(AAE)Z l2V; +nv,+- 2 2

~

F)V,

+ A2FVpA + 03(1) + .-.}4(0)

Since A and V are noncommuting, it is advantageous to write the expansion as an alternative of a minimal number of functions of Vand A These functions can be determined by expanding them as powers of A and V and identifying Eq. (70) with Eq. (69). The higher-order powers that are not uniquely determined will be chosen so as to present the functions in closed form. The simplest expression that expands Eq. (69) up to second order is the following:

4 ( ~= ) exp[fAVp(l + d)] exp(lAe) exp[$AV,(l

- d)]4(0)

(71)

where d

= d(x, Y ) =

c.0,Y)

-

(E/2)1/(&/2)

is the “potential eccentricity.” It is the relative deviation between the center of potential Z at the point (x,y) and the center of the slice 42. The wave function at the exit plane of the foil can then be calculated by repeated application of Eq. (71) in successive slices. This will be the basis for the real-space method. As explained in Section 111, the propagation part exp(lAe) can be calculated in reciprocal space and the phase-grating product in real space, both linked by fast Fourier transforms. However, it is possible to set up a new computational algorithm in which all the calculations are performed in real space.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

327

B. Compurison with Other Slice Methods In comparison with the usual multislice (MS) expression [Eq. (62)], the new expression [Eq. (7 l)] shows two corrections (a) Upper layer line effects, which account for the z dependence of the potential, are partly introduced by a kind of dipole expansion of the slice potential, using the concept of potential eccentricity 6. (A brief discussion will be given in Section V1I.F.) In the MS method, the slice potential is projected along z (the so-called projection approximation), which is a kind of monopole expansion that does not take into account the ULL effects within the slice. In order to account for ULL effects, the slices must be made sufficiently thin. However, in two cases the correction term 6 disappears so that the original MS expression becomes a second-order method, apart from the first and last slice (i) when the center of the slice can be chosen in the center of the potential; (ii) in the case of a perfect crystal, where the c parameter is small enough to be taken as slice thickness. Hence, the correction term - 6 of one slice is canceled by the term + 6 of the next slice. In this case the ULL effect is of the same magnitude as the scattering of one slice of unit cell thickness. (b) The first-order MS expression is transformed into second order by simply starting and ending with one-half of a phase grating, which is analogous to the transformation of the first-order rectangular integration rule into the second-order trapezium rule by shifting over one-half an integration step. In the case of a perfect crystal for which the periodicity along the beam direction is small enough to be taken as slice thickness, the second-order expression yields the same results as the multislice expression, apart from the first and last slice. This difference is expected to be large when thicker slices are used, as in the case of light scattering materials (e.g., oxides). Figure 6 shows a comparison between the accuracy obtained with the MS method and the new second-order method in the case of a light-scattering material (SiF,). From this it is clear that for medium foil thickness (below 10 nm) and when very accurate results are required, as is the case for the simulation of dark-field or weak-beam images, the accuracy of the secondorder method is one order of magnitude better for the same slice thickness. However, when normal accuracy is sufficient or when high-scattering materials, such as Au alloys, are involved for which the slice thickness needs to be much smaller, the accuracy of the second-order method is not much different from that of the first-order multislice. In cases where the unit cell periodicity along the beam direction is larger (e.g., > 1 nm) or in case of disorder, the

328

D. VAN DYCK

a

10

13

18

24 g M (nm-’)

270

I

540

I

1080

FIG.6. Accuracy Sz versus foil thickness for different values of the slice thickness E and number of beams. Computations are performed for SiF,[OOl] using 200-keV electrons. Full horizontal scale is 30 nm. Dashed line represents second-order method; solid line represents firstorder multislice method. For a definition of S2see Section VI.D.2. [From Van Dyck and Coene (1984).]

second-order expression may be much more suitable, but this effect is still under investigation. When ULL effects are neglected, i.e., the potential is assumed to be independent of z, the solution of the basic equation Eq. (35) can be written in closed form as

dJ(d= expCWA + V(X, Y)ldJ(O)

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

329

or, in shorthand notation,

4(4 = expCWA + V14(0)

(72) All slice methods are based on expansions of Eq. (72) within one slice. An interesting way of comparing different expansions makes use of the (A, V ) matrix as shown in Table I. Here the terms in the expansion of Eq. (72) are schematically arranged in increasing powers of V and increasing powers of A and the symbols have the meaning (A, V ) = (A2, V ) = (A, V z )=

A’&’ ~

2!

(AV

+ VA)

A2s3 7 (A’V + AVA + VA’)

7(AV’ + VAV + V’A)

A2E3

For the multislice expression, we have (see also Section V.E) exp(2sA) exp(AsV) = exp[Ae(A

+ V ) ] + $I’E~(AV

-

VA)

+ 03(A)

so that the first horizontal and vertical rows are properly included, but the second-order mixing term (A, V ) is only partially included. For the new expression [Eq. (71)] of second order, neglecting ULL effects (6 = 0), we have

exp(4hV) exp(AsA) exp(&V)

= exp[Ae(A

+ V ) ]+ 0 3 ( A )

Hence, the first horizontal and vertical rows of Table I are included as well as the first mixing term (A, V ) .The quality of an expansion is mainly determined by the importance of the omitted terms in Table I. However, the relative magnitude of these terms depends largely on the local potential in the foil. TABLE I CLASSIFICATION OF THE TERMS OF THE TAYLOR EXPANSION FOR EXP[AE(A + V ) ] I N POWERS OF B A N DV

1

A A2 A3

...

I

V

V2

v 3

...

I A A2 A3

V (A, V ) (A2, V )

V2 (A, V 2 )

...

v3

...

...

...

...

330

D. VAN DYCK

For instance, between the atoms, the potential is rather small and smooth and here the electron will not be scattered over large angles. Hence, in reciprocal space, the major diffracted beams will be close to the transmitted beam and their excitation errors will be small (Fig. 7a, region I). In a similar fashion, the term A, which is the real-space transformation of the diagonal matrix of the excitation errors, will be small so that AE g 1. In that case, the

ev

I

'1

Incident High Energy Electrons

L

[AJI It

I

i

z

z+ f

eA

m

Real-Space Electron propagatlon

FIG.7. (a) Regions of reciprocal lattice points, corresponding to increasing excitation errors (distances to the Ewald sphere). (b) Schematic representation of the local calculation of the propagation effect using a linear combination of the values of the wave function at adjacent points.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

331

leading terms in Table I are mainly in the first horizontal row, which is simply the phase-grating expression exp(AEV) = 1

+ LEV + ( 1 2 / 2 ! ) ~ 2+~ z. ..

i.e., for smooth potentials the phase-grating expression might be accurate. Closer to the atoms, the potential change becomes larger and so do the maximum beam angles and excitation errors (Fig. 7a, region 11). In this case the mixing terms in Table I become more important. Close to the atom cores, the potential change is larger and scatters the electrons over a wider angle, corresponding to a larger excitation error, so that the powers of A in Table I become increasingly important (Fig. 7a, region 111). It is clear that the atom cores provide the toughest challenge for a computational method. However, the closer to the atom core, the smaller the scattering area and the smaller the number of scattered electrons. Thus, paradoxically, most of the computational difficulties are imposed by the outermost diffracted beams that carry only a minor part of the total intensity. On the other hand, a correct calculation of the electron wave function at the atom core where the Coulomb potential is infinite requires an infinite number of beams, thus making calculation by any of the current computational procedures impossible. In the original proposal for the real-space method (Van Dyck, 1980) it was argued that, since Eq. (71) is a second-order expansion, it suffices to retain only the terms of Table I up to second order, so that the propagator, which is given by the first column of Table I, can be truncated as exp(1eA) z 1 + h A + $(AeA)’ (73) in order to cut the computational speed. However, as pointed out by Self (1 982), this approximation can lead to computational divergencies that stem from the fact that near the atomic core the electrons are scattered into beams with a large Bragg angle and a large excitation error. In principle, the intensity of these beams is damped by destructive interference caused by the rapid oscillation of eLEAfor large values of A. However, by truncating Eq. (73) after the second-order term, the beams are not necessarily damped but can even be artificially amplified. This effect can be elucidated using a simple analogy: The second-order expansion for

exp(-z) z 1 - z + $2 diverges for large z, whereas exp( - z) vanishes. In principle, these computational divergencies can be removed by using a very small slice thickness E, so that the second-order expansion converges even for the outermost beams. However, in this way the computing time becomes unnecessarily large, and the method is no longer competitive with conventional multislice calculations.

332

D. VAN DYCK

A more effective way of avoiding the divergencies without increasing the computing time is to use a better approximation for the propagation and/or by artificially eliminating the outermost beams from the diffraction process. In numerical calculations in reciprocal space, this problem is usually overcome by restricting the beams in the calculations to within a so-called dynamic aperture, centered around the transmitted beam (low-pass filter). The size of this dynamic aperture is determined by the criterion that the total neglected intensity must not exceed a certain tolerance and depends on the scattering power of the atoms and the thickness of the foil. In real space, this approximation corresponds to a smoothing of the wave function and the phase grating at the atom core by a kind of convolution with the Fourier transform of the dynamic aperture function. A similar, but smaller, smoothing effect occurs when using analytical approximations for the atom form factors (e.g., Doyle-Turner and Smith-Burge) or Debye-Waller temperature factors. A proper smoothing in real space enables us to increase the sampling distance and to decrease the calculation time. We will not go into the details of such a smoothing procedure. C . Numerical Procedure

For the computation of the wave function, the foil is divided into thin slices with thickness E, typically on the order of 0.1 to 1 nm, parallel to the entrance and exit faces of the foil. In each slice the wave function is calculated from the wave function in the previous slice by using Eq. (71):

4(x, Y, z + E ) = exp()lV,(x, Y ) C ~+ x e x P p p ( x , Y)C1

+, y)I} exp(lA~)

-

d(x, Y)l}4(X3 Y, z )

The wave function and the phase grating are sampled by a two-dimensional network of closely spaced points in the (x, y) plane of the slice. The phase grating can be directly calculated in real space by using the analytical Fourier transform of the Doyle-Turner expression for the atomic form factors (Doyle and Turner, 1968; Coene and Van Dyck, 1984a). Afterwards, the phasegrating function is smoothed in order to avoid computational divergencies near the atom core. The spacing between the sampling points is related to the fluctuations in the wave function and the potential and is on the order of 10 to 100 pm. Multiplication with the phase-grating factors expC(WVp(1 - 41

and

exPC(wyo(1

+ 611

is done by multiplying the value of the wave function in each mesh point with the corresponding value of the phase-grating function. The computation of the propagator elAa,which is a convolution operator, requires some special

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

333

attention. As shown in Section IV.D, the propagator describes the motion of the free electron between successive slices as a kind of complex random walk in which each point of the previous slice contributes to each point of the next slice. However, it is known that most of the high-energy electrons scatter forward within a cone with an apex on the order of the maximum Bragg angle (--1 1 ,oo l o rad) similar to the Takagi triangle (see, e.g., Subsection VII.B, Matsuhata et a!., 1984; Humphries and Spence, 1979; Spence et al., 1978). Hence, it can be proposed that the wave function after propagation through a slice may be evaluated from a linear combination of the values of the wave function at a small number of adjacent points in the previous slice, as shown in Fig. 7b. In practical calculations, the maximum diffraction angle between two slices is much larger than the maximum Bragg angle in order to account for the destructive interference of the outermost beams. By using this procedure for the calculation of the propagator, the total calculation time becomes proportional to N , the number of sampling points, which is of the same magnitude as the number of diffracting beams and smaller than the N log, N of the multislice calculations, especially for large N (large unit cells). In the original approach (Van Dyck, 1980), the linear combination was chosen so as to represent the propagator up to the second order in 1.However, as explained earlier, this algorithm can lead to computational difficulties at the atom cores and can only yield reliable results if the number of dynamical reflections (or sampling points) is increased and the slice thickness reduced. (Kilaas and Gronsky, 1983). Recently (Coene and Van Dyck, 1984b), a more systematic approach was introduced. For the sake of simplicity, we restrict ourselves to the onedimensional situation. In case of an orthogonal sampling of real space, the two-dimensional propagation operator can always be uncoupled as the product of two one-dimensional propagation operators. If the wave function after propagation from slice z to slice z + E is calculated from a linear combination of the values of the wave function before propagation in a limited number of sampling points (Fig. 7), 4 ( x n >z

+ E) = C a p 4 ( x n + ~ P

6Z>, = C a p + ( x n + p ,

Z)

(74)

P

with x, the sampling points and 6 the sampling distance, then the propagation in reciprocal space is described by the Fourier transform of Eq. (74), thus yielding

334

D. VAN DYCK

The coefficients up can be obtained by using a least-squares technique, i.e., by minimizing

where the integration is carried out over the reciprocal area S = ( - 1/6, and aM/aa, = o yields for the coefficients ap =

1

+ 1/6),

lS

2 f ( dexp( - 2ZiPW dg

i.e., the coefficients can be determined as the Fourier coefficients of the function f(g). In principle, the most natural choice for up is obtained by taking f ( g ) as the exact analytical propagation function, i.e., in case of normal electron incidence, f ( g ) = exp[ - i m g 2 / k ] (76) As pointed out by Self (1982), the coefficients ap obtained in this way are related to the Fresnel integrals C(x) and S(x) and decrease slowly with p . Hence, a large number of terms are required in Eq. (74) in order to obtain sufficient accuracy, so this algorithm is far from being competitive. The difficulty stems from the fact that, by discrete sampling of real space with a sampling interval 6, the reciprocal space is artificially made periodic in 216. Indeed, the Fourier approximation of the propagator [Eq. (75)] has this periodicity, whereas the exact function Eq. (76) does not, so that their behavior at the boundaries g = f 1/6 is entirely different, which causes the difficulties. A much better accuracy/speed compromise can be obtained by replacing the exact function, by a functionf(g) that differs from Eq. (16) only near the boundaries, where it has the desired continuity properties but where the intensity of the diffracted beams is negligible.

D. Analysis of Input Parameters 1 . Error Tests Used in Current Multislice Calculations

In multislice calculations, the choice of slice thickness and number of beams is generally connected with the normalization criterion, i.e., when absorption is neglected, the total intensity in the included diffracted beams at a certain thickness is compared to unity as a measure for the amount of electrons being lost by the numerical calculation. It is generally believed that 90%

CALCULATIONS I N HIGH-RESOLUTION ELECTRON MICROSCOPY

335

is a reliable threshold value, and if normalization is below this value, more beams and thinner slices are to be used. However, some authors (Lynch, 1971; O'Keefe, 1973) allow much lower normalization tests of about 70%, which is then justified by the assumption that even for such a diminishing accuracy the relative values of the diffracted beams are not drastically altered. Bursill et al. (1977, 1978) proposed a set of four tests on the calculations. First, a normalization test on the phase-grating function is carried out, and the preceding normalization test on the electron wave function is taken into account. The authors then considered an independent convergence test for an increasing number of beams and, separately, a decreasing slice thickness. The general validity of the normalization criterion has been critically questioned by Anstis (1977). In his paper, he proved that the multislice procedure converges to a normalized result for the slice thickness tending to zero, regardless of the number of included dynamic reflections. On the other hand, the analytical form of the slice method is normalized independent of the slice thickness. These statements show that the normalized criterion cannot be a satisfactory test for both the number of beams and the slice thickness, if considered independently. Shannon (1978) has presented a more sensitive but complicated test combining the slice thickness and the number of beams by checking some particular multislice results with eigenvalue calculations. 2. S2 Test Criterion

Since most of the dynamic calculations are used for high-resolution image simulations, a suitable error criterion should reflect the visibility of the calculation error in the final image. The final imaging procedure uses only a limited number of beams, selected by the objective aperture of the microscope (radius rA);these beams are situated close to the optical axis. The outermost reflections do not contribute directly to the HREM images, although they are more sensitive to change in the numerical parameters. The S2 test criterion (Coene and Van Dyck, 1984a) is defined as s2(E, g M ,

z, r A )

=

1

lgI5rA

Id)R(g, z ,

-

4 k z>12

where 4(g, z ) is the amplitude of the reflection g at a depth z calculated with a particular choice of input parameters (c, gM), where gM is the absolute value of the largest g in the calculations; and 4R(grz ) is the corresponding amplitude obtained with a reference calculation ( E very small, gM very large). As follows from Parseval's theorem, S2 also represents a kind of average quadratic deviation from the exact wave function in a discrete sampled real space, so that it can be considered as a measure of the contrast difference in the final image.

336

D. VAN DYCK

3. Study of the Input Parameters Coene and Van Dyck (1 984a,b) have calculated a large number of images for systems with high (Au,Mn) and with low (SiF,) scattering power and a range of choices for the input parameters. By visually comparing these images with the corresponding reference calculation it can be concluded that, in contrast to the usual normalization test, S 2 is a more faithful measure of the visible contrast error, whereby S2 = 1 x can be used as a safe criterion. Using this criterion, a rule of thumb for accurate calculations can be established where E is the slice thickness, k is the wave function, and 6 = 1/29, the sampling distance in real space. Equation (6) also ensures that the influence of pseudo-upper-layer reinforcements (Shannon, 1978; Lynch, 1971) is fully suppressed. In the literature, however, some authors have performed computations with input parameters that violate the compulsory rule [Eq. (77)] so that the reliability of such calculations is doubtful. In practice, the value 6 = 1/29, is mainly determined by the scattering power of the material. The slice thickness and g, are then obtained from Eq. (77). For instance, in the case of Au,Mn, gM is typically on the order of 50 nm- (6 = 20 pm), thus yielding a very small slice thickness E of about 20 pm (200-keV electrons). For SiF, on the other hand, g, zz 20 nm- ‘(6 % 50 pm) and E z 0.25 nm. It may seem surprising that the number of reflections that needs to be included in the dynamic calculations according to the S2 criterion clearly exceeds the number of observable reflections in the experimental diffraction pattern. However, before these outermost beams are extinguished by the destructive interference beyond a certain thickness, they may have been scattered into lower-order reflections so that their inclusion in the dynamical calculations is necessary. E . Discussion

The real-space method is mathematically sound and, in the limit, yields the “exact” solution of the high-energy diffraction equation. The forwardscattering character of the electrons can be exploited to decrease the computing time. Although the performance and accuracy are more critically dependent on the choice of the input parameters, such as number of beams and slice thickness, the computational algorithm can be tuned by using rules of thumb like Eq. (77). Since the difference of accuracy between 32- and 64-bit calculations is found to be smaller than the critical Sz value (Coene and Van Dyck (1984b), the optimized real-space calculations (the setup of the phase

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

337

grating not included) can be performed with 32-bit precision even in the case of the high-scattering material Au and for foil thicknesses up to 30 nm. However, 32-bit calculations are roughly twice as fast and require only half of the computer memory. In this way, calculation times on the order of 1 sec/slice (for 1000 sampling points) can be achieved on a 16-bit minicomputer, which is competitive when compared with other methods (see Section 111). Since the calculation time is proportional to the number of sampling points, i.e., to the sampled area, the method is very suitable for handling complicated systems. Indeed, the operator describing the propagation of the electron operates locally and is independent of the crystal structure. Hence, as opposed to fast Fourier techniques, all restrictions concerning the size of the unit have disappeared (see Subsection VILA). VII. RECENT DEVELOPMENTS, UNSOLVED PROBLEMS, AND PROSPECTS FOR THE FUTURE3

A . T h e Real-Space Patching Technique4

In the real-space method, very extensive structures, as required for calculating defect images, can be divided into a set of adjoining, nonoverlapping regions, called patches. The dynamic calculation for a thickness increment of one slice can be performed in each patch separately in the central computer memory. The other patches are stored on fast disk memory. The high-energy electron scattering, represented in a multiplication with the phase-grating expression, is a point-to-point operation and is implemented in a simple way. Since the electron propagation is calculated locally, it requires only some extra consideration at the edges of each patch where the information of the adjacent sampling points in the neighboring patches is needed. It should be noted that a small but not insignificant fraction of the computing time is required by the input/output of the patches from disk. This patching procedure is in principle also possible in FFT-multislice calculations. However, the convolutions calculated by using fast Fourier transforms do not operate locally, so that the input/output time becomes proportional to the product of the number of patches with the number of sampling points. Hence, the applications of a patching method in the FFTmultislice calculations is not so efficient. Another kind of patching technique was introduced by Olsen and Spence (1981). Here the dynamic diffraction is treated for the whole crystal thickness 3Van Dyck ( 1 984). 4Coene and Van Dyck (1984).

338

D. VAN DYCK

in each patch separately by using the multislice method of periodic continuation (see Subsection V1I.B). Then the obtained images are joined together. However, the periodic continuation artificially influences the image near the boundaries of each patch. Although this effect can be eliminated by introducing an overlap area between adjacent patches, the effective area of each patch is reduced, so that the total computing time increases, especially for thicker crystals. For the simulation of HREM images of a crystal defect, the FFTmultislice algorithm inevitably needs the periodic continuation of the defect, sometimes by artificially introducing returning defects (Wilson and Spargo, 1982). In real-space calculations, on the other hand, the defect area need not be embedded in an artificially repeated supercell, but can be surrounded by a perfect matrix, for which the wave function is calculated separately. Figure 8a shows an example of a very extensive structure, originating

1 1 1

1

z: 0

z=t

FIG.8. (a) High-resolution image of the overlap area of two coaxial Au,Mn variants. The images of the respective variants are visible at the extremes of the image, where the white dots correspond to the Mn columns. (b) Schematic drawing of the overlap area. [From Coene et al. (1985).]

339

340

D. VAN DYCK

from the overlap between two coaxial Au,Mn variants with the common axis along the electron beam. The compound Au,Mn is a column structure consisting of columns of Au and Mn. When the crystal is oriented with the columns parallel to the electron beam, the Mn columns often appear as white dots. Hence, it has been proposed [Van Tendeloo and Amelinckx (1978a,b)] that the white dots in the overlap area stem from overlapping Mn columns. The model used for image simulation consists of two overlapping wedges of both variants (Fig. 8b). A projection of the model structure along the beam axis is shown in Fig. 9. Simulations are performed for different specimen thicknesses, wedge angles, and focus values. Some of the results are shown in Fig. 10, from which it is clear that the white dots indeed correspond with

- 50 - 60

- 70 - 80 - 90 defocus(nm) FIG.10. Simulated images of the overlap area for different focus values. The wedge angle is 38" and the foil thickness 8.1 nm. The width of the overlap area is also indicated. [From Coene et al. (1985).]

CALCULATIONS I N HIGH-RESOLUTIONELECTRON MICROSCOPY

341

overlapping Mn columns. More work must still be done in order to elucidate the structure completely.

B. Considerations Concerning the Periodic-Continuation Method Apart from the real-space programs, the current computer programs for the simulation of high-resolution electron micrographs are devised for perfect crystals. Hence, the images of isolated defects are simulated by using the method of periodic continuation (Grinton and Cowley, 1971) in which the defects are repeated periodically, yielding a perfect crystal with a very large unit cell. Intuitively, this method is justified when the distance between adjacent defects is large enough to avoid overlap between their images. When the width of the defect image is large, as is the case for large defocus values, under conditions in which Fresnel-like fringes appear, the image simulations often require prohibitively large computer demands. Most of the computing time is required for the calculation of the amplitudes of the dynamically diffracted electron beams within the crystal. However, for thin crystals, it can be shown that the broadening of the defect images due to diffraction in the crystal is much smaller than that due to the electron transfer through the microscope, especially when large defocus values are used for imaging. Figure 1 l a shows the model system used for the calculations. It consists of a Au,Mn matrix in which periodic antiphase boundaries (APB) are introduced parallel to the incident beam. The distance between successive APBs is chosen as 1.277 nm, which leads to a nominal composition of Au,,Mn,; it corresponds to an actual one-dimensional periodic antiphase boundary structure found in the Au-Mn system. This particular model structure was selected in the first place because the scattering power is expected to be large and to provide a challenge for image simulation programs. The electrons crossing the APB plane can be visualized by calculating the electron density at the exit surface of the defective crystal by using the method of periodic continuation (Fig. 11b) and by subtracting from this the electron density calculated for the perfect crystal shifted along the APB. This image can then be displayed as a halftone picture (Fig. 1lc). It represents in fact an image of the deviation from the column approximation (Howie and Whelan, 1961). Figure 12 shows how the width of the APB image increases approximately linearly with increasing thickness. Similarly, it can be shown that in the electron microscope, the image widens further with defocus (Matsuhata et al., 1984). In both cases, the proportionality factor is about Since in most of the experimental situations the crystal thickness is smaller than the defocus value, it is profitable to perform the image simulation for isolated defects in

A.

342

D. VAN DYCK

(b)

(C)

FIG. 11. (a) Structure model, consisting of a perfect Au,Mn matrix with periodically introduced APB. (b) Electron density at the exit face of the foil (thickness, 8 nm). (c) Difference image using the method outlined in the text. [From Matsuhata et al. (1984).]

two stages so as to reduce the computing time and the required memory. For the simulation of the dynamic diffraction in the crystal, the artificial periodicity can be taken as small as possible, e.g., somewhat larger than of the crystal thickness. For the simulation of the transfer through the electron microscope, the artificial periodicity can be increased to a distance somewhat larger than of the largest defocus value. A similar reasoning holds mutatis mutandis also for the real-space method. It should be noted, however, that the technique cannot be applied for the calculation of the diffraction pattern. Indeed, as pointed out by Wilson and Spargo (1982), the simulation of the diffraction amplitudes requires special attention. In order to simulate the diffuse intensity scattered from an isolated defect by discrete sampling in reciprocal space, the scattering between successively repeated defects should be separated, which puts some restrictions on the nature of these defects and on their separation distance.

&

CALCULATIONS I N HIGH-RESOLUTION ELECTRON MICROSCOPY

343

FIG. 12. Width of the APB image as a function of thickness. Thicknesses are (from left to right) 8, 16,24, 32,40, and 48 nm. [From Matsuhata et al. (1984).]

C . Atomic Column Approximation Recently, a large amount of high-resolution electron microscopic work has been done on substitutional alloy systems with a column structure (e.g., Amelinckx, 1978; Van Sande et al., 1979; Schryvers et al., 1984; for a more complete reference list see Van Dyck et al., 1982). It has been shown by Van Dyck et al. (1982) that, under certain suitable experimental conditions, the high-resolution electron images of substitutional alloy systems can be interpreted directly in terms of the atomic columns, viewed along the beam direction. This feature is related to the fact that electron scattering is forward (see Subsection VII.B), and that for thin specimens, the wave function at the exit face reveals a local nature (Humphries and Spence, 1979; Spence et al., 1978; Bourret et al, 1983). This leads to the idea of the atomic column approximation in which the electron scattering is treated dynamically within the atomic columns but scattering between columns is neglected. Hence, for the image calculations in complicated alloy systems with a column structure, it suffices to perform the dynamic calculations for each type of atomic column only and to assemble the wave functions at the exit plane afterwards. This procedure can yield a drastic gain in computing time, especially for very complicated column structures (Matsuhata et al., 1983b). Figure 13 shows the intensity at the exit plane of the Au,Mn column structure for various crystal thicknesses, calculated with a full dynamic calculation and with the atomic column approximation (200 keV). Figure 14 shows the image after the electron microscopic transfer. From this it is clear that even for such heavily scattering material, the atomic column approximation holds to a thickness of somewhat less than 20 nm (Matsuhata et al., 1983b). It would be very interesting if the method could still be improved somewhat, so as to hold for thicknesses that cover most of the experimental situations. This work is still in progress.

344

D. VAN DYCK

FIG. 13. Electron density at the exit plane of a Au,Mn crystal, viewed along the atomic columns, for different foil thicknesses (a) using a complete dynamic calculation and (b) using the atomic column approximation.

D . Absorption, Inelastic Scattering, and Thermal DifSuse Scattering

The importance of absorption and inelastic scattering is generally underestimated and often ignored in simulation programs. However, in biological materials for instance, the mean-free-path for inelastic scattering is on the order of tens of nanometers so that the image simulation for thick crystals is probably only faithful if the effect of absorption is properly taken into account. In electron-diffraction calculations in reciprocal space, anomalous absorption is often introduced in a first approximation by using imaginary Fourier coefficients for the crystal potential (Yoshioka, 1957). In this way, the inelastically scattered electrons are eliminated from the diffraction process. The real-space method is particularly suitable for the investigation of absorption effects. By Fourier transforming the absorption coefficients into real space, an absorption potential is obtained that consists of bell-shaped imaginary absorption functions localized at the atom positions. The absorption

FIG.14. High-resolution image of a Au,Mn crystal as a function of foil thickness, calculated as for Fig. 13: (a) complete dynamic calculation and (b) atomic column approximation. (Cs = 1.1 nm, aperture = 7 nm-’, defocus = -45 nm).

CALCULATIONS I N HIGH-RESOLUTION ELECTRON MICROSCOPY

345

functions are composed of two major parts stemming from inner shell electron excitation and phonon excitation, respectively, whereas the plasmon contribution can be considered as constant over the crystal. The first part is dependent only on the atom type and can be tabulated, e.g., in the form of a sum of Gaussian functions. These values can be used for all real-space highresolution image calculations. The phonon contribution, however, which is sharply peaked, is also dependent on the structure of the crystal. Hence complete models for the absorption potential are only available for simple crystals, such as FCC metals (Humphries and Hirsch, 1968; Radi, 1970). In order to study the relative importance of the different inelastic processes, we performed dynamic calculations for Au[001] using different assumptions for the absorption potential (Matsuhata et al., 1983a). Figure 15 shows the calculated electron density at the Au atom core, where the sensitivity to absorption is maximal, as a function of crystal thickness. Calculations were performed using the model of Humphries-Hirsch with electron-phonon interaction (solid line), without electron-phonon interaction (dotted line), and without absorption (dashed line). It is interesting to note that, without absorption, the electron density at the atom core oscillates with crystal thickness with a periodicity of about 4 nm. This effect will be discussed in Subsection V1I.G. From Fig. 15 it is clear that calculations without absorption are only valid for small thicknesses and that the electron-phonon interaction constitutes the major part of the absorption. Unfortunately, an accurate calculation of the

20

10

0

0

2 0 nm

FIG. 15. Electron density at the atom core of Au [Ool] as a function of crystal thickness. Solid line uses the model of Humphries-Hirsch with electron-phonon interaction; dotted line uses the same model without electron-phonon interaction; and dashed line represents model without absorption.

346

D. VAN DYCK

electron-phonon absorption potential requires a knowledge of the lattice dynamics (phonons) in the crystal, which in turn requires a knowledge of the crystal structure. A single atom absorption function is no longer valid (Whelan, 1984). Hence, it might be intrinsically impossible to account properly for absorption in crystals with an unknown structure, which poses a serious problem. Furthermore, in electron image simulation, the inelastically scattered electrons are considered as lost for imaging. For thick crystals, it might also be worthwhile to account for the inelastically scattered electrons. However, since these electrons cannot be considered as coherent with the incident electrons, the computation can be a tremendous task. Although some theoretical progress has been made in the past (Rez, 1978; Serneels et al., 1980), there still remains a lot of work to be done (Howie, 1983). More experimental data, such as energy-filtered high-resolution images, would also be very welcome. Another problem that has not yet been treated adequately is the influence of the thermal motion of the atoms on the diffraction process. This effect is usually introduced by using a technique adapted from x-ray diffraction, i.e., by multiplying the atomic scattering factors with a Debye-Waller factor. In a sense, the atom is replaced by a kind of spatially averaged atom. In principle, the Debye-Waller factor can only be used to calculate the intensity of Bragg reflections when the diffraction is kinematic. However, at present, it is not clear whether the approximation also holds when the diffraction is highly dynamic. Furthermore, the diffuse intensity scattered in between the Bragg reflections cannot be treated. In principle, the thermal diffuse scattering (TDS) can be calculated by using a very large unit cell in which all the atoms may have different displacements. In order to account for the temporal averaging during exposure on the photographic plate, different images must be calculated and added. At present, such procedure is still limited by the computer demands.

E . The Reverse Problem: Direct Structure Retrieval Instead of using trial-and-error simulation techniques, which are only capable of testing small modifications of approximately known structures and are very time consuming, it would be much more interesting to search for methods that could extract the structural information directly from the electron micrographs. Such a direct method should proceed in two stages. First, the effects of the electron microscope transfer must be eliminated, i.e., we must return from the image to the exit plane of the specimen. Although much progress has been achieved in this field (e.g., Saxton, 1978; Kirkland, 1982; Saxton and Stobbs, 1984), especially for thin specimens, there still remains work to be done.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

347

In the second stage, the problem of deriving the crystal potential from the knowledge of the wave function at the exit plane of the crystal foil should be solved. Although it is sometimes argued that this problem cannot be solved mathematically due to questions of uniqueness, we show that in principle direct procedures are possible. Using the notation of Section IV and neglecting ULLs, the wave function at the exit plane of the crystal is given by 4(z)

= exPrMA

+ Vl 4 0

(78)

with $(O) = 1 for a normal incident plane wave. The problem now is how to derive the projected potential V, which is a two-dimensional real function, from the knowledge of 4(r), which is a two-dimensional complex function. If the crystal foil is sufficiently thin, so that the propagation effect can be neglected, the wave function can be expressed by the phase grating 4(z)

= exp(AzV)

(79)

from which V can be derived in a direct manner. If the crystal is thicker so that the phase grating does not hold, we start from the identity exp(AzV) = exp(AzV)

+ 4(z)

-

exp[Az(A

+ V)]

(80)

which can be used to construct a recursive algorithm yielding successive values V " , which converge toward V for n -+a, i.e., exp(AzV"+')

= exp(AzV")

+ 4(z)

-

exp[Az(A

+ V")]

(81)

with V o = 0. Note that V' is the result obtained by using the phase-grating expression. It is clear that if V" converges toward a limit, then Eq. (80) is satisfied. The limit thus presents the projected potential V if the solution of Eq. (78) is unique. In order to test this statement, we proceed as follows. First, the wave function 4 ( r ) is calculated with the RS method, starting from a known crystal potential V. In practice, we use the crystal SiF,, the projected model of which is depicted in Fig. 16. The recursive procedure [Eq. (Sl)] is started, and the successive results for V" are compared with V. Figure 17 shows the results for V' (dotted line), V 2 (dashed line), and V (solid line) taken along a trace that intersects the Si atoms (the origin is taken at the Si atoms and the two peaks are located at the F atoms). Crystal thicknesses are 1.0 and 1.6 nm, respectively. From this, it is clear that Eq. (78) converges to the exact solution for crystal thicknesses that exceed the phase-grating limit. This work is still in its initial stage. In the near future we must study other recursive algorithms instead of Eq. (78), which are expected to converge for larger crystal thicknesses, such as exp(lzV"+')

= exp( - AzA/2)$(z)

x exp[dz(A

-

exp( - AzA/2)

+ V " - ' ) ] + exp(AzV")

(82)

348

D. VAN DYCK

0

0

0

0

FIG. 16. Model of SiF,. Black circles represent Si atoms, and white circles represent F atoms; a=0.383 nm; c=0.540 nm.

It should be noted that in practical situations the crystal thickness z is not known and the wave function +(z) is only known to a certain resolution. These problems are still being studied. Another interesting approach that allows direct structure information is the so-called improved phase-grating approximation (Van Dyck, 1983). Substituting +(r) = expCWl (83) into the high-energy diffraction equation [Eq. (35)] yields, for normal incidence (in shorthand notation), ae/az =

n[(ve)z + A e + vl

(84)

FIG.17. Retrieval of the potential for SiF,. The abscissa represents the line intersecting the Si atoms (see Fig. 16.). Dotted line represents first iteration; dashed line represents second iteration; and solid line represents exact potential; t is the crystal thickness; (a) t = 1.0 nm and (b) t=1.6nm.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

349

where V and A are the gradient and Laplacian operators in the (x,y) plane, respectively, and 1 is defined as in Eq. (36). Equation (84) can also be written as an integral equation

O(d) = O(0)+ 1

Jod + Jod V dz

[(Ve)’

1

+ Ae] dz

where d is the foil thickness. When ULL effects are neglected, V is independent of z, so that B(d)

=

O(0)

+ 1Vd + 1

[(VO)’

Jod

+ AO] dz

(85)

Expansion of 0 in powers of 1 now yields with B(0) = 0 B(d) = 1Vd

+ (1’d2/2)AV+

(A3d3/3)(VV)’ + (d3d3/6)A2V+ O(A4)) (86)

which can be approximated up to second order as

e(d) =

+(i~/2)p

(87)

where A V = p is proportional to the projected charge density in the (x,y ) plane and hence, from Eq. (83),

4 = exp(Al/d + (1’d2/2)p)

(88)

An approximation of Eq. (88) can also be derived by using the Bloch wave formalism (Pirouz, 1981), so that in this way a direct link can be made between Bloch-wave and real-space theory. Since from Eq. (36) 1 is an imaginary constant, the electron density at the exit face of the crystal is given by

14(d)12 = exp(L2d2p)= 1 + L2d2p

(89)

Whereas for a pure phase grating the electron density remains unity, there is now a term, proportional to the projected charge density and the square of the crystal thickness, which is positive at the atom core and negative in the electron cloud. That is, if 141’ is visualized as a halftone picture, such as in Fig. 14, the atom core is imaged as white peak surrounded by a dark ring. This effect was observed by Ishizuka and Uyeda (1977) by using multislice calculations on Cu phthalocyanine, where it was attributed to a lens effect of the atoms. From these calculations, it can be estimated that in such material Eq. (88) holds to about 15 nm for 500 keV electrons, which is far better than the pure phase grating. Equation (88) or (89) can also be used to find an approximate solution to the reverse problem. Indeed from the phase of the wave function at the exit plane we know dV, and from the amplitude we know d’p, from which, by

350

D. VAN DYCK

solving Poisson’s equation, we obtain V. Hence, the projected potential is obtained in two ways and the ratio yields the crystal thickness. As a test for the validity of Eq. (88), this ratio must be constant over the whole exit face. Deviations are expected to start from the atom cores on increasing crystal thickness.

F . Upper Layer Lines and Beam Tilt Most of the algorithms currently in use for image simulation are based on the projection approximation in which ULLs or HOLZs are neglected (Spence, 1981). Although Lynch (1971) argued that ULLs are important for Aurlll], it is generally believed that the projection approximation is valid when the incident beam is nearly parallel to a zone axis and when the periodicity c along the beam direction is not too large. This is also clear from Subsection VLB, where it is shown that the ULL effect is of the same order of magnitude as the scattering of a slice of unit cell thickness. However, there are a number of experimental situations in which ULL effects can be important. First, when the periodicity c is large (Nihoul and Cesari, 1984) or when disorder is present along the c axis. If the crystal planes are perpendicular to the zone axis, we can use either the usual multislice simulation programs, provided c is divided in several slices, or the real-space method corrected for potential eccentricity [Eq. (71)], which allows larger slice thicknesses to be used. However, a problem arises when the crystal planes are not perpendicular to the zone axis as is the case for monoclinic crystals viewed along the monoclinic axis. Here the two-dimensional unit cell of the projected potential can be small if the projection extends over the whole c, but may be infinitely large for each of the subdividing slices (Wood, 1982). Hence, although the three-dimensional unit cell may have relatively small dimensions, the use of the multislice programs will require the use of the periodic continuation technique. In this case the use of the real-space method can be much more appropriate. Another experimental situation where ULL effects can be important is the imaging of surface structures by using dynamically forbidden reflections (e.g., Yagi, 1984). When the c axis is not too large, the real-space formula can be used, since the surface structures are reflected in the potential eccentricity 6 of the first and last slice. Another very important experimental situation in which ULL effects can also be important arises when the incident beam is inclined from the zone axis so that the Ewald sphere cuts the higher-order Laue zones, as is the case of convergent electron beam diffraction, STEM, and high-resolution electron microscopy with a convergent incident beam (open condensor). Although the multislice method (Ishizuka, 1982) and the real-space method can in principle be adapted for beam tilt, the current programs are

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

351

not adequate for larger tilt angles. If the various incident beams of the illumination core are incoherent, different images must be calculated for each of these directions of incidence, and the images must be superimposed afterwards. This procedure, however, is very time consuming and prohibitive. Perturbation methods (e.g., Von Hugo, et al., 1984) are promising, but more theoretical work must be done.

G. Application to General Quantum-Mechanical Problems As already stated by Van Dyck (1980), and Gratias and Portier (1983), the high-energy diffraction equation for normal incidence is of the same type as the time-dependent Schrodinger equation in two dimensions, provided the depth z is replaced by the time t. This can easily be demonstrated by the following example. Imagine a column structure such as Au,Mn. If ULLs are neglected, the potential is only dependent on (x,y), and, on moving through the crystal, the electron sees atom “tubes.” If we calculate the wave function at successive crystal depths, we obtain a time sequence of the motion of an electron in a twodimensional assembly of atom potentials. The electron density at various depths in Au, Mn is shown in Fig. 18 (200 keV). This figure can be compared with Fig. 13. From this, it is clear that at the atom core of Au the intensity

FIG.18. Electron density at the exit face of Au,Mn (cf. Fig. 13) for different thicknesses. The Mn positions are located at the dots in the corners. The Au positions are located at the four dots in the center of the squares.

352

D. VAN DYCK

increases and decreases periodically with the depth with a periodicity of approximately 4 nm, which is in agreement with the oscillations observed in Fig. 15. This phenomenon can be explained intuitively as due to an oscillator-like behavior of the electron in the electrostatic potential of the atom. The periodicity of 4 nm corresponds with a time period of about 2 x lo-’’ and a frequency of about 5 x 1OI6 Hz. A similar behavior is observed at the Mn atom cores but, due to the smaller electrostatic potential, the periodicity is much larger. It would be interesting to study the relation between this phenomenon and the interaction between Bloch waves as well as the correspondence with atomic and molecular orbitals [as investigated by Buxton et al., (1978)l or Bloch wave channeling (Humphries, 1980). From the foregoing, it is clear that the slice methods, and in particular the real-space method, which is optimized for solving the basic equation [Eq. (35)] by depth slicing, can also be used to solve the time-dependent Schrodinger equation by time slicing, i.e., the slice methods are in fact numerical approximations to the Feynman path integrals (Feynman and Hibbs, 1965; Van Dyck, 1975). Hence, the real-space method might be put into good use for the solution of more general quantum-mechanical problems, but this requires further investigation. H . Conclusion and Prospects for the Future For the present type of image simulations in crystals, the current multislice programs are still adequate if the unit cell is not too large. When large unit cells, defects, disorder, and upper layer lines are involved, the real-space method, combined with the patching technique (Subsection VILA), seems to be more appropriate. Probably the simulation technique can be extended to other areas, such as the diffuse scattering from thermal atom motion or substitutionally disordered systems. However, in order to make image simulation more reliable, a more quantitative comparison between experimental and simulated images, combined with a least-squares refinement, is highly desirable. For this purpose, a speeding up of the computational algorithms is very welcome. In the real-space method some extra speeding up can be expected from better sampling schemes that take into account the local nature of the potential and probably also from the improved phase grating approach [Eq. (85)]. In this respect, a dedicated image computer can be put into very good use either for the treatment of the digitized experimental images or for the calculation of the simulated images. Here the disposal of an array processor is also advantageous. More studies in the near future must be concentrated toward the effect of beam tilt, for the treatment of convergent beam electron diffraction, and toward inelastic scattering. Also, a real-space description of the image forma-

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

353

tion in the electron microscope can yield some benefits (Marks, 1984). It is our belief however, that, on a long term, the direct structure retrieval methods, combined with on-line image capturing, filtering, and processing devices, are the most promising.

ACKNOWLEDGMENT The author is indebted to Dr. Coene, Dr. Matsuhata, Prof. R. Serneels, Prof. J. Van Landuyt, Dr. G. Van Tendeloo, and Prof. S. Amelinckx for many fruitful discussions and for the use of photographic materials.

REFERENCES Alpress, J. G., Hewat, E. A., Moodie, A. F., and Sanders, J. V. (1972). Acta Crystallogr. Sec. A 28, 528. Amelinckx, S. (1978). Chem. Scr. 14, 197. Anstis, G. R. (1977). Acta Crystallogr. See. A 33, 844. Berry, M. V. (1971). J. Phys. C4,697. Bethe, H. A. (1928). Ann. Phys. 87,55. Bourret, A., Thibault-Desseaux, J., DAnterroches, C., Penisson, J. M., and De Crecy, A. (1983). J. Microsc. 129. 337. Bursill, L. A., and Wilson, A. R. (1977). Acta Crystallogr. See. A 33,672. Bursill, L. A,, and Wood, G. J. (1978). Philos. Mag. A 38,673. Buxton, B. F., Loveluck, J. E., and Steeds, J. W. (1978). Philos. Mag. A 3, 259. Coene, W., and Van Dyck, D. (1984a). Ultramicroscopy. 15,41. Coene, W., and Van Dyck, D. (1984b). Ultramicroscopy. 15,287. Coene, W., Van Dyck, D., Van Tendeloo, G. and Van Landuyt, J. (1985) Phil. Mag. A. In press. Cowley, J. M., and Moodie, A. F. (1951). Proc. Phys. Soc. B 70,486. Cowley, J. M., and Moodie, A. F. (1957). Acta Crystallogr. 10,609. Doyle, P. A,, and Turner, P. S. (1968). Acta Crystallogr. Sec. A 24,390. Feynman, R. P., and Hibbs, A. R. (1965). “Quantum Mechanics and Path Integrals,” McGrawHill, New York. Fields, P. M., and Cowley, J. M. (1978). Acta Crystallogr. See. A 34, 103. Fujimoto, F. (1959). J. Phys. Soc. Jpn. 14,1558. Fujiwara, K. (1959). J. Phys. Soc. Jpn. 14, 1513. Fujiwara, K. (1961). J. Phys. Soc. Jpn. 16,2226. Gevers, R., Serneels, R., and David, M. (1974). Phys. Status Solidi B 66,471. Goodman, P., and Moodie,A. F. (1974). Acta Crystallogr. See. A 30,280. Gratias, D., and Portier, R. (1983). Acta Crystallogr. Sec. A 39, 576. Grinton, G. R., and Cowley, J. M. (1971). Optik 34,221. Howie, A., (1978). in “Diffraction and Imaging Techniques in Material Science,” Vol 2 (S. Amelinckx. R. Gevers, and J. Van Landuyt, eds.), p. 457. North-Holland Publ., Amsterdam. Howie, A. (1983). J. Microsc. 129,239. Howie, A., and Basinsky, Z. S. (1968). Philos. Mag. 17 1039.

354

D. VAN DYCK

Howie, A., and Whelan, M. J. (1961). Proc. R. SOC.(London) Ser. A 263,217. Humphries, C. J. (1980). Proc. 6th Int. ConJ HVEM Antwerp, (P. Brederoo and J. Van Linduyt, eds.) p. 68. Humphries, C. J., and Hirsch, P. B. (1968). Philos. Mag. 18, 115. Humphries, C. J., and Spence, J. C. H. (1979), Proc. 37th E M S A Meet. Baton Rouge, Louisiana, p. 554. Claitor’s Publishing Division. Ishizuka, K. (1982). Acra Crystallogr. 138,773. Ishizuka, K., and Uyeda, N. (1977). Acta Crystallogr. See. A 33,740. Kilaas, R.,and Gronsky, R. (19x3). Ultramicroscopy 11,289. Kirkland, E. J. (1982). Ul/ramicroscopy 9,45,65. Lynch, D. F. (1971). Acta Crystallogr. See. A 27,399. Lynch, D. F. (1982). Proc. Austr. Conf: Electron Microsc., 7th. Lynch, D. F., and OKeefe, M. A. (1972). Acta Crystalogr. See. A 28,536. Marks, L. D. (1984). Ultramicroscopy 12,237. Matsuhata, H., Van Dyck, D. and Coene, W. (1983a). Proc. Joinr Meet. Belg. Ger. SOC.Electron Microsc. Antwerp, p. 135. Matsuhata, H., Van Dyck, D., and Coene, W. (1983b), Proc. Joint Meet. Belg. Ger. Soc. Electron Microsc., Antwerp, p. 136. Matsuhata, H., Van Dyck, D., Van Landuyt, J., and Amelinckx, S. (1984). Ultramicroscopy 13, 343.

Moodie, A. F. (1981). In “Fifty Years of Electron Diffraction,” (P. Goodman, ed.), p. 327. D. Reidel, Dordrecht. Niehrs, H., and Wagner, F. N. (1955). 2. Phys. 143,285. Nihoul, G., and Cesari, C. (1984). Phys. Status Solid; A 81.87. OKeefe, M. A. (1973). Acta Crystallogr. See. A 29,389. Olsen, A., and Spence, J. C. H. (1981). Philos. Mag. A 43, 945. Pirouz, P. (1981). Acta Crystallogr. See. A 37,465. Radi, G. (1970). Acta Crystallogr. Sec. A 26.41. Rez, P., (1978). Inst. Phys. ConJ Ser. 41,61. Rez, P. (1980). Proc. Annu. Meet., Electron Microsc, SOC.Am. p. 180. Sanders, J. V., and Goodman, P. (1981). In “Fifty Years of Electron Diffraction,” (P. Goodman, ed.), p. 281. D. Reidel, Dordrecht. Saxton, W. 0. (1978). In “Advances in Electronics and Electron Physics,” Supplement 10 Academic Press, New York. Saxton, W. 0.and Stobbs, W. M. (1984). Proc. EUREM 8, Budapest 1,287. Scherzer, 0. (1949). J. Appl. Phys. 20,20. Schiske, P. (1950)). Ph.D. Thesis, University of Vienna. Schryvers, D., Van Landuyt, J., and Amelinckx, S. (1984). Proc. EUREM 8, Budapest I , 891. Self, P. G. (1982). J. Microsc. 127,293. Self, P G., OKeefe, M. A., Buseck, P. R., Spargo, A. E. C. (1983). Ultramicroscopy 11,35. Serneels, R., Haantjens, D., and Gevers, R. (1980). Pilos. Mag. A 42, 1. Shannon, M. D. (1978). Inst. Phys. Con$ Ser. 41,41. Skarnulis, A. J., Wild, D. L., Anstis, G. R.,Humphries, C. J., and Spence, J. C. H. (1981). Inst. Phys. ConJ Ser. 61,347. Spence, J. C. H. (1978). Acra Crystallogr. Sec. A 34, 112. Spence, J. C. H. (1981). “Experimental High-Resolution Electron Microscopy.” Clarendon Press, Oxford. Spence, J. C. H., and Whelan, M. J. (1975). Acta Crystallogr. See. A 31, 242. Spence, J. C. H., OKeefe, M., and Iilima, S. (1978). Phillos. Mag. A 38,463. Sturkey, L. (1957). Acta Crystallogr. 10, 858.

CALCULATIONS IN HIGH-RESOLUTION ELECTRON MICROSCOPY

355

Sturkey, L. (1962). Proc. Phys. Soc. 80, 321. Tournarie, M. (1960). Bull. Soc. Franc. Miner. Cryst. 83, 179. Tournarie, M. (1961). C. R. Acad. Sci. 252, 1961. Tournarie, M. (1962). J. Phys. Soc. Japan 17, Suppl. B. 11. 98. Van Dyck, D. ( 1975) Phys. Status Solidi B 72,32 1. Van Dyck, D. (1976) Phys. S/atus Solidi B 77,301. Van Dyck, D. (1978a). In “Diffraction and Imaging Techniques in Material Science,” Vol. I , (S. Amelinckx, R. Gevers, and J. Van Landuyt, eds.), p. 355. North-Holland Publ., Amsterdam. Van Dyck, D. (1978b). Proc. Inter. Congr. Electron Microsc., 9th, Toronto, Vol. 1, (J. M. Sturgess, ed.) p. 196. Van Dyck, D. (1979). Phys. Status. Solidi A 52,283. Van Dyck, D. (1980).J. Microsc. 119, 141. Van Dyck, D. (1983). J. Micrsoc. 132, 31. Van Dyck, D. (1984). Proc. EUREM 8, Budapest 1,261. Van Dyck, D., and Coene, W. (1984). Ultramicroscopy. 15.29. Van Dyck, D., DeRidder, R., and DeSitter, J. (1980). J. Comp. Appl. Math. 6(1), 83. Van Dyck, D., Van Tendeloo, G., and Amelinckx, S. (1982). Ulrrumicroscopy 10,263. Van Sande, M., Van Tendeloo, G., Amelinckx, S., and Airo, P. (1979). Phys. Sfatus Solidi A 54, 499. Van Hugo, D., Kohl, H., and Rose, H. (1984). Proc. EUREM8, Budapest 1,135. Van Tendeloo, G., and Amelinckx, S. (1978a). Phys. Slat. Solidi A 47, 555. Van Tendeloo, G., and Amelinckx, S. (1978b). Phys. Stat. Solidi A 49, 337. Whelan, M. J. (1984). Proc. EUREM 8, Budapest 1,429. Wilson, A. R., and Spargo, A. E. (1982). Philos. Mag. A 46, 435. Wood, G. J. (1982). Private communication. Yagi, K. (1984). Proc. I n / . Congr. Crysrullogr. 13th. p. C2. Yoshioka, H. (1957). J. Phys. Soc. Jpn. 12,618.