Computation of dynamical wave amplitudes

Colliex, C.; Cowley, J. M.; Dudarev, S. L.; Fink, M.; Gjønnes, J.; Hilderbrandt, R.; Howie, A.; Lynch, D. F.; Peng, L. M.; Ren, G.; Ross, A. W.; Smith, V. H. Jr; Spence, J. C. H.; Steeds, J. W.; Wang, J.; Whelan, M. J.; Zvyagin, B. B.

doi:10.1107/97809553602060000593

International
Tables for
Crystallography
Volume C
Mathematical, physical and chemical tables
Edited by E. Prince

pdf | chapter contents | chapter index | related articles

International Tables for Crystallography (2006). Vol. C. ch. 4.3, pp. 414-416

Section 4.3.6. Computation of dynamical wave amplitudes

C. Colliex,^a J. M. Cowley,^b ^‡ S. L. Dudarev,^c M. Fink,^d J. Gjønnes,^e R. Hilderbrandt,^f A. Howie,^g D. F. Lynch,^h L. M. Peng,ⁱ G. Ren,^j A. W. Ross,^d V. H. Smith Jr,^k J. C. H. Spence,^l J. W. Steeds,^m J. Wang,^k M. J. Whelan^c and B. B. Zvyaginⁿ ^‡

^a Laboratoire Aimé Cotton, CNRS, Campus d'Orsay, Bâtiment 505, F-91405 Orsay CEDEX, France,^bDepartment of Physics and Astronomy, Arizona State University, Tempe, AZ 85287-1504, USA,^cDepartment of Materials, University of Oxford, Parks Road, Oxford OX1 3PH, England,^dDepartment of Physics, The University of Texas at Austin, Austin, TX 78712, USA,^eDepartment of Physics, University of Oslo, PO Box 1048, Blindern, N-0316 Oslo, Norway,^fChemistry Division, Room 1055, The National Science Foundation, 4201 Wilson Blvd, Arlington, VA 22230, USA,^gCavendish Laboratory, Madingley Road, Cambridge CB3 0HE, England,^hCSIRO Division of Materials Science & Technology, Private Bag 33, Rosebank MDC, Clayton, Victoria 3169, Australia,ⁱDepartment of Electronics, Peking University, Beijing 100817, People's Republic of China,^jBeijing Laboratory of Electron Microscopy, Chinese Academy of Sciences, PO Box 2724, Beijing 100080, People's Republic of China,^kDepartment of Chemistry, Queen's University, Kingston, Ontario, Canada K7L 3N6,^lDepartment of Physics, Arizona State University, Tempe, AZ 85287, USA,^mH. H. Wills Physics Laboratory, University of Bristol, Tyndall Avenue, Bristol BS8 1TL, England, and ⁿInstitute of Ore Mineralogy, Akad. Nauk Russia, Staromonetny 35, 109017 Moscow, Russia

4.3.6. Computation of dynamical wave amplitudes

| top | pdf |

4.3.6.1. The multislice method

| top | pdf |

The calculation of very large numbers of diffracted orders, i.e. more than 100 and often several thousand, requires the multislice procedure. This occurs because, for N diffracted orders, the multislice procedure involves the manipulation of arrays of size N, whereas the scattering matrix or the eigenvalue procedures involve manipulation of arrays of size N by N.

The simplest form of the multislice procedure presumes that the specimen is a parallel-sided plate. The surface normal is usually taken to be the z axis and the crystal structure axes are often chosen or transformed such that the c axis is parallel to z and the a and b axes are in the xy plane. This can often lead to rather unconventional choices for the unit-cell parameters. The maximum tilt of the incident beam from the surface normal is restricted to be of the order of 0.1 rad. For the calculation of wave amplitudes for larger tilts, the structure must be reprojected down an axis close to the incident-beam direction. For simple calculations, other crystal shapes are generally treated by the column approximation, that is the crystal is presumed to consist of columns parallel to the z axis, each column of different height and tilt in order to approximate the desired shape and variation of orientation.

