Iterative Phase Retrieval Algorithms for Scanning Transmission Electron Microscopy

We start by considering a 4D dataset composed of diffraction intensities $I(\boldsymbol{R},\boldsymbol{k})$, where $\boldsymbol{R}$ and $\boldsymbol{k}$ denote the real-space probe position and reciprocal-space spatial frequency, respectively. We form a vectorial virtual image using a first moment detector $\boldsymbol{k}$:

where $\alpha$ denotes each of the two Cartesian directions $x$ and $y$ and the virtual CoM image, $\mathcal{I}_{\alpha}^{V}(\boldsymbol{R})$, is proportional to the gradient of the electrostatic sample potential $\partial_{\alpha}V(\boldsymbol{R})$. The electrostatic sample potential can be reconstructed by Fourier-integrating eq. 1 iteratively (Frankot, Chellappa, 1988):

where $\mathcal{F}$ is the two-dimensional Fourier transform operator, $\nabla$ is the gradient operator which can be computed numerically using a centered finite-difference, $\beta$ is the step size, $k$ is the spatial frequency magnitude, $\Re\{\cdot\}$ denotes the real-part of the expression and ^′ indicates the next iteration. Note that here, and in what follows, we express the phase-shift imparted by projected potentials in units of radians. We define a self-consistent error metric using eq. 2, by comparing the current estimate of the finite-difference gradient against the experimental CoM:

To avoid periodic artifacts arising from the use of Fourier transforms, we pad and mask the reconstructed object at every iteration (Savitzky et al., 2021). Finally, we implement a simple backtracking algorithm where the step-size is automatically halved if the error increases relative to the previous iteration and the current update step is rejected.

Due to the helical path accelerated electrons take along the optic axis, as well as the orientation of the detector in the STEM column relative to the beam raster grid, there is typically a rotational offset between the real- and reciprocal-space coordinate systems. Since the gradient of the electrostatic potential is expected to be a conservative vector field, we can solve for this rotation by minimizing the curl or maximizing the divergence of the vectorial CoM virtual image as a function of rotational offset (Savitzky et al., 2021; Zachman et al., 2022). These techniques have a $180^{\circ}$ ambiguity in the relative rotation which can be solved by requiring the phase-shift to be positive everywhere except over vacuum (approximately true for most samples), or by performing an independent “shadow-image” calibration using a defosued ronchigram (Savitzky et al., 2021). While solving for the relative rotation angle is only somewhat important for iterative DPC, ptychographic reconstructions are especially sensitive to this value (see section 4).

Figure 2 illustrates our iterative DPC approach with a simulated double-walled carbon nanotube (CNT) structure (fig. 2b). We introduce a -17^∘ rotational offset between real-space and reciprocal-space and transpose the two axes of the diffraction intensities. This transpose is common to many experimental set-ups depending on the detector readout and STEM scanning direction, and can be solved for using with the same curl minimization approach. The rotationally-aligned and transposed CoM images are shown in fig. 2c-d. The curl minimization algorithm recovers the proper alignment (fig. 2e), as seen by the CoM images’ contrast being aligned with the $\hat{x}$- and $\hat{y}$-axes respectively. The convergence profile of our DPC reconstruction is shown in fig. 2f, illustrating a large drop off in error after the first few iterations. Comparing the virtual dark-field image in fig. 2g and the iterative DPC reconstruction in fig. 2h from the same dataset, we can observe how much more information-rich the iterative DPC reconstruction is for this simulated dataset.

There are several advantages to using iterative DPC for recovering sample-induced phase-shifts. First, compared to other phase-retrieval techniques it is relatively straightforward and computationally-inexpensive, so it can be used to quickly characterize the phase-shift imparted by the sample. Moreover, similar to other phase-contrast techniques, it is dose-efficient and enables linear contrast (see section 5), although additional considerations are required for quantitative analysis. However, iterative DPC reconstructions only solve for the phase of the sample and cannot deconvolve the probe’s wavefunction. Therefore, any aberrations in the probe, including defocus, will affect the transfer of information in the reconstruction. Lastly, the pixel-size of the reconstruction is limited by the real-space step size of the STEM probe and the absolute resolution limit is set by twice the convergence semi-angle (see section 5).

We start by forming a virtual image from each pixel within the bright-field disk. Each of these virtual images will be offset from the optical axis and can be computationally aligned to the center pixel using the geometric parallax operator. In short, the operator corrects the geometric intercepts of each virtual image to be properly focused by aligning the virtual images. In our iterative algorithm, we typically start with a heavily-binned dataset to enable alignment of low-dose datasets and progressively decrease the amount of binning, using the previous iteration’s shifts as a starting point for cross-correlating the virtual images. Figure 3c shows the relative shifts of the virtual images with a high degree of binning, and as we iteratively decrease the level of binning, we recover a finely-sampled vector field used to align the simulated CNT dataset (fig. 3d). The high signal-to-noise aligned bright-field image is show in fig. 3e.

The cross-correlation vector shifts $\boldsymbol{w}(\boldsymbol{k})$ can be expressed using the Cartesian derivatives of the aberration surface

where $\chi(\boldsymbol{k})$ can be written as a function of spatial frequency:

