2.8. prismatique.stem
For specifying STEM systems, STEM simulation parameters, and running STEM simulations.
The PRISM algorithm can be used to perform STEM simulations as an alternative to the multislice algorithm. Here we discuss very briefly aspects of both the PRISM and multislice algorithms, although it is recommended that the user see Ref. [Ophus1] for more details on the former and Ref. [Kirkland1] for an exposition on the latter.
As discussed in the documentation for the class
prismatique.thermal.Params, the intensity pattern for a given probe in
a STEM experiment or a given beam in a HRTEM experiment depends on the state
operator \(\hat{\rho}_{t}\) of a transmitted beam electron, which is assumed
to be a weighted sum of pure states
where \(\left|\psi_{t}\left(\delta_{f};\mathbf{u}_{1},\ldots,\mathbf{u}_{N}; \boldsymbol{\delta}_{\beta}\right)\right\rangle\) being the state vector of a transmitted beam electron for a perfectly coherent beam and a sample in a frozen atomic configuration \(\left\{ \mathbf{u}_{j}\right\} _{j=1}^{N}\), with the \(\delta_{f}\) and \(\boldsymbol{\delta}_{\beta}\) implicitly specifying the defocus and beam tilt respectively. Note that the weighted sum is over the defocus, beam tilt, and \(\left\{ \mathbf{u}_{j}\right\} _{j=1}^{N}\). Each \(\left|\psi_{t}\left(\delta_{f};\mathbf{u}_{1},\ldots,\mathbf{u}_{N}; \boldsymbol{\delta}_{\beta}\right)\right\rangle\) is assumed to have been evolved from an initial pure state \(\left|\psi_{i}\left(\delta_{f};\boldsymbol{\delta}_{\beta}\right) \right\rangle\) that is determined by the conditions of the corresponding perfectly coherent beam at incidence with the sample. Each \(\left|\psi_{i}\left(\delta_{f};\boldsymbol{\delta}_{\beta}\right) \right\rangle\) is evolved using the fast electron approximation. In this approximation, the full beam electron wavefunction \(\Psi\left(\mathbf{r} \right)\), for a given perfectly coherent beam and frozen atomic configuration, is assumed to be of the following factored form [Kirkland1]:
where
\(k=\left|\mathbf{k}\right|=1/\lambda\) is the magnitude of the electron’s total wavevector with \(\lambda\) being the beam electron wavelength; \(e^{2\pi ikz}\) is a plane wave travelling along the \(z\)-axis [i.e. the optic axis]; \(\psi\left(\mathbf{r}\right)\) satisfies
with
\(\hbar\) being Planck’s constant, \(e\) being the elementary charge, \(m_{e}\) being the rest mass of an electron, and \(V\left(\mathbf{r}\right)\) being the total potential due to the sample, which is treated classically; and
with \(\left|x,y\right\rangle\) being the electron transverse position eigenvector. As discussed in the documentation for the class :prismatique.discretization.Params, the total potential \(V\left(\mathbf{r}\right)\) is replaced by an effective potential \(\tilde{V}\left(\mathbf{r}\right)\), which is obtained by coarse-graining \(V\left(\mathbf{r}\right)\) along the \(z\)-axis:
where \(N_{\text{slices}}\) is a positive integer;
with
\(\Delta Z\) is the \(z\)-dimension of the sample’s supercell in units
of length. See the documentation for the class
prismatique.discretization.Params for a discussion on sample supercells
and the calculation of \(\tilde{V}_{n}\left(x,y\right)\). Note that
\(z=z_{0}\) and \(z=z_{N_{\text{slices}}}\) are the
\(z\)-coordinates of the entrance and exit surfaces of the sample’s
supercell respectively. This coarse-graining approximation introduces a global
error that is linear in the slice thickness \(\delta z\). Upon replacing
\(V\left(\mathbf{r}\right)\) with \(\tilde{V}\left(\mathbf{r}\right)\),
Eq. (2.8.4) can be solved using the following iterative
scheme [Kirkland1]:
where
with the entrance surface of the sample being located at \(z=0\) and \(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) being the incident beam wavefunction;;
with the exit surface of the sample being located at \(z=N_{\text{slices}}\delta z=\Delta Z\) and \(\psi^{\left(\text{exit}\right)}\left(x,y\right)\) being the so-called exit wavefunction; the notation \(\hat{O}\left\{ g\right\} \left(q_{x},q_{y}\right)\) denotes an operator \(\hat{O}\) transforming a function \(g\) into a space with coordinates \(q_{x}\) and \(q_{y}\) [e.g. \(\hat{O}\) could be a two-dimensional Fourier transform]; the notation \(\left(\hat{O}_{1}\circ\hat{O}_{2}\right)\left\{ g\right\} \left(q_{x},q_{y}\right)\) denotes the composition of two operators \(\hat{O}_{1}\) and \(\hat{O}_{2}\) transforming a function \(g\) into a space with coordinates \(q_{1}\) and \(q_{2}\) with \(\hat{O}_{2}\) being applied first;
with
and
The operator \(\hat{C}_{n;\varphi}\) expresses the transformation that the wavefunction \(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) undergoes when a beam/probe electron passes through a 2D potential \(\tilde{V}_{n}\left(x,y\right)\) perdendicularly, i.e. the wavefunction acquires a phase shift \(\varphi\tilde{V}_{n}\left(x,y\right)\) while keeping the amplitude the same. The operator \(\hat{K}_{\delta z^{\prime}}\) expresses the transformation that the wavefunction \(\psi_{x_{p},y_{p}}\left(x,y,z\right)\) undergoes when a beam/probe electron propagates through free space in the \(z\)-direction by a distance \(\delta z^{\prime}\). In this discussion, we will refer to \(\hat{K}_{\delta z^{\prime}}\) as the free-space propagator. If \(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) is the incident wavefunction of a perfect coherent probe in a STEM experiment, then Eqs. (2.8.10)-(2.8.18) essentially yield the multislice algorithm for STEM simulations.
The key observation involved in the PRISM algorithm is that the incident wavefunction of a perfect coherent probe in a STEM experiment [Eq. (2.8.11)] can be expanded in terms of plane waves as [Ophus1]:
where the \(\phi_{k_{x},k_{y}}\left(x^{\prime},y^{\prime}\right)\) are the plane wave basis functions:
\(\alpha_{k_{x},k_{y}}\left(x_{p},y_{p};\theta_{x},\theta_{y}\right)\) are the expansion coefficients:
with \(\left(x_{p},y_{p}\right)\) denoting the coordinates in the
\(xy\)-plane of the probe position, \(\theta_{x}\) and
\(\theta_{y}\) being the probe tilt angles along the \(x\)- and
\(y\)-axes respectively, \(h\) being the height [i.e. the length along
the \(z\)-axis] of the sample’s supercell,
\(\chi\left(k_{x},k_{y}\right)\) being the phase deviation due to coherent
lens aberrations, and \(\zeta\left(k_{x},k_{y}\right)\) being an aperture
function. See the documentation for the class
embeam.coherent.Aberration for a definition of
\(\chi\left(k_{x},k_{y}\right)\), and Ref. [DaCosta1] for the expression of
\(\zeta\left(k_{x},k_{y}\right)\) implemented in prismatic. One
important property of the aperture function, as defined in Ref. [DaCosta1], is
that for some \(k_{\text{max}}>0\) we have
At this stage, it is important to point out that both real-space and momentum/Fourier-space need to be discretized in order to handle wavefunctions numerically. The smallest possible Fourier-space pixel sizes in the \(k_{x}\)- and \(k_{y}\)-directions, which we denote here as \(\Delta k_{x}\) and \(\Delta k_{y}\) respectively, are determined by the \(x\)- and \(y\)-dimensions of the sample’s supercell:
where \(\Delta X\) and \(\Delta Y\) are the \(x\)- and \(y\)-dimensions of the sample’s supercell in units of length. It is worth noting at this point that our system is subject to the following periodic boundary conditions:
In prismatique/prismatic, for pixels of size \(\left(\Delta k_{x},
\Delta k_{y}\right)\), the \(k_{x}\)- and \(k_{y}\)-dimensions of
Fourier-space in units of said pixels is equivalent to half the \(x\)- and
\(y\)-dimensions of the sample’s supercell in units of pixels. The reason
that the former has half the dimensions of the latter is because prismatic
removes all momenta beyond the Nyquist limit, which is half of the original
Fourier-space pixel set, in order to prevent aliasing.
In discretized Fourier-space with pixels of size \(\left(\Delta k_{x}, \Delta k_{y}\right)\), we rewrite \(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) as
where \(N_x\) and \(N_y\) are the \(x\)- and \(y\)-dimensions of the sample’s supercell in units of pixels.
The PRISM algorithm involves expanding \(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) in terms of plane waves in a manner similar to that in Eq. (2.8.26):
where \(f_{x}\) and \(f_{y}\) are positive integers called the \(x\)
and \(y\) interpolation factors introduced in the documentation for the
class prismatique.discretization.Params, and
are what we refer to as the “sample supercell’s reduced xy-dimensions in units
of pixels”, also introduced in the documentation for the class
prismatique.discretization.Params. Note that for \(f_{x}+f_{y}>2\),
Eq. (2.8.27) approximates
\(\psi^{\left(\text{incident}\right)}\left(x,y\right)\) further as compared
to Eq. (2.8.26). Hence, choosing interpolation factors that
satisfy \(f_{x}+f_{y}>2\) will result in a loss of simulation accuracy since
we are discarding some plane wave basis functions, however the upside is an
improvement in simulation wall time. Further below we discuss how wall time
depends on \(f_{x}\) and \(f_{y}\).
By discarding plane wave basis functions according to \(f_x\) and \(f_y\), and implementing the anti-alias procedure described above, we are effectively working in a different discretized Fourier-space that is a subset of the original discretized Fourier-space: the new Fourier-space pixel sizes in the \(k_{x}\)- and \(k_{y}\)-directions are respectively
and the \(k_{x}\)- and \(k_{y}\)-dimensions of this new discretized Fourier-space in units of pixels are respectively
The pixelated STEM detector that is simulated in prismatic shares the same
pixel sizes and dimensions as this new discretized Fourier-space. Noting that a
beam electron’s scattering angles with the \(x\)- and \(y\)-axes,
\(\theta_x\) and \(\theta_y\), are related to the electron’s momentum
by:
it follows from Eqs. (2.8.23), and (2.8.28)-(2.8.30) that the angular field of view of the aforementioned pixelated STEM detector in the \(x\)- and \(y\)-directions are
Equations (2.8.29)-(2.8.32) are also applicable to the case where the multislice algorithm is used, i.e. we are still operating in this new discretized Fourier-space, however the interpolation factors are trivially set to \(f_x=f_y=1\).
In the PRISM algorithm, we use Eq. (2.8.10) again but this time with
and
for each \(\left(m_{x},m_{y}\right)\) pair in the double sum of Eq. (2.8.27). Once we have calculated all the necessary elements of \(S_{m_{x},m_{y}}\left(x,y\right)\), we can then use it to calculate the exit wavefunctions:
with \(S_{m_{x},m_{y}}\left(x,y\right)\) being the so-called \(S\)-matrix or scattering matrix. In other words, the plane wave basis functions are transformed individually according to Eq. (2.8.10), and then added together using the same expansion coefficients \(\alpha_{m_{x}f_{x}\Delta k_{x},m_{y}f_{x}\Delta k_{y}}\left(x_{p},y_{p};\theta_{x},\theta_{y}\right)\) in accordance with Eq. (2.8.27) to get the exit wavefunction \(\psi^{\left(\text{exit}\right)}\left(x,y\right)\). In the PRISM algorithm, a slight modification to Eq. (2.8.35) where instead we calculate
with
i.e. \(\tilde{S}_{m_{x},m_{y};x_{p},y_{p}}\left(x,y\right)\) is obtained by
cropping \(S_{m_{x},m_{y}}\left(x,y\right)\) to a rectangular region that is
centered about \(\left(x_{p},y_{p}\right)\), and has \(x\)- and
\(y\)-dimensions in units of length smaller than those of the sample’s
supercell by factors of \(4f_x\) and \(4f_y\) respectively. Similarly
to discarding plane wave basis functions, the cropping of
\(S_{m_{x},m_{y}}\left(x,y\right)\) will result in a loss of simulation
accuracy, especially for significant beam spread, however this loss comes with
faster wall times. To minimize the loss of simulation accuracy, prismatic
offers the option to refocus the exit wavefunction upon calculating it as
described above. To refocus the exit wavefunction, the free-space propagator
\(\hat{K}_{\delta z^{\prime}}\) [Eq. (2.8.13)] is applied to the
\(S\)-matrix [DaCosta1]:
where the propagation displacement \(\delta z^{\prime}\) is determined
automatically by prismatic such that the spread of the refocused exit
wavefunction has been reduced, thus reducing the loss of information due to
cropping. Since the free-space propagator imparts only a phase shift to the
Fourier components of the exit wavefunction, the corresponding diffraction
pattern — where only intensity measurements are recorded — is unaffected
[DaCosta1].
The main advantage of the PRISM algorithm is that \(S_{m_{x},m_{y}}\left(x,y\right)\) can be recycled to calculate the exit wavefunction for each probe position \(\left(x_{p},y_{p}\right)\). If we assume that \(N_x=N_y\gg1\), then we have [Ophus1]:
where \(t_{\text{Multi}}\) and \(t_{\text{PRISM}}\) are the wall times for the multislice and PRISM algorithms respectively, \(N_{\text{probes}}\) is the number of probe positions considered, and \(N_{\text{PW}}\) is the number of plane waves considered in our expansion of Eq. (2.8.27). If \(N_{\text{probes}}\gg8f_{x}f_{y}N_{\text{slices}}\), which is often the case, then the performance of the PRISM algorithm will scale like:
As mentioned above, the state operator \(\hat{\rho}_{t}\) of a transmitted
beam electron is assumed to be a weighted sum of pure states
\(\hat{\rho}_{t}\left(\delta_{f};\mathbf{u}_{1},\ldots,\mathbf{u}_{N};
\boldsymbol{\delta}_{\beta}\right)\) over defocus, beam tilt, and frozen atomic
configurations. In prismatique, the weighted sum is actually only over
defocus and frozen atomic configurations, not beam tilt, for STEM
simulations.
Modules
For specifying the output parameters for STEM simulations. |
|
For running STEM simulations. |
|
For specifying simulation parameters related to the modelling of STEM systems. |