1. Introduction

Ptychography recovers the phase shift of scattered radiation within a probed specimen from a series of diffraction patterns [1]. It evolved due to improvements in computing power, detector hardware and algorithm development from a technique with relatively limited practical application to a mature technique that is popular in light microscopy [2], the X-ray community [3–5] and electron microscopy, where it has recently reached sub-angstrom resolution [6]. Its principle is based on illuminating the specimen with a spatially limited probe that is shifted such that an overlap between illuminated areas provides information redundancy.

One of the main benefits of ptychography is that spatial resolution is not limited by the imaging optics [1]. However, in practice, it comes at the price of a challenging task for the reconstruction algorithm to solve the phase problem. While the algorithm must recover the phase, the quality of the reconstruction depends strongly on the precise knowledge of the shape [7–10], the coherence [11,12], and the positions [13–15] of the probe. If these parameters are not initially known, they are required to be retrieved by the algorithm during the reconstruction process. Initially, ptychography relied on the projection approximation, which restricted the method to relatively thin samples [1,8]. To account for multiple scattering in thicker samples, a multislice method can be adopted [16–22]. In cases where the noise level in the diffraction data is too high and/or the oversampling limit is not reached, ptychography reconstruction algorithms tend to either converge to the wrong solution or are too unstable to reach convergence [23].

A variety of ptychographic reconstruction algorithms exist, utilizing either iterative [9] or direct inversion [24] to solve for the phase of the object, and sometimes for the probe shape and/or positions [18,25]. These iterative algorithms have introduced gradient based optimization [26] (including the use of finite difference derivatives in the optimization of the beam positions [27]), multislice reconstructions [17,18] as well as regularization [17,26,28], but these innovations have not been previously combined.

We introduce regularized optimization for ptychography (ROP), which is a multislice-capable conjugate gradient descent reconstruction algorithm that can reconstruct the phase while simultaneously retrieving the probe shape and positions. ROP is an extension to the inversion of dynamical electron scattering (IDES) algorithm, which has been successfully demonstrated to achieve inversion from simulated high resolution transmission electron microscope images, diffraction data and experimental optical Fourier ptychography data [11,16,17,29]. The reconstruction algorithm adapts the design of an Artificial Neural Network (ANN), a model very well established in the machine learning community, to efficiently calculate the gradients and the reconstructed phase, which can be regularized during optimization to ensure sparseness.

In this paper it is shown that the conditioning of the reconstruction problem, characterized by the oversampling ratio [30], gets worse for desirable experimental conditions such as smaller observed diffraction patterns and larger step sizes, which can increase recording speed, field of view, or reduce the amount of acquired data, and more convergent electron beams, which improve resolution. A smaller beam support which improves computation speed also worsens the conditioning.

The power of using ROP is demonstrated by relaxing the experimental requirements and significantly reducing the oversampling ratio to even below unity, where the problem becomes underdetermined. It is thus shown that the incoherent resolution can be reached from just the intensities in the central diffraction discs, due to the regularization better conditioning the phase problem. Furthermore, only partially known experimental conditions such as the beam shape and the beam positions can be treated, as can relatively thick specimens due to the incorporation of the multislice algorithm. These improvements will enable the application of ptychography at higher frame rates and scanning speeds, more convergent probes, and larger fields of view. The reduced data load per recording will furthermore facilitate a fast computation of the reconstruction.

The paper is organized as follows. In Section 2. the image formation in electron ptychography and its impact on the inverse problem are treated. Furthermore, an overview of the settings for the simulations and the experiments is given. In Section 3. the probe and position correction are demonstrated on simulated ${\textrm {Mo}\textrm {S}_2}$ data. Then, the effect of regularization on underdetermined data is investigated, and reconstructions from experimental ${\textrm {Mo}\textrm {S}_2}$ and ${\textrm {Nb}_3\textrm {Cl}_8}$ data are presented. Finally, the conclusions are drawn in Section 4.

2. Methods

2.1 Image formation in electron ptychography

We consider the multi-slice method for the scattering process of an electron beam that passes through a sample. In this framework, the Schrödinger equation has a solution that is expressed for an electron wavefunction that evolves in a specimen which has been sectioned into Z multiple slices. The wavefunction $\psi ^{z}(\vec {r})$ is alternatively transmitted through a slice zand propagated to the next slice $z+1$ by the following equation:

The main variables of the ptychography set-up are illustrated in Fig. 1. The beam is scanned with a step size of $\Delta x$, and has a real-space support of width w, divided into m pixels of width d. This translates into an $m\times m$ diffraction pattern in reciprocal space with pixel size $\delta = 1/w$. The experimentally recorded diffraction patterns are $n \times n$ pixels wide, with $n \leq m$, and the variable ${\theta _{\textrm {obs}}}$ is their width. Furthermore, ${\theta _{\textrm {con}}}$ stands for the beam’s convergence semi-angle. Reconstruction aims for a resolution r, corresponding to the angular frequency ${\theta _{\textrm {res}}}$. As detailed in Appendix A, when single scattering is dominant, ${\theta _{\textrm {res}}}$ gets smeared out between ${\theta _{\textrm {res}}} - {\theta _{\textrm {con}}}$ and ${\theta _{\textrm {res}}} + {\theta _{\textrm {con}}}$. This implies that ptychography is able to achieve the incoherent resolution limit of $2 {\theta _{\textrm {con}}}$ from just the intensities within the central disc.

 

Fig. 1. Overview of the parameters of the ptychography set-up, in real and reciprocal space.

Download Full Size | PDF

Retrieving the electrostatic object potential V from the measurements is a mathematical problem with a number of knowns, $N_{\textrm {k}}$, given by the number of pixels in the measurements, and a number of unknowns, $N_{\textrm {u}}$, determined by double the number of pixels in the reconstructed complex potential. The ratio $N_{\textrm {k}}/N_{\textrm {u}}$ is the so-called oversampling ratio [30], for which in Appendix A an expression is given. In the limiting case of a large scan area, the incoherent resolution r, and thin specimens, the expression reduces to

2.2 Inverse problem

The following derivations are carried out for a single beam position and the associated diffraction pattern. The index p is thus dropped for the time being, but will be reintroduced from Eq. (15) onward to treat the interdependence. Furthermore, the variable $\sigma$ is absorbed into V for brevity.