Here, $\lambda$ is the electron wavelength, $C_{m,n}^{x/y}$ are the coefficients of the two orthogonal aberrations of radial order $m+1$ and angular order $n$ in units of Ångströms, and $\mathrm{atan2}\left(k_{y},k_{x}\right)$ is the arctangent function which returns the correct sign of the aberration axis in all quadrants. When $n=0$ the aberrations are radially symmetric (e.g. defocus, spherical aberration, etc.) and no $C_{m,0}^{y}$ term is necessary. Equations 4 and 5 form a linear system of equations, which can be used to solve for the aberration coefficients (Cowley, 1979; Lupini et al., 2010, 2016).

It is instructive to investigate the aberration coefficients fit explicitly for the lowest-order aberrations, namely defocus and astigmatism ($m_{\mathrm{max}}=1$). Let the initial pixel positions, $\boldsymbol{v}(\boldsymbol{k})$, and measured vector shifts, $\boldsymbol{w}(\boldsymbol{k})$, be denoted by the real-valued arrays $V\in\mathbb{R}^{G\times 2}$ and $W\in\mathbb{R}^{G\times 2}$ respectively, where $G$ is the number of bright-field pixels. The linear system of equations can then be expressed using the operator $H\in\mathbb{R}^{2\times 2}$:

which can be estimated using least-squares and further decomposed into rotational $U$ and radial components $P$:

These can further be used to recover estimates for the $\theta$ (see section 2), defocus $C_{1,0}$, and astigmatism $C_{1,2}^{x/y}$ coefficients:

We can use our aberration coefficient estimates to apply a contrast transfer function (CTF) correction to the reconstruction (fig. 3f) using the following phase-flipping filter in Fourier space (Lupini et al., 2016):

where $\chi^{\mathrm{even}}(\boldsymbol{k})$/$\chi^{\mathrm{odd}}(\boldsymbol{k})$ are even/odd components of eq. 5 and $\operatorname{sign}\left[\cdot\right]$ returns the sign of its argument. The complex-valued phase ramp imparted by the odd contribution ensures the filter remains conjugate symmetric, returning a real-valued signal.

Finally, we note that the virtual images can be upsampled after alignment using kernel-density estimation of the cross-correlation vector field (Ophus et al., 2016a), creating a high signal-to-noise reconstruction (fig. 3g). The final upsampled reconstruction in fig. 3g illustrates the structure of the simulated CNT dataset. This is a dose-efficient technique since it uses all the scattered information inside the bright-field disk.

More complex deconvolution of signal in the bright-field disk is possible through scattering matrix reconstruction, which has been used for 3D-reconstructions of thicker samples (Brown et al., 2022). Similarly, a larger convergence angle will improve the depth-of-field, and the parallax operator can be applied to the upsampled reconstructions to study the out of plane structure of materials (Terzoudis-Lumsden et al., 2023).

Overall, parallax imaging is a dose-efficient and relatively computationally-inexpensive approach for recovering the structure of weak-phase specimens in a defocused dataset. Unlike iterative DPC, the resolution of the upsampled reconstruction is not limited by the scan step-size, but instead is equal to twice the convergence semi-angle similarly to other phase contrast STEM methods (Ophus et al., 2016c). For low-dose experiments with weakly scattering samples, parallax additionally provides a good estimate for low-order aberrations, which may be used as input for subsequent ptychographic reconstructions. Due to the low-computational overhead of parallax reconstructions, they can also be incredibly beneficial for checking microscope alignments and sample conditions during data acquisition.

In this section, we set up the inverse problem for the simplest ptychographic formulation and adjust the models as appropriate in the following sections to capture more complex scattering physics. While as we saw in sections 2 and 3, typical 4D-STEM experimental diffraction intensities are acquired from a 2D raster scan, this is not strictly required and complicates ptychographic notation. Instead, in the following sections we consider a dataset, $I(\boldsymbol{R}_{m},\boldsymbol{k})=I_{m}(\boldsymbol{k})$, of diffraction intensity measurements collected from any $M$ probe positions, $\boldsymbol{R}_{m}$. Single-slice ptychography assumes the exit-wave, $\psi_{m}(\boldsymbol{r})$, is given by a simple multiplication of the complex-valued incoming-wave and the complex-valued sample transmission function, henceforth referred to as the probe, $\mathcal{P}(\boldsymbol{r})$, and object, $\mathcal{O}(\boldsymbol{r}-\boldsymbol{R}_{m})$, respectively, propagated to the farfield detector:

where $\Pi_{O}$ is referred to as the "overlap projection" and $J_{m}(\boldsymbol{k})$ is the estimated beam intensity at probe position $\boldsymbol{R}_{m}$. The inverse problem is to find probe and object operators which minimize the distance between the modelled and measured diffraction intensities given by an appropriate error metric such as:

From the many existing algorithms to solve this inverse problem (Fienup, 1982; Rodenburg, Faulkner, 2004; Thibault et al., 2008; Maiden, Rodenburg, 2009; Maiden et al., 2017), we have implemented stochastic gradient descent and proximal gradient solvers due to their robustness and intuitive geometric interpretation, respectively. The first step is recognizing that the Euclidean projection of the model estimate onto the measured data, i.e. the minimal modification to the current estimate of the complex-valued exit-wave to match the measured intensities, is given by replacing its Fourier amplitude with the square-root of the measured intensity while retaining its phase (Fienup, 1982):

where $\Pi_{F}$ is referred to as the “Fourier projection.” This allows us to define the appropriate inverse operators for the stochastic gradient descent and proximal gradient methods respectively:

In order to efficiently solve the system of eqs. 16, 17, 18, 19, 20, 21 and 22 the object, probe, and exit-wave functions introduced in section 4.1 need to be stored numerically as evenly sampled arrays. This introduces constraints during data acquisition and reconstruction which we outline below and ultimately control the numerical accuracy, computational cost, and resolution.

The acquired diffraction intensities are stored in the real-valued array, $I\in\mathbb{R}^{M\times Q\times Q}$, for an array of probe positions, $R\in\mathbb{R}^{M\times 2}$. Note that when applying a rotation/transpose to align real and reciprocal space coordinate systems, we transform the probe positions and not the diffraction intensities, in order to avoid resampling the measured data (see sections 2 and 3).

The reciprocal-space pixel size, $\Delta k$, defines the real-space extent, $e=1/\Delta k$, of the complex-valued probe array, $P\in\mathbb{C}^{Q\times Q}$, where $\Delta k$ has to be selected such that $e$ is large enough to comfortably fit the probe without wrap-around artifacts (see section 10.3). The largest acquired spatial frequency, $2k_{\mathrm{max}}=\Delta k*Q$, defines the reconstruction real-space pixel size, $\Delta x=1/(2k_{\mathrm{max}})$, which together with the rotated/transposed probe positions define the object field-of-view in pixels.

The probe can be initialized with an estimate of the aberration coefficients and probe-forming aperture, obtained using microscope metadata for the probe semiangle or, better, from direct measurement of the probe under vacuum. The object array is then padded by as much as half the probe extent to define the size of the complex-valued array, $O\in\mathbb{C}^{S\times T}$, such that $S,T\geq Q$ (see fig. 1d).

In evaluating expressions of the form $\psi_{m}(\boldsymbol{r})=\mathcal{P}(\boldsymbol{r})\mathcal{O}(\boldsymbol{r}-\boldsymbol{R}_{m})$, $Q\times Q$ region-of-interest patches are extracted from the larger $S\times T$ object array centered around $\boldsymbol{R}_{m}$, rounded to the nearest pixel, multiplied by the probe array, Fourier-shifted by any remaining sub-pixel offsets, and stored in the complex-valued exit-waves array, $\psi\in\mathbb{C}^{M\times Q\times Q}$. Finally, note that the diffraction intensities, probe, and exit-waves are all internally stored as top-left corner-centered arrays to enable the unambiguous centering of even and odd dimension datasets and compatibility with fast Fourier transform libraries.

Consider the same double-walled CNT dataset reconstructed in sections 2 and 3, simulated using $150\textup{\AA}$ defocus and a $1\textup{\AA}$ scan step-size. Figure 4 shows a single-slice ptychographic reconstruction in comparison to the simulated ground truth electrostatic potential. Compared with the DPC and parallax reconstructions in figs. 2 and 3, the improvement ptychography achieves over these methods, by enabling super-resolution set by the largest acquired spatial frequency, is apparent. While $k_{\mathrm{max}}$ determines the resolution limit of the reconstruction, the extent of real-space probe-overlap and the dynamic range of the diffraction intensities control the redundancy in the dataset and ultimately the quality of the reconstruction. Due to the high dose and large $k_{\mathrm{max}}$ used here, ptychography performs better than the other phase retrieval techniques considered here.

The phase contrast CTF is defined as the ratio between the measured output signal (fig. 5a, black) and the simulated phase signal (fig. 5a, gray), typically measured as a function of experimental parameters such as wave aberration coefficients or the probe convergence semiangle $k_{\mathrm{probe}}$. CTF values range from -1 to +1, producing negative values when the measured signal is inversely proportional to the input. Non-unity CTF values represent systematic phase errors inherent to the reconstruction method (fig. 5a, magenta). Assuming sign-flips due to negative CTF regions are corrected, the residual systematic errors are given by:

where $N$ is the electron fluence per probe.

The phase error $\Delta\phi^{\mathrm{noise}}$ due to finite electron fluence is estimated by taking the difference of the absolute magnitude of a reconstruction with (fig. 5a, cyan) and without (fig. 5a, black) noise. The dependence of this phase error on the simulation’s signal magnitude $\phi_{0}$ and total electron dose $N$ can be removed by calculating the normalized random phase error:

We define the normalized total phase error, equal to the systematic and noise errors combined in quadrature as:

We use this normalized error metric to estimate which phase contrast STEM method produces the lowest phase error, shown in fig. 5b-c as ternary “phase diagrams” for spatial frequencies below the diffraction limit and convergence-angle respectively. The model predicts that ptychography will produce lower phase errors for higher electron fluences, while DPC and parallax will produce lower errors at lower electron fluences for low and high defocus values, respectively. This is because at a high electron fluence systematic errors dominate, while at lower fluence the shot noise error dominates. Dashed lines have been overlaid to delineate the “phase-boundaries” between each technique, as a guide to the eye. For spatial frequencies beyond the diffraction limit, ptychography produces the lowest total error under all conditions as expected. It should however be noted that at very low electron fluences, the shot noise error in ptychography is high-enough to introduce larger phase errors than the null hypothesis of a zero phase estimate, suggesting a smaller collection angle which is subsequently zero-padded to a desired $k_{\mathrm{max}}$ might be preferable (Sha et al., 2023). Both of these observations suggest that ptychography reconstructions require a minimum electron fluence for high accuracy. Ptychography also continues to benefit from increasing electron fluence, whereas the error for both DPC and parallax phase measurements asymptotically approaches a constant non-zero value when the accuracy is limited by the systematic errors.

