The motivation for the development of gradient-flipping and phase prediction assisted flux-preserving full-resolution wave reconstruction (GPFRWR) is that, in most TEM investigations of specimens that do not generate magnetic and electrostatic fringing fields, the phase in vacuum regions within the field of view has a constant value. The slab geometry of many TEM samples may also result in approximately flat regions of phase in large parts of the specimen, at least at medium spatial resolution. In other words, the gradient of the phase is often quite sparse, especially when excluding high spatial frequencies.

In its original application, the charge flipping algorithm for solving crystal structures from X-ray diffraction data^{30} is very effective at finding a sparse solution in the charge density domain by flipping the signs of small values while retaining values above a given threshold and enforcing consistency with measured diffraction intensities. After demonstrating its feasibility for removing low spatial frequency noise in TIE reconstructions^{23.24}we adapted the principle to non-linear inline electron holography by inserting a phase-modifying procedure every few iterations (eg, every third iteration) in an iterative reconstruction algorithm (the FRWR algorithm^{25.26}), flipping the signs of small values of each of the two components of the gradient of the phase and reducing their amplitudes, yielding a modified gradient ({overrightarrow{G}}^{^{prime}}left(overrightarrow{r}right)). This procedure was implemented simply by multiplying these values by a scaling factor β that was slightly larger than − 1 (eg, β = − 0.97). This operation was only performed within the field of view defined by the experimental data. The size of the array defining the reconstructed phase was larger than the field of view, in order to be able to accommodate non-periodic boundary conditions. Whereas the larger array has periodic boundary conditions, the array corresponding to the field of view of the experimental data can have any boundary because it lies within the typically 1.5 to 2 times wider reconstruction array that has periodic boundary conditions^{16}. Once the small gradients have been flipped, the modified phase is obtained by integrating the modified gradient ({overrightarrow{G}}^{^{prime}}left(overrightarrow{r}right)) The Fourier transform of the modified phase ({phi }^{^{prime}}left(overrightarrow{q}right)) is obtained from ({overrightarrow{G}}^{^{prime}}left(overrightarrow{r}right)) by the following operation:

$${phi }^{^{prime}}left(overrightarrow{q}right)=phi left(overrightarrow{q}right)left[1-mathit{exp}left(-{r}_{c}^{2}{q}^{2}right)right]+FTright[overrightarrow{nabla }cdot {overrightarrow{G}}^{^{prime}}left(overrightarrow{r}right)right]frac{mathit{exp}left(-{r}_{c}^{2}{q}^{2}right)}{-{q}^{2}},$$

(3)

where *r*_{c} is the length scale below which the flipping of the gradient has very little effect, ie this value can be chosen to apply gradient flipping only to spatial frequencies for which the iterative reconstruction algorithm would otherwise converge only very slowly. Dividing by – q^{two} effectively implements an inverse Laplace operator. At *what*= 0, we multiply by 0 instead of dividing by q^{two}. This approach can be justified by the argument that the absolute phase is not a well-defined physical quantity. Multiplication by 0 at *what*= 0 will cause the mean of the reconstructed phase to be set to 0. After reconstruction, an offset that corresponds to the mean of the phase-in vacuum can be subtracted. This subtraction was applied to the phase maps shown in this paper.

Using the above expression, gradient flipping affects mostly spatial frequencies in the phase that are larger than *r*_{c} by using a Gaussian taper to transition between the ranges *r* > *r*_{c} and*r*< *r*_{c}. Since the iterative part of the focal series reconstruction algorithm reconstructs primarily high spatial frequencies in the phase and requires many iterations to have an effect on lower spatial frequencies^{16}, Eq. (3) ensures that gradient flipping does not affect the convergence of the iterative reconstruction algorithm significantly. It may even help to speed up convergence, especially in regions where large areas of the phase are flat (*eg*for nanoparticles on homogeneous supports or if the field of view contains large areas of empty space).

In order to reduce the detrimental effect of incoherent scattering (eg, electron–phonon and electron-plasmon interactions or the excitation of core electrons) on the interpretability of experimental TEM images using elastic scattering theory, our algorithm is based on the following flux-preserving imaging model, which is modified to include an inconsistent background:

$${I}_{Delta f}left(overrightarrow{r}right)=left[{left|{Psi }_{Delta f}left(overrightarrow{r}right)right|}^{2}+{I}_{incoherent}left(overrightarrow{r}right)right]otimes {E}_{s,Delta f}left(overrightarrow{r}right),$$

(4)