The discrepancy between simulation and experiment of a single diffraction pattern is quantified by an error metric $\mathcal {E}$. In this paper, we define three different metrics:

Derivative based methods allow a free choice of error metric. When the noise is purely Poissonian, as is the case with direct detectors [32], the log-likelihood $\mathcal {E}_3$ is well-suited, while in [11], it was shown how $\mathcal {E}_1$ was advantageous in the presence of small modelling errors.

In addition to the error metric, we introduce a regularization term $\mathcal {R}$ that acts as a sparsity constraint on the object. Among different alternatives [17,26,28], we express that atomic resolution images are sparse in the pixel basisNew, as the atoms are separated by regions of zero or low intensity, and thus define $\mathcal {R}$ by the $\ell _1$-norm of the projected potential:

 

Fig. 2. Schematic of the network model. The core ANN is divided into sub-networks by the number of probe positions. For each sub-network, the impinging wave $\psi ^{z}$ is multiplied with a transmission function $t^{z}$. The next layer of edges and nodes encodes a real-space convolution with the Fresnel propagator, resulting in $\psi ^{z+1}$. This is repeated until the exit wave $\psi ^{Z+1}$ is produced and the intensity is taken in the far-field.

Download Full Size | PDF

The derivative of the error w.r.t. the atomic potential at each slice is given by:

The derivative of the loss $\ell$ with respect to the incoming probe, $\psi ^{0} = \psi ^{0}_{re} + i \psi ^{0}_{im}$, is calculated with a translated object instead of a translated beam to ensure that the result is centered on the origin. Following Appendix C we have:

The forward propagation of the electron wave and the subsequent backpropagation of the loss are computed independently for each diffraction pattern p. The updates for atomic potential V and incoming beam $\psi ^{0}$ are then performed with a loss function and gradients that are respectively formed from the entire batch of diffraction patterns P,

The projected potential of the object, the probe shape and the probe positions are retrieved through a non-linear conjugate-gradient method, making use of Polak-Ribière search-directions d that depend on the batched gradients and the search-directions of the previous iteration [34–36]:

More details on the technical implementation are given in Appendix D.

2.3 Settings for simulated and experimental data

To investigate the performance of the algorithm and the regularization effects on the reconstruction, we conducted one simulation and one experiment on both a ${\textrm {Mo}\textrm {S}_2}$ and a ${\textrm {Nb}_3\textrm {Cl}_8}$ specimen. The first simulation tested for successful convergence of the algorithm when either probe or positions were distorted and the result was supported by an experiment. For the simulation, the ptychographic data of a ${\textrm {Mo}\textrm {S}_2}$ monolayer was generated using the code described in [37] and with parameters that were mostly chosen to match those used in the experiment [6]. This data, as well as all following simulated data, can be found in Dataset 1 [38]. This experiment, also investigating a ${\textrm {Mo}\textrm {S}_2}$ monolayer, was performed on the Titan Themis 300 microscope, operated with an acceleration voltage of 80 kV ($\lambda = 4.18$ pm) and a convergence semi-angle ${\theta _{\textrm {con}}}$ of 21.4 mrad. The EMPAD direct electron detector [32] recorded $124 \times 124$ pixel diffraction patterns and while a $87 \times 51$ probe position scan with a step size $\Delta x$ of 0.021 nm was used in the experiment, a scan of $100 \times 100$ probe positions was chosen in the simulation.

Our second simulation and experiment aimed at improving the resolution of the reconstruction for limited data from a ${\textrm {Nb}_3\textrm {Cl}_8}$ specimen. A resolution corresponding to twice the beam convergence semi-angle, ${\theta _{\textrm {res}}} = 2 {\theta _{\textrm {con}}}$, was aimed for. During reconstruction the angular frequencies were calculated out to ${\theta _{\textrm {cal}}} = 3 {\theta _{\textrm {con}}}$. In this case, a measure of spatial resolution is the minimum resolvable distance between two atoms.

${\textrm {Nb}_3\textrm {Cl}_8}$ is a trigonal system with lattice parameters $a = b = 6.831$ Å and $c = 13.75$ Å. In the $[001]$ orientation Nb-dumbbells with an atom spacing of 0.67 Å can be observed. This is shown in Fig. 5(c). At an acceleration voltage of 120 kV, equivalent to $\lambda = 3.35$ pm, the dumbbell spacing corresponds to a scattering angle of 50 mrad.

We simulated ptychographic data using the code described in [37], with a semi-convergence angle, ${\theta _{\textrm {con}}}$, set to 28 mrad so that when $\rho = 2$ the angular frequency of the dumbbell falls well below ${\theta _{\textrm {res}}} = \rho {\theta _{\textrm {con}}} = 56$ mrad. The number of pixels in the beam support of the simulation was chosen to be four times higher than needed for the reconstruction to avoid the so-called “inverse crime” fallacy [39]. Because this results in a four-times smaller pixel size in reciprocal space, a four-by-four binning of the simulated diffraction patterns was done before reconstruction. The number of slices in the multislice simulation was set to 35 with a thickness of 0.1338 nm per slice, thus covering the crystal’s vertical extent of $3c = 4.125$ nm. The scan pattern followed a Halton sequence to avoid artifacts in the reconstruction that could originate from the translation symmetry of a grid scan pattern which has been dubbed “raster grid pathology” [23]. A total of 6,600 points within a disc of diameter 4.19 nm were probed, resulting in an average distance between neighboring beam positions of $\Delta x = 0.0457$ nm.

The desired resolution corresponded to $2 {\theta _{\textrm {con}}} = 56$ mrad, which is therefore sufficient to resolve the dumbbells at 50 mrad. The diffraction patterns only contained the central disc, i.e. ${\theta _{\textrm {obs}}} = {\theta _{\textrm {con}}}$. The binning resulted in a beam support of $w = 0.232$ nm wide and diffraction patterns of $n = 4$ pixels, and $m = 20$. The resulting oversampling ratio was 0.47, or underdetermined by a factor of over two.

In addition, we generated simulated data similar to the aforementioned data, differing only in the number of beam positions by a factor of four, totalling 26,400 positions. The diffraction patterns were subsequently binned down by a factor of two, resulting in an oversampling ratio of 6.6. The reconstruction from this strongly overdetermined and noiseless data acted as ground truth in the evaluation of the other reconstructions. The settings of the simulation and reconstruction are detailed in Table 1.