To investigate the information-rich ternary “phase-diagrams” shown in fig. 5b-c, we elaborate on each error component in eq. 25. For a DPC measurement, we expect the CTF to be given by the normalized intensity probe autocorrelation function:

where $A(\boldsymbol{k})$ is the probe-forming aperture, and we only consider the defocus $C_{1}$ contribution to the aberration surface. This expression is plotted in fig. 5d as a function of normalized scattering vector $k/k_{\mathrm{probe}}$ and defocus. Figure 5d clearly indicates that DPC measurements should be performed as close to zero defocus as possible, and that contrast flipping is possible at large defocus values.

For a parallax measurement, the CTF may be written as the aperture autocorrelation function multiplied by the bright-field TEM CTF (Saxton, 2000; Meyer et al., 2002, 2004; Lupini et al., 2010), equal to:

This expression is plotted in fig. 5e. The aperture autocorrelation is the ideal signal transfer for a converged probe measurement, and is also known as the “double-overlap function” (Yang et al., 2015). The bright-field CTF produces concentric rings referred to as “Thon rings”, comprising bands of alternating sign with decreasing spacing with increasing defocus (Saxton, 2000; Meyer et al., 2002, 2004; Lupini et al., 2010). Parallax reconstructions require some defocus to produce any contrast, but note that too much defocus may reduce the measured signal due to coherence limitations. Recovering the sample’s electrostatic potential from parallax measurements additionally requires “CTF-correction”, where the CTF is estimated so that spatial frequencies with negative CTF can be sign-flipped as in section 3.2.

Our numerical tests of the CTF for the GD ptychographic reconstruction algorithm shown in fig. 5f indicate that:

i.e. the algorithm correctly recovers the phase and amplitude of all spatial frequencies when reconstructing a white noise object, even those beyond the diffraction limit. Note that at the lowest spatial frequencies we estimate CTF values slightly below 1 due to slow convergence speeds. Reconstructions of the lowest spatial frequencies in phase contrast TEM have been shown to require a large number of iterations to converge (Ophus, Ewalds, 2012). Our testing suggests that given enough iterations, all spatial frequencies can be fully recovered in a noise-free reconstruction.

The normalized random phase error for DPC and parallax are analytically derived to be given by:

Figure 5g-h plots the inverse of these expressions, respectively, while fig. 5i shows the numerical estimate for ptychography. We plot the inverse errors so that the regions with a higher signal-to-noise ratio (SNR) can be easily identified. Overall, ptychography shows increased phase error relative to both DPC and parallax. This is especially prevalent at the lowest spatial frequencies, indicating that ptychography may produce more noise in the lowest spatial frequencies than DPC. Note this is not the case with parallax imaging as its CTF approaches zero at the lowest spatial frequencies. Interestingly, some spatial frequencies show reduced random error for ptychographic reconstructions using a defocused probe, with the pattern directly following the CTF calculated for parallax. However, these “Thon rings” of increased SNR are considerably sharper than those of parallax or BF TEM, roughly equal to the parallax CTF to the fourth power.

Figure 5j-l show the estimated ratio of signal to random errors as $|\mathrm{CTF}|/\Delta\phi^{\mathrm{rand}}$ for DPC, parallax, and GD ptychography, with line traces shown in fig. 5m-o respectively. The DPC SNR shown in fig. 5j and m approximately recovers the profile estimated using a different method (Yang et al., 2015). DPC produces the highest signal-to-noise close to $k=k_{\mathrm{probe}}$ and at zero defocus, falling off at higher and lower spatial frequencies.

The parallax SNR shown in fig. 5k and n approaches the ideal converged probe signal transfer as the defocus is increased, though with roughly half of the spatial frequencies transferred with rapid oscillations. As mentioned above, care must also be taken in choosing a defocus high enough to produce significant signal transfer without introducing coherence limitations.

The GD ptychography SNR shown in fig. 5l and o is much flatter than that of either DPC or parallax. Ptychography is also the only method which can recover signal beyond the diffraction limit, $k>2k_{\mathrm{probe}}$. Additionally, ptychography does not produce zero SNR at any spatial frequency, showing that it models the signal transfer in a much more accurate manner. This model also indicates that using a defocused probe for ptychographic imaging can somewhat improve the signal transfer relative to imaging in-focus, at least for low electron fluences, though this is a relatively minor difference. However, note that the overall ptychography SNR is approximately half of the mean value for DPC and parallax below the diffraction limit.

Finally, we discuss limitations of our CTF and phase error models. First, we warn readers against over-interpreting any specific numbers given by these models, but rather try to understand the general trends. In particular, real experimental samples are never equivalent to white noise pure phase objects. Atomic potentials fall off as $k^{-2}$ for high spatial frequencies, and thus we expect reaching super-resolution in GD ptychography will require much larger electron fluences than those given above. There may also be strong correlations in the electrostatic potential signals, particularly for crystalline Bragg samples, which may influence the resulting measurements. We have also not explicitly considered direct, non-iterative ptychographic methods such as single sideband ptychography, though we expect it to fall roughly inline with the DPC analysis considered above.

Each of the methods we have considered have other potential sources of error. For example DPC relies on accurate rotation calibration (section 2.2). Parallax requires accurate cross-correlation estimates between all of the virtual images (section 3.2), which can be challenging at low electron fluence, though this can partly compensated using fiducial markers. Ptychography requires accurate calibration of all of the experimental parameters, and ptychographic reconstructions may not converge if these calibrations are not sufficiently accurate (section 4.2). Ptychographic reconstructions can also fail if the data is not recorded with sufficient redundancy.