The numerical procedure involves calculation of the transmission function through a thin slice, calculation of the vacuum propagation between centres of neighbouring slices, followed by evaluation in a computer of the iterated equation $[u_n(h,k)=p_n\{p_{n-1}\ldots p_3[\,p_2(p_1q_1*q_2)*q_3]*\ldots *q_n\}\eqno (4.3.6.1)]$ in order to obtain the scattered wavefunction, [u_n(h,k)] , emitted from slice n, i.e. for crystal thickness $[H=\Delta z_1+]$ $[\Delta z_2+\ldots+\Delta z_n]$ ; the symbol * indicates the operation `convolution' defined by $[f_1(x)*f_2(x)=\textstyle\int\limits^\infty_{-\infty}f_1(w)f_2(x-w)\,{\rm d} w,]$ and $[p_n=\exp\big(\!-i2\pi\Delta z_n(\lambda/2)\{[h(h-h'')/a^2]+[k(k-k'')/b^2]\}\big)]$ is the propagation function in the small-angle approximation between slice n − 1 and slice n over the slice spacing $[\Delta z_n]$ . For simplicity, the equation is given for orthogonal axes and h′′, k′′ are the usually non-integral intercepts of the Laue circle on the reciprocal-space axes in units of (1/a), (1/b). The excitation errors, $[\zeta(h,k)]$ , can be evaluated using $[\zeta(h,k)=-(\lambda/2)\{[h(h-h'')/a^2]+[k(k-k'')/b^2]\}. \eqno (4.3.6.2)]$ The transmission function for slice n is $[q_n(h,k)=F\{\exp[i\sigma\varphi_n(x,y)\Delta z_n]\}, \eqno (4.3.6.3)]$ where F denotes Fourier transformation from real to reciprocal space, and $[\varphi_n(x,y)\Delta z_n= {^p}\varphi(x,y)=\textstyle\int\limits^{z_{n-1}+\Delta z_n}_{z_{n-1}} \varphi(x,y,z)\,{\rm d} z]$ and $[\sigma={\pi\over W\lambda}\,{2\over{1+(1-\beta^2)^{1/2}}}]$ $[\beta={v\over c},]$ where W is the beam voltage, v is the relativistic velocity of the electron, c is the velocity of light, and λ is the relativistic wavelength of the electron.

The operation * in (4.3.6.1) is most effectively carried out for large N by the use of the convolution theorem of Fourier transformations. This efficiency presumes that there is available an efficient fast-Fourier-transform subroutine that is suitable for crystallographic computing, that is, that contains the usual crystallographic normalization factors and that can deal with a range of values for h, k that go from negative to positive. Then, $[u_n(h,k)=F\{F^{-1}[u_{n-1}(h,k)]F^{-1}[q_n(h,k)]\}, \eqno (4.3.6.4)]$ where F denotes $[u(h,k)={1\over n_x n_y}\; \sum^{n_x}_{x=1} \sum^{n_y}_{y=1} U(x,y)\exp\left\{2\pi i\left[\displaystyle{hx\over n_x},{ky\over n_y}\right]\right\}]$ and $[F^{-1}]$ denotes $[U(x,y) ={\sum^{n_h}_{h=-n_h}} \sum^{n_k}_{k=-n_k} u(h,k)\exp\left\{2\pi i\left[\displaystyle{hx\over n_x},{ky\over n_y}\right]\right\},]$ where [n_h=(n_x/2)-1] , [n_k=(n_y/2)-1] , and [n_x,n_y] are the sampling intervals in the unit cell. The array sizes used in the calculations of the Fourier transforms are commonly powers of 2 as is required by many fast Fourier subroutines. The array for [u_n(h,k)] is usually defined over the central portion of the reserved computer array space in order to avoid oscillation in the Fourier transforms (Gibbs instability). It is usual to carry out a $[64\times64]$ beam calculation in an array of $[128\times128]$ , hence the critical timing interval in a multislice calculation is that interval taken by a fast Fourier transform for 4N coefficients. If the number of beams, N, is such that there is still appreciable intensity being scattered outside the calculation aperture, then it is usually necessary to impose a circular aperture on the calculation in order to prevent the symmetry of the calculation aperture imposing itself on the calculated wavefunction. This is most conveniently achieved by setting all p(h, k) coefficients outside the desired circular aperture to zero.

It is clear that the iterative procedure of (4.3.6.1) means that care must be taken to avoid accumulation of error due to the precision of representation of numbers in the computer that is to be used. Practical experience indicates that a precision of nine significant figures (decimal) is more than adequate for most calculations. A precision of six to seven (decimal) figures (a common 32-bit floating-point representation) is only barely satisfactory. A computer that uses one of the common 64-bit representations (12 to 16 significant figures) is satisfactory even for the largest calculations currently contemplated.

The choice of slice thickness depends upon the maximum value of the projected potential within a slice and upon the validity of separation of the calculation into transmission function and propagation function. The second criterion is not severe and in practice sets an upper limit to slice thickness of about 10 Å. The first criterion depends upon the atomic number of atoms in the trial structure. In practice, the slice thickness will be too large if two atoms of medium to heavy atomic weight $[(Z\ge30)]$ are projected onto one another. It is not necessary to take slices less than one atomic diameter for calculations for fast electron (acceleration voltages greater than 50 keV) diffraction or microscopy. If the trial structure is such that the symmetry of the diffraction pattern is not strongly dependent upon the structure of the crystal parallel to the slice normal, then the slices may be all identical and there is no requirement to have a slice thickness related to the periodicity of the structure parallel to the surface normal. This is called the `no upper-layer-line' approximation . If the upper-layer lines are important, then the slice thickness will need to be a discrete fraction of the c axis, and the contents of each slice will need to reflect the actual atomic contents of each slice. Hence, if there were four slices per unit cell, then there would need to be four distinct q(h, k), each taken in the appropriate order as the multislice operation proceeds in thickness.