For the experiment, a 4 nm thick ${\textrm {Nb}_3\textrm {Cl}_8}$-crystal was produced through exfoliation and transfer on a SiN TEM grid with holes. It was investigated on the Titan Themis 300 microscope, with an acceleration voltage of 120 kV ($\lambda = 3.35$ pm) and a convergence semi-angle ${\theta _{\textrm {con}}}$ of 24 mrad. $124 \times 124$ pixel diffraction patterns were recorded on the EMPAD direct electron detector, with a $64 \times 64$ scan and a scan step size $\Delta x$ of 0.043 nm.

From these raw data, three reconstructions were calculated, the first from the unaltered data, and then two from reduced data sets. For the first reduced data set (the second reconstruction) the diffraction patterns were binned $3 \times 3$, and only the central $12 \times 12$ pixels were retained, thus yielding diffraction patterns out to ${\theta _{\textrm {obs}}} = 36.8$ mrad, and a beam support of $w = 0.571$ nm. Since ${\theta _{\textrm {obs}}} + {\theta _{\textrm {con}}} = 60.8$ mrad encompasses the dumbbell signal at 50 mrad, resolved dumbbells were expected. The second reduced data set (the third reconstruction) is obtained by a $6 \times 6$ binning, retaining only the central $14 \times 14$ pixels. This yielded a beam support of $w = 0.285$ nm and diffraction patterns out to ${\theta _{\textrm {obs}}} = 82.3$ mrad, well beyond the farthest extent of the dumbbell signal at 50 mrad $+ {\theta _{\textrm {con}}} = 74$ mrad, and hence resolved dumbbells were expected. The three data sets had an oversampling ratio of 70, 5.6 and 7.5, respectively. All settings are summarized in Table 1.

3. Results

3.1 Probe and position correction on simulated MoS$_2$ data

We investigated the performance of the update functions (20) and (21) on the simulated MoS$_2$ data described in 2.3. For all the tests, optimization was done for 400 iterations, regularization has been neglected, the initial atomic potential was assumed to be vacuum and since a single layer MoS$_2$ sample exhibits weak phase difference, the depth of the ROP model has been set to one slice.

For the analysis of the probe shape optimization, a false incoming probe that had been deviated in defocus, relative to the correct probe, was read in. In Fig. 3, we apply the optimization for incoming probes that deviate in defocus by up to 12 nm. The correct probe positions were used, leaving the algorithm to alternate only between the object and probe shape update function. Comparing different optimization strategies, a sub-epoch of 5 iterations for the object update function and 10 iterations for the probe shape update function achieved the lowest root mean squared error (RMSE) between the probe shape of the updated probe and the correct probe. We further chose the error metric $\mathcal {E}_1$ and the initial step size $\alpha _0$ for the object update function was set to 1, while the initial step size $\beta _0$ for the probe shape update function was set to 1e-5. Using the aforementioned parameter settings, Fig. 3(a), shows the intensity line profile of a probe with a defocus of 10 nm that was taken as the initial guess (Probe A), the updated final probe (Probe B) and the correct probe that was used to generate the data (Probe C). The RMSE of the probe shape between Probe B and Probe C was 2.05e-6. Figure 3(b) shows for all the test cases the RMSE of the probe shape between the probe at each iteration of the optimization process and Probe C.

 

Fig. 3. a) Comparison of the intensity profile between the correct incoming probe, a false probe (i.e. defocused by 10nm) and the updated probe. b) RMSE between the correct incoming probe and differently strong deviated incoming probes determined by different defocus values.

Download Full Size | PDF

In our second analysis, we focused on the probe position optimization. Here, the reconstruction algorithm started the optimization with probe positions (Positions A) that were randomly deviated from the correct positions (Positions C) and finally converged to the updated probe positions (Positions B). The performance was then evaluated by taking the RMSE between the Positions at every iteration and Positions C. In this analysis, the correct probe shape was used and therefore only update functions (19) and (21) were alternated. Sub-epochs of 3 iterations for the object update function and 4 for the probe position update function gave the best result. The initial step size $\alpha _0$ for the object update function was set to 1e 3 and the initial step size $\gamma _0$ for the probe position update function was set to 1e 5. The optimal error metric appeared to be $\mathcal {E}_2$, showing a much lower final RMSE than $\mathcal {E}_1$. Figure 4(a) shows the Positions A for a subsection of the scan grid and an averaged deviation of $0.54\Delta x$. They are compared to the Positions B and the Positions C. The trajectories they form during the optimization process are indicated by a dashed blue line. In Fig. 4(b), we show the optimization for cases with initial mean deviation starting from $0.54 \Delta x$ up to $1.92 \Delta x$ and compared the RMSE during the optimization process.

 

Fig. 4. a) Update map of probe positions that have been deviated on average by $0.54 \Delta x$ (with $\Delta x = 0.021\textrm {nm}$) as they converge to the true probe positions. b) RMSE between the correct probe positions and differently strong randomly deviated probe positions.

Download Full Size | PDF

3.2 Effect of regularization on under-determined data

In this section, we analyzed the influence of the regularization $\mathcal {R}$ on the reconstructed object potential. In particular, we concentrated on the simulated and underdetermined ${\textrm {Nb}_3\textrm {Cl}_8}$ data set, with $4\times 4$ diffraction patterns and an oversampling ratio of 0.47 as summarized in the column “Sim. Rec.” in Table 1. We created 7 data sets that differed by the electron count per pixel in the central disc, further denoted as electron count. We increased the electron count from 100 to 1000 and restricted the analysis to only update function (19). The optimization has been limited to 100 iterations and again the initial atomic potential was chosen to be vacuum. The thickness of the ${\textrm {Nb}_3\textrm {Cl}_8}$ specimen has been covered in our model by using 7 slices, each separated by 0.669 nm and equally updated according to Eq. (11). The initial step size $\alpha _0$ for the object update function was set to 1 and we selected the error metric $\mathcal {E}_2$.