We have also ignored the influence of amplitude variations in the object due to electrons scattered outside the range of the detector, or the effect of very strong scattering. Thicker samples too will violate the projected potential assumption which is implicit in DPC, parallax, and single-slice ptychographic reconstructions. Despite these limitations, we hope the above analysis will give some insight into the transfer of information of STEM phase retrieval techniques.

We highlight the strengths and limitations of STEM phase retrieval techniques for materials characterizations using a defocused experimental dataset of gold nanoparticles on amorphous carbon film. Figure 6 illustrates virtual dark-field imaging as-well as different reconstructions using DPC, parallax, and ptychography. These images are different reconstructions from the same experimental 4D-STEM dataset. Annular dark-field and DPC (fig. 6a-b) are limited both by aberrations in the incident beam and real-space scan sampling. Thus for a defocused dataset with a 2Å sampling, these methods provide low-resolution images, where it is not possible to resolve any atomic features. Dark-field imaging, which has monotonic contrast with the projected sample potential, does accurately capture the outline of the gold nanoparticle. By contrast DPC, which has poor signal-to-noise for low-spatial frequencies (section 5), is not able to accurately describe the contrast difference between the gold nanoparticle and the amorphous carbon support.

Figure 6c highlights the advantages of parallax imaging for defocused datasets such as this one. Unlike dark-field and DPC, parallax reconstructions can correctly account for the aberrations of the incoming probe, computationally bringing the dataset in-focus (section 3). In combination with subpixel accurate cross-correlation shifts and kernel-density estimation upsampling, this enables the recovery of high-resolution features, such as gold Bragg peaks.

For this dataset, the highest resolution reconstruction is enabled using GD ptychography. Figure 6d-e show the reconstructed projected potential after 5 and 150 iterations, respectively. As discussed above, low spatial frequencies are particularly slow to converge in ptychography, so the largest difference between fig. 6d-e is that the thickness contrast from the gold nanoparticle has not converged yet. The improved convergence is further apparent by the appearance of “Thon rings” in the Fourier transform of the image. The remaining improvements are enabled by the ability to enforce various object, probe, and position constraints (sections 14.1, 14.2 and 14.3). First, we have applied position correction to remedy the scan distortion visible in the dark-field image. Second, the probe was constrained at every iteration by replacing its amplitude with a vacuum probe aperture measurement while smoothing its phase to a low-order aberration surface expansion. Finally, we also regularized the object by forcing it to be a positive potential object and applying a small amount of low-pass Fourier filtering to reduce the high-frequency noise. Figure 6 highlights the advantages of ptychography for high-resolution characterization. Importantly, our reconstruction model is flexible, allowing one to introduce various object, probe, and position regularization to arrive at a physical solution.

The advantages of ptychography for materials-science samples are further explored by characterizing a few-layer hexagonal boron nitride (hBN) sample containing vacancy defects (Meyer et al., 2009; Tran et al., 2016). Figure 7a shows a dark-field STEM image of the few-layer hBN sample. While it is possible to see lattice contrast and the large triangular vacancy defect cluster in the sample, it is difficult to resolve individual atoms, including any single point defects in the sample. The Fourier transform illustrates how limited in resolution the dark-field image is at this convergence angle, voltage, electron dose, and sampling.

Figure 7b-d compare the dark-field image to ptychographic reconstructions of a similar field of view of the same sample with the same microscope alignment. In addition to the signal-to-noise benefits of ptychography, the maximum scattering angle determines the resolution. Here, we pad the data in diffraction-space to reconstruct high-frequency information that would otherwise be aliased in the reconstruction. We highlight the differences between fig. 7b, a reconstruction without any regularization applied during the reconstruction, and fig. 7c, an optimized and regularized reconstruction(see  section 14.1). Comparing these two images, the systematic low-frequency error is reduced in the optimized reconstruction, as we would expect from prior knowledge of the sample that the contrast across the field of view be nearly uniform. Moreover, the atomic columns are sharper, which is evidenced in both the real-space image and the absence of “streaking” in the Fourier transform due to scan noise. Both the low and high spatial frequencies can be further improved in the reconstruction using a mixed-state model with 2 probes (fig. 7d), which will be further discussed in section 11. Figure 7d shows that the ptychographic reconstruction reaches an information limit of 29 pm (3.4 Å^-1), which is much higher than the resolution achieved with dark-field imaging with these parameters, making this a “super-resolution" reconstruction.

Traditionally, imaging of biological samples has been performed with parallel illumination TEM or with stained specimens to help guard against beam-induced damage. However, recent work has highlighted the possibilities for STEM phase retrieval to characterize vitrified hydrated proteins, organelles and whole cells (Kucukoglu et al., 2024; Seifer et al., 2024; Yu et al., 2022, 2024). In particular, STEM phase retrieval methods produce high contrast micrographs, which have been proposed to be especially beneficial for very small proteins or thick biological assemblies (Rez et al., 2016; Seifer et al., 2024). In this section, we highlight how the dose-efficiency of the phase retrieval methods we have outlined above, which have predominately been used in the physical sciences, can be applied to study biological systems.