The multislice procedure has two checks that can be readily performed during a calculation. The first is applied to the transmission function, q(h, k), and involves the evaluation of a unitarity test by calculation of $[\textstyle\sum\limits_{h'} \sum\limits_{k'}\,q(h',k')q^*(h+h',k+k')=\delta(h,k) \eqno (4.3.6.5)]$ for all h, k, where q* denotes the complex conjugate of q, and δ(h,k) is the Kronecker delta function. The second test can be applied to any calculation for which no phenomenological absorption potential has been used in the evaluation of the q(h, k). In that case, the sum of intensities of all beams at the final thickness should be no less than 0.9, the incident intensity being taken as 1.0. A value of this sum that is less than 0.9 indicates that the number of beams, N, has been insufficient. In some rare cases, the sum can be greater than 1.0; this is usually an indication that the number of beams has been allowed to come very close to the array size used in the convolution procedure. This last result does not occur if the convolution is carried out directly rather than by use of fast-Fourier-transform methods.

A more complete discussion of the multislice procedure can be obtained from Cowley (1975) and Goodman & Moodie (1974). These references are not exhaustive, but rather an indication of particularly useful articles for the novice in this subject.

4.3.6.2. The Bloch-wave method

| top | pdf |

Bloch waves, familiar in solid-state valence-band theory, arise as the basic wave solutions for a periodic structure. They are thus always implicit and often explicit in dynamical diffraction calculations, whether applied in perfect crystals, in almost perfect crystals with slowly varying defect strain fields or in more general structures that (see Subsection 4.3.6.1) can always, for computations, be treated by periodic continuation.

The Schrödinger wave equation in a periodic structure, $[\nabla^2\psi+4\pi^2\bigg[\chi^2+\textstyle\sum\limits_{\bf g} \,U_{\bf g} \exp (2\pi i{\bf g}\cdot{\bf r})\bigg]\psi=0, \eqno (4.3.6.6)]$ can be applied to high-energy, relativistic electron diffraction, taking $[\chi=\lambda^{-1}]$ as the relativistically corrected electron wave number (see Subsection 4.3.1.4). The Fourier coefficients in the expression for the periodic potential are defined at reciprocal-lattice points g by the expression $[U_{\bf g}=U^*_{-{\bf g}}={m\over m_0} {\exp(-M_{\bf g}) \over\pi\Omega} \sum_j \, f_j[\sin(\theta_{\bf g})/\lambda]\exp(-2\pi i{\bf g}\cdot{\bf r}_j), \eqno (4.3.6.7)]$ where [f_j] is the Born scattering amplitude (see Subsection 4.3.1.2) of the jth atom at position $[{\bf r}_j]$ in the unit cell of volume $[\Omega]$ and $[M_{\bf g}]$ is the Debye–Waller factor.