In Fig. 5 we demonstrate the influence of the regularization on the reconstruction quality with a particular focus on the data set that has an electron count of 316. Figure 5(a) shows that for each of the 7 data sets, a different regularization parameter $\mu$ is required. The higher the electron count in the data, i.e. the less noise present, the less regularization is needed for an optimal reconstruction. Here, we estimated the optimal $\mu$ based on the RMSE between the final reconstruction of each test case and a reconstruction from the over-determined data set, acting as a ground truth. Tested regularization parameters covered a range from 5.18 to 2.68e2 and the optimal $\mu$ was determined as the closest sampling point to the local minimum of a 3rd-order polynomial function that was fitted to a sub-region of the test range. To increase the accuracy of our fit function, we sampled more densely the range close to the local minimum. Figure 5(b) exemplifies the parameter estimation based on the RMSE for the data set with an electron counts of 316. The parameter resulting in the lowest RMSE was found to be $\mu = 3.20$e1. Figure 5(c) shows the reconstruction from the ground truth data that is compared to the reconstructions of the test cases. Figure 5(d) shows a strongly under-regularized reconstruction, corresponding to $\mu = 5.18$, Fig. 5(e) shows an optimally regularized reconstruction, corresponding to $\mu = 3.20$e1 and Fig. 5(f) shows a strongly over-regularized reconstruction, corresponding to $\mu = 2.68$e2.

 

Fig. 5. Evaluation of the optimal regularization for (a) the 7 simulated data sets with different noise levels of a Nb$_3$Cl$_8$ specimen, based on (b) the RMSE between the respective regularized reconstruction and (c) the ground truth reconstruction. For the case of an electron count of 316, marked with an arrow in (a), the object potential is reconstructed d) under-regularized with $\mu = 5.18$, e) optimally regularized with $\mu = 3.20$e1 and f) over-regularized with $\mu = 2.68$e2.

Download Full Size | PDF

Although $\mu$ was determined through a sweep, other approaches have been considered in [40,41]. Furthermore, in Appendix B a similar analysis was carried out showing reconstructions from underdetermined data with a total generalized variation regularization [41,42].

3.3 Reconstructions from experimental MoS$_2$ and Nb$_3$Cl$_8$ data

Figure 6 shows the result from the MoS$_{2}$, focusing on the different parts of the algorithm. The settings are shown in Tables 1 and 2, where $\mu$ is the regularization parameter, $\alpha _0$ is the initial object update step size, $\beta _0$ is the initial probe update step size, and sub ep. $\alpha$ and sub ep. $\beta$ are the sub epoch lengths of each in iterations, respectively. Figure 6(a) shows the resulting object when only the object is updated, and the rest of the update-able options (probe positions and probe shape) are left at their defaults. This results in an object with a large background ramp of illumination from top to bottom, which mostly obscures the atomic potentials. However, fitting and removing this background plane of potential does reveal underlying atomic potentials. A possible reason for the phase ramp in the background is a sub-pixel error in positioning the diffraction patterns [43]. Figure 6(b) shows the result of then adding just the probe update, which helps to remove the ramp of illumination, but does not accurately reconstruct the expected spherically symmetric atomic potentials. Figure 6(c) shows the result when a previously optimized probe, from 6(b), is passed in as an initial argument to the reconstruction, resulting in much more symmetric atomic potentials. Figure 6(d) shows the same as 6(c), but now also allowing for the positions to be updated. This results in even more circularly symmetric atoms, particularly in the row between the arrows in 6(c), in which the atoms are elongated. The change in atom shape can be seen in the subpanels for 6(c) and 6(d), in which the average shape of the atom in the rows marked by the white arrows is shown, along with the best fit ellipse to the half maximum intensities of the class average. A clear difference in ellipse major and minor axis length can be seen in 6(c), with a much less noticeable difference in 6(d). The updated positions are shown in 6(e), showing the ability to compensate for global inaccuracies in probe position, due to imperfect rotation and pixel size scaling calibrations, as well as local inaccuracies in positions in the region corresponding to the region between arrows in 6(c) with the elongated atoms. The corresponding region to the area between the white arrows in 6(c) and 6(d) is between the black arrows in 6(e), where a clear shift in positions can be seen, possibly due to microscope instabilities. Figure 6(f) finally shows the reconstructed probe at the sample’s exit used for 6(c) and 6(d). Here we show that only through the optimization of probe positions, probe shape, and object itself do we obtain the highest quality reconstructions using ROP.

 

Fig. 6. MoS$_2$ object potential after (a) only updating the object, (b) updating the object and probe from an initial default, (c) updating the object and probe, starting with an optimized probe, but with the original grid of probe positions, and (d) with updated probe positions. Notice the elongation of the atoms between the arrows in (c), but not in (d). Inset in both (c) and (d) are the class averages of the atoms between the white arrows, overlaid with the best fit ellipses to the atoms’ half maximum intensities. (e) The updated probe positions with black arrows corresponding to the white arrows in (c) and (d), and (f) the updated probe used for the final reconstruction, where amplitude to the 1/4 power is represented by intensity, and phase is represented by the HSV colorwheel.

Download Full Size | PDF

Figure 7 shows the result from using ROP on data acquired from a 4 nm slab of Nb$_{3}$Cl$_{8}$, and the benefits of using the multislice approach. The parameters used can be seen in Tables 1 and 3, where $\Delta z_{\textrm {slice}}$ is the propagation distance between slices. In all cases, updated probe positions were used from a prior reconstruction in order to focus on the multislice capabilities of the algorithm. Figure 7(a) shows the result of reconstructing the phase object with only one slice in the beam direction. This ignores all aspects of multiple scattering, and results in the dumbbells being unresolved. In the unit cell class average [44] (lower right), the location of the dumbbells are resolved, but not the spacing between the atoms. Figure 7(b) shows the same result, but this time having 6 slices, with a fixed propagation distance of 0.66875 nm, corresponding to half the vertical lattice parameter, thereby matching the sample’s total thickness of 4 nm. The dumbbells are sharper, resolvable in individual cases, and can be seen in the class average. Figure 7(c) and 7(d) both show data reduction techniques, in which the acquired diffraction patterns were binned and cropped around the central beam. In Fig. 7(c), in which diffraction patterns are binned by 3 and then cropped to $12\times 12$ pixels, the dumbbells are still visible in both the class average and object itself. However, in Fig. 7(d), where the patterns are first binned by 6 and then the central $14\times 14$ pixels were cropped, the dumbbells again are no longer resolvable, although like in 7(a), their positions can be found. Fig. 7(e) shows the Fourier Transform of 7(c), with 6-fold symmetric reflections around the dotted circle indicating that the desired resolution to resolve the dumbbells was achieved. In Fig. 7(f), the Fourier Transform of 7(d), 6-fold symmetric spots are only observed corresponding to a maximum resolution of 0.84 Å when our reduced settings using sixfold binning are employed. Both Fourier Transforms were averaged using a six-fold symmetry operation and weighted with a $r^{0.4}$ strength ramp to emphasize the signal from the noise. As Fig. 7(c) and 7(e) show, the potential for much faster data acquisition as well as reconstruction exists due to less pixels being required, when appropriately chosen acquisition settings are used.

 

