A Practical Guide to Scanning and Transmission Electron Microscopy Simulations

Colin Ophus; Stephanie M Ribet; Georgios Varnavides; Jacob Madsen; Toma Susi

Computational Article

Contents

Here we briefly outline the mathematical concepts required to run STEM simulations.

Scalars, Vectors, Tensors, and Arrays¶

You are probably very familiar with the real numbers, which are values which measure a continuous quantity. These numbers can be integers such as 64, fractions such as $5/7$ , decimals such as -0.827 (a subset of fractions), and irrational numbers such as $\sqrt{2}$ . These values are scalars, meaning there is only one number associated with each value, making them one-dimensional (1D).

By contrast, a vector is a value which cannot be represented by a single number. Examples include three-dimensional (3D) position coordinates $(x,y,z)$ , magnetic moments, or a set of three Euler angles that describe an arbitrary rotation in 3D space. Vectors are always one-dimensional. In the same way a vector generalizes the concept of a scalar, a tensor generalizes both vectors and scalars by extending a value to any number of dimensions. You are also likely implicitly familiar with image arrays, which are 2D tensors for a greyscale image, and 3D for a color image.

Complex Numbers¶

To model scattering in electron microscopy, we work with a special number system referred to as complex numbers. A complex number is a length-2 vector, where the two values are called the real and imaginary values. The real part of each number is the kind of real number as described as above. The imaginary part is also a real number, but multiplied by the imaginary unit $i$ . This unit is a number which satisfies the equation

i^2 = -1.

(1.1)

Examples of complex numbers include $3+5i$ , $4/5-(9/7)i$ , $-0.23+4.33i$ , and $1+\sqrt{5}i$ . Figure 1.1 shows examples of complex numbers. These are enormously powerful, and essential when performing calclations using quantum mechanics.

Complex numbers. Every complex value consists of a real and an imaginary component, which can be used to compute the magnitude and phase of each value. — Figure 1.1:**Complex numbers.** Every complex value consists of a real and an imaginary component, which can be used to compute the magnitude and phase of each value.

Canvas(footer_visible=False, header_visible=False, layout=Layout(height='360px', width='400px'), resizable=Fal… — Figure 1.1:**Complex numbers.** Every complex value consists of a real and an imaginary component, which can be used to compute the magnitude and phase of each value.

In addition to the real and imaginary parts, we can define two other useful quantities. The magnitude of a complex number $z=a+ib$ is defined by its Euclidean length

|z| = \sqrt{a^2 + b^2}.

(1.2)

This quantity is also sometimes referred to as the magnitude or the absolute value. The second quantity we define is the complex argument, defined as the angle ϕ from the real axis to the complex vector. Note that we usually want the signed angle, meaning we use which quadrant the complex value $(a,b)$ is located in to determine the signed angle using the expression

\phi = \arctan2(b, a).

(1.3)

In programming, functions which measure this value are often named $\arg(z)$ or $\rm{angle}(z)$ . One of the most powerful identities using complex numbers is Euler’s formula

\exp(i z) = \cos(z) + i \sin(z).

(1.4)

This expression links complex numbers to the trigonometric functions. This relationship foreshadows a key upcoming concept, namely that the trigonometric oscillations describe waves, where the physical concept of the complex angle will be mapped onto relative phase shifts of electron waves. We will therefore usually refer to the complex argument as the phase of a given complex value.

As we will see, even though computers internally store complex numbers using real and imaginary components as $z=(a,b)$ , it is often more intuitive for us to represent these numbers using their magnitude and phase, i.e. as $|z| \exp{(i \phi)}$ .

We also sometimes take the conjugate of complex numbers, which means that we reverse the sign of the imaginary part. This operation is defined for a complex number $z=a+i b$ as

z^* = a - i b

(1.5)

or for $z=|z| \exp{(i \phi)}$ as

z^* = |z|\exp{(-i \phi)}.

(1.6)

1D Fourier Transform¶

The Fourier transform is an integral transformation that decomposes a function into its constituent spatial frequencies. This means that we write a function as a series of orthogonal sine waves with different frequencies / periodicities. A set of functions are orthogonal if we can multiply any pair of them and integrate over their domain to get a result of zero. The key property conveyed by this orthogonality is that the set of sine waves output by the Fourier Transform form a complete basis, and thus can represent any arbitrary function.

We can map a function $f(x)$ over the $x$ coordinates to its Fourier transform $F(\omega)$ over the ω coordinates using the expression

F(\omega) = \int_{-\infty}^{\infty} f(x) \exp\left( -i \omega x \right) dx,

(1.7)

and we can recover the original function with the inverse transform

f(x) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) \exp\left( -i \omega x \right) d\omega.

