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

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.

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}$pasos comparados con el problema de diagonalización de la solución de onda de Bloch que escala como${\displaystyle N^{2}}$dónde${\displaystyle N}$es el número de átomos en el sistema. (Véase la Tabla 1 para comparar el tiempo de cálculo).

El paso más importante para realizar un cálculo multicapa es configurar la celda unitaria y determinar el grosor adecuado de la capa. En general, la celda unitaria utilizada para simular imágenes será diferente de la que define la estructura cristalina de un material en particular. Esto se debe principalmente a los efectos de aliasing que se producen por errores de envoltura en los cálculos de la transformada rápida de Fourier (FFT). La necesidad de añadir "relleno" adicional a la celda unitaria ha dado lugar a la denominación de "supercelda", y la necesidad de añadir estos píxeles adicionales a la celda unitaria básica conlleva un coste computacional.

Para ilustrar el efecto de elegir un espesor de corte demasiado delgado, consideremos un ejemplo sencillo. El propagador de Fresnel describe la propagación de ondas electrónicas en la dirección z (la dirección del haz incidente) en un sólido:

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

Dónde${\displaystyle \mathbf {u} }$es la coordenada de la red recíproca , z es la profundidad en la muestra y${\displaystyle \lambda }$es la longitud de onda de la onda electrónica (relacionada con el vector de onda por la relación${\displaystyle k=2\pi /\lambda }$). En el caso de la aproximación de ángulo pequeño (${\displaystyle \theta \sim }$100 mRad) podemos aproximar el desplazamiento de fase como${\displaystyle \Delta z}$. Para 100 mRad el error${\displaystyle d-S}$es del orden del 0,5% desde${\displaystyle \cos(0.1)=0.995}$Para ángulos pequeños, esta aproximación se mantiene independientemente de la cantidad de cortes que haya, aunque elegir un${\displaystyle \Delta z}$Un valor mayor que el parámetro de red (o la mitad del parámetro de red en el caso de las perovskitas ) para una simulación de múltiples capas daría como resultado la ausencia de átomos que deberían estar en el potencial del cristal.

Otras cuestiones prácticas son cómo incluir eficazmente efectos como la dispersión inelástica y difusa, las excitaciones cuantizadas (por ejemplo, plasmones, fonones, excitones), etc. Existía un código que tomaba en consideración estos aspectos mediante un enfoque de función de coherencia ^{[ 12 ]} llamado Yet Another Multislice (YAMS), pero dicho código ya no está disponible ni para su descarga ni para su compra.

Existen varios paquetes de software para realizar simulaciones multicorte de imágenes. Entre ellos se encuentran NCEMSS, NUMIS, MacTempas y Kirkland. Existen otros programas, pero lamentablemente muchos no han recibido mantenimiento (por ejemplo, SHRLI81 de Mike O'Keefe del Laboratorio Nacional Lawrence Berkeley y Cerius2 de Accerlys). En la Tabla 2 se presenta una breve cronología de los códigos multicorte, aunque esta lista no es exhaustiva.