Fig. 7. Nb$_3$Cl$_8$ object potential updating the probe and object with (a) the full data set, and one slice, (c) a 3$\times$3 binned data set, with the center 12$\times$12 pixels cropped and a six slice multislice reconstruction, (b) the full dataset with a six slice multislice reconstruction, and (d) binned by 6, with center 14$\times$14 pixels cropped, also with six slices. Inset in the top right corner are representative diffraction patterns, and in the bottom right corner, class averages of the objects’ unit cell. (e) is the Fourier Transform of (c), with the dashed black circle drawn corresponding to the dumbbell resolution of 0.67 Å, and (f) is the Fourier Transform of (d), with white circles drawn around two representative spots corresponding to a 0.84 Å resolution.

Download Full Size | PDF

4. Conclusions

In this paper we present the regularized optimization for ptychography (ROP) algorithm that uses derivative-based non-linear conjugate gradients optimization with Polak-Ribière search directions for the inversion. Three error metrics were used: the sum of absolute differences, the sum of squared differences and the negative log-likelihood for Poisson noise.

Through incorporation of the multislice formalism, multiple scattering was accounted for, resulting in an improved lateral resolution as demonstrated by its necessity in resolving the Nb-dumbbells in 4 nm thick ${\textrm {Nb}_3\textrm {Cl}_8}$ with and without data reduction.

Data compromised by partially known experimental parameters, such as beam shape and positions, was successfully treated by optimizing said parameters along with the object. Relatively large errors of up to 10 nm in defocus and up to 0.024 nm or 1.15 times the step size in beam position were shown to be recoverable.

It was shown how the experimentally favorable settings of large step sizes and convergence angles, and small beam support and width of the recorded diffraction patterns, worsen the oversampling ratio and make the inverse problem less well-conditioned. Regularization made reconstruction possible nonetheless, even from noisy simulated data that was underdetermined by a factor of 0.47. The experimental data had its overdetermination significantly reduced, from a factor of 70 down to 5.6, while still obtaining the desired resolution, and it was shown that higher noise-levels required a stronger regularization to stabilize the reconstruction.

These improvements relax the experimental requirements and hence are likely to extend the applicability of ptychography to thicker samples, and higher frame rates and scanning speeds, more convergent probes, and larger fields of view. The reduced data load per recording furthermore facilitates a fast transfer of the experimental data and computation of the reconstruction. This will help solidify the applicability of ptychography as a standard experimental technique for electron microscopy and optical setups, and increase its applicability to a broader class of microscopes and detectors.

A. Oversampling ratio for ptychography

An expression for the oversampling ratio [30] $N_{\textrm {k}}/N_{\textrm {u}}$ is derived. It depends on the experimental settings needed to attain a certain resolution r, corresponding to an angular frequency of ${\theta _{\textrm {res}}}$.

When single scattering is dominant the angular frequency ${\theta _{\textrm {res}}}$ gets smeared out between ${\theta _{\textrm {res}}} - {\theta _{\textrm {con}}}$ and ${\theta _{\textrm {res}}} + {\theta _{\textrm {con}}}$. Attaining a resolution r hence requires measuring diffraction patterns out to an angular frequency ${\theta _{\textrm {obs}}}$ between these two values. For the same reason, in order for the forward model to accurately capture features of dimension r, the diffraction patterns need to be simulated out to at least ${\theta _{\textrm {cal}}} = {\theta _{\textrm {res}}} + {\theta _{\textrm {con}}}$. In the case of thin specimens, like the two-dimensional materials treated in this paper, negligible scattering is expected to exceed ${\theta _{\textrm {res}}} + {\theta _{\textrm {con}}}$, and calculations beyond this limit are not needed.

The resolution r is matched to the incoherent resolution, attained in high angle annular dark field scanning electron microscopy, by setting ${\theta _{\textrm {res}}}$ equal to $2 {\theta _{\textrm {con}}}$, and thus ${\theta _{\textrm {con}}} < {\theta _{\textrm {obs}}} \leq 3 {\theta _{\textrm {con}}}$. The lower bound implies that ptychography is able to achieve the incoherent resolution limit from just the intensities in the central disc.

Within the calculated diffraction patterns, the region surrounding the cut-off disc, i.e. the inner disc with diameter $2m/3$, is zeroed to avoid wrap-around artifacts [31]. The angular frequency corresponding to the cut-off disc’s edge is the maximum angle calculated and is denoted ${\theta _{\textrm {cal}}}$. During the reconstruction process, the calculations are fitted to observations that extend out to ${\theta _{\textrm {obs}}}$, with generally ${\theta _{\textrm {obs}}} \leq {\theta _{\textrm {cal}}}$. The width in pixels of the calculated diffraction patterns, m, is $3 {\theta _{\textrm {cal}}} / (\lambda \delta )$, and the equivalent width of the fitting region, n, is given by $2 {\theta _{\textrm {obs}}} / (\lambda \delta )$. The pixel width in reciprocal space, $\delta$, is $1/w$. The pixel width in real space is $d = \lambda / (3 {\theta _{\textrm {cal}}})$.

The number of knowns and unknowns, $N_{\textrm {k}}$ and $N_{\textrm {u}}$, is expressed as,

$$\begin{aligned} N_{\textrm{k}} &= n^{2} \left( \frac{s}{\Delta x} + 1 \right)^{2},\\ N_{\textrm{u}} &= 2 \left( \frac{ s + w }{d} \right)^{2}. \end{aligned}$$

The term of 1 for $N_{\textrm {k}}$ prevents fencepost errors, and $\Delta x$ is the scan step. The factor of 2 in the expression for $N_{\textrm {u}}$ accounts for the complex nature of the reconstruction, and s is the width of the scanned area. This yields,