(1.8)

1D Discrete Fourier Transform¶

The Fourier transform is very useful for many mathematical analyses. However, we cannot use it in this form for this article, because unlike continuous variables such as position $x$ and frequency ω, our experimental images and simulated electron waves are sampled on discrete grids. We instead use the discrete Fourier transform, which transforms a finite-length function that has been sampled on a uniformly-spaced grid $f(r)$ into its Fourier transform $F(k)$ using the expression

\begin{aligned} F(k) &= \mathscr{F}_{r \rightarrow k}\{ f(r) \} \\ &= \sum_{n=0}^{N-1} f_n \exp\left( -2 i \pi k n / N \right), \end{aligned}

(1.9)

where $N$ is the number of samples for both functions, and $f_n$ is the $n$ ’th sample of $f(x)$ . The inverse DTF is defined as

\begin{aligned} f(r) &= \mathscr{F}_{k \rightarrow r}\{ F(k) \} \\ &= \frac{1}{N} \sum_{k=0}^{N-1} F_k \exp\left( 2 i \pi k n / N \right), \end{aligned}

(1.10)

where $F_k$ is the $k$ ’th sample of $F(k)$ . Both $f(x)$ and $F(k)$ can be complex-valued. If either sequence is purely real-valued, the other will be conjugate symmetric, which means that if reverse the sequence and take the complex conjugate, we will recover the same sequence.

Figure 1.2 shows an interactive version of the 1D discrete Fourier transform. Try switching between the different preset functions. You can also change the amplitude of $F(k)$ by clicking and moving the points up and down, or the phase by moving left and right. It can be difficult at first to gain intution for how changing the values of a given spatial frequency $F_k$ will change the output, but working in Fourier space will become second nature to you when you routinely perform STEM simulations.

Figure 1.2:**Interactive discrete 1D Fourier transform.**