Figure 8 plots the parallax and ptychographic reconstructions of a sample composed of virus-like particles (VLP) embedded in vitreous ice. For low electron dose experiments such as these, it can be challenging to properly converge the parallax reconstruction which relies on the accurate alignment of virtual bright-field images. This can, however, be achieved by running initial iterations at high real-space binning values, which are then used to seed lower bin alignments iteratively. Due to the sub-pixel accuracy of the calculated cross-correlation shifts, we use KDE upsampling (section 3.2) to obtain a high-resolution phase reconstruction of the VLPs. While the reconstruction in Figure 8a shows some low-spatial frequency artifacts, these can be removed with filtering (Figure 8c).

For the ptychographic reconstruction, we optimize parameters to produce the ptychographic potential image in Figure 8b. For physical-sciences samples we often constrain the object to be a positive potential object. However, for biological samples plunge-frozen in vitreous ice, we do not use the positivity constraint due to the outline of small negative phase shifts due to the gap between the protein and water molecules (Shang, Sigworth, 2012). This may lead to a noisier reconstruction, but allows us to recover the structure of weakly scattering objects that may fall into the noise floor and otherwise be cut off with a strong positivity constraint. Similar probe constraints are used for biological and materials samples. After reconstruction, the ptychographic reconstruction can also be filtered to improve contrast (Figure 8c), which provides an image with similar resolution to parallax imaging as highlighted by the Fourier transforms of the filtered reconstructions in Figure 8e-f.

It can be challenging to compare the resolution of the parallax and ptychography in fig. 8 from 2D micrographs alone. While the ptychographic potential image has higher sampling than the parallax reconstruction, the cone-weighted Fourier transforms highlight that there is little transfer of information beyond the semi-convergence angle in parallax. Further work is required to compare the trade-offs between these two techniques for characterizing biological structures (Kucukoglu et al., 2024).

There are two ways to increase probe-overlap, while maintaining the same real- and reciprocal-space sampling: i) decreasing the scan step-size between adjacent probes, or ii) increasing the size of the real-space probe, e.g. by defocusing further or using a smaller probe convergence angle. The former strategy is the most straight-forward, but necessitates larger acquisition times and dataset sizes. The latter strategy can be quite effective however is limited both in using larger defocus, which requires larger region-of-interest probe arrays to avoid wrap-around artifacts, and smaller convergence angles, which change the effective transfer of information.

Figure 9 illustrates the dependence on reconstruction quality by starting with the same dataset as fig. 4 (shown on the top-right tile) and decreasing the electron dose moving vertically down and decreasing the probe overlap moving horizontally to the left. While the reconstruction quality deteriorates along both these axes, it is interesting to note the difference in the two failure modes. Decreasing the probe overlap at high dose (top row moving left) results in high-frequency artifacts, while decreasing the electron dose at large overlap (rightmost column moving down) results in substantial loss of atomic features. It should also be noted that achieving both high electron-dose and large probe overlap is often experimentally challenging for dose-sensitive samples.

Electron ptychography is a high-dimensional, non-convex inverse problem over the complex numbers. As such, rigorous proofs of convergence and comparisons between different algorithms are few and far between. In this section, we instead comment on the use of batching in the stochastic gradient descent solver, and various well-known projection-set algorithms based on empirical observations.

First, note that when the batch size in eqs. 16, 17 and 18 is equal to a single probe position, $\mathcal{M}=1$, the stochastic gradient descent algorithm reduces to a regularized version of the well-known extended ptychographic iterative engine (e-PIE) (Maiden, Rodenburg, 2009; Maiden et al., 2017), and when $\mathcal{M}=M$ the algorithm is no longer stochastic. For simulated datasets with high electron dose, we find that using a vanishing normalization parameter, $\alpha\to 0$, achieves the lowest error. For low-dose experimental datasets, we almost always observe better behavior as $\alpha\to 1$. While in general we find that using as many probe positions as the available memory permits, $\mathcal{M}\to M$, appears to be both more accurate as well as computationally more performant, we note that using fewer probe positions per iteration allows the algorithm to converge the probe faster. This is particularly important for mixed-state ptychography, which we elaborate further in section 11.

The scalar parameters $\left(a,b,c\right)$ introduced in eqs. 19, 17 and 18 describe a whole family of proximal gradient algorithms, with specific values encoding well-known named algorithms such as alternating projections (AP, $a=0,b=1,c=1$) (Neumann von, 1950), difference-map (DM, $a=-1,b=1,c=2$) (Elser, 2003), relaxed averaged alternating reflections (RAAR, $a=1-2\gamma,b=\gamma,c=2$) (Luke, 2004), and relax reflect reflect (RRR, $a=-\gamma,b=\gamma,c=2$)(Elser, 2017). We make the following observations, summarized in fig. 10:

1. i.)

   The AP (fig. 10b) or “error-reduction” algorithm is very similar to the gradient descent algorithm when $\mathcal{M}=M$ and $\beta=1$ (fig. 10a).

2. ii.)

   While the DM algorithm (fig. 10c), tends to diverge for most datasets we have investigated, variants of the algorithm with relaxation parameters namely RAAR (fig. 10d) and RRR (fig. 10e) show better convergence behaviour.

3. iii.)

   The space of parameters defined by $\left(a,b,c\right)$ lends itself well to dataset-specific hyper-parameter tuning. Figure 10g uses Bayesian optimization with gaussian processes to arrive at a set of parameters with better convergence properties $\left(a=-0.34,b=0.47,c=1.1\right)$. This optimization is further discussed in section 14.4.