(22)$$ \frac{N_{\textrm{k}}}{N_{\textrm{u}}} = \frac{2}{9} \frac{\theta_{\textrm{obs}}^{2}}{\theta_{\textrm{cal}}^{2}} \frac{w^{2}}{\Delta x^{2}} \left( \frac{ 1 + \Delta x / s }{ 1 + w/s } \right)^{2}, $$

(23)$$\overrightarrow {{s\to\infty}} \frac{2}{9} \frac{\theta_{\textrm{obs}}^{2}}{\theta_{\textrm{cal}}^{2}} \frac{ w^{2} }{ \Delta x^{2} }, $$

with the second equation describing a typical wide-field scan. For thin specimens and the incoherent resolution limit, ${\theta _{\textrm {cal}}}$ is substituted with $3 {\theta _{\textrm {con}}}$ and (23) reduces further to (3).

A.1. Guidelines

From (3) guidelines for experimental settings are deduced. Expressing that ${N_{\textrm {k}}}/ {N_{\textrm {u}}} > 1$ leads to

(24)$$w > 2.1 \Delta x, \textrm{ for } {\theta_{\textrm{obs}}} = 3 {\theta_{\textrm{con}}}, \textrm{ and} $$

(25)$$w > 6.4 \Delta x, \textrm{ for } {\theta_{\textrm{obs}}} = {\theta_{\textrm{con}}}. $$

The approximate beam width is given by the Rayleigh criterion as $0.61 \lambda / {\theta _{\textrm {con}}}$, and practical experience with both the experimental and simulated data sets in this paper, has indicated that $\Delta x$ should be considerably lower than this value, for example two-thirds or a half.

To ensure an acceptable forward simulation, w must be wide enough to encompass a sizeable fraction of the beam intensity. At $w \simeq (k+0.24) \lambda / {\theta _{\textrm {con}}}$ the support encloses the first k minima of a focused diffraction limited probe. It is hence recommended to choose w of the order $2 \lambda / {\theta _{\textrm {con}}}$ or above.

B. Total generalized variation regularization

When the object can be approximated as piecewise linear, its Laplacian is sparse and a total generalized variation (TGV) regularization can be persued. The regularization term in (8) is replaced by the $\ell _1$-norm of the approximate Laplacian, calculated as a convolution with the kernel

(26)$$\begin{pmatrix} 0 &1 &0 \\ 1 &-4 &1 \\ 0 &1 &0 \end{pmatrix}. $$

In Fig. 8, results are shown for a reconstruction with TGV-regularization with identical settings as in Section 3.2, for an electron count of 316. The RMSE in Fig. 8(a) again points to an optimum $\mu$-value, this time at 3.2, and the reconstruction at this $\mu$-value displays clearly resolved dumbbells in Fig. 8(b).

 

Fig. 8. Evaluation of the optimal regularization parameter for TGV. (a): The RMSE for the regularized reconstructions showing an optimum for $\mu = 3.2$, and (b) the reconstruction at the optimal $\mu$-value.

Download Full Size | PDF

C. Derivatives

Here, the derivative of the loss function $\ell$ with respect to the incoming beam shape and beam positions as given by Eqs. (12), (13) and (14) are derived.

For the incoming beam $\psi ^{0}$, which is a complex quantity, the derivative can be expressed for its real $\psi ^{0}_{re}$ and imaginary part $\psi ^{0}_{im}$, respectively. According to the chain rule for Wirtinger derivatives [45], we have:

(27)$$\frac{ \partial \ell }{\partial \psi^{0}_{re}} = \frac{\partial \ell}{\partial \psi^{0} t^{0}} \frac{\partial \psi^{0} t^{0}}{\partial \psi^{0}_{re}} + \frac{\partial \ell}{\partial (\psi^{0} t^{0})^{*}} \frac{\partial (\psi^{0} t^{0})^{*}}{\partial \psi^{0}_{re}},$$

where $*$ being the complex conjugate, $\partial \ell / \partial \psi ^{0} t^{0}$ obtainted from the backpropagation and $\partial \ell / \partial (\psi ^{0} t^{0})^{*} = ( \partial \ell / \partial \psi ^{0} t^{0})^{*}$. The equivalent expression for the imaginary part is omitted. Equation (12) then follows from $\partial \psi ^{0} t^{0} / \partial \psi ^{0}_{re} = t^{0}$ and $\partial (\psi ^{0} t^{0})^{*} / \partial \psi ^{0}_{re} = t^{0*}$ and Eq. (13) follows from $\partial \psi ^{0} t^{0} / \partial \psi ^{0}_{im} = it^{0}$ and $\partial (\psi ^{0} t^{0})^{*} / \partial \psi ^{0}_{im} = it^{0*}$.

For the derivative of the beam positions we keep the notation light, by only treating the one-dimensional case, the extension to two dimensions being trivial. The first element of ${\vec {r}}$ and ${\vec {R}}$ are denoted x and X, respectively. The probe positioned in the origin is denoted as $\psi ^{0}$, and in position X as $\psi ^{X}$, with $\psi ^{X}(x) = \psi ^{0}( x - X )$. Following the chain rule:

(28)$$\frac{ \partial \ell }{\partial X} = \sum_k \frac{\partial \ell}{\partial \psi^{X}_k} \frac{\partial \psi_k^{X}}{\partial X} + \frac{\partial \ell}{ \partial \psi^{X*}_k} \frac{\partial \psi_k^{X*}}{\partial X},$$

with $\partial \ell / \partial \psi _k^{X}$ provided by the backpropagation. Equation (14) is then obtained with $\partial \psi ^{X} / \partial X = -\partial \psi ^{X} / \partial x$ and $\partial \ell / \partial \psi _k^{X*} = ( \partial \ell / \partial \psi _k^{X} )^{*}$.

D. Algorithm implementation

All the reconstructions presented were generated with a ROP implementation that utilizes the parallel hardware architecture of a NVIDIA V100 GPU, allowing a simultaneous propagation through the network for a batch of probe positions. For a batch size of 100, convergence of the algorithm for the tasks described in section 3.1 ranged between 3 to $4\textrm {h}$ and for every task in section 3.2 roughly $5\textrm {min}$. Experimental results in section 3.3 were obtained with a batch size of 87 for the reconstructions of the ${\textrm {Mo}\textrm {S}_2}$ data and a batch size of 10 for the reconstructions of the ${\textrm {Nb}_3\textrm {Cl}_8}$ data, with an average time of $15\textrm {min}$ and $40\textrm {min}$, respectively.

