Deriving the Multislice Algorithm

This is a short description on the theory of the multislice simulation method for electron scattering. For a more complete derivation including detailed theory, see Advanced Computing Computing in Electron Microscopy by EJ Kirkland[Kir10].

In electron microscopy, the energy of the incident electron waves (20–1000 keV) are generally much larger than the electrostatic potential wells of the atoms, which we will assume to only provide small perturbations on the forward motion of the electrons. Hence, we write the wave function, \(\psi_\mathrm{full}\), of the propagating electrons as a slowly varying plane wave along the optical axis, \(z\), with a modulation in \(x\) and \(y\) given by \(\psi(x,y)\)

\[ \psi_\mathrm{full}(x,y,z) = \psi(x,y)\exp(2\pi iz/\lambda) , \]

where \((x,y,z)\) are the three real-space coordinates and \(\lambda\) is the de Broglie wavelength of the electrons. Substituting this into the Schrödinger equation we obtain

\[ -\frac{\hbar^2}{2m} \left[\nabla_{xy}^2 + \frac{\partial^2}{\partial z^2} + \frac{4\pi i}{\lambda}\frac{\partial}{\partial z} + \frac{2meV(\vec{r})}{\hbar^2} \right] \psi(\vec{r}) = 0 \qquad \nabla_{xy}^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} . \]

In the high energy approximation, we assume that the wavefunction varies slowly in the \(z\)-direction compared to the potential and that the wavelength is small, thus

\[ \left| \frac{\partial^2 \psi}{\partial z^2} \right| \ll \left| \frac{1}{\lambda} \frac{\partial \psi}{\partial z} \right| . \]

Hence, the Schrödinger equation simplifies to a first order differential equation in \(z\)

\[ \frac{\partial \psi(\vec{r})}{\partial z} = \left[\frac{i\lambda}{4\pi} \nabla_{xy}^2 + i \sigma V(\vec{r}) \right] \psi(\vec{r}) , \]

where \(\sigma=2\pi me\lambda/h^2\) is the interaction parameter. This equation is integrated numerically by slicing the potential into thin slices, such that the influence of each slice can be approximated as a simple phase shift of the wave function. The wave function is propagated between slices as a small angle outgoing wave (Fresnel diffraction). The transmission and propagation across a single slice can be written

\[ \psi(x, y, z + \Delta z) = p(x,y,\Delta z) * [t(r) \psi(\vec{r})] + \mathcal{O}(\Delta z^2) , \]

where \(*\) represents a convolution. The transmission function, \(t(r)\), for the portion of the potential between \(z\) and \(z+\Delta z\) is