Finally, we note that in all the experimental datasets we have encountered so far - we find stochastic GD reconstructions to be superior to proximal gradientmethods, mostly due to the flexibility afforded by the choice of batch-size.

The 2D raster-scan pattern diffraction intensities are typically acquired on during 4D-STEM experiments results in a periodic “gridding” artifact, also called raster scan pathology (Fannjiang, 2019; Fannjiang, Chen, 2020).

First, note that this arises only in the case of “blind” ptychography – i.e. when both the probe and the object need to be solved simultaneously. It follows that one of the ways to avoid the gridding artifacts is using perfect knowledge of the illuminating probe and only solve for the reconstructed object. While this is not practical for experimental datasets, it can be approximated by recording a high-quality vacuum probe measurement to use as the aperture, and estimating the aberration coefficients using parallax as discussed in section 3. While this can indeed help converge the reconstruction, in particular for low-dose datasets, we note that any misalignment of the vacuum probe measurement and the experimental datasets will result in model-mismatch which can also result in gridding.

Gridding can be suppressed with increased data redundancy. As discussed in section 10.1, this can be achieved using either fine scan-step sampling or larger probe sizes, e.g. by introducing defocus. As discussed above, we warn the reader about the latter strategy as the windowed object patch size, set by the detector pixels and maximum scattering angle, needs to be large enough to comfortably fit the long probe tails. Periodic boundary conditions imposed by the windowed Fourier transforms result in self-interaction of the probe which again introduce model-mismatch and ultimately gridding artifacts.

Gridding can also be partly ameliorated by using small batch-sizes, which allow the probe to update more frequently. This gives the algorithm an opportunity to self-correct model-mismatch and climb out of local minima. Taking this a step further, mixed-state ptychography (section 11) can be used to include a sacrificial low-intensity probe to absorb these artifacts.

Object regularization includes ensuring the reconstructed object is smooth, using Gaussian or Butterworth filters in real- and reciprocal-space respectively, as well as total variation denoising. These smoothing filters can be applied both in the plane of the object and along the beam-direction for multi-slice and joint ptychographic tomography. For accurate recovery of slices, multi-slice ptychography in particular requires regularization along the beam-direction (Chen et al., 2021). The strongest such regularization is that all slices have identical phase contributions. This strict assumption can be relaxed by instead enforcing first-order total variation denoising, allowing for non-identical slices which ideally approach scattering from the substrate/vacuum at large distances from the sample.

Separating the phase and amplitude components of the complex-valued reconstructed object is a challenge in electron ptychography. For weakly-scattering objects (or equivalently sufficiently thin slices of strongly-scattering objects in multi-slice ptychography), the amplitude of the reconstructed object should vary negligibly from unity. As such, we can reduce the solution-space dimensionality by a factor of two by employing the pure-phase approximation and constraining the amplitude to unity at each iteration. Taking this a step further, the object itself can be stored as the real-valued potential directly, $\mathcal{O}(\boldsymbol{r})=\mathrm{e}^{\mathrm{i}V(\boldsymbol{r})}$. In addition to the memory efficiencies this affords, it further allows us to impose additional regularization constraints on the potential directly. Notably, atomic potentials are, by nature, non-negative and thus we can enforce positivity by clipping negative values to zero. Moreover, positivity can be combined with shrinkage filters, whereby a constant phase-shift is subtracted from the current estimate at each iteration prior to zero-clipping, to promote atomicity by forcing the background to zero. We note however that this constraint is rather aggressive, especially in early iterations, and can result in artificially suppressing scattered signal.

Similarly, the current probe estimate can be constrained using physical regularizations, such as a fixed probe aperture, a smooth phase aberration surface, and center-of-mass fixing. As mentioned in section 4.2, the probe aperture can be initialized using either a vacuum measurement or an estimate of the probe semi-angle, i.e. a perfect circular aperture. In general, a high signal-to-noise vacuum measurement of the probe aperture is important for probe regularization and calibration. This can be acquired by averaging several diffraction intensity measurements over a vacuum region with the same experimental parameters as the 4D-STEM dataset, albeit with a smaller step-size to avoid descan artifacts.

A common artifact in ptychographic reconstructions is probe-object mixing (fig. 14a), where object features from the scattered diffraction intensities appear in the Fourier amplitude of the probe. This can be ameliorated by replacing the Fourier amplitude of the probe estimate using the square-root of a vacuum measurement, effectively reducing the solution-space dimensionality by a factor of two. If such a measurement is not available (and cannot be estimated from a vacuum region inside the field of view), then the Fourier amplitude of the probe estimate can be fitted at every iteration to an angularly-varying sigmoid function to suppress intensity outside the probe aperture and ensure a constant top-hat intensity inside the aperture (fig. 14c-e).

Moreover, the probe aberration surface is expected to be smoothly-varying. This is not guaranteed a-priori (fig. 14b), and can instead be imposed as an iterative constraint. An elegant way of enforcing this is to fit the aberration surface $\chi(\boldsymbol{k})$ to a linear combination of low-order aberrations (Ophus et al., 2016d), by restricting the maximum radial and angular orders of the expansion in eq. 5. Figure 14f-g show the fitted aberration surface and some basis set examples respectively. Finally, we note that for this to work, the aberration surface needs to be phase-unwrapped first which can introduce additional sources of error. If an accurate initial estimate of the probe aberrations exists, e.g. using the parallax aberration coefficients fitting introduced in section 3.2, these can be removed prior to smoothing and re-introduced after to avoid having to phase-unwrap the aberration surface.