Table 1. Settings for simulation, and reconstruction from simulated and experimental measurements of $\textrm {Mo}\textrm {S}_2$ and $\textrm {Nb}_3\textrm {Cl}_8$. The acceleration voltage is 80 kV ($\lambda = 4.18$ pm) and 120 kV ($\lambda = 3.35$ pm), respectively. To treat the non-square scan patterns, the width s of the scan area is approximated as the width of a square of equal area. $N_{\textrm {k}}/N_{\textrm {u}}$ has been calculated from (22).

View Table | View all tables in this article

Table 2. Additional ROP settings for the $\textrm {Mo}\textrm {S}_2$ experimental reconstructions, corresponding to panels in Figure 6.

View Table | View all tables in this article

Table 3. Additional ROP settings for the $\textrm {Nb}_3\textrm {Cl}_8$ experimental reconstructions, corresponding to panels in Figure 7.

View Table | View all tables in this article

Funding

Deutsche Forschungsgemeinschaft (182087777 - SFB 951, BR 5095/2-1); Defense Advanced Research Projects Agency (TEE-D18AC00009); National Science Foundation (DMR-1539918, DMR-1719875).

Acknowledgments

The authors thank I. El Baggari and L.F. Kourkoutis from Cornell University for the ${\textrm {Nb}_3\textrm {Cl}_8}$ sample preparation, and B. Haas from the Humboldt Universität zu Berlin for his help with the template matching.

Z.C. and D.A.M. are supported by PARADIM, a National Science Foundation (NSF) Materials Innovation Platform program. This work made use of the Cornell Center for Materials Research facility supported by NSF.

Experimental data acquired at Cornell is available from D.A. Muller (david.a.muller@cornell.edu) on request. Code developed at Humboldt Universität zu Berlin is available for scientific purposes from the first author (schlozma@hu-berlin.de) upon request.

Disclosures

The authors declare no conflicts of interest.

References

1. J. M. Rodenburg, “Ptychography and related diffractive imaging methods,” Adv. Imaging Electron Phys. 150, 87–184 (2008). [CrossRef]  

2. J. Rodenburg, A. Hurst, and A. Cullis, “Transmission microscopy without lenses for objects of unlimited size,” Ultramicroscopy 107(2-3), 227–231 (2007). [CrossRef]  

3. J. Rodenburg, A. Hurst, A. Cullis, B. Dobson, F. Pfeiffer, O. Bunk, C. David, K. Jefimovs, and I. Johnson, “Hard-x-ray lensless imaging of extended objects,” Phys. Rev. Lett. 98(3), 034801 (2007). [CrossRef]  

4. M. Holler, M. Guizar-Sicairos, E. H. Tsai, R. Dinapoli, E. Müller, O. Bunk, J. Raabe, and G. Aeppli, “High-resolution non-destructive three-dimensional imaging of integrated circuits,” Nature 543(7645), 402–406 (2017). [CrossRef]  

5. M. J. Cherukara, T. Zhou, Y. Nashed, P. Enfedaque, A. Hexemer, R. J. Harder, and M. V. Holt, “AI-enabled high-resolution scanning coherent diffraction imaging,” Appl. Phys. Lett. 117(4), 044103 (2020). [CrossRef]  

6. Y. Jiang, Z. Chen, Y. Han, P. Deb, H. Gao, S. Xie, P. Purohit, M. W. Tate, J. Park, S. M. Gruner, V. Elser, and D. A. Muller, “Electron ptychography of 2D materials to deep sub-ångström resolution,” Nature 559(7714), 343–349 (2018). [CrossRef]  

7. M. Guizar-Sicairos and J. R. Fienup, “Phase retrieval with transverse translation diversity: a nonlinear optimization approach,” Opt. Express 16(10), 7264–7278 (2008). [CrossRef]  

8. P. Thibault, M. Dierolf, A. Menzel, O. Bunk, C. David, and F. Pfeiffer, “High-resolution scanning x-ray diffraction microscopy,” Science 321(5887), 379–382 (2008). [CrossRef]  

9. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109(10), 1256–1262 (2009). [CrossRef]  

10. M. Kahnt, J. Becher, D. Brückner, Y. Fam, T. Sheppard, T. Weissenberger, F. Wittwer, J.-D. Grunwaldt, W. Schwieger, and C. G. Schroer, “Coupled ptychography and tomography algorithm improves reconstruction of experimental data,” Optica 6(10), 1282–1289 (2019). [CrossRef]  

11. X. Jiang, W. Van den Broek, and C. T. Koch, “Inverse dynamical photon scattering (IDPS): an artificial neural network based algorithm for three-dimensional quantitative imaging in optical microscopy,” Opt. Express 24(7), 7006–7018 (2016). [CrossRef]  

12. Z. Chen, M. Odstrcil, Y. Jiang, Y. Han, M.-H. Chiu, L.-J. Li, and D. A. Muller, “Mixed-state electron ptychography enables sub-angstrom resolution imaging with picometer precision at low dose,” Nat. Commun. 11(1), 2994 (2020). [CrossRef]  

13. A. Maiden, M. Humphry, M. Sarahan, B. Kraus, and J. Rodenburg, “An annealing algorithm to correct positioning errors in ptychography,” Ultramicroscopy 120, 64–72 (2012). [CrossRef]  

14. F. Zhang, I. Peterson, J. Vila-Comamala, A. Diaz, F. Berenguer, R. Bean, B. Chen, A. Menzel, I. K. Robinson, and J. M. Rodenburg, “Translation position determination in ptychographic coherent diffraction imaging,” Opt. Express 21(11), 13592–13606 (2013). [CrossRef]  

15. A. Tripathi, I. McNulty, and O. G. Shpyrko, “Ptychographic overlap constraint errors and the limits of their numerical recovery using conjugate gradient descent methods,” Opt. Express 22(2), 1452–1466 (2014). [CrossRef]  