\[ t(x,y) = \exp\left[i\sigma \int_z^{z+\Delta z} V(x,y,z) dz'\right] , \]

and the conventional paraxial Fresnel propagator \(p(x, y, \Delta z)\) (first-order term of the Taylor expansion, see below) is

\[ p(x,y,\Delta z) = \frac{1}{i \lambda \Delta z}\exp\left[\frac{i\pi}{\lambda \Delta z}(x^2+y^2)\right] . \]

The wave at the exit plane of the specimen is obtained by sequentially propagating and transmitting the wave function starting with an assumed input wave.

Convolutions can be performed efficiently by utilizing the Fast Fourier Transform (FFT). The implemented form of Eq. the single slice propagation is

\[ \psi_{z + \Delta z}(x,y) = p(x,y) * \left[t(\vec{r}) \psi(\vec{r}) \right] = \mathcal{F}^{-1}\{P(k_x, k_y, \Delta z) \ \mathcal{F}[t(\vec{r}) \psi(\vec{r})] \} , \]

where \(P\) is the Fourier transform of the fresnel propagator. The computational cost for the FFT scales as \(N \log(N)\) with the number of samples \(N\).

Higher-order Fourier- and real-space formulations

1. Relativistic electron wavelength

A fast electron accelerated through a potential \(V\) acquires kinetic energy \(E = eV\). Its de Broglie wavelength must account for relativistic mass increase. Starting from the relativistic energy–momentum relation \(E_\mathrm{tot}^2 = (pc)^2 + (m_0 c^2)^2\) with total energy \(E_\mathrm{tot} = m_0 c^2 + eV\), the momentum is

\[p = \frac{1}{c}\sqrt{E_\mathrm{tot}^2 - (m_0 c^2)^2} = \sqrt{2m_0 eV\!\left(1 + \frac{eV}{2m_0 c^2}\right)}.\]

The de Broglie relation \(\lambda = h/p\) then gives the relativistically corrected wavelength:

\[\lambda = \frac{h}{\sqrt{2m_0 eV\!\left(1 + \dfrac{eV}{2m_0 c^2}\right)}}.\]

At 200 keV the correction factor \(1 + eV/(2m_0 c^2) \approx 1.20\), shortening \(\lambda\) by ~10% relative to the non-relativistic value. At 10 keV it is \(\approx 1.01\) and essentially negligible.

2. Helmholtz equation

In free space (no crystal potential), the electron wavefunction satisfies the stationary Schrödinger equation, which reduces to the Helmholtz equation:

\[\left(\nabla^2 + k_0^2\right)\psi(\mathbf{r}) = 0,\]

where the total wavenumber is \(k_0 = 1/\lambda\) in the crystallographic convention (\(k = 1/\lambda\) rather than the physics convention \(k = 2\pi/\lambda\), so plane waves are \(e^{i2\pi\mathbf{k}\cdot\mathbf{r}}\)). Writing the Laplacian in Cartesian coordinates with \(z\) along the optical axis:

\[\frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} + \frac{\partial^2 \psi}{\partial z^2} + \frac{1}{\lambda^2}\,\psi = 0.\]

3. Fourier decomposition in the transverse plane

The multislice method treats the specimen as a stack of thin slices perpendicular to \(z\). Within each slice the potential is projected onto a 2D plane (the projection approximation), and between slices the electron propagates through free space. Since the simulation cell is periodic in \(x\) and \(y\), we expand \(\psi\) in a discrete 2D Fourier series over transverse spatial frequencies \(\mathbf{k}_\perp = (k_x, k_y)\):

\[\psi(\mathbf{r}) = \sum_{\mathbf{k}_\perp} \hat{\psi}(\mathbf{k}_\perp, z)\; e^{i2\pi(k_x x + k_y y)}.\]

Each Fourier component represents a plane wave tilted at angle \(\theta \approx \lambda k_\perp\) to the optical axis (small-angle regime). Substituting into the Helmholtz equation, the transverse derivatives act on the exponential giving \(-4\pi^2(k_x^2 + k_y^2)\), while the \(z\) derivative acts on \(\hat{\psi}\):

\[\frac{\partial^2 \hat{\psi}}{\partial z^2} + 4\pi^2\!\left(\frac{1}{\lambda^2} - k_x^2 - k_y^2\right)\hat{\psi} = 0.\]

This is a 1D harmonic equation in \(z\) for each transverse frequency. The forward-propagating solution is \(\hat{\psi}(\mathbf{k}_\perp, z) = \hat{\psi}(\mathbf{k}_\perp, 0)\,e^{i2\pi k_z z}\), where the longitudinal wavenumber \(k_z\) satisfies the dispersion relation:

\[k_x^2 + k_y^2 + k_z^2 = \frac{1}{\lambda^2}.\]

Solving for \(k_z\) (positive root for forward propagation):

\[k_z = \sqrt{\frac{1}{\lambda^2} - k_\perp^2} = \frac{1}{\lambda}\sqrt{1 - \lambda^2 k_\perp^2}.\]

The dispersion relation encodes the geometry: tilted plane waves (large \(k_\perp\)) have smaller \(k_z\) and advance more slowly along \(z\), accumulating a phase lag relative to the axial wave (\(k_\perp = 0\), \(k_z = 1/\lambda\)). This differential phase accumulation is exactly what the propagator must encode.

4. Exact free-space propagator (Fourier space)

Propagation through a vacuum slab of thickness \(\Delta z\) multiplies each Fourier component by \(e^{i2\pi k_z \Delta z}\). In the multislice formalism we work with the reduced wavefunction \(\tilde{\psi}\) that has the fast axial oscillation \(e^{i2\pi z/\lambda}\) factored out. The reduced propagator is:

\[P(\mathbf{k}_\perp) = e^{i2\pi(k_z - 1/\lambda)\,\Delta z} = \exp\!\left[\frac{i2\pi\Delta z}{\lambda}\left(\sqrt{1 - \lambda^2 k_\perp^2} - 1\right)\right].\]

Defining \(x \equiv \lambda^2 k_\perp^2 = \lambda^2(k_x^2 + k_y^2)\), the phase argument is:

\[\boxed{\varphi_\mathrm{exact} = \frac{2\pi\Delta z}{\lambda}\left(\sqrt{1 - x} - 1\right)}\]

This is the exact free-space propagator phase. It is valid for all scattering angles and naturally handles both propagating (\(x < 1\)) and evanescent (\(x > 1\)) components.

4a. Evanescent modes

When \(k_\perp > k_0 = 1/\lambda\) (i.e. \(x > 1\)), \(k_z\) becomes purely imaginary and the wave decays exponentially:

\[\left.\sqrt{1 - x}\;\right|_{x>1} = i\sqrt{x - 1},\]

\[P = \underbrace{e^{-2\pi\Delta z\sqrt{x-1}/\lambda}}_{\text{exponential decay}}\;\cdot\;\underbrace{e^{-i2\pi\Delta z/\lambda}}_{\text{constant phase}}.\]

These modes carry no energy and decay rapidly. They are absent in the order 1 and 2 approximations, which implicitly assume \(x \ll 1\). In practice the antialiasing aperture zeros most evanescent components, but the exact propagator handles them correctly regardless.

4b. Numerical stability

For small \(x\) (paraxial regime, \(k_\perp \ll 1/\lambda\)), evaluating \(\sqrt{1-x} - 1\) directly subtracts two nearly equal numbers, losing significant digits. Multiplying numerator and denominator by the conjugate:

\[\sqrt{1-x} - 1 = \frac{(\sqrt{1-x})^2 - 1^2}{\sqrt{1-x} + 1} = \frac{-x}{\sqrt{1-x} + 1},\]

\[\boxed{\varphi_\mathrm{exact} = \frac{2\pi\Delta z}{\lambda}\cdot\frac{-x}{\sqrt{1-x} + 1}}\]

This is the form used in the implementation. For \(x \ll 1\) the denominator is \(\approx 2\) and the expression smoothly gives \(\varphi \approx -\pi\Delta z\lambda k_\perp^2\) without cancellation. For \(x\) near 1 it remains exact.

5. Taylor approximations of the Fourier-space propagator

The binomial series gives:

\[\sqrt{1-x} = 1 - \frac{x}{2} - \frac{x^2}{8} - \frac{x^3}{16} - \cdots\]

Order 1 (Fresnel / paraxial)

Keeping only the first correction term:

\[\varphi_1 = -\frac{2\pi\Delta z}{\lambda}\cdot\frac{x}{2} = -\pi\Delta z\,\lambda\,k_\perp^2\]

\[\boxed{P_1 = \exp\!\left(-i\pi\Delta z\,\lambda\,k_\perp^2\right)}\]

Since \(k_\perp^2 = k_x^2 + k_y^2\) appears linearly in the exponent, this factorizes as \(P_1 = e^{-i\pi\Delta z\lambda k_x^2}\cdot e^{-i\pi\Delta z\lambda k_y^2}\) — the standard separable Fresnel propagator.

Order 2

Keeping two correction terms:

\[\varphi_2 = -\pi\Delta z\,\lambda\,k_\perp^2 - \frac{\pi\Delta z\,\lambda^3}{4}\,k_\perp^4\]

\[\boxed{P_2 = P_1 \cdot \exp\!\left(-\frac{i\pi\Delta z\,\lambda^3}{4}\,k_\perp^4\right)}\]

Note that \(k_\perp^4 = (k_x^2 + k_y^2)^2 = k_x^4 + 2k_x^2 k_y^2 + k_y^4\). The cross-term \(2k_x^2 k_y^2\) is essential — the separable form \(k_x^4 + k_y^4\) is incorrect.

Truncation error scaling

The per-slice phase error of the order-\(n\) approximation is dominated by the first dropped term. For order 1:

\[\Delta\varphi_1 \approx \frac{\pi\Delta z\,\lambda^3}{4}\,k_\perp^4 \;\propto\; \lambda^3.\]

For order 2:

\[\Delta\varphi_2 \approx \frac{\pi\Delta z\,\lambda^5}{8}\,k_\perp^6 \;\propto\; \lambda^5.\]

The steeper \(\lambda\) scaling explains why order 2 converges to exact much faster with increasing beam energy — the \(\lambda^5\) vs \(\lambda^3\) dependence gives roughly two additional orders of magnitude suppression at 200 keV compared to 10 keV.

6. The exact multislice equation

Sections 2–5 treated free-space propagation only. Inside the specimen, the wavefunction satisfies the Schrödinger equation with a crystal potential \(V(\mathbf{r})\). In the high-energy (forward-scattering) approximation, the envelope \(\tilde{\psi}\) obeys:

\[i\frac{\partial\tilde{\psi}}{\partial z} = \hat{H}\,\tilde{\psi}, \qquad \hat{H} = \underbrace{\frac{\nabla_\perp^2}{4\pi k_0}}_{\hat{K}} + \underbrace{\sigma V(x,y,z)}_{\hat{V}},\]

where \(\sigma = 2\pi m e\lambda / h^2\) is the interaction parameter (relativistically corrected) and \(k_0 = 1/\lambda\). Under the projection approximation (potential constant within each slice of thickness \(\Delta z\)), the exact solution for one slice is:

\[\boxed{\psi(z + \Delta z) = \exp\!\left[i\Delta z\left(\hat{K} + \hat{V}\right)\right]\psi(z)}\]

This is the exact multislice equation. All multislice methods approximate it; they differ in how they handle the exponential of the combined operator.

7. Two independent error sources

Error source 1: split-step (Trotter) factorization

The operators \(\hat{K}\) (kinetic, acts in Fourier space) and \(\hat{V}\) (potential, acts in real space) do not commute: \([\hat{K}, \hat{V}] \neq 0\). The Fourier-space multislice method factorizes the combined exponential into separate steps:

\[e^{i\Delta z(\hat{K} + \hat{V})} \approx e^{i\Delta z\hat{K}} \cdot e^{i\Delta z\hat{V}}.\]

By the Baker–Campbell–Hausdorff formula, this Lie–Trotter splitting introduces an error:

\[e^{i\Delta z\hat{K}} \cdot e^{i\Delta z\hat{V}} = \exp\!\left[i\Delta z(\hat{K} + \hat{V}) - \frac{\Delta z^2}{2}[\hat{K}, \hat{V}] + \cdots\right].\]

The leading error term \(O(\Delta z^2)\) is present at every slice and can only be reduced by making \(\Delta z\) smaller.

Error source 2: Taylor truncation of the propagator

Within \(e^{i\Delta z\hat{K}}\), the exact kinetic phase \(\sqrt{1 - \lambda^2 k_\perp^2} - 1\) is approximated by a low-order Taylor polynomial, as described in Section 5.

8. Fourier-space method (PR #298)

The Fourier-space multislice uses the split-step approach: transmit in real space (multiply by the transmission function \(T = e^{i\sigma V_z}\)), then propagate in Fourier space (multiply by \(P(\mathbf{k}_\perp)\)):

\[\tilde{\psi}(z + \Delta z) = \mathcal{F}^{-1}\!\left[P(\mathbf{k}_\perp) \cdot \mathcal{F}\!\left[T \cdot \tilde{\psi}(z)\right]\right].\]

The three propagator options are:

Propagator	Phase	Error source 2
Order 1	\(-\pi\Delta z\lambda k_\perp^2\)	\(O(\lambda^3 k_\perp^4)\)
Order 2	\(-\pi\Delta z\lambda k_\perp^2 (1 - \frac{\lambda^2}{4}k_\perp^2)\)	\(O(\lambda^5 k_\perp^6)\)
Exact	\(\frac{2\pi\Delta z}{\lambda}\!\left(\sqrt{1-\lambda^2 k_\perp^2}-1\right)\)	none

The exact propagator (PR #298) eliminates error source 2 but retains error source 1 — the split-step factorization. It is fast because the propagator array is built once per unique slice geometry and reused as a simple Fourier-space multiply.

9. Real-space method (PR #236)

The real-space multislice takes a fundamentally different approach. Instead of factorizing \(e^{i\Delta z(\hat{K}+\hat{V})}\) into separate transmission and propagation steps, it Taylor-expands the exponential of the combined operator directly:

\[e^{i\Delta z\hat{H}}\psi = \sum_{n=0}^{N} \frac{(i\Delta z)^n}{n!}\,\hat{H}^n\,\psi, \qquad \hat{H} = \hat{K} + \hat{V}.\]

Since \(\hat{H}\) is never factorized, the powers \(\hat{H}^n\) naturally contain all cross-terms between \(\hat{K}\) and \(\hat{V}\):

\[\hat{H}^2 = \hat{K}^2 + \hat{K}\hat{V} + \hat{V}\hat{K} + \hat{V}^2.\]

The commutator \([\hat{K}, \hat{V}]\) is automatically accounted for to all orders by the Taylor series. There is no split-step error.

The kinetic operator \(\hat{K}\) involves \(\sqrt{1 - \lambda^2 k_\perp^2}\), which in real space becomes an infinite series of Laplacian powers. The order parameter controls how many terms are kept, and the expansion_scope parameter controls where the truncation is applied.

Propagator scope (expansion_scope="propagator")

Only the kinetic operator \(\hat{K}\) is expanded to order \(n\) in powers of \(\nabla_\perp^2\):

\[\hat{K}_n = \sum_{j=1}^{n} c_j\!\left(\frac{\nabla_\perp^2}{4\pi k_0}\right)^{\!j}, \qquad \hat{H} = \hat{K}_n + \hat{V}.\]

The coefficients \(c_j\) are the same Taylor coefficients as in the Fourier-space orders. At order 1, \(\hat{K}_1 = \nabla_\perp^2/(4\pi k_0)\), which is the Fresnel (paraxial) kinetic operator. At higher orders, additional powers of the Laplacian are included. The potential \(\hat{V}\) enters the operator once, but the outer exponential series mixes \(\hat{K}_n\) and \(\hat{V}\) through the repeated application of \(\hat{H}\).

The Laplacian is computed in real space using centered finite-difference stencils of configurable accuracy (default 6th order), avoiding the need for FFTs.

Full scope (expansion_scope="full")

The full conventional operator \(\hat{K}_1 + \hat{V}\) is applied repeatedly at each order of the inner expansion:

\[\hat{H}_\mathrm{full}^{(n)} = \sum_{j=1}^{n} c_j\,(\hat{K}_1 + \hat{V})^j.\]

The higher-order terms already contain kinetic–potential cross-terms within the operator itself, not just in the outer exponential. At order 1 this reduces to the same operator as the propagator scope. At higher orders, the potential enters the repeated Laplacian applications, so even the representation of \(\hat{H}\) is more complete.

This mode also enables computation of backscattered wave components by evaluating the correction term that arises from the difference between forward- and backward-propagating solutions at the slice boundary.

10. Relationship between the methods

The two approaches can be understood as addressing different layers of approximation in the multislice algorithm:

Method	Propagator truncation	Split-step error	Backscatter
Fourier, order 1	\(O(\lambda^3)\)	yes	no
Fourier, order 2	\(O(\lambda^5)\)	yes	no
Fourier, exact	none	yes	no
Real-space, propagator scope, order \(n\)	\(O(\lambda^{2n+1})\)	none	no
Real-space, full scope, order \(n\)	mixed to order \(n\)	none	yes

The exact Fourier-space propagator is the Fourier-space equivalent of the real-space method with expansion_scope="propagator" at order \(n \to \infty\). They would agree exactly in the limit of infinitely thin slices, where the Trotter splitting error vanishes. For any finite slice thickness, the real-space method is inherently more accurate because it exponentiates the combined operator rather than factorizing it.

The fully-corrected real-space method (expansion_scope="full") goes further still: it mixes kinetic and potential contributions in the operator expansion itself, and can additionally capture backscattered wave components that are absent from all Fourier-space formulations.

The tradeoff is computational cost. The Fourier-space method builds the propagator array once per unique slice geometry (a function of thickness, grid, sampling, and energy) and reuses it as a single Fourier-space multiply at every slice for every beam position. The real-space method must iterate a finite-difference Taylor series to convergence at every slice for every beam position, with each iteration requiring a stencil convolution over the full grid.