Multislice

The multislice algorithm^[1] is a method for the simulation of the elastic scattering of an electron beam with matter, including all multiple scattering effects. The method is reviewed in the book by John M. Cowley,^[2] and also the work by Ishizuka.^[3] The algorithm is used in the simulation of high resolution transmission electron microscopy (HREM) micrographs, and serves as a useful tool for analyzing experimental images.^[4] This article describes some relevant background information, the theoretical basis of the technique, approximations used, and several software packages that implement this technique. Some of the advantages and limitations of the technique and important considerations that need to be taken into account are described.

The multislice method has found wide application in electron microscopy and crystallography. The mapping from a crystal structure to its image or electron diffraction pattern is relatively well understood and documented. However, the reverse mapping from electron micrograph images to the crystal structure is generally more complicated. The fact that the images are two-dimensional projections of three-dimensional crystal structure makes it tedious to compare these projections to all plausible crystal structures. Hence, the use of numerical techniques in simulating results for different crystal structure is integral to the field of electron microscopy and crystallography. Several software packages exist to simulate electron micrographs.

There are two widely used simulation techniques that exist in literature: the Bloch wave method,^[5] derived from Hans Bethe's original theoretical treatment,^[6] and the multislice method. This article focuses on the multislice method for simulation of dynamical diffraction, including multiple elastic scattering effects. Most of the packages that exist implement the multislice algorithm along with Fourier analysis to incorporate electron lens aberration effects to determine electron microscope image and address aspects such as phase contrast and diffraction contrast. For electron microscope samples in the form of a thin crystalline slab in the transmission geometry, the aim of these software packages is to provide a map of the crystal potential, however this inversion process is greatly complicated by the presence of multiple elastic scattering.

The first description of what is now known as the multislice theory was given in the classic paper by Cowley and Moodie.^[1] In this work, the authors describe scattering of electrons using a physical optics approach without invoking quantum mechanical arguments. Many other derivations of these iterative equations have since been given using alternative methods, such as Greens functions, differential equations, scattering matrices or path integral methods, see for instance the book by Lianmao Peng, Sergei Dudarev and Michael Whelan.^[7]

A summary of the development of a computer algorithm from the multislice theory of Cowley and Moodie for numerical computation was reported by Goodman and Moodie.^[8] They also discussed in detail the relationship of the multislice to the other formulations. Specifically, using Zassenhaus's theorem, this paper gives the mathematical path from multislice to 1. Schrödinger equation, 2. Darwin's differential equations, widely used for diffraction contrast Transmission electron microscopy (TEM) image simulations - the Howie-Whelan equations,^[9] 3. Sturkey's scattering matrix method.^[10] 4. the free-space propagation case, 5. The phase grating approximation, 6. A new "thick-phase grating" approximation, which has never been used, 7. Moodie's polynomial expression for multiple scattering, 8. The Feynman path-integral formulation, and 9. relationship of multislice to the Born series. The relationship between algorithms is summarized in Section 5.11 of Spence (2013),^[11] (see Figure 5.9).

The form of multislice algorithm presented here has been adapted from Peng, Dudarev and Whelan 2003.^[7] The multislice algorithm is an approach to solving the Schrödinger equation:

${\displaystyle {\begin{aligned}-{\frac {\hbar ^{2}}{2m}}{\frac {\partial ^{2}\Psi (x,t)}{\partial x^{2}}}+V(x,t)\Psi (x,t)&=E\Psi (x,t)\end{aligned}}}$

In 1957, Cowley and Moodie showed that the Schrödinger equation can be solved analytically to evaluate the amplitudes of diffracted beams.^[1] Subsequently, the effects of dynamical diffraction can be calculated and the resulting simulated image will exhibit good similarities with the actual image taken from a microscope under dynamical conditions. Furthermore, the multislice algorithm does not make any assumption about the periodicity of the structure and can thus be used to simulate HREM images of aperiodic systems as well.

The following section will include a mathematical formulation of the multislice algorithm. The Schrödinger equation can also be represented in the form of incident and scattered wave as:

${\displaystyle {\begin{aligned}\Psi ({\mathbf {r} })&=\Psi _{0}({\mathbf {r} })+\int {G({\mathbf {r,r'} })V({\mathbf {r'} })\Psi ({\mathbf {r'} })d{\mathbf {r'} }}\end{aligned}}}$

where ${\displaystyle G(\mathbf {r,r'} )}$ is the Green's function that represents the amplitude of the electron wave function at a point ${\displaystyle \mathbf {r} }$ due to a source at point ${\displaystyle \mathbf {r'} }$.

Hence for an incident plane wave of the form ${\displaystyle \Psi (r)=\exp(i\mathbf {k\cdot r} )}$ the Schrödinger equation can be written as

${\displaystyle {\begin{aligned}\Psi ({\mathbf {r} })=\exp(i{\mathbf {k\cdot r} })-{\frac {m}{2\pi \hbar ^{2}}}\int {\frac {\exp(ik\cdot {\mathbf {|r-r'|} })}{\mathbf {|r-r'|} }}V({\mathbf {r'} })\Psi ({\mathbf {r'} })dr'.\end{aligned}}}$

(1)

We then choose the coordinate axis in such a way that the incident beam hits the sample at (0,0,0) in the ${\displaystyle {\hat {z}}}$-direction, i.e., ${\textstyle \mathbf {k} =(0,0,k)}$. Now we consider a wave-function ${\displaystyle \Psi (r)=\phi (\mathbf {r} )\exp(i\mathbf {k\cdot r} )}$ with a modulation function ${\displaystyle \phi ({\mathbf {r} })}$ for the amplitude. Equation (1) becomes then an equation for the modulation function, i.e.,

${\displaystyle {\begin{aligned}\phi ({\mathbf {r} })&=1-{\frac {m}{2\pi \hbar ^{2}}}\int {{\frac {\exp[ik|{\mathbf {r-r'} }|-i{\mathbf {k} }\cdot ({\mathbf {r-r'} })]}{|{\mathbf {r-r'} }|}}V({\mathbf {r'} )\phi ({\mathbf {r'} })}dr'}\end{aligned}}}$.

Now we make substitutions with regards to the coordinate system we have adhered, i.e.,

${\displaystyle {\begin{aligned}{\mathbf {k} }\cdot ({\mathbf {r-r'} })&=k(z-z')\\|{\mathbf {r-r'} }|&\approx (z-z')+({\mathbf {X-X'} })^{2}/{2(z-z')}\end{aligned}}}$

where ${\displaystyle {\boldsymbol {X}}={\begin{pmatrix}x\\y\end{pmatrix}}}$.

Thus

${\displaystyle {\begin{aligned}\phi ({\mathbf {r} })=1-i{\frac {\pi }{E\lambda }}\int \int \limits _{z'=-\infty }^{z'=z}V({\mathbf {X'} },z')\phi ({\mathbf {X'} },z'){\frac {1}{i\lambda (z-z')}}\exp \left(ik{\frac {|{\mathbf {X-X'} }|^{2}}{2(z-z')}}\right)d{\mathbf {X'} }dz'\end{aligned}}}$,

where ${\displaystyle \lambda =2\pi /k}$ is the wavelength of the electrons with energy ${\displaystyle E=\hbar ^{2}k^{2}/{2m}}$ and ${\displaystyle {\begin{aligned}\sigma =\pi /E\lambda \end{aligned}}}$ is the interaction constant. So far we have set up the mathematical formulation of wave mechanics without addressing the scattering in the material. Further we need to address the transverse spread, which is done in terms of the Fresnel propagation function

${\displaystyle {\begin{aligned}p({\mathbf {X} },z)={\frac {1}{iz\lambda }}\exp \left(ik{\frac {{\mathbf {X} }^{2}}{2z}}\right)\end{aligned}}}$.

The thickness of each slice over which the iteration is performed is usually small and as a result within a slice the potential field can be approximated to be constant ${\displaystyle V({\mathbf {X'} },z)}$. Subsequently, the modulation function can be represented as:

${\displaystyle {\begin{aligned}\phi ({\mathbf {X} },z_{n+1})=\int p({\mathbf {X} }-{\mathbf {X'} },z_{n+1}-z_{n})\phi ({\mathbf {X} },z_{n})\exp \left(-i\sigma \int \limits _{z_{n}}^{z_{n+1}}V({\mathbf {X'} },z')dz'\right)dX'\end{aligned}}}$

We can therefore represent the modulation function in the next slice

${\displaystyle {\begin{aligned}\phi _{n+1}=\phi ({\mathbf {X} },z_{n+1})=[q_{n}\phi _{n}]*p_{n}\end{aligned}}}$

where, * represents convolution, ${\displaystyle p_{n}=p({\mathbf {X} },z_{n+1}-z_{n})}$ and ${\displaystyle q_{n}({\mathbf {X} })}$ defines the transmission function of the slice.

${\displaystyle {\begin{aligned}q_{n}({\mathbf {X} })=\exp\{-i\sigma \int \limits _{z_{n}}^{z_{n+1}}V({\mathbf {X} },z')dz'\}\end{aligned}}}$

Hence, the iterative application of the aforementioned procedure will provide a full interpretation of the sample in context. Further, it should be reiterated that no assumptions have been made on the periodicity of the sample apart from assuming that the potential ${\displaystyle V(\mathbf {X} ,z)}$ is uniform within the slice. As a result, it is evident that this method in principle will work for any system. However, for aperiodic systems in which the potential will vary rapidly along the beam direction, the slice thickness has to be significantly small and hence will result in higher computational expense.

Table 1 - Computational efficiency of Discrete Fourier Transform compared to Fast Fourier Transform

Data Points	${\displaystyle \mathbf {log_{2}} }$N	Discrete FT	Fast FT	Ratio
64	6	4,096	384	10.7
128	7	16,384	896	18.3
256	8	65,536	2,048	32
512	9	262,144	4,608	56.9
1,024	10	1,048,576	10,240	102.4
2,048	11	4,194,304	22,528	186.2

The basic premise is to calculate diffraction from each layer of atoms using fast Fourier transforms (FFT) and multiplying each by a phase grating term. The wave is then multiplied by a propagator, inverse Fourier transformed, multiplied by a phase grating term yet again, and the process is repeated. The use of FFTs allows a significant computational advantage over the Bloch Wave method in particular, since the FFT algorithm involves ${\displaystyle N\log N}$ steps compared to the diagonalization problem of the Bloch wave solution which scales as ${\displaystyle N^{2}}$ where ${\displaystyle N}$ is the number of atoms in the system. (See Table 1 for comparison of computational time).

The most important step in performing a multislice calculation is setting up the unit cell and determining an appropriate slice thickness. In general, the unit cell used for simulating images will be different from the unit cell that defines the crystal structure of a particular material. The primary reason for this due to aliasing effects which occur due to wraparound errors in FFT calculations. The requirement is to add additional “padding” to the unit cell has earned the nomenclature “super cell” and the requirement to add these additional pixels to the basic unit cell comes at a computational price.

To illustrate the effect of choosing a slice thickness that is too thin, consider a simple example. The Fresnel propagator describes the propagation of electron waves in the z direction (the direction of the incident beam) in a solid:

${\displaystyle {\tilde {\phi }}(\mathbf {u} ,z)={\tilde {\phi }}(\mathbf {u} ,z=0)\exp(\pi i\lambda \mathbf {u} ^{2}z)}$

Where ${\displaystyle \mathbf {u} }$ is the reciprocal lattice coordinate, z is the depth in the sample, and ${\displaystyle \lambda }$ is the wavelength of the electron wave (related to the wave vector by the relation ${\displaystyle k=2\pi /\lambda }$). In the case of the small-angle approximation (${\displaystyle \theta \sim }$ 100 mRad) we can approximate the phase shift as ${\displaystyle \Delta z}$. For 100 mRad the error ${\displaystyle d-S}$ is on the order of 0.5% since ${\displaystyle \cos(0.1)=0.995}$. For small angles this approximation holds regardless of how many slices there are, although choosing a ${\displaystyle \Delta z}$ greater than the lattice parameter (or half the lattice parameter in the case of perovskites) for a multislice simulation would result in missing atoms that should be in the crystal potential.

Additional practical concerns are how to effectively include effects such as inelastic and diffuse scattering, quantized excitations (e.g. plasmons, phonons, excitons), etc. There was one code that took these things into consideration through a coherence function approach ^[12] called Yet Another Multislice (YAMS), but the code is no longer available either for download or purchase.

There are several software packages available to perform multislice simulations of images. Among these is NCEMSS, NUMIS, MacTempas, and Kirkland . Other programs exist but unfortunately many have not been maintained (e.g. SHRLI81 by Mike O’Keefe of Lawrence Berkeley National Lab and Cerius2 of Accerlys). A brief chronology of multislice codes is given in Table 2, although this is by no means exhaustive.

This software is developed by Earl Kirkland of Cornell University. This code is freely available as an interactive Java applet and as standalone code written in C/C++. The Java applet is ideal for a quick introduction and simulations under a basic incoherent linear imaging approximation. The ACEM code accompanies an excellent text of the same name by Kirkland which describes the background theory and computational techniques for simulating electron micrographs (including multislice) in detail. The main C/C++ routines use a command line interface (CLI) for automated batching of many simulation. The ACEM package also includes a graphical user interface that is more appropriate for beginners. The atomic scattering factors in ACEM are accurately characterized by a 12-parameter fit of Gaussians and Lorentzians to relativistic Hartree–Fock calculations.

This package was released from the National Center for High Resolution Electron Microscopy. This program uses a mouse-drive graphical user interface and is written by Roar Kilaas and Mike O’Keefe of Lawrence Berkeley National Laboratory. While the code is no longer developed, the program is available through the Electron Direct Methods (EDM) package written by Laurence D. Marks of Northwestern University. Debye-Waller factors can be included in as a parameter to account for diffuse scattering, although the accuracy is unclear (i.e. a good guess of the Debye-Waller factor is needed).

The Northwestern University Multislice and Imaging System (NUMIS) is a package is written by Laurence Marks of Northwestern University. It uses a command line interface (CLI) and is based on UNIX. A structure file must be provided as input in order to run use this code, which makes it ideal for advanced users. The NUMIS multislice programs use the conventional multislice algorithm by calculating the wavefunction of electrons at the bottom of a crystal and simulating the image taking into account various instrument-specific parameters including ${\displaystyle C_{s}}$ and convergence. This program is good to use if one already has structure files for a material that have been used in other calculations (for example, Density Functional Theory). These structure files can be used to general X-Ray structure factors which are then used as input for the PTBV routine in NUMIS. Microscope parameters can be changed through the MICROVB routine.

This software is specifically developed to run in Mac OS X by Roar Kilaas of Lawrence Berkeley National Laboratory. It is designed to have a user-friendly user interface and has been well-maintained relative to many other codes (last update May 2013). It is available (for a fee) from here.

This is a software for multislice simulation was written in FORTRAN 77 by J. M. Zuo, while he was a postdoc research fellow at Arizona State University under the guidance of John C. H. Spence. The source code was published in the book of Electron Microdiffraction.^[13] A comparison between multislice and Bloch wave simulations for ZnTe was also published in the book. A separate comparison between several multislice algorithms at the year of 2000 was reported.^[14]

The Quantitative TEM/STEM (QSTEM) simulations software package was written by Christopher Koch of Humboldt University of Berlin in Germany. Allows simulation of HAADF, ADF, ABF-STEM, as well as conventional TEM and CBED. The executable and source code are available as a free download on the Koch group website.

This is a code written by Vincenzo Grillo of the Institute for Nanoscience (CNR) in Italy. This code is essentially a graphical frontend to the multislice code written by Kirkland, with more additional features. These include tools to generate complex crystalline structures, simulate HAADF images and model the STEM probe, as well as modeling of strain in materials. Tools for image analysis (e.g. GPA) and filtering are also available. The code is updated quite often with new features and a user mailing list is maintained. Freely available on their website.

Multi-slice image simulations for high-resolution scanning and coherent imaging transmission electron microscopy written by Juri Barthel from the Ernst Ruska-Centre at the Jülich Research Centre. The software comprises a graphical user interface version for direct visualization of STEM image calculations, as well as a bundle of command-line modules for more comprehensive calculation tasks. The programs have been written using Visual C++, Fortran 90, and Perl. Executable binaries for Microsoft Windows 32-bit and 64-bit operating systems are available for free from the website.

OpenCL accelerated multislice software written by Adam Dyson and Jonathan Peters from University of Warwick. clTEM is under development as of October 2019.

The code cudaEM is a multi-GPU enabled code based on CUDA for multislice simulations developed by the group of Stephen Pennycook.