Prismatic 2.0-Simulation Software for Scanning and High Resolution Transmission Electron Microscopy (STEM and HRTEM)

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


Introduction
Transmission electron microscopy (TEM) is heavily used in both materials science and biological studies of materials on the nanoscale, due to its high spatial resolution and the flexibility of operating modes (Egerton et al. 2005). TEM experiments can be performed using plane wave illumination, where users can either record the far field intensity as a diffraction pattern (Williams and Carter 2009) or by forming an image of the electron wave after it has been transmitted through the sample, often referred to as high-resolution transmission electron microscopy (HRTEM) (Buseck et al. 1989). Alternatively, the electron beam can be focused into a small spot and scanned over the sample surface, which is referred to as scanning transmission electron microscopy (STEM) (Pennycook and Nellist 2011). The introduction of spherical aberration correctors in the past few decades enable the formation of a finer probe in STEM (Batson et al. 2002) and point-spread function in HRTEM ). Aberration-corrected TEM and STEM have greatly facilitated many atomic resolution experiments, including imaging single-layer graphene sheets (Gass et al. 2008, Robertson andWarner 2013), elemental mapping (Kothleitner et al. 2014), atomic electron tomography (Yang et al. 2017), vibrational spectroscopy (Venkatraman et al. 2019), observation of polar skyrmions (Das et al. 2019), and many others.
With the widespread adoption of charge-coupled device cameras (Krivanek and Mooney 1993), and later direct electron detectors (MacLaren et al. 2020), both STEM and TEM are now fully digital sciences. Augmenting STEM and TEM experiments with modern data science methods holds enormous promise for future experimentation (Spurgeon et al. 2020). One of the most data-intensive families of STEM experiments is the use of fast direct electron detectors to collect thousands or even millions of 2D images of the diffracted probes over a 2D grid of probe positions, often referred to as a four dimensional-STEM (4D-STEM) dataset (Ophus 2019). One of the key ingredients for developing and implementing new computational methods is the ability to perform forward simulations of the various STEM and TEM experimental methods (LeBeau et al. 2008, Zhang et al. 2020. HRTEM simulations and STEM simulations for the 2D, 3D, or 4D experiments described above are typically performed using the multislice method Moodie 1957, Ishizuka andUyeda 1977). However, these simulations have traditionally required very long compute times, due to the fact that a separate multislice simulation must be performed for each new probe position (Kirkland 2020).
In a previous study, Ophus (2017) introduced a new algorithm, named plane-wave reciprocal-space interpolated scattering matrix (PRISM) for simulation of STEM experiments. The PRISM algorithm involves calculating and storing a compact scattering matrix operator which can be rapidly applied to each of the probe wavefunctions to model their propagation through the sample. This algorithm can potentially increase the speed of STEM image simulations by multiple orders of magnitude. The PRISM algorithm was first implemented into the original Prismatic simulation code by . Since then, PRISM has been adapted for double-channeling STEM-EELS simulation by Brown et al. (2019); it has been separately implemented for image simulation by Brown et al. (2020) and Madsen and Susi (2021).
Prismatic is intended for use as a fully-featured TEM/STEM simulation software for electron microscopy, for diverse use cases such as experimental validation, database generation, or teaching.
It also serves as the reference implementation for the aforementioned PRISM algorithm. Prismatic is an open-source and cross-platform software package that can be easily installed, easily used, and that comes with a GUI. The software is primarily written in C++ with CUDA modules for GPU acceleration to take advantage of available computing and HPC resources and is readily integrated into other scientific open-source software for microscopy applications such as py4DSTEM .
In this paper, we present Prismatic version 2.0, a software package for image and diffraction space simulations of electron scattering for both STEM and HRTEM. Prismatic v2.0 has many new basic features, such as performing HRTEM simulations, increased support for arbitrary aberrations, support for arbitrary STEM scan patterns, focal series simulations, and enhanced support for generalized input and output. We have also added several improvements to enhance the accuracy and speed of simulations, including a new approach for the correct sub-slicing of 3D atomic potentials with sub-pixel shifting, a refocusing approach to the scattering matrix calculation for PRISM simulations that increases accuracy for thicker samples and delocalized probes, and various post-processing methods for coherence or shot-noise limitations. We have added many upgrades to the code, including new pipelines for compiling, a unit-testing suite, an overhauled GUI system (previewed in Fig. 1), and precompiled Docker containers for ready-to-use installations of the command-line and pyprismatic interfaces. In this paper, we explain these methods and additions in detail, as well as demonstrate several new applications for Prismatic and some uses of the Python bindings in pyprismatic.