16. W. Van den Broek and C. T. Koch, “Method for retrieval of the three-dimensional object potential by inversion of dynamical electron scattering,” Phys. Rev. Lett. 109(24), 245502 (2012). [CrossRef]  

17. W. Van den Broek and C. T. Koch, “General framework for quantitative three-dimensional reconstruction from arbitrary detection geometries in TEM,” Phys. Rev. B 87(18), 184108 (2013). [CrossRef]  

18. A. M. Maiden, M. J. Humphry, and J. Rodenburg, “Ptychographic transmission microscopy in three dimensions using a multi-slice approach,” J. Opt. Soc. Am. A 29(8), 1606–1614 (2012). [CrossRef]  

19. E. H. Tsai, I. Usov, A. Diaz, A. Menzel, and M. Guizar-Sicairos, “X-ray ptychography with extended depth of field,” Opt. Express 24(25), 29089–29108 (2016). [CrossRef]  

20. P. Li and A. Maiden, “Multi-slice ptychographic tomography,” Sci. Rep. 8(1), 2049 (2018). [CrossRef]  

21. M. Gilles, Y. Nashed, M. Du, C. Jacobsen, and S. Wild, “3D x-ray imaging of continuous objects beyond the depth of focus limit,” Optica 5(9), 1078–1086 (2018). [CrossRef]  

22. V. Nikitin, S. Aslan, Y. Yao, T. Biçer, S. Leyffer, R. Mokso, and D. Gürsoy, “Photon-limited ptychography of 3D objects via Bayesian reconstruction,” OSA Continuum 2(10), 2948–2968 (2019). [CrossRef]  

23. P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy 109(4), 338–343 (2009). [CrossRef]  

24. H. Yang, R. Rutte, L. Jones, M. Simson, R. Sagawa, H. Ryll, M. Huth, T. Pennycook, M. Green, H. Soltau, Y. Kondo, B. Davis, and P. Nellist, “Simultaneous atomic-resolution electron ptychography and Z-contrast imaging of light and heavy elements in complex nanostructures,” Nat. Commun. 7(1), 12532 (2016). [CrossRef]  

25. A. M. Maiden, M. J. Humphry, F. Zhang, and J. M. Rodenburg, “Superresolution imaging via ptychography,” J. Opt. Soc. Am. A 28(4), 604–612 (2011). [CrossRef]  

26. P. Thibault and M. Guizar-Sicairos, “Maximum-likelihood refinement for coherent diffractive imaging,” New J. Phys. 14(6), 063004 (2012). [CrossRef]  

27. P. Dwivedi, A. Konijnenberg, S. Pereira, and H. Urbach, “Lateral position correction in ptychography using the gradient of intensity patterns,” Ultramicroscopy 192, 29–36 (2018). [CrossRef]  

28. V. Katkovnik and J. Astola, “Sparse ptychographical coherent diffractive imaging from noisy measurements,” J. Opt. Soc. Am. A 30(3), 367–379 (2013). [CrossRef]  

29. C. T. Koch and W. Van den Broek, “Measuring three-dimensional positions of atoms to the highest accuracy with electrons,” C. R. Phys. 15(2-3), 119–125 (2014). [CrossRef]  

30. J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15(6), 1662–1669 (1998). [CrossRef]  

31. E. J. Kirkland, Advanced Computing in Electron Microscopy (Springer, 2010).

32. M. W. Tate, P. Purohit, D. Chamberlain, K. X. Nguyen, R. Hovden, C. S. Chang, P. Deb, E. Turgut, J. T. Heron, D. G. Schlom, D. C. Ralph, G. D. Fuchs, K. S. Shanks, H. T. Philipp, D. A. Muller, and S. M. Gruner, “High dynamic range pixel array detector for scanning transmission electron microscopy,” Microsc. Microanal. 22(1), 237–249 (2016). [CrossRef]  

33. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations Trends Mach. Learn. 3(1), 1–122 (2010). [CrossRef]  

34. E. Polak and G. Ribière, “Note sur la convergence de méthodes de directions conjuguées,” R.I.R.O. 3(16), 35–43 (1969). [CrossRef]  

35. B. T. Polyak, “The conjugate gradient method in extremal problems,” USSR Comput. Math. Math. Phys. 9(4), 94–112 (1969). [CrossRef]  

36. J. Nocedal and S. Wright, Numerical Optimization (Springer, 2006).

37. W. Van den Broek, X. Jiang, and C. Koch, “FDES, a GPU-based multislice algorithm with increased efficiency of the computation of the projected potential,” Ultramicroscopy 158, 89–97 (2015). [CrossRef]  

38. M. Schloz, T. C. Pekin, W. Van den Broek, and C. T. Koch, “Overcoming information reduced data and experimentally uncertain parameters in ptychography with regularized optimization,” Zenodo (2020), https://doi.org/10.5281/zenodo.3924495.

39. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory (Springer, 1992).

40. E. van den Berg and M. P. Friedlander, “Probing the Pareto frontier for basis pursuit solutions,” SIAM J. Sci. Comput. 31(2), 890–912 (2009). [CrossRef]  

41. W. Van den Broek, B. W. Reed, A. Béché, A. Velazco, J. Verbeeck, and C. T. Koch, “Various compressed sensing setups evaluated against Shannon sampling under constraint of constant illumination,” IEEE Trans. Comput. Imaging 5(3), 502–514 (2019). [CrossRef]  

42. K. Papafitsoros and C. B. Schönlieb, “A combined first and second order variational approach for image reconstruction,” J. Math. Imaging Vis. 48(2), 308–338 (2014). [CrossRef]  

43. M. Guizar-Sicairos, A. Diaz, M. Holler, M. S. Lucas, A. Menzel, R. A. Wepf, and O. Bunk, “Phase tomography from x-ray coherent diffractive imaging projections,” Opt. Express 19(22), 21345–21357 (2011). [CrossRef]  

44. J.-M. Zuo, A. B. Shah, H. Kim, Y. Meng, W. Gao, and J.-L. Rouvière, “Lattice and strain analysis of atomic resolution Z-contrast images based on template matching,” Ultramicroscopy 136, 50–60 (2014). [CrossRef]  

45. W. Wirtinger, “Zur formalen Theorie der Funktionen von mehr komplexen Veränderlichen,” Math. Annalen 97(1), 357–375 (1927). [CrossRef]