The simplest solution to (4.3.6.6) is a single Bloch wave, consisting of a linear combination of plane-wave beams coupled by Bragg reflection. $[\psi({\bf r})=b({\bf k},{\bf r})=\textstyle\sum\limits_h C_{\bf h}\exp[2\pi i({\bf k}+{\bf h})\cdot{\bf r}]. \eqno (4.3.6.8)]$ In practice, only a limited number of terms N, corresponding to the most strongly excited Bragg beams, is included in (4.3.6.8). Substitution in (4.3.6.6) then yields N simultaneous equations for the wave amplitudes $[C_{\bf g}.]$ $[[\chi^2+U_0-({\bf k}+{\bf g})^2]C_{\bf g}+\textstyle\sum\limits_{\bf g'\ne0}U_{\bf g'}C_{\bf g-{\bf g}'}=0. \eqno (4.3.6.9)]$ Usually, χ and the two tangential components [k_x] and [k_y] are fixed by matching to the incident wave at the crystal entrance surface. [k_z] then emerges as a root of the determinant of coefficients appearing in (4.3.6.9).

Numerical solution of (4.3.6.9) is considerably simplified (Hirsch, Howie, Nicholson, Pashley & Whelan, 1977a) in cases of transmission high-energy electron diffraction where all the important reciprocal-lattice points lie in the zero-order Laue zone [g_z=0] and $[\chi^2\gg |U_{\bf g}|]$ . The equations then reduce to a standard matrix eigenvalue problem (for which efficient subroutines are widely available): $[\textstyle\sum\limits_{\bf h} M_{\bf gh}C_{\bf h}=\gamma C_{\bf g}, \eqno (4.3.6.10)]$ where $[M_{\bf gh}=U_{\bf g-h}/2\chi+s_{\bf g}\delta_{\bf gh}]$ and $[s_{\bf g}=[k^2-({\bf k}+{\bf g})^2]/2\chi]$ is the distance, measured in the z direction, of the reciprocal-lattice point g from the Ewald sphere.

There will in general be N distinct eigenvalues $[\gamma=k_z-\chi_z]$ corresponding to N possible values $[k^{(\,j)}_z]$ , $[j=1,2,\cdots N]$ , each with its eigenfunction defined by N wave amplitudes $[C^{(\,j)}_0, C^{(\,j)}_{\bf g},\ldots, C^{(\,j)}_{\bf h}]$ . The waves are normalized and orthogonal so that $[\textstyle\sum\limits_{\bf g} C^{(\,j)}_{\bf g}{^*} C^{(l)}_{\bf g}=\delta_{jl};\quad\textstyle \sum\limits_jC_{\bf g}^{(\,j)*}C^{(l)}_{\bf h}=\delta_{\bf gh}. \eqno (4.3.6.11)]$ In simple transmission geometry, the complete solution for the total coherent wavefunction $[\psi({\bf r})]$ is $[\psi({\bf r})=\textstyle\sum\limits_j \psi^{(\,j)}\exp[-2\pi q^{(\,j)}z]\sum\limits_{\bf g}C^{(\,j)}_{\bf g}\exp [2\pi i(\chi+{\bf g})\cdot {\bf r}]. \eqno (4.3.6.12)]$ Inelastic and thermal-diffuse-scattering processes cause anomalous absorption effects whereby the amplitude of each component Bloch wave decays with depth z in the crystal from its initial value $[\psi^{(\,j)}=C^{(\,j)*}_0]$ . The decay constant is computed using an imaginary optical potential $[iU'({\bf r})]$ with Fourier coefficients $[iU'_{\bf g}=iU'^{*}_{-{\bf g}}]$ (for further details of these see Humphreys & Hirsch, 1968, and Subsection 4.3.1.5 and Section 4.3.2). $[q^{(\,j)} = {m\over h^2\chi_z} \sum_{\bf g,h} C^{(\,j)*}_{\bf g} U'_{\bf h}C^{(\,j)}_{\bf g-h}. \eqno (4.3.6.13)]$ The Bloch-wave, matrix-diagonalization method has been extended to include reciprocal-lattice points in higher-order Laue zones (Jones, Rackham & Steeds, 1977) and, using pseudopotential scattering amplitudes, to the case of low-energy electrons (Pendry, 1974).

The Bloch-wave picture may be compared with other variants of dynamical diffraction theory, which, like the multislice method (Subsection 4.3.6.1), for example, employ plane waves whose amplitudes vary with position in real space and are determined by numerical integration of first-order coupled differential equations. For cases with $[N\lt50]$ beams in perfect crystals or in crystals containing localized defects such as stacking faults or small point-defect clusters, the Bloch-wave method offers many advantages, particularly in thicker crystals with t > 1000 Å. For high-resolution image calculations in thin crystals where the periodic continuation process may lead to several hundred diffracted beams, the multislice method is more efficient. For cases of defects with extended strain fields or crystals illuminated at oblique incidence, coupled plane-wave integrations along columns in real space (Howie & Basinski, 1968) can be the most efficient method.

The general advantage of the Bloch-wave method, however, is the picture it affords of wave propagation and scattering in both perfect and imperfect crystals. For this purpose, solutions of equations (4.3.6.9) allow dispersion surfaces to be plotted in k space, covering with several sheets j all the wave points $[{\bf k}^{(\,j)}]$ for a given energy E. Thickness fringes and other interference effects then arise because of interference between waves excited at different points $[{\bf k}^{(\,j)}]$ . The average current flow at each point is normal to the dispersion surface and anomalous-absorption effects can be understood in terms of the distribution of Bloch-wave current within the unit cell. Detailed study of these effects, and the behaviour of dispersion surfaces as a function of energy, yields accurate data on scattering amplitudes via the critical-voltage effect (see Section 4.3.7). Static crystal defects induce elastic scattering transitions $[{\bf k}^{(\,j)}\rightarrow{\bf k}^{(l)}]$ on sheets of the same dispersion surface. Transitions between points on dispersion surfaces of different energies occur because of thermal diffuse scattering, generation of electronic excitations or the emission of radiation by the fast electron. The Bloch-wave picture and the dispersion surface are central to any description of these phenomena. For further information and references, the reader may find it helpful to consult Section 5.2.10 of Volume B (IT B, 2001).

References

Cowley, J. M. (1975). Diffraction physics. New York: North-Holland.Google Scholar

Goodman, P. & Moodie, A. F. (1974). Numerical evaluation of N-beam wave-functions in electron scattering by the multislice method. Acta Cryst. A30, 280–290.Google Scholar

Hirsch, P. B., Howie, A., Nicholson, R. B., Pashley, D. W. & Whelan, M. J. (1977a). Electron microscopy of thin crystals. New York: Krieger.Google Scholar

Howie, A. & Basinski, Z. S. (1968). Approximations of the dynamical theory of diffraction contrast. Philos. Mag. 17, 1039–1063.Google Scholar

Humphreys, C. J. & Hirsch, P. B. (1968). Absorption parameters in electron diffraction theory. Philos. Mag. 18, 115–122.Google Scholar

International Tables for Crystallography (2001). Vol. B, 2nd ed. Dordrecht: Kluwer Academic Publishers.Google Scholar

Jones, P. M., Rackham, G. M. & Steeds, J. W. (1977). High-order Laue zone electron diffraction. Proc. R. Soc. London Ser. A, 354, 192–222.Google Scholar

Pendry, J. B. (1974). Low energy electron diffraction. New York: Academic Press.Google Scholar

International Tables for Crystallography (2006). Vol. C. ch. 4.3, pp. 414-416