where*Ψ*_{Δf} ( *r*) is the coherent electron wave function defocused by an amount Δ*F*in the plane of the detector, *AND*_{s, Δf}(*r* ) is the inverse Fourier transform of the spatial coherence envelope^{fifteen} and *Yo*_{inconsistent} (*r* ) is the intensity distribution of the incoherent background that is determined in addition to the electron wave function *Ψ*_{0} ( *r*) at the exit surface of the specimen. Since only a single inconsistent image is assumed, this residual image is taken to be the same for all images in the focal series. The inconsistent intensity is strictly positive. At each iteration, after the average of the remaining difference between the forward simulated image and the experimental image has been added to the incoherent intensity array, this array is multiplied by a damping factor, which is set to 0.98 by default but can be changed when calling the reconstruction program. Thus, only a fraction of the minimum difference between simulated and experimental images in each pixel is assigned to this inconsistent intensity array, ensuring that only image counts observed that cannot be attributed to the coherent imaging process are considered for assignment to the inconsistent background.

The FRWR algorithm^{25.26} accelerates the retrieval of low spatial frequency phase information by using a ‘phase prediction’ mechanism. In addition to directly minimizing the difference between simulated and experimental images in a manner that is somewhat similar to conventional Gerchberg-Saxton-type focal series reconstruction algorithms.^{25}, the FRWR algorithm also explicitly updates the phase. This is done by first computing a phase update that is motivated by the transport of intensity equation (TIE) as follows:

$$begin{array}{c}Delta varphi left(overrightarrow{r}right)={varphi }^{exp}left(overrightarrow{r}right)-{varphi } sim}left(overrightarrow{r}right)\ =frac{2pi }{lambda }frac{{nabla }^{-1}{I}_{0}^{ -1}left(overrightarrow{r}right){nabla }^{-1}left({I}_{+Delta f}^{mathrm{exp}}left(overrightarrow{ r}right)-{I}_{-Delta f}^{mathrm{exp}}left(overrightarrow{r}right)right)}{2Delta f}-frac{2 pi }{lambda }frac{{nabla }^{-1}{I}_{0}^{-1}left(overrightarrow{r}right){nabla }^{-1 }left({I}_{+Delta f}^{sim}left(overrightarrow{r}right)-{I}_{-Delta f}^{sim}left(overrightarrow{ r}right)right)}{2Delta f}\ =frac{2pi }{lambda }frac{{nabla }^{-1}{I}_{0}^{ -1}left(overrightarrow{r}right){nabla }^{-1}left({I}_{+Delta f}^{mathrm{exp}}left(overrightarrow{ r}right)-{I}_{+Delta f}^{sim}left(overrightarrow{r}right)right)}{2Delta f}-frac{2pi }{ lambda }frac{{nabla }^{-1}{I}_{0}^{-1}left(overrightarrow{r}right){nabla }^{-1}left( {I}_{-Delt a f}^{mathrm{exp}}left(overrightarrow{r}right)-{I}_{-Delta f}^{sim}left(overrightarrow{r}right)right) }{2Delta f}end{array}$$

(5)

*ie*a TIE-like update of the phase that takes into account all*N*plans of focus can be computed using the expression

$$Delta varphi (overrightarrow{r})=frac{2pi }{Nlambda }{sum }_{n=1}^{N}frac{{nabla }^{- 1}{I}_{0}^{-1}(overrightarrow{r}){nabla }^{-1}left({I}_{Delta {f}_{n}}^{ exp}(overrightarrow{r})-{I}_{Delta {f}_{n}}^{sim}(overrightarrow{r})right)}{Delta {f}_{n} }.$$

(6)

Before adding this phase update to the phase of the current estimate of the retrieved wave function, its magnitude is limited to a maximum value*ϕ*_{max}the so-called phase prediction threshold (PPT), by using exponential compression for high values of Δ*ϕ*(*r*) to avoid a sharp threshold:

$$varphi left(overrightarrow{r}right)= left{begin{array}{ll}signleft(varphi left(overrightarrow{r}right)right)frac {{varphi }_{max}}{2}left[2-expleft(frac{frac{{varphi }_{max}}{2}-left|varphi left(overrightarrow{r}right)right|}{frac{{varphi }_{max}}{2}}right)right]& if left|varphi left(overrightarrow{r}right)right|ge frac{{varphi }_{max}}{2}\ varphi left(overrightarrow{r} right)& if left|varphi left(overrightarrow{r}right)right|< frac{{varphi }_{max}}{2}end{array}right.$$

The use of this limit is important because inconsistent contributions to the signal and noise prevent a perfect match from being obtained between experimental and simulated intensities, causing the expression in parentheses in Eq. (6) to practically never reach zero. The continual addition of a non-zero update to the phase would therefore result in diverging behavior of the reconstruction algorithm.