The Multislice Method
We describe the electron beam using a complex wavefunction ψ(x, y, z), where x, y, z are the real space coordinate system. When considering only forward scattering, we can reduce the problem to the evolution of a 2D wavefunction ψ(r) along the optical axis z, where r = x, y . This evolution is described by the paraxial Schrödinger equation for fast electrons (Van Dyck 1985) where i is the imaginary constant and V (r) is the electrostatic potential, usually corresponding to the sample. The relativistically-corrected electron-matter interaction constant σ is given by where γ, m e , q e , λ, and h are the relativistic correction factor, the electron mass, the electron charge, the relativistically corrected electron wavelength, and the Planck constant respectively.
The two operators on the right-hand side of equation (1) do not commute so a widely utilized numerical approach to its solution is a split-step method first derived by Cowley and Moodie (1957). For small changes in z, Eq.
(1) can be solved in two steps, taking first only the ∇ 2 term and then the V term into account. First, we divide up the sample into a series of N slices, V n (r), which are 2D arrays that integrate the electrostatic potential contained in a given slice of thickness ∆z, given by The solution to Eq. (1) taking into consideration just V on the right-hand side is where ψ 0 (r) is the input wavefunction and T (r) is referred to as the transmission function. This is equivalent to the so-called "phase object" approximation which holds for samples thin enough to ignore the effects of thickness. For the next part of the split-step solution, where we assume V n (r) = 0, the operator ∇ 2 can be efficiently applied in Fourier space (Ishizuka andUyeda 1977, Kirkland 2020), giving The term e iλ∆z|k| 2 , referred to as the propagation operator, uses the 2D reciprocal space coordinates k = (k x , k y ), and the forward and inverse 2D Fourier transform operations denoted byF r→k andF −1 k→r respectively. Equations (4) and (5) are alternately applied to calculate the final wavefunction after interacting with the sample, which is typically referred to as the exit wave. This numerical solution is called the "multislice method" (Cowley and Moodie 1957). It requires N transmission operations and N − 1 propagation operations, and is the most common simulation algorithm for modeling TEM experiments (Kirkland 2020).
Alternatively the operator-product in Eq. (6) can be encapsulated as single matrix equation (Sturkey 1962) where we have opted to start with the probe in reciprocal space on the right hand side with the exit-surface wave function on the left hand side being in real space. The concept of the scattering matrix is common in quantum mechanics for calculating scattering behavior of electrons and other charged particles (Weinberg 1995) and is typically purely in reciprocal space. For our purposes it will be more convenient to use the reciprocal space to real space formulation in Eq. (7) since the STEM probe is compact in reciprocal space and the PRISM algorithm will involve cropping the exit wave in real space. The rows or columns of S can be changed from real to reciprocal space or vice versa with the operation of the appropriateF r→k orF −1 k→r applied either to the left or right side of S.