Because ptychographic algorithms rely on overlap between adjacent probe positions, accurate probe positions are required for object convergence (Maiden et al., 2012a). The probe positions estimate can be refined during the ptychographic reconstruction as well as constrained to fit a linear drift distortion model. Drift in STEM data acquisition causes the real probe positions to deviate from the ideal raster positions (Ophus et al., 2016b). In general, faster 4D-STEM data acquisition times, which are determined by detector read-out speeds, step size, and the required field of view, lead to fewer artifacts in the positions of the probes.

During ptychographic reconstruction, we can refine our probe positions estimate using a gradient step on the intensity measurements in parallel (Dwivedi et al., 2018). This results in a displacement vector for each probe position, which can be further constrained by fitting the displacement vectors to a six-parameter (x-scaling, y-scaling, shear, rotation, x-translation, y-translation) global affine transformation.

In order to successfully perform a ptychographic reconstruction of experimental data, it is necessary that certain experimental parameters be known to high precision, particularly the scan step size, the reciprocal pixel size of the diffraction patterns, and the real-space/reciprocal-space rotation angle. Other parameters, such as the defocus and other aberrations are refined during reconstruction, but a better guess still affects the convergence. In addition, reconstruction parameters such as the step size, batch size, and regularization strengths, also impact the quality of the reconstruction. In order to find the optimum values of both experimental and reconstruction parameters that give the highest quality reconstruction of a dataset, we have implemented an optimizer using the Bayesian optimization (BO) with Gaussian processes (GP) algorithm. BO-GP is commonly used in the tuning of hyperparameters for machine learning (Bergstra et al., 2011), and has been successfully applied to ptychography both for parameter tuning (Cao et al., 2022) and directly as a procedure for iterative phase retrieval (Pelz et al., 2017). The BO-GP algorithm models the function to be minimized as a mixture of Gaussian processes, and chooses evaluation points through a combination of exploration of unsampled regions of the parameter space and exploitation of local minima of the GP model. Our implementation provides a simple interface for specifying optimization parameters, and the optimization is performed by the BO-GP routine provided by scikit-optimize (Head et al., 2021). Various reconstruction quality metrics can be used as the objective function of the optimizer, including the final data error of the reconstruction as well as the contrast, entropy, or total variation of the reconstructed object.

While the ideal experimental parameters depend strongly on the research question, the microscope, and the specimen, there are a few suggestions that may be helpful for most phase retrieval experiments. First, we recommend starting by simulating a proxy of the system to test the range of experimental parameters that may be possible with a given set-up. For example, the number of pixels in available detectors, combined with the convergence semi-angle limit the maximum defocus that is possible. When recording diffraction data, it is important the detector can sample the diffraction intensities finely enough to recover features of the specimen in the ronchigram, which further set constraints on the maximum possible defocus value.

In designing an experiment, it is also important to consider realistic limitations of the microscope and specimen. For example, although the methods introduced above can provide some degree of probe-position refinement, the reconstruction will be severely limited with substantial scan distortion. Therefore, the number of probe positions in a scan should be determined based on stability, magnification, and detector speed. Similarly, while ptychography and parallax can correct for aberrations in the incident beam, when possible it is better to align the microscope as well as possible before starting the experiment, to remove higher-order residual aberrations. There can be other sources of artifacts in experimental data that are not usually accounted for ptychographic models, such as incoherence of the source, post-specimen aberrations introduced by the beam through a spectrometer, and dead detector pixels. All of these should be minimized when possible. This may require using a smaller beam current to minimize incoherence, aligning the spectrometer, dark references, and preprocessing the diffraction data. Preprocessing may include weighted-binning to account for dead pixels, hot- and cold-pixel filtering, and gain correction.

Further, calibration scans can be helpful for initializing parameters for phase retrieval experiments. These scans may include experiments on samples, such as the gold nanoparticles shown in fig. 6, which provide a large phase shift, used to calibrate the rotation between real and reciprocal space. As described above, it is important to use the right range of defocus values for ptychography experiments. However, the nominal microscope defocus value and true defocus value (and sign) might not agree. Therefore, acquiring calibration scans at a few defocus values can solve for the relationship between indicated and true defocus. Calibrations may be similar day-to-day on the same microscope especially for similar alignment conditions. Therefore, these measurements are most important for initial phase retrieval experiments on a microscope.

As highlighted a few times in this manuscript, vacuum scans are invaluable for ptychography experiments. For initializing the probe for ptychography and possibly constraining it during the reconstruction, a vacuum scan should produce a high signal-to-noise image of the aperture. Ideally this is done with the same diffraction shift, so the vacuum and data are centered on the detector. A vacuum scan does not requre as many images as a full 4D-STEM scan, since all the scan positions are averaged together to create one 2D vacuum probe. Instead it may be better to take fewer scan positions at a higher dose, if allowed by the detector. Although descan can be accounted for computationally after the experiment, it is better to minimize this artifact in vacuum measurements by taking the scan at as high magnification as possible, or keeping the STEM probe stationary while recording these images. It may also be helpful to take a vacuum scan with a very large camera length, so the probe is evenly illuminating the entire detector. This flat-field reference can be used for gain correction of the 4D-data. Lastly, when encountering any unexpected challenges on the microscope, it is always a good idea to consider recruiting a few additional rubber ducks to help you debug the issue.