Canvas(footer_visible=False, header_visible=False, layout=Layout(height='370px', width='675px'), resizable=Fal… — Figure 1.2:**Interactive discrete 1D Fourier transform.**

2D Fourier Transforms¶

When simulating transmission electron microscopy images, diffraction patterns, and other wavefunctions, we are typically working with two-dimensional (2D) arrays. Thus we need to expand our previous definition to work in two dimensions. However, because this transformation is orthogonal over different Cartesian coordinates, we can compute a 2D discrete Fourier transform by simply sequentially multiplying the two transformations, giving

\begin{aligned} F(k_x, k_y) &= \mathscr{F}_{x \rightarrow k_x}\{ \mathscr{F}_{y \rightarrow k_y}\{ f(x,y) \} \} \\ &= \sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f_{nm} \exp\left( -2 i \pi k_x m / M \right) \exp\left( -2 i \pi k_y n / N \right), \end{aligned}

(1.11)

where $f_{nm}$ is the pixel at indices $(m,n)$ for the function $f(x,y)$ , which has $(M,N)$ total number of pixels in the $x$ and $y$ directions respectively. We can simplify this notation by combining the two sequential operations into one operator, giving the forward and inverse transforms

\begin{aligned} F(\bm{k}) &= \mathscr{F}_{r \rightarrow k}\{ f(\bm{r}) \} \\ f(\bm{r}) &= \mathscr{F}_{k \rightarrow r}\{ F(\bm{k}) \}, \end{aligned}

(1.12)

where $\bm{r} = (x,y)$ and $\bm{k} = (k_x,k_y)$ are the 2D coordinate systems for real and diffraction space respectively. As in the 1D case, in practice we use the fast Fourier transform (FFT) algorithm to compute these transforms with small computational overhead. Figure 1.3 shows an interactive version of the forward and and inverse 2D FFT.

Interactive discrete 2D Fourier transform. Left panel shows diffraction space, right shows real space upsampled by 4x.
Both images have been FFT-shifted, with origins shown as a red + symbol.
Amplitude is shown as the pixel values, while phase is shown as a periodic colorwheel. — Figure 1.3:**Interactive discrete 2D Fourier transform.** Left panel shows diffraction space, right shows real space upsampled by 4x. Both images have been FFT-shifted, with origins shown as a red `+` symbol. Amplitude is shown as the pixel values, while phase is shown as a periodic colorwheel.

On the left panel of Figure 1.3, you can click to swap individual pixels between 0 and 1 in Fourier space, while observing the response on the real space function. On the right panel, you can click and drag to shift the function by a distance $\Delta \bm{r}$ , which is applied by using the Fourier shift theorem as

f(\bm{r} - \Delta \bm{r}) = \mathscr{F}_{k \rightarrow r}\{ \mathscr{F}_{r \rightarrow k}\{ f(\bm{r}) \} \exp( -2i \pi \bm{k} \, \cdot \Delta \bm{r} ) \},

(1.13)

where $\exp( -2i \pi \bm{k} \, \cdot \Delta \bm{r} )$ is the expression corresponding to a plane wave. If you try dragging the real space function in Figure 1.3 a short distance, you will see the characteristic coloring of a plane wave, a periodically repeating rainbow function.

A deep understanding of the 2D Fourier transform is essential to both simulate and understand TEM experiments. Therefore in Table 1.1 we provide a table of Fourier transform pairs which are important for TEM. Note that for simplicity we have omitted some details such as some of the numerical prefactors in front of each function.

Fourier transform pairs¶

Table 1.1:2D Fourier transform pairs important in TEM simulation and analysis. Note that these transform pairs assume that both $\bm{r}$ and $\bm{k}$ are 2D, i.e. $\bm{r}=(x,y)$ and $\bm{k}=(k_x,k_y)$ .

	Diffraction space	Real space
Name	Function	Function	Name
Delta Dirac	$\delta(\bm{k})$	1	plane wave along optic axis
Shifted Dirac function	$\delta(\bm{k} - \bm{k}_0)$	$\exp(-2i \pi \, \bm{r} \cdot \bm{k}_0)$	plane wave tilted to $\bm{k}_0$
diffraction spot pair	$\delta(\bm{k} - \bm{k}_0) + \delta(\bm{k} + \bm{k}_0)$	$\cos(-2 \pi \, \bm{r} \cdot \bm{k}_0)$	cosine wave
plane wave mult.	$F(\bm{k}) \exp(-2i \pi \bm{k} \cdot \Delta \bm{r})$	$f(\bm{r} - \Delta \bm{r})$	Shifted function
circular aperture	$u(\|\bm{k}\| - k_{\rm{max}})$	$\frac{J_1(2 \pi k_{\rm{max}} \|\bm{r}\| )}{\pi k_{\rm{max}} \|r\|}$	Airy disk function
defocus operator	$\exp(-i\pi \lambda \|\bm{k}\|^2 C_1)$	$\exp \left( i \frac{\pi \|\bm{r}\|^2}{\lambda C_1} \right)$	Fresnel quadratic phase factor
Gaussian envelope	$\exp \left( -\frac{\|\bm{k}\|^2}{2 \sigma^2} \right)$	$\exp \left( -2 \pi^2 \sigma^2 \|\bm{r}\|^2 \right)$	Gaussian envelope
exponential decay	$\exp(-\lambda \|\bm{k}\|)$	$\frac{\lambda^3}{(\lambda^2 + \|\bm{r}\|^2)^{3/2}}$	3D Lorentzian function
Lorentzian function	$\frac{\gamma^2}{\gamma^2 + \|\bm{k}\|^2}$	$K_0(2 \pi \gamma \|\bm{r}\|)$	Modified Bessel func., $2^{\rm{nd}}$ kind

Additional transform pairs can be found in the Wikipedia table of important Fourier transforms.

Guide to STEM Simulations

Scientific Python Code