Calculation of Projected Potentials
The calculation of the atomic scattering potential V n (r) is one of the most crucial aspects of scattering simulations. The discretization of the scattering potential limits the accuracy and total amount of information that can be transmitted and propagated; any artifacts created in the calculation of V n (r) can lead to nonphysical behavior in the elastic scattering calculations. Many implementations of the scattering potential utilize the isolated atom approach, where the total potential is the sum of potentials of independent, isolated atoms such that J o u r n a l P r e -p r o o f where i is the index of each atom. The isolated atom approach, since each atom of a given species is independent and equivalent, can take advantage of a precomputed look-up table. Generally, this approach is computationally cheap and simple to implement, and correctly accounts for the dominant aspect of nuclear scattering but inherently does not capture bonding effects. Analytical scattering potentials do not exist for elements other than hydrogen; the rest of the atomic species are typically parameterized based on single atom scattering factors determined through Hartree-Fock calculations, like in the parameterization by Kirkland (2020) where a 0 is the Bohr radius and a i , b i , c i , and d i are fitted parameters. This parameterization can also be integrated analytically along the beam direction so that the contribution of an atom to its two-dimensional slice can be calculated directly where K 0 is the modified Bessel function.
However, 2D scattering potential calculations can not perfectly capture the 3D scattering of real atoms, since atoms must first be assigned to discrete projected potential slices before integration-an explicitly 3D integration of the potential into the different potential slices can more accurately capture this 3D scattering (Lobato and Van Dyck 2015). Simple isolated atom implementations are also often pixel-limited in the planar directions, such that the atomic center can only exist on discrete pixels. This can fail to capture the subtlety of the real positions of atoms, especially when considering thermal vibration effects. Sub-pixel accuracy, however, can be costly to implement and hard to achieve. One possible strategy is to forgo the use of a look-up table and instead integrate each atom directly on sub-pixel grids, which can can become prohibitive for especially large cells. Alternatively, one could apply some sort of sub-pixel shift in the transverse direction, which can introduce artifacts into the scattering potential if implemented poorly. Alternative 3D integration methods for computing atomic potentials have been given by Lobato and Van Dyck (2015) and Madsen and Susi (2021). Atomic potentials can also be directly evaluated electron densities as calulated using Density Functional Theory calculations, as by Madsen and Susi (2021 In a conventional multislice simulation, each STEM probe position requires independent evaluation of the propagation through the potential as described by Eqs. In PRISM, instead of directly forming the incident probe wavefunction and propagating the wavefunction through the projected potential via the multislice algorithm, we instead calculate and store the S-matrix in Eq. (7) for a basis set of incident plane waves. The S-matrix can be reused for any number of STEM probes once it has been calculated, and thus, we can trade the upfront computation of S for the much greater acceleration in the output calculation stage. Further speed ups for a modest loss of simulation accuracy can be achieved through coarsening of the input planewave basis (by only calculating every f th plane wave from the original simulation grid; f is known as the the PRISM interpolation factor) and reciprocal-space interpolation of the output wave function (achieved by real-space cropping of the output wave to a square region measuring 1/f of the simulation grid). The reciprocalspace interpolation grants the user a large amount of control of the speed and accuracy of the simulation, which could otherwise only be done in the multislice algorithm by tuning the resolution of the projected potential or resampling of the diffraction patterns. For a more in-depth discussion of the PRISM algorithm, we refer readers to Ophus et al. (2016a)

Numerical Calculation of 3D Potentials
We address the issues sometimes encountered with projection of the atomic potentials into single 2D slices by implementing 3D potential integration with subpixel accuracy of the atomic positions by use of a pre-

Lor
Au at (5, 5, 5.125) Å  calculated look-up table, which is sampled more finely than the multislice slice sampling in the propagation direction, alongside sub-pixel accuracy Fourier shifting of the potentials in the plane perpendicular to propagation. Our implementation is designed to improve upon simple 2D integration in a manner that is robust in various simulation conditions, is easy to understand and tune, and does not sacrifice a significant amount of computation. We first calculate the atomic potential for each unique species in the system on a local pixel grid R = (x r , y r , z r ), using the three-dimensional parameterization described in Eq. (9), where x r , y r are coordinates in the plane perpendicular to propagation and at the same resolution as the final potential field, and z r is the coordinate along the propagation axis at some integer subsampling N z of the final slice thickness t. Equation (9) is evaluated out to maximum radius r max , which along with the real space potential samplings r x , r y , and t/N z determine the dimensions of the grid R. The minimal coordinate of the atomic potential is R = (r x , r y , 0) which avoids the discontinuity otherwise present at R = 0. To prevent sharp steps in the potential at the extremities of r max , we subtract the value of V (r) at the smallest extent of the atom in any of x r , y r , and J o u r n a l P r e -p r o o f z r and then clamp any negative values to zero such that where . This procedure is sketched in Fig. 2(a). The potential for each unique species is then stored in a look-up table after a 2D Fourier transform in x and y such that since the subpixel shifts will be applied in Fourier space.
To achieve sub-pixel positioning of the ith atom, to position (x i , y i ), in the plane perpendicular to the beam propagation direction we apply a Fourier shift to the transformed potential V (k, z) based on the difference between the rounded pixel position of the atom round(x i /x y ) and its ideal fractional where B(k) is a soft bandwidth limit given by applied to the potential in reciprocal space, and k max is the maximum reciprocal space coordinate of the grid on which the lookup tables is constructed. Application of the bandwidth limit B(k) limits oscillations in real space. This applied shift moves the atom to the correct place in real space without creating significant artifacts in the potential in Fourier space (where the propagation calculation occurs) below the Nyquist limit. We note that a coarse real space sampling (i.e., large values for r x and r y ) can inject oscillations into the tails of the potential. This can result in unphysical negative values of the scattering potential in real space. V shifted (k, z) is then Fourier transformed back to real space upon which it is again bandlimited such that where r band = 1.0, ( xr max(xr) ) 2 + ( yr max(yr) ) 2 ≤ 1.0 0.0, otherwise . (16) This final potential is then added into the final cell array. Along the propagation direction, the superresolution z values of V(r) are simply binned to their closest rounding slice. We found that interpolation of the potential along the propagation direction caused little to no increase in accuracy at the expense of artifact introduction, and thus, was not implemented.

STEM Probe Formation
The most common probe configuration for STEM experiments is given by a circular aperture in the condenser plane. The use of fast Fourier transforms (FFTs) to implement the multislice solution to Eq. (1) and to set up STEM probes enforces periodic boundary conditions which results in interference artifacts as the tails of the STEM probe approaches the edges of the simulation grid. We find that implementing these STEM probes in Prismatic by using a soft aperture minimizes these artifacts. Fig. 3 shows a comparison between STEM probes defined using hard and soft aperture edges. A hard aperture probe is defined by the function where k probe defines the maximum scattering vector included in the STEM probe, and is proportional to the maximum scattering angle α probe by the expression α probe = λ k probe . The soft aperture probe used in Prismatic is defined by where ∆k are the 2D pixel sizes in Fourier space, is the Hadamard (element-wise) product, and || · || 2 is the 2-norm. The term in the denominator of Eq. (18) is equal to the Fourier space pixel size ∆k multiplied by the Fourier space coordinate k, which depends on the local orientation of k if the two pixel sizes are not equal. If these two pixel sizes are equal (typical of square simulation cells), Eq. (18) reduces to Ψ 0 (k) = min max k probe − |k| ∆k + 1 2 , 0 , 1 .
Comparison of hard-aperture and softaperture STEM probes. Initial STEM probe defined with (a) a hard aperture edge, and (b) a soft aperture edge. Initial STEM probes for PRISM with fx = fy = 3, for (c) a hard aperture edge, and (d) a soft aperture edge. The soft edge minimizes artifacts resulting from the periodic boundary conditions inherent in all STEM simuations, PRISM simulations are colored inside the cropping box in real space.
full size PRISM simulation. These cases are fairly similar, though the rings corresponding to the probe tails become inaccurate after the third ring due to interference with probe tails on the opposite side of the grid resulting from the inherent periodic boundary conditions.
In most cases this refinement will make only a minor improvement to simulation accuracy, but for smaller simulation cells in multislice, or smaller cropping boxes in PRISM simulations, using the correct aperture becomes more important. Fig. 3c and d show hard and soft apertures respectively for a PRISM f = 3 simulation. In the real space image of the STEM probe, the regions outside of the cropping box are shown in a grey color scale. The soft aperture produces two accurate probe tail rings, while the hard aperture is inaccurate for all probe tails. The hard aperture in Fig. 3c even produces an asymmetric center lobe of the STEM probe. These examples are at a deliberately lower resolution and smaller grid real-space size than for a typical STEM simulation to more clearly demonstrate the importance of using an accurate soft aperture function for defining the initial STEM probes. To our knowledge, this is a novel feature in Prismatic.

Anti-Aliasing of the Transmission and Propagation Operators
As both the multislice and PRISM algorithms rely on propagation through the calculation of successive discrete FFTs, we must prepare arrays and calculations in a way that prevents aliasing artifacts. Aliasing artifacts arise when signals of different frequencies become indistinguishable from each other due to the discrete sampling of the signals. To prevent aliasing of the propagation wave function Ψ, we multiply the array with a binary anti-aliasing mask which removes pixels above the Nyquist limit (1/2 the extent of the array) from the calculation at each propagation and transmission step. This method causes an intensity loss of the parts of the wave-function that are scattered to high angles due to the anti-aliasing aperture, and may cause some inaccuracies in electrons that might multiply-scatter back into low angles. Applying an anti-aliasing filter at the Nyquist limit completely prevents highly scattered electrons from "wrapping around" the wavefunction and scattering through the periodic boundary. An alternative approach for anti-aliasing is to apply anti-aliasing filters to both the propagated wave function and the projected potential slices at 2/3rds the extent of the arrays (Kirkland 2020, Lobato and Van Dyck 2015). We note for clarity that in our implementation, the projected potential is not band-limited. One should set the real space sampling of the simulation such that the error converges in the scattered regions of interest such that the application of the anti-aliasing aperture does not affect the interpretation of the simulation results.

HRTEM Simulations
TEM simulation can be performed using the exact same methods as STEM multislice by replacing the incident wave function with a plane wave, i.e, where k tilt is the tilt of the beam. Since the Smatrix calculation in the PRISM method forms the transmitted STEM probe using a scattered plane-wave basis set, the implementation can be readily used for HRTEM simulations.
In Prismatic, HRTEM simulations are now implemented through utilization of the S-matrix infrastructure, with added functionality for more granular control of the beam selection for propagating sets of plane waves simultaneously with the multislice algorithm.
Specific tilt ranges and selections can be controlled both through radial and rectangular mask generation, and through use of the PRISM interpolation factor (which coarsens the grid of tilts available), for example, to test different beam tilts or to model angular coherence. HRTEM simulations in Prismatic are focused to the middle of the sample cell instead of the entrance surface as implicitly done for STEM simulations. In simulation, this is accomplished through application of a defocus propagator to each beam in the S-matrix after its standard calculation. The final output is formally the 3D complex=valued S-matrix, with unique plane wave tilts along the beam direction, which can be saved either as the complex wavefunction or as integrated intensities.

Aberrations
Coherent aberrations of the probe-forming lens in STEM and the image-forming lens in HRTEM simulations can be modeled through the inclusion of an aberration function χ when calculating the wave function in Fourier space such that where Ψ 0 is the unaberrated wavefunction. χ can be expressed through a variety of appropriate basis sets (Krivanek 1994, Thust et al. 1996; in Prismatic, which now includes functionality to apply arbitrary sets of aberrations, we employ a unitless basis set similar to the one described in Ophus et al. (2016b). For aberrations with radial order m and azimuthal order n, the aberration function χ( k) can be fully described for coherent aberrations as where λ is the electron wavelength, C mag m,n and C ang m,n are dimensionless coefficients describing the aberration magnitude in rads and azimuthal phase of an aberration, and atan2 is the 2-argument arctangent function which returns the polar angle of (k x , k y ) in the correct quadrant. The dimensionless aberration coefficients are related to the convention in Ref. (Kirkland 2011) which have units of length by C mag m,n = 2πĈ mag m,n /(mλ). With this basis set, it is relatively simple to describe aberrations and provide input data in the form of delimited text files. In Figure 4, we showcase aberration phase plates at 300kV up to 2 inverse angstroms, for spherical aberration, 3-fold astigmatism, axial coma, and an arbitrary mixture of the above with defocus (C 2,0 ) and 4-fold astigmatism (C 4,4 ).

FIG. 4.
Example aberration phase plates. From topleft to bottom-right, spherical aberration, 3-fold astigmatism, axial coma, and an arbitrary mixture of aberrations at 300kV. Borders represent a cutoff of 2 inverse angstroms.

Refocusing of the Scattering Matrix
An inherent limitation of the PRISM algorithm lies within the interaction of the propagated probe and the cropping box applied to the scattering matrix. As the reciprocal-space sampling of the plane-wave basis set coarsens (by increasing the interpolation factor f ), the cropping box applied to the propagated plane-waves becomes smaller. At high interpolation factors, much of the information carried by the propagated wavefunction can become destroyed or otherwise compromised as the beam spreads beyond the edges of the cropping box.  (5) sample running through entire cell. After S-matrix is propagated through the full cell, it can be refocused to any plane before intensity is calculated. STEM probe intensity column shows the radius containing different fractions of the probe as a function of defocus. Phase of all complex waves is shown using a cyclic color scale, similarly to Fig. 4. We note that these simulations are chosen with an arbitrary sample size and beam energy, and this schematic is intended to generally depict the behavior of a propagated scattering matrix undergoing refocusing.
The strong interaction of the beam at high interpolation factors thus severely limits the accuracy of the PRISM method. This becomes a more prominent effect when simulating thick samples, due to the natural broadening of the converged beam as it propagates through the sample, or when simulating STEM probes with large defocus values which are already large in extent. To overcome this limitation of the PRISM method while still retaining the computational speed-up achieved with large interpolation factors, we propose scattering matrix refocusing, an algorithmic technique which can be used to boost the accuracy of PRISM simulations with significant beam spread.
The scattering matrix S, when calculated by the PRISM algorithm, represents a basis of propagated plane waves which are inherently in focus at the entrance surface of the sample. To reduce the severity of the converged probe interaction with the cropping box once the probe is calculated from S, we apply a free-space propagator to each beam within S to propagate the exit wave to a plane where the intensity distribution is generally known a-priori to be more compact, and is applied in Fourier space by application of the convolution theorem. S refocused can then be used to calculate the exit wavefunction without any further alterations to the image formation algorithm. Since the propagator in Eq. (22) imparts only a phase shift to each beam, the intensity of the exit wavefunction in the diffraction pattern, where STEM intensity measurements are recorded, is unaffected.
One way to interpret the mechanism of matrix of refocusing is as shifting of the optical plane at which the cropping box is applied. In the original PRISM algorithm, the cropping box is applied at the exit surface of the sample. With refocusing, the cropping box can instead be thought to be applied at the point along the optic axis where the propagated probe is most converged-therefore, the probe has minimal opportunity to interact with the cropping box and errors that may have otherwise been incurred by application of the cropping box itself are reduced.
The process is illustrated in Fig. 5 where the refocusing procedure is applied to beams focused on the entrance, Fig. 5(a), and exit surfaces, Fig. 5(b), of different samples. The electrostatic potential V (r) of different fictitious samples is shown in the leftmost panel of the figure and includes samples with no atoms, atoms at the top surface, an atom mid-sample, an atom at the bottom surface and a single column of atoms. The simulated complex wavefunction of a STEM probe is plotted for each of the samples in the next column. For the vacuum sample the converged probe is seen to spread from its cross-over point in Fig. 5(a) and condense to its cross-over point in Fig. 5(b). For the cases where single atoms are introduced at the top, middle and bottom of the sample we see electron scattering to high angles emanating from the positions of these atoms. For the case of a column of atoms we see beam channeling behavior characteristic of STEM (Hovden et al. 2012) in Fig. 5(a) where the beam couples to the column at its cross-over point atop the sample. This channeling behavior is diminished for the equivalent panel in Fig. 5(b) since the beam is not brought to cross over until the exit surface of the sample. The next 3 columns in the figure show how the constituent plane waves that form the basis of the Smatrix interact with the specimen. In simulation we can "refocus" the exit surface wave function to any arbitrary plane within the specimen by freespace propagation of the exit wave as described previously and a depth section of the exit wave propagated to different planes is plotted in the second to last column of Fig. 5 for each of the different cases. The final column plots the radially integrated intensity of each wavefunction. We aim to find the plane where the wave-function is most compact and apply the PRISM cropping box after propagating the exit wave to this plane. For all but one of the cases in Fig. 5(a) this plane is close to the entrance surface of the specimen -i.e. the crossover point of the probe -for the case of the column of atoms channeling of the beam down the column means that the exit surface of the specimen is the plane where the exit wave function is most compact. In Fig. 5 the beam is always more compact toward the exit surface of the specimen. These results suggest that the original crossover point of the incident STEM probe is generally the plane where the beam is most compact after propagation even after interaction with the sample except for the case where strong channeling of the beam by a column of atoms occurs. This principle can be used to guide selection of the optimal plane for S-matrix refocusing and subsequent cropping.

Post-processing
A crucial aspect of bridging the interpretation gap between "perfect" simulated images to those acquired in experiment is to process the simulated image in ways that model experimental distortions and noise conditions. These corrections are becoming increasingly important given the advent of high-throughput workflows using machine learning, such as automated analysis (Zhang et al. 2020) and denoising (Vincent et al. 2021) with deep neural networks, which require large amounts training data with latent space distributions similar (or, ideally, equivalent) to that of the target application data. In this release of Prismatic, we have implemented a small series of post-processing routines into the pyprismatic package to supplement the standard simulation routine so that coherence effects, dosage effects, and aberration effects as previously described can be easily applied to simulated data.

Coherence Effects
An implicit assumption in image simulation is that electrons emanate from an ideal, infinitesimal point source, whereas in reality, imaging electrons are emitted from a finite area of the surface of the electron emitter with varying energies and thus some degree of quantum partial spatial and temporal incoherence. Spatial coherence effects for modern electron microscopes, equipped with field emission gun sources such that the distribution of incident angles is relatively small, can be accounted for through integration of independently propagated wave functions by assuming spatial coherence Allen et al. (2004). The detected wave function can then be described as where k is the wave-vector nominally in the direction of propagation, p(k β ) is a probability density function describing the distribution of incident directions, and χ is the aberration function given by Eq. (20) which models the aberrations of the post-specimen HRTEM objective lens. Including the effects of spatial coherence in HRTEM simulation, then, can be easily achieved by simulating a set of tilted plane waves up to the spread of incident directions and performing incoherent averaging of their intensities once aberrations effects have been included.
For STEM simulatons, we can apply spatial coherence effects by considering finite source-sizes through a convolution model. While the precise distribution describing source-size broadening effects can be quite complex (as shown when measured through holography (Verbeeck et al. 2012)), source-size effects on the detected image can be well approximated by the convolution of the final image intensity with a blurring kernel such that where I is the image intensity and K is a blurring kernel.
Prismatic has implemented source-size effects for postprocessing of STEM images through convolution with standard kernels; kernel functions for Gaussian and Cauchy kernels are included, which can adequately represent the effects of source-size blurring for most uses (Verbeeck et al. 2012). Convolution is performed with unit kernel normalization such that the intensity distribution maintains its physical definition as a probability density function for detection of scattered electron. Source-size blurring by means of convolution can be applied independently and in any order with any other incoherent averaging procedure, such as the averaging over frozen phonon configurations.
Chromatic aberration (C c ) is a temporal coherence effect that occurs as a result of the non-uniform energy of imaging electrons as emitted by the electron gun. Lower energy electrons are deflected more strongly by a magnetic lens than higher energy electrons, resulting in an energy-dependent focal point of the microscope lenses. C c becomes the dominant aberration when C s corrective lenses are introduced and becomes more important when low energy beams are used for imaging, as shown by the C c limited resolution where E 0 is beam energy, ∆E is energy spread, β is the angle of collection, and C c is the chromatic aberration coefficient. Ignoring any issues in instrument power supply, the standard deviation of defocus ∆z can be obtained through the following relationship where ∆E is the energy spread in the electron gun (Williams and Carter 2009). Assuming the spread of incident beam energies is small, temporal coherence effects can be analogously included as an integration over the distribution of focal planes, as done for incident angles in Eq. 23. Thus, chromatic aberration effects can be accounted for in simulated results through incoherent averaging of propagated wavefunctions at a series of focal planes.

Shot Noise
Modern STEM experiments are often noise limited by the maximum dose a sample can tolerate (Egerton 2013). As a series of discrete electron scattering events, dose-limited noise profiles are well modeled by Poisson distributions. Applying Poisson noise to simulated images (which exist as evaluations of quantum-mechanical probabilities) is thus one of the most important post-processing steps for comparing simulated images and experimental images.
Prismatic implements Poisson noise application for HRTEM images with electron dose measured in counts per area; for STEM images, dose is applied in counts per probe. For both, intensities are scaled by the relevant measure of electron dose and then pixel counts are sampled from a Poisson distribution using the scaled intensities as independent variances. To avoid degradation of the electron distribution statistics, application of Poisson noise is performed as the final step of post-processing for any simulated images.

Results and Discussion
The previously described multislice and PRISM algorithms, as well as HRTEM simulation, S-matrix refocusing, simulation series, and arbitrary probe positions, have been implemented in the newest release of Prismatic, alongside utility features such as data importing and an updated GUI. Brief discussion of these implementations will follow, to serve the ongoing need of documenting open-source scientific software. Relevant case studies to the microscopy community will be shown in conjunction with each major feature that is discussed.

Running a STEM Imaging Simulation
One of the primary uses of Prismatic is to perform imaging simulations for STEM experiments. Fig. 6 shows one such example, simulating the interaction of a STEM probe with a bilayer of WS 2 with a twist angle of 7.34 • . This simulation was performed using pyprismatic with these settings: In this simulation, we have stored the outputs of two common STEM imaging modalities using additional "save" parameters. The first output is the expectation value of the diffracted STEM probe's momentum, shown in Fig. 6b for the k x and k y directions. The changes in the probe's "center-of-mass" (COM) momentum k CoM (r) as measured in the diffraction plane, as it is scanned across atomic sites is clearly visible as an initial rise and then fall for each site. These signals can be directly measured in a STEM experiment, either by segmented detectors (Shibata et al. 2010) or from pixelated detectors (Müller-Caspary et al. 2019, Ophus 2019) as in the simulation here. These measurements are usually intended for differential phase contrast (DPC) measurements. There are various ways to reconstruct the phase shift of the sample from these COM-DPC signal channels, including the iterative method described in  and shown in Fig. 6c.
The second output is for monolithic annular detectors, which are commonly used for annular bright field or dark field imaging. For maximum flexibility, Prismatic integrates the diffracted probe intensity using virtual detectors shaped in concentric rings. These rings are finely sampled (the default bin width is 1 mrad), so that the user can generate many different annular detector configurations from a single simulation. For example, in this simulation a bright field image can be generated by summing the intensity outputs from 0 to 21 mrads. The contrast of the different bins are shown as vertical slices in Fig. 6d.  Fig. 7 shows 4D-STEM outputs from a simulation of the same WS 2 sample shown in Fig. 6. The projected potentials and probe positions are shown in Fig. 7a and b, while the diffracted probe intensities are shown in Fig. 7c. These probe positions span 3 atomic sites, and the DPC signal at these sites manifests as a shift of the average moment towards these sites. Additionally, when some portion of the the STEM probe overlaps with these sites, a significant number of electrons are scattered outside of the initial probe's angular range. The simulation settings of this example are nearly identical to the previous simulation, except for these changes: Note that in Fig. 7a, the potentials have been shifted by half of the total cell dimensions in x and y, in order to place the probe positions in the center of the figure. No discontinuity is visible at the cell boundaries, due to the periodicity of the input cell and the simulation.

Running an HRTEM Imaging Simulation
Performing HRTEM simulations in pyprismatic is similarly simple, with notation as follows for a plane wave simulation with complex valued output (before integration into intensity): m.algorithm = "t" m.saveComplexOutputWave = True We include example Jupyter notebooks for other imaging scenarios in the Prismatic repository.

Comparative Thermal Convergence of Potential Parameterization Methods
The previous version of Prismatic used 2D lookup tables for the projected atomic potentials. In this release, we have implemented 3D lookup tables for the atomic potentials with subpixel shifting of the atomic sites as previously described above. The 3D potential integration scheme with subpixel shifting can improve the accuracy of a scattering simulation by offering a more realistic representation of a sample in terms of its projected potential. This is accomplished through more accurate accounting of atomic positions along the beam direction via 3D integration and perpendicular to the beam direction through subpixel shifting. Calculation of the projected potential is thus much less sensitive to the input simulation parameters such as slice thickness and projected potential resolution. In a 2D integration scheme, incorrect setting of these parameters can introduce artifacts; for example, atoms can jump unstably across slices and from pixel to pixel in the transverse direction under thermal perturbations.
Here, we examine further how the two projected potential integration techniques compare under frozen phonon convergence tests. Since these schemes generate fundamentally different projected potentials at a given level of resolution, we can expect them to display qualitatively different convergence behavior over a range of temperature conditions.
To investigate the comparative convergence behavior of the projected potential integration schemes with respect to number of frozen phonons, we simulated both STEM and HRTEM images of twisted bilayer WS 2 at a 7.34 • moiré angle. The convergence of scattering simulations with respect to the number of frozen phonon configurations for two-dimensional materials requires more frozen phonons than for a thicker sample because fewer unique atomic configurations can be sampled in the beam direction for any given frozen phonon configuration. Some results studying the precise thermal convergence of scattering simulations include those by Aarholt et al. (2020). The moiré cells were simulated at 80kV with a real space potential input resolution of 0.075Å, a 2Å slice thickness, and a sampling factor of 40 for 3D potential integration. 128 unique frozen phonon configurations were generated for each simulation, with increasing Debye-Waller factors from 0.01 to 0.16 in multiplicative steps of 2 representing five distinct temperatures for the sample; results for each frozen phonon were saved independently. STEM simulations utilized the PRISM algorithm at an interpolation factor of 2 and a probe step size of 25 pm over a 16x reduced window of the input cell.
Convergence results are shown in Figure 8. Convergence was measured at each of these five represented temperatures for HAADF STEM at with detector angles between 60 and 120 mrad; for ABF STEM with detector angles between 10 and 30 mrad; and for HRTEM with a single untilted plane wave. To measure convergence, we used a k-fold cross-validation based scheme. For each number N of frozen phonons averaged, we generate 64 unique subsets of N configurations from the set of 128 total configurations. We then form the incoherently averaged image for each subset, measure the pixel-wise standard deviation of the intensities between the different subsets, and finally measure the mean over the field-ofview of the standard deviations. At no point are results from the differing projected potential integration schemes directly compared.
Here, we note that the convergence rate with respect to frozen phonons is largely similar across both simulation modes and temperatures. Overall cross-validation errors for both schemes drop as temperature drops. This is expected, as there is less variance in the overall atomic positions of the sample. At higher temperatures, the 3D potential with subpixel shifting scheme has comparatively higher cross-validation error, which drops below the 2D potential scheme without subpixel shifting as the temperature decreases and progressively drops faster. Crossover temperatures for cross-validation error between the two schemes occur at different points for HAADF STEM as compared to ABF STEM and HRTEM. The crossover point ultimately results from the difference in convergence behavior between the two schemes and can be understood through a sampling argument-it is likely primarily a result of the subpixel shifting. Since the 2D potential scheme must place atoms on top of discrete pixels, it is essentially sampling a small finite set of configurations for each atom, while the 3D potential with subpixel shifting samples the true continuous distribution. At high temperatures, it is easy for the 2D potential scheme to saturate a statistically significant sampling of this finite set, while at low temperatures, it becomes difficult to perturb the atom far enough to generate different projected potential slices. In contrast, the 3D potential with subpixel shifting scheme can appropriately sample its distribution of frozen phonon configurations at all temperatures, leading to improved convergence behavior at low temperatures and slightly worse convergence at high temperatures. Again, we note that the representation of the sample given by the 3D potential with subpixel shifting integration scheme is fundamentally more physically realistic and consistent against simulation parameters than 2D potentials without subpixel shifting. While the new scheme is more computationally expensive, the importance of this difference in physicality should guide the usage of the technique in simulations, especially since, for large scale simulations, the computation of the projected potential is still an overall negligible cost as compared to the computation of wavefunction propagation.

Comparison of Soft and Hard Probe Apertures in STEM Diffraction
The differences in aliasing behavior between soft and hard STEM probe apertures can be hard to predict for an actual simulation. It is important to remember FIG. 9. STEM diffraction simulation of SrTiO3. Comparative diffraction simulations for STEM diffraction of SrTiO3 at 80kV and with a 4mrad probe semiangle for multislice and PRISM simulations with both soft and hard probe apertures. Simulation of a 12x12x16 unit cell block of SrTiO3 with 64 frozen phonons; input probe is placed above the center of a Sr site. Color scale is set to have equivalent contrast between multislice and PRISM simulations. that aliasing effects may not always arise as ringing artifacts in simulations. In Figure 9, we demonstrate the differences in simulation behavior between soft and hard aperture probes for STEM diffraction of a SrTiO 3 crystal at 80kV with a 4mrad probe, for both multislice and PRISM simulations. The PRISM simulations were run with an interpolation factor of 2. All simulations were run with 3D potentials, a 0.1Å pixel size, and 64 frozen phonons. For the multislice simulations, there is no discernible difference in the intensity or contrast of the diffracted disks between the soft apertures. For the PRISM simulations, however, there is a significant difference in the contrast (and apparent shape) of the disks, and it is seen that the soft aperture PRISM simulation much more accurately matches that of the multislice simulation. While not every use case or simulation may show such drastic differences, the soft aperture probe will overall improve the accuracy of any simulation and is more likely to be a crucial component in improving the quality of inexpensive simulations, with coarser samplings and smaller input simulation cells. The effect of a finite source-size upon the simulated STEM images of twisted bilayer WS 2 at 100kV are shown in Figure 10, panels (a) and (b) for the as-simulated HAADF image and the HAADF image convolved with a Gaussian kernel with a FWHM of 80pm, respectively. The effect of Poisson noise is demonstrated in Figure 10 panel (c) for a dose rate of 5000 counts per probe.

Limited Coherence Effects in STEM Simulations
In Prismatic, it is now easy to account for coherence effects caused by chromatic aberration by performing defocus series simulations. As a case study, we present convergence results for chromatic aberration effects against the number of defocus planes simulated. STEM simulations of twisted bilayer WS 2 at a 7.34 • moiré angle were performed at 20kV and 100kV over a defocus range of ±150Å and ±30Å, respectively, with 65 focal planes each, a potential resolution 0.1Å and a probe step size of 20pm. We note for clarity that each defocus series was performed with a single execution of the simulation

FIG. 11.
Convergence series of chromaticallyaberrated STEM simulations. Results for BF (red) and HAADF (blue) STEM for twisted bilayer WS2 at 20kV (solid) and 100kV (dashed) with defocus spreads of 50Å and 10 A, respectively. Error is measured as the root-mean square error of total probe intensity against the fully averaged STEM image.
code, and that Prismatic provides new output formats that allow for easy data access for such simulation series.
The PRISM method was used with an interpolation factor of 2; we note that PRISM is a doubly advantageous technique for defocus series simulations even at an interpolation factor of 1 as the scattering matrix can be reused for each defocus value, whereas in multislice, each probe for each defocus must be simulated independently. For convergence, we measure the pixelwise root-meansquare error of a series of images averaged against differing numbers of defocus planes. In computing the defocus spread ∆z, we used a fixed energy spread for both operating voltages; the defocus spreads we used were 50Å and 10Å for beam energies of 20kV and 100kV, respectively. Results for these two simulations are shown in Figure 11. As expected, the fraction of probe error at 100kV in both bright field and HAADF detectors is very low in all cases; at 20kV, where the energy spread is more substantial, the error is more significant when including less than four additional defocus planes (five planes in total).

Limited Coherence Effects in HRTEM Simulations
With the reuse of the S-matrix infrastructure, in Prismatic, we have easy access to tilt-series HRTEM simulations. Tilt series can be used to model coherence effects, for example, or experimental procedures like precession electron diffraction. Here, we will showcase simulation results taking advantage of this infrastructure to model spatial coherence effects that can arise in HRTEM experiments. We continue with WS 2 as a model material. This time, to ease the image interpretation, we simulate a monolayer WS 2 cell that is approximately 317Å by 183Å in size. Simulations were performed at 300kV with an input real-space resolution of 0.075Å in each direction; eight total simulations were run, each with four values of defocus-2.6 nm, 9.8 nm, 31.4 nm, and 56.9 nm relative to the middle of the monolayerto represent different contrast mechanisms and over two ranges of output tilts, ± 3mrad and ±15 mrad. For each set of simulated tilts, we perform incoherent averaging of the plane-waves with Gaussian weights up to ±3σ to emulate spatial coherence effects, where σ is the standard deviation of the beam tilt. We similarly include temporal coherence effects by incoherently averaging each tilt series over a range focal planes with Gaussian weights for a focal spread with 10Å standard deviation for each simulation. Results are shown in Figure 12 for coherence levels σ of 0.0, 0.3, 1.0, and the unrealistically poor 5.0 mrad. As shown, emulating spatial coherence effects in this manner captures well the loss of contrast at high defocus. For small σ a large input cell is needed to achieve proper sampling of tilts as the resolution in Fourier space is inversely proportional to the sample size. Therefore, modeling these effects is more computationally intensive, particularly with regards to memory, than a typical HRTEM simulation; on a workstation with a 16 core Intel Xeon Gold 6130 processor, 250 GB system RAM and a Nvidia Quadro P5000 GPU (16 GB RAM), these simulations took approximately 1 minute to perform for each frozen phonon, as compared to the sub-10 second simulation time for single plane-wave imaging.

Matrix Refocusing of STEM Simulations with Large Defocus
Focal series experiments are common in studies that involve 3D reconstruction, as they provide varying sets of information with respect to the sample phase. In order to validate these methods, such as through-focal tomography (Hovden et al. 2014 arbitrary off-zone axis. We compare the multislice algorithm, the PRISM algorithm with an interpolation factor of 2, and the PRISM algorithm with the same interpolation and with refocusing. The simulation was performed at 200 kV with a semi-convergence angle of 21 mrad and a real-space atomic potential resolution of 0.1Å. Results comparing the accuracy of PRISM with and without refocusing to the results of the multislice algorithm are shown in Figure 13 for two STEM probes propagated through the center and through the edge of the nanoparticle over a defocus range of ± 250 nm (measured relative to the center of the nanoparticle). Starting from the center of the nanoparticle, the PRISM algorithm without refocusing begins to accumulate significant amount of error, particularly within the center disk, as the applied defocus reaches large values (in this simulation, around ± 100 nm). In contrast, PRISM with refocusing shows much more stable behavior through this defocus range within the center disk and accumulates only a small amount of error outside the beam convergence angle. At the side of the nanoparticle, both simulations show relatively consistent behavior when measuring error in binned virtual detector. There is a notable error in both methods at the the edge of the center disk, which can be mostly attributed to bin assignment and precision artifacts when integrating the full pattern to the virtual detector in the simulation. When comparing against the actual 2D CBED pattern before integration onto the virtual detector, we can clearly see a great discrepancy between results using PRISM with and without refocusing. Visual comparison of the features within the center disk shows that, at large defoci, information and structure are almost completely lost without refocusing, whereas applying refocusing to the scattering matrix maintains much more of the information at high defocus values. This may be important, for example, in generating CBED images for use in applications such as machine learning, where visual similarity may be more important than integrated measures of intensity. Finally, in spite of the results presented here, we note that application of the refocusing technique is highly specific to certain imaging conditions and samples of interest and is not a broadly applicable technique to PRISM simulations overall, and that the use of refocusing may be guided by investigating its performance on the metrics of interest for one's own application.

Conclusions
In this manuscript, we have presented version 2.0 of our Prismatic code for simulation of transmission electron microscopy experiments. This version has added various features, including 3D potential integration, HRTEM plane wave simulations, refocusing of the S-matrix for STEM simulations with large probe defocus values, more accurately sampled STEM probes, and unit tests for more reliable updates of the code. In addition we have also added several post-processing methods, including simulation of coherence limits generated by source size and energy spread effects, and shot noise calculations using Poisson random distributions of electron counts. With these additions, Prismatic is more useful to researchers who are performing STEM and TEM simulations. Comparison of error for STEM simulations using the PRISM algorithm with and without matrix refocusing at an interpolation factor of 2, as measured against the equivalent multislice simulation. (a) Projected potential for a spherical Au nanoparticle with 5 nm diameter, tilted off-axis with center-probe and side-probe positions marked (blue and green, respectively). (b) CBED patterns for multislice (left), PRISM (middle), and PRISM with refocusing (right); panels are split between in-focus (left side) and 250nm overfocused (right side). (c) Heatmaps of absolute error of PRISM (left) and PRISM with refocusing (right) as measured against multislice over a range of ±250nm of defocus and over scattering angles from 0-60mrad. Top panels in (b) and (c) represent probes propagated through the center of the nanoparticle, while bottom panels represent probes propagated along the side of the nanoparticle. We note the horizontal error lines visible in (c) are largely a result of coordinate artifacts that arise when comparing PRISM and multislice results; error results presently shown in (c) include slight blurring by a Gaussian kernel for purposes of visual clarity.

Code Availability
The source code for Prismatic is available at our Github repository. For other downloads, walkthroughs, and information, please visit prism-em.com. Simulations ran using pyprismatic to produce figures in the manuscript can be found in the source code repositories. Other scripts and data used can be made available upon request. For purposes of manuscript review, the most current development version of the source repository can be found at the author's development repository. A full updated software release is planned to be timed with publication of this manuscript.