© EDP Sciences, 2025

1 Introduction

Within recent years, scanning transmission electron microscopy (STEM) has evolved into an attractive tool for the investigation of beam-sensitive objects such as viruses [1,2], 2D materials [3–5], zeolites [6–10], Lirich oxides [11–13], polymers [14], perovskites [15–17] or metal-organic frameworks (MOF) [18,19]. When imaging such fragile materials, the damage following the transfer of energy from interacting electrons [20], such as e.g. knock-on displacement of atoms [21,22], heating [23] or radiolysis [24], imposes a critical electron dose [25] beyond which the specimen structure is lost. This critical dose then constitutes the main experimental limitation, thus in practice determining the best achievable resolution [26,27].

More generally, the prevalence of beam damage requires both a re-evaluation of the maximum electron dose to be invested and an improvement in detector quantum efficiency (DQE). The latter was fulfilled by the introduction of direct electron detectors (DED) [28–38] surpassing the capacities of conventional scintillator cameras [39], including in terms of their modulation transfer function (MTF) [40–42]. The gain in recording speed, allowed by faster electronics, furthermore enabled the acquisition of a convergent beam electron diffraction (CBED) pattern at each scan position [43], a technique often referred to as momentum-resolved STEM (MR-STEM) [44] or 4D-STEM [45]. More recently, event-driven detection [46], based on the Timepix [29,32,36,47] technology, permitted the extension of this technique to sub-microsecond single-pattern acquisition times [48–50], thus reaching the same speed as conventional STEM.

The subsequent knowledge on the far-field intensity distribution furthermore enables the use of a class of computational imaging methods known as ptychography [51–53] for the measurement of the projected electrostatic potential of the illuminated object, in the form of a phase shift map. Those methods consist in the correlative use of a series of coherent scattering experiments, in which a redundancy of imprinted specimen information permits the retrieval of a common illuminated object. They can be thus be considered as an extension of the well-established coherent diffractive imaging (CDI) technique [54–59] to the situation where multiple independent recordings are employed and where, at least in the basic case, no prior information is available.

Among ptychographic methods based on the focused-probe geometry, iterative phase retrieval [60] approaches have recently met some success e.g. with biological specimens [1,61–64]. Those approaches consist in the probe position-dependent simulation of the elastic scattering of the incident electrons, thus leading to a repeated update of the multiplicative transmission function $T\left({{\displaystyle \overrightarrow{r}}}_0\right)$ used to represent the specimen, given the error made against the experimental recordings and while cycling through the corresponding scan positions. This update is performed a number of times for each complete cycle, depending on the chosen batch size, and usually follows one of the several variants [65–72] of the Gerchberg-Saxton algorithm [73–74] for sequential projections or is given by the gradient of a specific loss function [75–78], i.e. the maximum likelihood approach. The process may also include a regularization term [71,79–83] or be based on a parameterization strategy [83–86].

Due to the wide range of parameters available, encompassing e.g. the coupled loss and regularization functions, the update strength, the batch size, the initial guess on the reconstructed object or the possible use of a supplementary momentum term [71], iterative methods possess a high degree of flexibility. On the other hand, while a specific choice of parameter set may permit a degree of adaptation to particular cases, for instance with regards to the noise model [78,87–89], this also implies the need for a complex tuning step to achieve numerical convergence [90]. Different results may otherwise be obtained through separate reconstruction processes or algorithms, hence creating reproducibility issues. Achieving convergence may furthermore prove more challenging in the low-dose case [76,91], where the exploited far-field patterns are underdetermined, independently of the dose-efficiency demonstrated by the converged result in itself. Finally, iterative ptychography remains numerically intensive and often requires advanced computation capacities to avoid exceedingly long processing times [92–96].

For those reasons, there is still a high interest in pursuing work on analytical solutions [97–99] which, since they lead to method-unique results through direct and well-understood imaging processes, arguably constitute more reproducible measurement approaches. In particular, as they are also fast, their application in a wide range of conditions or for large collections of specimens can be streamlined, hence making those methods especially useful for challenging experimental cases. In-line treatment options permitting live imaging [100,101] have been reported as well, while this remains challenging in the framework of an iterative process [102]. Analytical ptychography has moreover demonstrated a high dose-efficiency [103–106], including with a sensitivity to light elements [107–110], and was successfully applied to a variety of beam-sensitive objects [9–14,16] in recent years.

In this publication, the fundamental capacities of analytical ptychography methods to image a beam-sensitive specimen are explored in conditions of very low electron dose. In particular, interest is taken in the resolution achieved for different numerical apertures, in the dose-dependent precision of the measurements and in the fundamental frequency transfer capacity of different approaches. The methods investigated are the Wigner distribution deconvolution (WDD) [98,111–113], integrated center of mass (iCoM) [114,115] imaging as well as the sideband integration (SBI) [99,103,104,116] method, sometimes referred to as single sideband (SSB) [106]. Benchmarking is continued with an overfocused probe [117,118] and an adapted process permitting the direct correction of known aberrations. This specific recording geometry has recently attracted interests for the imaging of beam-sensitive objects and bears similarity with the original idea of ref. [119]. After an initial review of the theory in the fully coherent case, practical implementation is demonstrated through the newly introduced scan-frequency partitioning algorithm (SFPA), permitting a straightforward parallelization and offering high flexibility in the size and pixel resolution of the reconstruction window. All demonstrations made here are based on MR-STEM simulations, hence allowing direct control over the illumination parameters and the dose, while ensuring sparsity in the electron counts. Two model objects are employed: monolayer MoS₂ and ice-embedded apoferritin [110].

2 Theory of analytical ptychography and new highly parallelizable implementation

2.1 Coherent and elastic interaction under the phase object approximation

2.1.1 Transmission function

In its conventional form [97,98], analytical ptychography makes use of the phase object approximation (POA) [120]. In this context, the imaged material is considered thin enough so that no variation of wave amplitude occurs within it, thus making the scattering-induced phase shift additive along the propagation axis. For thicker objects, the applicability of the POA is limited due to the role of near-field propagation, leading to a finite depth of focus [121,122] and dynamical diffraction effects such as channeling [123].

Formally, the elastic interaction of the electron probe $P\left({{\displaystyle \overrightarrow{r}}}_0\right)$ with the specimen is then modeled by a multiplication with a transmission function $T\left({{\displaystyle \overrightarrow{r}}}_0\right)$, defined for each real-space position ${{\displaystyle \overrightarrow{r}}}_0$ in the specimen plane by

$$T\left({{\displaystyle \overrightarrow{r}}}_0\right)=e^{i\sigma \mu \left({{\displaystyle \overrightarrow{r}}}_0\right)},$$(1)

where $\mu \left({{\displaystyle \overrightarrow{r}}}_0\right)$ is the projected electrostatic potential of the specimen, i.e. the integral of the three-dimensional potential along the propagation axis. The interaction parameter σ, expressed in V⁻¹ ⋅ m⁻¹, is given by

$$\sigma =\frac{2\pi e}{hc}\frac{m{c}^2+eU}{\sqrt{eU\left(2m{c}^2+eU\right)}},$$(2)

with e is the elementary charge, m the electron rest mass, h the Planck constant and c the speed of light. The product of σ with $\mu \left({{\displaystyle \overrightarrow{r}}}_0\right)$ thus represents the local phase shift imposed to the electron wavefunction by the specimen and is typically given in radians.

Importantly, due to the dependence of σ on U, the acceleration voltage affects this phase shift in a non-linear manner, independently of the specimen itself. An empirical absorption term may also be added to $\mu \left({{\displaystyle \overrightarrow{r}}}_0\right)$, as an imaginary number, to improve the agreement with experimental results [120,124], e.g. by accounting for amplitude variations in a computationally retrieved transmission function. Typically, this term is related to inelastic scattering [125,126], which otherwise leads to a diffuse component in the far-field [127–130], and specimen vibrations [131–133]. In this work, it is left out for simplicity, i.e. the interaction is assumed to not affect the coherence of the electron beam.

2.1.2 Convergent illumination

Continuing, given a fully coherent illumination, the electron probe $P\left({{\displaystyle \overrightarrow{r}}}_0\right)$ is found equal to

$$P\left({{\displaystyle \overrightarrow{r}}}_0\right)={ℱ}^{-1}\left[A\left({{\displaystyle \overrightarrow{q}}}_0\right)e^{-i\chi \left({{\displaystyle \overrightarrow{q}}}_0\right)}\right]\left({{\displaystyle \overrightarrow{r}}}_0\right),$$(3)

with $\chi \left({{\displaystyle \overrightarrow{q}}}_0\right)$ the geometrical aberration function and $A\left({{\displaystyle \overrightarrow{q}}}_0\right)$ representing the aperture in the focal plane of the probe-forming lens, being equal to 1 for $\left|\right|{{\displaystyle \overrightarrow{q}}}_0\left|\right|<q_A$ and 0 otherwise. The quantity q_A = sin(α)/λ introduced here is the reciprocal space cut-off imposed by the aperture, with α the semi-convergence angle and λ the relativistically corrected wavelength [134]. Noteworthily, the aperture function actually used in the numerical implementation is further normalized, as being representative of the wavefunction in the focal plane.

The notations ℱ and ℱ⁻¹ respectively refer to a Fourier transform and an inverse Fourier transform, given by

$$\begin{array}{ll}{\displaystyle \tilde{\phi }}\left({\displaystyle \overrightarrow{v}}\right)=ℱ\left[\phi \left({\displaystyle \overrightarrow{u}}\right)\right]\left({\displaystyle \overrightarrow{v}}\right)\hfill & ={\displaystyle {\displaystyle \sum }_{{\displaystyle \overrightarrow{u}}}}{e}^{-i2\pi {\displaystyle \overrightarrow{v}}\cdot {\displaystyle \overrightarrow{u}}}\phi \left({\displaystyle \overrightarrow{u}}\right)\hfill \\ \phi \left({\displaystyle \overrightarrow{u}}\right)={ℱ}^{-1}\left[{\displaystyle \tilde{\phi }}\left({\displaystyle \overrightarrow{v}}\right)\right]\left({\displaystyle \overrightarrow{u}}\right)\hfill & ={\displaystyle {\displaystyle \sum }_{{\displaystyle \overrightarrow{v}}}}{e}^{i2\pi {\displaystyle \overrightarrow{v}}\cdot {\displaystyle \overrightarrow{u}}}{\displaystyle \tilde{\phi }}\left({\displaystyle \overrightarrow{v}}\right),\hfill \end{array}$$(4)

with Fourier normalization left implicit.

2.1.3 Scattered and measurable intensity

At a given scan position ${{\displaystyle \overrightarrow{r}}}_s$, a localized exit wave ${\Psi }_{{{\displaystyle \overrightarrow{r}}}_s}\left({{\displaystyle \overrightarrow{r}}}_0\right)$ is formed which is given by

$${\Psi }_{{{\displaystyle \overrightarrow{r}}}_s}\left({{\displaystyle \overrightarrow{r}}}_0\right)=P\left({{\displaystyle \overrightarrow{r}}}_0-{{\displaystyle \overrightarrow{r}}}_s\right)T\left({{\displaystyle \overrightarrow{r}}}_0\right),$$(5)

hence the diffracted intensity $I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)$ accessible in the far-field, with ${\displaystyle \overrightarrow{q}}$ a scattering vector, is

$$I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)=|ℱ\left[{\Psi }_{{{\displaystyle \overrightarrow{r}}}_s}\left({{\displaystyle \overrightarrow{r}}}_0\right)\right]\left({\displaystyle \overrightarrow{q}}\right){|}^2.$$(6)

As such, this quantity is interpretable as a probability distribution among the locations, at the concerned optical plane ${\displaystyle \overrightarrow{q}}$, for the collapse of the electron wave.

Finally, the intensity $I_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ measured by the camera, across detector space ${{\displaystyle \overrightarrow{q}}}_d$, includes a possible point spread effect represented by the MTF $M\left({{\displaystyle \overrightarrow{r}}}_d\right)$, ${{\displaystyle \overrightarrow{r}}}_d$ being the reciprocal dimension of ${{\displaystyle \overrightarrow{q}}}_d$. This leads to

$$\begin{array}{ll}{I}_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill & =ℱ\left[M\left({{\displaystyle \overrightarrow{r}}}_d\right){ℱ}^{-1}\left[I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)\right]\left({{\displaystyle \overrightarrow{r}}}_d\right)\right]\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill \\ \hfill & ={\displaystyle \tilde{M}}\left({{\displaystyle \overrightarrow{q}}}_d\right){\otimes }_{{{\displaystyle \overrightarrow{q}}}_d}{I}_{{{\displaystyle \overrightarrow{r}}}_s}\left({{\displaystyle \overrightarrow{q}}}_d\right).\hfill \end{array}$$(7)

Given that ${\displaystyle \tilde{M}}\left({{\displaystyle \overrightarrow{q}}}_d\right)=ℱ\left[M\left({{\displaystyle \overrightarrow{r}}}_d\right)\right]\left({{\displaystyle \overrightarrow{q}}}_d\right)$ is a real quantity, if the MTF-induced information spread effect is isotropic, then ${\displaystyle \tilde{M}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ and $M\left({{\displaystyle \overrightarrow{r}}}_d\right)$ are both real and point-symmetric. This assumption is implicit in the rest of this work.

2.2 Wigner distribution formalism

2.2.1 Scattering data reformulation

Given the prior recording of $I_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ through an MR-STEM experiment, a Fourier transform with respect to the scan coordinates ${{\displaystyle \overrightarrow{r}}}_s$ towards an arbitrarily sampled spatial frequency space ${\displaystyle \overrightarrow{Q}}$ leads to a complex distribution ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$. As long as real-space is sampled finely enough by the scan points ${{\displaystyle \overrightarrow{r}}}_s$, i.e. under the condition of sufficient overlap ratio ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ between successively illuminated areas [135], ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ can then be interpreted as a map of the specimen-dependent ${\displaystyle \overrightarrow{Q}}$-responses attributed to the camera pixels ${{\displaystyle \overrightarrow{q}}}_d$. In particular, each scattering vector ${\displaystyle \overrightarrow{q}}$ in the far-field is assimilated to a single conventional TEM image by arguments of reciprocity [136,137]. More details on the redundancy condition and area overlap are provided in Appendix A, including with a mathematical criterion.

The distribution ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ is found equal to

$$\begin{array}{ll}{{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill & ={\displaystyle \tilde{M}}\left({{\displaystyle \overrightarrow{q}}}_d\right){\otimes }_{{{\displaystyle \overrightarrow{q}}}_d}\left({\displaystyle \tilde{P}}\left({{\displaystyle \overrightarrow{q}}}_d\right){{\displaystyle \tilde{P}}}^{*}\left({{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right)\right)\hfill \\ \hfill & {\otimes }_{{{\displaystyle \overrightarrow{q}}}_d}\left({\displaystyle \tilde{T}}\left({{\displaystyle \overrightarrow{q}}}_d\right){{\displaystyle \tilde{T}}}^{*}\left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}}\right)\right)\hfill \\ \hfill & =ℱ\left[M\left({{\displaystyle \overrightarrow{r}}}_d\right)\Gamma \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right){\rm Y}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)\right]\left({{\displaystyle \overrightarrow{q}}}_d\right),\hfill \end{array}$$(8)

where $\Gamma \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)$ and ${\rm Y}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)$ are Wigner distributions [138], i.e. autocorrelations of the probe and of the transmission function.

Formally, they are given by

$$\begin{array}{ll}{\rm Y}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)\hfill & =ℱ\left[T\left({{\displaystyle \overrightarrow{r}}}_0+{{\displaystyle \overrightarrow{r}}}_d\right)T^{*}\left({{\displaystyle \overrightarrow{r}}}_0\right)\right]\left({\displaystyle \overrightarrow{Q}}\right)\hfill \\ \hfill & ={ℱ}^{-1}\left[{\displaystyle \tilde{T}}\left({{\displaystyle \overrightarrow{q}}}_0\right){{\displaystyle \tilde{T}}}^{*}\left({{\displaystyle \overrightarrow{q}}}_0-{\displaystyle \overrightarrow{Q}}\right)\right]\left({{\displaystyle \overrightarrow{r}}}_d\right),\hfill \end{array}$$(9)

and

$$\Gamma \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)={ℱ}^{-1}\left[A\left({{\displaystyle \overrightarrow{q}}}_0\right)A\left({{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)\theta \left({{\displaystyle \overrightarrow{q}}}_0;{{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)\right]\left({{\displaystyle \overrightarrow{r}}}_d\right).$$(10)

The function $\theta \left({{\displaystyle \overrightarrow{q}}}_0;{{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)$ encompasses the imperfections in the illumination and is equal to

$$\theta \left({{\displaystyle \overrightarrow{q}}}_0\text{'};{{\displaystyle \overrightarrow{q}}}_0\right)=e^{-i\left(\chi \left({{\displaystyle \overrightarrow{q}}}_0\text{'}\right)-\chi \left({{\displaystyle \overrightarrow{q}}}_0\right)\right)}.$$(11)

As such, the insertion of this term in a ptychographic processing allows correcting for geometrical aberrations.

A subsequent inverse Fourier transform from the camera dimensions ${{\displaystyle \overrightarrow{q}}}_d$ to an arbitrary set of reciprocal real-space coordinates ${\displaystyle \overrightarrow{R}}$ leads to $J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)$, a new complex four-dimensional distribution, equal to the product of the Wigner distributions with the MTF, i.e.

$$J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)=M\left({\displaystyle \overrightarrow{R}}\right)\Gamma \left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right){\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right).$$(12)

2.2.2 Wigner distribution deconvolution and direct extraction of the transmission function

The WDD method for analytical ptychography [98–113] thus first consists in the recovery of ${\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ through

$${\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)\equiv \frac{M\left({\displaystyle \overrightarrow{R}}\right){\Gamma }^{*}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)}{ϵ+|M\left({\displaystyle \overrightarrow{R}}\right)\Gamma \left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right){|}^2},$$(13)

where ϵ is a small number introduced to avoid divisions by zero, i.e. the actual deconvolution is done via a Wiener filter process [139]. Noteworthily, a careful choice of the Wiener parameter ϵ is also important to avoid an amplification of the noise propagated from detector space [106].

As a second step in the processing, the summation of ${\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ along ${\displaystyle \overrightarrow{R}}$ permits the recovery of the transmission function by

$$\begin{array}{ll}f\left({\displaystyle \overrightarrow{Q}}\right)\hfill & ={\displaystyle {\displaystyle \sum }_{{\displaystyle \overrightarrow{R}}}}{\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)={\displaystyle \tilde{T}}\left({\displaystyle \overrightarrow{0}}\right){{\displaystyle \tilde{T}}}^{*}\left(-{\displaystyle \overrightarrow{Q}}\right)\hfill \\ T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)\hfill & ={\left({ℱ}^{-1}\left[\frac{f\left({\displaystyle \overrightarrow{Q}}\right)}{\sqrt{f\left({\displaystyle \overrightarrow{0}}\right)}}\right]\left({\displaystyle \overrightarrow{r}}\right)\right)}^{*},\hfill \end{array}$$(14)

where $f\left({\displaystyle \overrightarrow{Q}}\right)$ is introduced as an intermediary result. With the deconvolution done, $\sigma {\mu }^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$, i.e. the measurement of the phase shift map $\sigma \mu \left({\displaystyle \overrightarrow{r}}\right)$ by WDD ptychography, can be obtained from $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$ by extracting its angle. In most practical cases, including those presented here, the values obtained remain small enough to avoid phase discontinuities, thus making an unwrapping process unnecessary.

2.2.3 Zero-frequency component

As illustrated by equation (14), the measurement is performed such that $\text{arg}\left[f\left({\displaystyle \overrightarrow{0}}\right)\right]=0$. Consequently, $f\left({\displaystyle \overrightarrow{0}}\right)$ is a real number, which implies that the mean of the ptychographic phase in the reconstruction window, i.e. its DC component, remains inaccessible. This is consistent with the fact that phase, as a mathematical abstraction rather than a significant physical quantity, remains unmeasurable unless compared to a reference, e.g. by wave interference like in the case of off-axis electron holography [140,141].

In particular, in the coherent and elastic interaction regime, only relative local phase shifts created within the illuminated patch, thus requiring gradients in the specimen-induced phase shift map, can lead to measurable changes in the momentum distribution [142]. As a side-note, the normalization of $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$ by the real constant $\sqrt{f\left({\displaystyle \overrightarrow{0}}\right)}$ only modifies its amplitude, and not its angle. It thus does not affect the measurement of the projected potential itself.

2.2.4 Resolution limit

Continuing, since $\Gamma \left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)=0$ for $\left|\right|{\displaystyle \overrightarrow{Q}}\left|\right|\ge 2q_A$, owing to formula (10), it appears at first sight as though the best resolution achievable by the WDD approach is equal to half the conventional Abbe criterion δr_Abbe = 0.5/q_A. In reality, super-resolution, i.e. the transfer of frequencies exceeding the 2q_A diffraction limit [74,98,143–146], is still possible based on the so-called stepping out method [98,113].

This is nevertheless done at a high cost in dose [147], as it is specifically dependent on the availability of dark field electrons. Since this publication focuses particularly on the dose-efficiency of analytical ptychography, this aspect is left for future work. The WDD method can otherwise make use of the dark field electrons outside of the stepping out paradigm, as apparent in the equations. Moreover, and as long as the interaction can still be faithfully described using the POA, geometrical aberrations are corrected through the introduction of the term $\theta \left({{\displaystyle \overrightarrow{q}}}_0;{{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)$, defined in equation (11).

2.3 Sideband formalism for a weak scatterer

2.3.1 Simplification of Wigner distribution formalism via a first order Taylor expansion

The special case where the specimen can be considered weakly scattering, in addition to fulfilling the POA, occurs when the range of phase shift covered by the measurable $\sigma \mu \left({{\displaystyle \overrightarrow{r}}}_0\right)$, i.e. accounting for the reduction due to resolution limit, is significantly smaller than 1. Equivalently to the well-known small-angle approximation, the transmission function may be then replaced by the following first order Taylor expansion

$$T\left({{\displaystyle \overrightarrow{r}}}_0\right)\approx 1+i\sigma \mu \left({{\displaystyle \overrightarrow{r}}}_0\right).$$(15)

This condition is usually referred to as the weak phase object approximation (WPOA) [120].

It then follows that

$$\begin{array}{ll}T\left({{\displaystyle \overrightarrow{r}}}_0+{{\displaystyle \overrightarrow{r}}}_d\right)T^{*}\left({{\displaystyle \overrightarrow{r}}}_0\right)\approx 1\hfill & +i\sigma \mu \left({{\displaystyle \overrightarrow{r}}}_0+{{\displaystyle \overrightarrow{r}}}_d\right)\hfill \\ \hfill & -i\sigma \mu \left({{\displaystyle \overrightarrow{r}}}_0\right)\hfill \\ \hfill & +{\sigma }^2\mu \left({{\displaystyle \overrightarrow{r}}}_0+{{\displaystyle \overrightarrow{r}}}_d\right)\mu \left({{\displaystyle \overrightarrow{r}}}_0\right).\hfill \end{array}$$(16)

As the condition of a weakly scattering object also implies that ${\sigma }^2\mu \left({{\displaystyle \overrightarrow{r}}}_0+{{\displaystyle \overrightarrow{r}}}_d\right)\mu \left({{\displaystyle \overrightarrow{r}}}_0\right)<<1$, equation (8) leads to

$${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)\approx B\left({{\displaystyle \overrightarrow{q}}}_d\right)\delta \left({\displaystyle \overrightarrow{Q}}\right)+i\sigma \omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right){\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right).$$(17)

Notably, this is also justified by the realness of $\mu \left({{\displaystyle \overrightarrow{r}}}_0\right)$, which implies that ${{\displaystyle \tilde{\mu }}}^{*}\left(-{\displaystyle \overrightarrow{Q}}\right)={\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right)$.

The functions $B\left({{\displaystyle \overrightarrow{q}}}_d\right)$ and $\omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$ are given by

$$\begin{array}{ll}B\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill & =ℱ\left[M\left({{\displaystyle \overrightarrow{r}}}_d\right){ℱ}^{-1}\left[A^2\left({{\displaystyle \overrightarrow{q}}}_0\right)\right]\left({{\displaystyle \overrightarrow{r}}}_d\right)\right]\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill \\ \omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)\hfill & =ℱ\left[M\left({{\displaystyle \overrightarrow{r}}}_d\right)\Gamma \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)\left(e^{i2\pi {\displaystyle \overrightarrow{Q}}\cdot {{\displaystyle \overrightarrow{r}}}_d}-1\right)\right]\left({{\displaystyle \overrightarrow{q}}}_d\right).\hfill \end{array}$$(18)

${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ thus consists of a zero-frequency term, associated to the unscattered portion of the electron beam, and a two-sideband term resulting from the Fourier transform of the distribution $\Gamma \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{r}}}_d\right)$.

2.3.2 Deconvolutive extraction of the potential

Equation (17) serves as a basis for the sideband method of analytical ptychography [99,103,104,116], in this work referred to as SBI. It can thus be understood as a special case of the Wigner distribution approach, applicable when the object fulfills the WPOA. Here, a deconvolutive form is used, similarly to e.g. reference [116]. For clarity, it will be referred to as SBI-D in the rest of this text, while the conventional summative form [103,104] will be referred to as SBI-S.

The SBI-D process thus consists in performing

$$\begin{array}{ll}g\left({\displaystyle \overrightarrow{Q}}\right)\hfill & =\frac{1}{i\Omega }{\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_d}}\frac{{\omega }^{*}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right){{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)}{ϵ+|\omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right){|}^2}\hfill \\ {\mu }^{SBI}\left({\displaystyle \overrightarrow{r}}\right)\hfill & =\frac{1}{\sigma }{ℱ}^{-1}\left[g\left({\displaystyle \overrightarrow{Q}}\ne {\displaystyle \overrightarrow{0}}\right)\right]\left({\displaystyle \overrightarrow{r}}\right).\hfill \end{array}$$(19)

$g\left({\displaystyle \overrightarrow{Q}}\right)$ is an intermediary result and $\Omega ={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_d}}{A}^2\left({{\displaystyle \overrightarrow{q}}}_d\right)$ is introduced for normalization, i.e. the projected potential is obtained by calculating a mean among the scattering coordinates $\parallel {{\displaystyle \overrightarrow{q}}}_d\parallel <q_A$, post-division by $\omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$.

Importantly, like for the WDD case, the DC component is not recoverable, since $\omega \left({\displaystyle \overrightarrow{0}};{{\displaystyle \overrightarrow{q}}}_d\right)=0$, as is shown by formula (18). The inclusion of $\theta \left({{\displaystyle \overrightarrow{q}}}_0;{{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)$ also permits the correction of aberrations, at least as long as the WPOA is fulfilled.

2.3.3 Summative extraction of the potential

If the influence of the MTF is neglected, i.e. $M\left({{\displaystyle \overrightarrow{r}}}_d\right)=1$, then equation (17) becomes

$$\begin{array}{ll}{{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)\approx \hfill & A^2\left({{\displaystyle \overrightarrow{q}}}_d\right)\delta \left({\displaystyle \overrightarrow{Q}}\right)\hfill \\ \hfill & +i\sigma A\left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}}\right)A\left({{\displaystyle \overrightarrow{q}}}_d\right)\theta \left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{Q}}\right){\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right)\hfill \\ \hfill & -i\sigma A\left({{\displaystyle \overrightarrow{q}}}_d\right)A\left({{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right)\theta \left({{\displaystyle \overrightarrow{q}}}_d;{{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right){\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right).\hfill \end{array}$$(20)

Hence, ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ can be described as a superposition of two sidebands terms with a zero-frequency component.

In practice, this sideband-like geometry means that, upon visualizing the values across the ${{\displaystyle \overrightarrow{q}}}_d$-dimensions, for a given specimen frequency ${\displaystyle \overrightarrow{Q}}$ and as long as $\chi \left({{\displaystyle \overrightarrow{q}}}_0\right)=0$, the double overlap area will be homogeneously equal to $\pm \sigma {\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right)$. This constitutes the basis of the conventional SSB workflow [103,104] and provides an opportunity for straightforward aberration correction [107,148].

In this context, the SBI-S process consists in performing a summation within the double overlaps, while excluding triple overlaps where the terms cancel out. Formally, it consists in

$$\begin{array}{ll}\hfill & g\text{'}\left({\displaystyle \overrightarrow{Q}}\right)=\frac{1}{i}{\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_d}}\left({\beta }^{+}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)-{\beta }^{-}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)\right){{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill \\ \hfill & \zeta \left({\displaystyle \overrightarrow{r}}\right){\otimes }_{{\displaystyle \overrightarrow{r}}}{\mu }^{SBI}\left({\displaystyle \overrightarrow{r}}\right)=\frac{1}{\sigma }{ℱ}^{-1}\left[g\text{'}\left({\displaystyle \overrightarrow{Q}}\right)\right]\left({\displaystyle \overrightarrow{r}}\right),\hfill \end{array}$$(21)

where the ${\beta }^{\pm }\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$ terms are given by

$$\begin{array}{ll}{\beta }^{+}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)\hfill & =\frac{A\left({{\displaystyle \overrightarrow{q}}}_d\right)A\left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}}\right)\left(1-A\left({{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right)\right)}{\theta \left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{Q}}\right)}\hfill \\ {\beta }^{-}\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)\hfill & =\frac{A\left({{\displaystyle \overrightarrow{q}}}_d\right)A\left({{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right)\left(1-A\left({{\displaystyle \overrightarrow{q}}}_d-{\displaystyle \overrightarrow{Q}}\right)\right)}{\theta \left({{\displaystyle \overrightarrow{q}}}_d;{{\displaystyle \overrightarrow{q}}}_d+{\displaystyle \overrightarrow{Q}}\right)}.\hfill \end{array}$$(22)

Each thus aims to access one of the double overlap areas within the ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ distribution, with a phase shift term inserted to compensate aberrations.

2.3.4 Phase contrast transfer function

The contrast transfer function (CTF) ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$, depicted in Figure 1, is introduced due to the summation over ${\beta }^{\pm }\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$ and is given by

$${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_0}}A\left({{\displaystyle \overrightarrow{q}}}_0\right)A\left({{\displaystyle \overrightarrow{q}}}_0-{\displaystyle \overrightarrow{Q}}\right)\left(1-A\left({{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)\right).$$(23)

As such, it is equal to the surface of the double overlap region, across camera space ${{\displaystyle \overrightarrow{q}}}_d$, corresponding to each spatial frequency ${\displaystyle \overrightarrow{Q}}$. It is interesting to note that, while this CTF is peaked at intermediary frequencies, i.e. close to q_A, it decays for both higher and lower frequencies.

Here, it should be furthermore highlighted that, while the CTF of SBI-S is explicit, SBI-D possesses the same fundamental characteristics with regards to frequency transfer, as shown by the dependence of $\omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$ on ${\displaystyle \overrightarrow{Q}}$. In essence, it is not the choice between the summative or the deconvolutive forms that leads to the frequency weighting, but rather the assumption of a weakly scattering object in itself.

In particular, if equation (20) is fulfilled, and if no geometrical aberrations are present, the parts of the distribution ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ found outside the double overlap areas are expected to carry no useful information on the specimen, and thus to contain only noise. ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ then reflects the information content in the scattering data itself, being equal to the proportion of available scattering vectors ${\displaystyle \overrightarrow{q}}$ that are useful to recover a specific frequency component ${\displaystyle \overrightarrow{Q}}$ of the specimen, i.e. it constitutes the phase contrast transfer function (PCTF) of the experiment in the sense of e.g. reference [104]. When the WPOA is fulfilled, it is thus expected to intrinsically apply to all STEM-based phase retrieval methods, irrespective of whether their formulation assumes a weakly scattering object in the first place.

As such, the SBI method, which consists in a treatment based explicitly on equation (17), permits to exclude all pixels outside double overlaps, thus in principle minimizing the total noise in the real-space result. As explained in details in reference [149], the SBI-based treatment of counts in the detector space ${{\displaystyle \overrightarrow{q}}}_d$, which follow Poisson statistics [150], then leads to a predictable noise level added to the frequency spectrum of the recovered object, given by the square root of the PCTF.

In this context, the option to deconvolve the reconstructed phase shift $\sigma \mu \left({\displaystyle \overrightarrow{r}}\right)$ post-process with ${ℱ}^{-1}\left[\sqrt{{\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)}\right]\left({\displaystyle \overrightarrow{r}}\right)$ has been proposed as an effective noise normalization strategy [106,149], i.e. rendering the noise level homogeneous across spatial frequencies. Deconvolving by the complete $\zeta \left({\displaystyle \overrightarrow{r}}\right)$, which in the conventional SSB workflow [103,104] is equivalent to averaging ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ within the double overlap areas instead of performing a summation, may otherwise permit to homogenize frequency transfer. This is nevertheless only practical when the dose is high enough, as the amplitude of the frequency components to be amplified may then be below the noise level.

2.3.5 Applicability of the weak scatterer approximation

It should be understood that the specific situation where $\sigma \mu \left({\displaystyle \overrightarrow{r}}\right)$ possesses the low value range of a weak phase object only occurs in a handful of cases. This may not only be due to excessive atomic potentials or material thicknesses, but also because lower acceleration voltages U imply non-linearly higher values for the interaction parameter σ, as shown by equation (2). In the case where the illuminated object is not a weak scatterer, as noted e.g. in reference [106,151], the SBI process still imposes a frequency-wise attenuation following ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$, due to the forms of $\omega \left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$ and ${\beta }^{\pm }\left({\displaystyle \overrightarrow{Q}};{{\displaystyle \overrightarrow{q}}}_d\right)$. The CTF is then method-induced rather than reflective of the PCTF of the experiment itself.

In this situation, the underlying sideband-like geometry in the distribution ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ also cannot be expected to occur strictly, i.e. the values taken by ${{\displaystyle \overrightarrow{q}}}_d$-coordinates within the ${\displaystyle \overrightarrow{Q}}$-wise double overlap areas may not be homogeneously equal to $\pm \sigma {\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right)$ anymore and exploitable information may be present outside as well. On that second aspect, it is noteworthy that, under the more general POA, non-zero coordinates of ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ include both the triple overlap areas and the dark field. In contrast, the scattering of electrons outside the primary beam is not possible in the framework of the WPOA, as directly noticeable in equations (17) and (20). Whereas this is not the case for the iCoM and WDD methods, the SBI-S and SBI-D processes are thus unable to exploit dark field electrons.

As a side-note, the CTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ leads to artificial image features, e.g. negative halos around atomic sites [152]. If the specimen is not weakly scattering, such artificial features are not expected to occur through methods based only on the POA, hence demonstrating the possible violation of the WPOA upon comparison of results.

2.3.6 Interests of the deconvolutive approach

Finally, when comparing the two forms of sideband ptychography, SBI-D presents a few advantages compared to the already established SBI-S approach. First, the MTF $M\left({{\displaystyle \overrightarrow{r}}}_d\right)$ can be explicitly included. Second, SBI-D permits the use of an arbitrarily shaped aperture [116] where the selection of specific overlapping regions would be less obvious, and thus the straightforward application to specific phase plate designs [153–155].

2.4 Center of mass imaging

2.4.1 Scan position-wise average momentum transfer

The center of mass (CoM) ${⟨{\displaystyle \overrightarrow{q}}⟩}_{{{\displaystyle \overrightarrow{r}}}_s}$ of the scan position-dependent CBED patterns constitutes a measurement of the average momentum transfer between the scattered electrons and the specimen. Under the POA, the CoM is linearly related to the local gradient of the projected potential [114]. Formally, this means that

$$\begin{array}{ll}{⟨{\displaystyle \overrightarrow{q}}⟩}_{{{\displaystyle \overrightarrow{r}}}_s}\hfill & ={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_d}}{{\displaystyle \overrightarrow{q}}}_d{I}_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill \\ \hfill & =\gamma \left({{\displaystyle \overrightarrow{r}}}_s\right){\otimes }_{{{\displaystyle \overrightarrow{r}}}_s}\left(\frac{\sigma }{2\pi }{{\displaystyle \widehat{{\displaystyle \overrightarrow{▽}}}}}_{{{\displaystyle \overrightarrow{r}}}_s}\mu \left({{\displaystyle \overrightarrow{r}}}_s\right)\right),\hfill \end{array}$$(24)

with ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ a CTF given by

$${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{q}}}_0}}A\left({{\displaystyle \overrightarrow{q}}}_0\right)A\left({{\displaystyle \overrightarrow{q}}}_0-{\displaystyle \overrightarrow{Q}}\right).$$(25)

${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ is depicted in Figure 2. This CTF is peaked at low frequencies and smoothly decays as a function of ${\displaystyle \overrightarrow{Q}}$, thus indicating difficulties in transferring higher frequencies.

In the absence of a fast DED to perform an MR-STEM experiment, the average momentum transfer is conventionally approximated by using quadrants of a segmented detector [156,157], a technique usually referred to as differential phase contrast (DPC) in relation to historical references [158,159]. In that context, CoM imaging can be understood as a more accurate approach to measuring the DPC signal [44], in particular considering that the use of segmented detectors leads to a non-isotropic CTF [160]. As a side-note, another existing detector paradigm consists in a position-sensitive non-pixelated design [161].

2.4.2 Fourier-based extraction of the potential

Following the measurement of ${⟨{\displaystyle \overrightarrow{q}}⟩}_{{{\displaystyle \overrightarrow{r}}}_s}$, an extraction of the projected potential can be performed through a simple Fourier integration scheme, as conventionally used e.g. for the integrated DPC (iDPC) [115,162] counterpart to iCoM. This consists in

$$\begin{array}{ll}h\left({\displaystyle \overrightarrow{Q}}\right)\hfill & =\frac{{\displaystyle \overrightarrow{Q}}\cdot ℱ\left[⟨{\displaystyle \overrightarrow{q}}⟩\left({{\displaystyle \overrightarrow{r}}}_s\right)\right]\left({\displaystyle \overrightarrow{Q}}\right)}{i\left(ϵ+\left|\right|{\displaystyle \overrightarrow{Q}}|{|}^2\right)}\hfill \\ \gamma \left({\displaystyle \overrightarrow{r}}\right){\otimes }_{{\displaystyle \overrightarrow{r}}}{\mu }^{iCoM}\left({\displaystyle \overrightarrow{r}}\right)\hfill & =\frac{1}{\sigma }{ℱ}^{-1}\left[h\left({\displaystyle \overrightarrow{Q}}\right)\right]\left({\displaystyle \overrightarrow{r}}\right),\hfill \end{array}$$(26)

with $h\left({\displaystyle \overrightarrow{Q}}\right)$ an intermediary result. Importantly, just like for the WDD and SBI methods, the DC component is inaccessible, as shown by the scalar product with ${\displaystyle \overrightarrow{Q}}$.

Importantly, in contrast to analytical ptychography, aberration correction does not seem straightforward with CoM imaging. The dependence on focus and thickness has nevertheless been investigated in recent years [163–167], in particular with the objective of maintaining an interpretable contrast when the POA is not strictly fulfilled anymore.

2.4.3 Optical transfer function

Continuing, in contrast to the PCTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$, which is applicable in the situation where the WPOA is fulfilled and is then reflective of the information content of the scattering data itself, ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ is fully process-induced and is derivable in the more general context of the POA. In particular, it is equal to the Fourier transform of the unaberrated probe intensity [114] and, as such, can be termed as an optical transfer function (OTF) in the sense of [168]. Beyond that, the immediate consequence of this OTF is the higher weighting of low-frequency features, compared to the rest of the object spectrum, being then attenuated.

As a result, the iCoM imaging mode is prone to low-frequency artefacts [162], which may constitute a limit to the use in the low-dose case [169]. This is nevertheless not problematic for many of the common applications of DPC and CoM consisting in the imaging of long-range features, such as e.g. charge density gradients [170], magnetic domain structures [171], large proteins [2], interfaces between materials [172], particle shapes [173], skyrmions [174] or stray electrostatic fields [175].

In principle, and like in the SBI case, it should furthermore be possible to compensate this effect by directly deconvolving the real-space result with $\gamma \left({\displaystyle \overrightarrow{r}}\right)$, though the limitation is then whether the frequency components to be amplified have been brought below the noise level. Hence, such a solution is not practical at low doses. In the specific situation where iCoM imaging is employed on a weak phase object, both the PCTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ and the OTF ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ can be expected to apply.

2.5 Scan-frequency partitioning algorithm

Algorithm 1 SFPA

1: Choose imaging method

2: Partition ${{\displaystyle \overrightarrow{r}}}_s$-coordinates in packets $P_{{{\displaystyle \overrightarrow{r}}}_s}$

3: Define ${\displaystyle \overrightarrow{Q}}$-grid

4: Partition ${\displaystyle \overrightarrow{Q}}$-coordinates in domains $D_{{\displaystyle \overrightarrow{Q}}}$

5: if imaging method is WDD then

6: Initialize intermediary result as $f\left({\displaystyle \overrightarrow{Q}}\right)=0$

7: else if imaging method is SBI-D then

8: Initialize intermediary result as $g\left({\displaystyle \overrightarrow{Q}}\right)=0$

9: else if imaging method is SBI-S then

10: Initialize intermediary result as $g\text{'}\left({\displaystyle \overrightarrow{Q}}\right)=0$

11: else if imaging method is iCoM then

12: Initialize intermediary result as $h\left({\displaystyle \overrightarrow{Q}}\right)=0$

13: Distribute $P_{{{\displaystyle \overrightarrow{r}}}_s}/D_{{\displaystyle \overrightarrow{Q}}}$ couples asynchronously

14: for each packet $P_{{{\displaystyle \overrightarrow{r}}}_s}$ do

15: for each domain $D_{{\displaystyle \overrightarrow{Q}}}$ do

16: if imaging method is WDD then

17: Calculate $J_{{\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}}^{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{R}}\right)$

18: Calculate $f_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

19 Add $f_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ to $f\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

20: else

21: Calculate ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}}^{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$

22: if imaging method is SBI-D then

23:Calculate $g_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

24: Add $g_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ to $g\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

25: else if imaging method is SBI-S then

26: Calculate $g{\text{'}}_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

27: Add $g{\text{'}}_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ to $g\text{'}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

28: else if imaging method is iCoM then

29: Calculate $h_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

30: Add $h_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ to $h\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$

31: if imaging method is WDD then

32: Divide $f\left({\displaystyle \overrightarrow{Q}}\right)$ by $\sqrt{f\left({\displaystyle \overrightarrow{0}}\right)}$

33: Inverse Fourier transform along ${\displaystyle \overrightarrow{Q}}$

34: Transmission function is measured

35: Extract angle of transmission function

36: Phase shift map is measured

37: else

38: Inverse Fourier transform along ${\displaystyle \overrightarrow{Q}}$

39: Phase shift map is measured

2.5.1 Motivations and numerical basis

One of the main limiting factor for the practical implementation of analytical ptychography is the necessity to first load the full dataset in computer memory, in order to perform a collective treatment consisting in a succession of fast Fourier transforms (FFT) and deconvolution/summation steps. Such a process requires a large available memory and makes e.g. GPU implementation difficult. This publication thus proposes a new scan-frequency partitioning algorithm, i.e. the SFPA solution mentioned in the introduction, which constitutes a straightforward, memory-limited and parallelizable implementation of the WDD, SBI and iCoM methods. As explained in more details below, this approach also relaxes sampling conditions that would normally be imposed by the scan grid.

The basis for the algorithm is the replacement of the FFT leading from the scan dimension ${{\displaystyle \overrightarrow{r}}}_s$ to the spatial frequencies ${\displaystyle \overrightarrow{Q}}$ with an explicit, term-by-term, summation, e.g. following the Einstein notation. This is conventionally referred to as the einsum algorithm, as included e.g. in several Python packages. A similar explicit construction of the Fourier series was used for instance in ref. [100] for live processing. The formal procedure is described in the following paragraphs, and is otherwise provided in Algorithm 1. Noteworthily, under the current implementation developed for this work, the PyTorch package [176] was chosen for its capacities in straightforward GPU-based programming.

2.5.2 Partitioning of calculation steps

The SFPA encompasses two distinct levels of partitioning among the calculation steps needed for the complete process. A first one is ensured by cutting the complete four-dimensional dataset $I_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ in packets of scan positions $P_{{{\displaystyle \overrightarrow{r}}}_s}$, each containing a user-defined number of arbitrarily chosen coordinates ${{\displaystyle \overrightarrow{r}}}_s$. The packets are treated individually, in particular the einsum-based Fourier transform, itself done for specific spatial frequency domains $D_{{\displaystyle \overrightarrow{Q}}}$. This then represents the second partitioning introduced in the algorithm.

The complete set of calculations is thus divided in a number of single independent operations, each involving a specific $P_{{{\displaystyle \overrightarrow{r}}}_s}/D_{{\displaystyle \overrightarrow{Q}}}$ couple. Those operations yield partial Fourier transformed datasets given by

$$\begin{array}{ll}{{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}}^{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill & ={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{r}}}_s\in P_{{{\displaystyle \overrightarrow{r}}}_s}}}{e}^{-i2\pi {\displaystyle \overrightarrow{Q}}\cdot {{\displaystyle \overrightarrow{r}}}_s}{I}_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)\hfill \\ J_{{\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}}^{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{R}}\right)\hfill & ={\displaystyle {\displaystyle \sum }_{{{\displaystyle \overrightarrow{r}}}_s\in P_{{{\displaystyle \overrightarrow{r}}}_s}}}{e}^{-i2\pi {\displaystyle \overrightarrow{Q}}\cdot {{\displaystyle \overrightarrow{r}}}_s}{ℱ}^{-1}\left[I_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)\right]\left({\displaystyle \overrightarrow{R}}\right),\hfill \end{array}$$(27)

depending on the type of reconstruction performed, i.e. $J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)$ is the input of a WDD process while ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ is needed for SBI and iCoM. Note that in equation (27), the term ℱ⁻¹ indicates an inverse Fourier transform done over the camera space, for the few CBED patterns in $P_{{{\displaystyle \overrightarrow{r}}}_s}$.

Each partial Fourier transformed dataset is used for one of the following calculations

$$\begin{array}{ll}f\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill & \leftarrow f\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)+f_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill \\ g\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill & \leftarrow g\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)+g_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill \\ g\text{'}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill & \leftarrow g\text{'}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)+g{\text{'}}_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill \\ h\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)\hfill & \leftarrow h\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)+h_{P_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right),\hfill \end{array}$$(28)

where the packet-specific intermediary result $f_{D_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ is obtained through equations (13) and (14), $g_{D_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ through equation (19), $g{\text{'}}_{D_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ through equation (21) and $h_{D_{{{\displaystyle \overrightarrow{r}}}_s}}\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ through equations (24) and (26). Performing the same process for the entirety of the ${{\displaystyle \overrightarrow{r}}}_s$-to-${\displaystyle \overrightarrow{Q}}$ components of the Fourier series finally yields the full reconstruction result.

Noteworthily, in the case of the WDD method, the complete, four-dimensional, Wigner distribution ${\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ is not explicitly retrieved. Instead, in the implementation described by Algorithm 1, each $P_{{{\displaystyle \overrightarrow{r}}}_s}/D_{{\displaystyle \overrightarrow{Q}}}$ couple leads to an increment of $f\left({\displaystyle \overrightarrow{Q}}\in D_{{\displaystyle \overrightarrow{Q}}}\right)$ directly. The same einsum-based Fourier transform strategy could nevertheless be used for this purpose, i.e. without an immediate summation step across ${\displaystyle \overrightarrow{R}}$, straightforwardly as well.

2.5.3 Opportunity for parallelization

A first interest of the scan-frequency partitioning algorithm is its low need in active memory, since the size of the packets $P_{{{\displaystyle \overrightarrow{r}}}_s}$ and domains $D_{{\displaystyle \overrightarrow{Q}}}$ are chosen by the user directly. This in turn permits to adapt the process to the computer memory available, including as part of a straightforward implementation on a GPU, e.g. involving specialized Python-based procedures [176]. Furthermore, since the treatment of each individual $P_{{{\displaystyle \overrightarrow{r}}}_s}/D_{{\displaystyle \overrightarrow{Q}}}$ couple is independent of all others, parallelization is possible along both the ${{\displaystyle \overrightarrow{r}}}_s$ and ${\displaystyle \overrightarrow{Q}}$ dimensions. In comparison, the implementation reported in reference [100] only allowed it along ${{\displaystyle \overrightarrow{r}}}_s$, though it was already enough for live processing using a computer with sufficient performance.

Given the low memory requirement of a single $P_{{{\displaystyle \overrightarrow{r}}}_s}/D_{{\displaystyle \overrightarrow{Q}}}$ calculation, such a parallelization strategy is in principle implementable on a wider range of devices, including low-end. Though extensive numerical benchmarking was left for future work, it should be noted that avoiding the two-dimensional ${{\displaystyle \overrightarrow{r}}}_s$-to-${\displaystyle \overrightarrow{Q}}$ FFT can be expected to lead to an increment in the numerical complexity of the complete process. Specifically, it then goes from the typical O (N_s;x ⋅ N_s;y ⋅ log(N_s;x ⋅ N_s;y)) to $O\left(N_{s;x}\cdot N_{s;y}\cdot N_{{\displaystyle \overrightarrow{Q}}}\right)$, with N_s;x/N_s;y the number of positions along the two scan axes and $N_{{\displaystyle \overrightarrow{Q}}}$ the total number of frequencies ${\displaystyle \overrightarrow{Q}}$ used. This number is equal to

$$N_{{\displaystyle \overrightarrow{Q}}}=4\pi S_{rec}{q_A}^2,$$(29)

where π (2q_A) ² is the reconstructed frequency surface and S_rec is the reconstructed real-space surface. Specifically, the discretized frequencies ${\displaystyle \overrightarrow{Q}}$ are distributed within a disk of radius 2q_A, with a pixel density determined by the real-space extent of the reconstruction window.

Note that, in order to avoid aliasing and periodicity artefacts, S_rec has to be made sufficiently larger than the scanned surface S, with a portion of it maintained at a value of zero in the real-space reconstruction window. The same is true in the frequency space window [177]. Overall, the preparation of the reconstructed object in both ${\displaystyle \overrightarrow{Q}}$- and ${\displaystyle \overrightarrow{r}}$-space follows the conventional approach used, for instance, in typical multislice simulations, as is described e.g. in reference [178].

2.5.4 Decorrelation of scan and reconstruction grids

Perhaps most importantly, and as implied by equation (29), the employment of the einsum algorithm permits the explicit decorrelation of the scan and frequency dimensions. As such, the real-space reconstruction grid, and thus the actual choice of ${\displaystyle \overrightarrow{Q}}$-coordinates for which the result is calculated, is prepared independently of the scan grid, and the formal contribution of each given ${{\displaystyle \overrightarrow{r}}}_s$-coordinate to a single arbitrary frequency ${\displaystyle \overrightarrow{Q}}$ may be determined separately from all others. In contrast, in the conventional full FFT solution [103,113], as well as in reference [100], each scan point equates one pixel in the result. The SFPA however permits a calculation at frequencies exceeding e.g. the maximum that would then be allowed by the finite scan interval. In that context, it becomes possible, for instance, to perform a reconstruction given a strongly defocused probe and less scan positions [118], while conserving an appropriately resolved reconstruction window in real-space.

At first sight, this development may be understood as breaking the Nyquist-Shannon sampling theorem. The possibility of retrieving information beyond the Nyquist frequency of the scan grid should however be seen as resulting from the usage of the information available along the ${{\displaystyle \overrightarrow{q}}}_d/{\displaystyle \overrightarrow{R}}$ dimensions. In particular, in defocused conditions, the amount of details contained in the far-field intensity is greater as it then consists in a shadow image of the specimen [179]. As such, including the appropriate aberration function in $\theta \left({{\displaystyle \overrightarrow{q}}}_0;{{\displaystyle \overrightarrow{q}}}_0+{\displaystyle \overrightarrow{Q}}\right)$ permits to correctly translate this information back into the result and extend the area, in the real-space reconstruction window, that may be informed by a single scan position. Beyond this, since such a measurement geometry is already commonly used in combination with iterative approaches [117,118], where a similar decorrelation of the scan and reconstruction pixels is implicit and where the same scattering data is used, this fundamental ability of analytical ptychography is expected.

Another interest of using an arbitrary set of reconstruction frequencies is the facilitated implementation of high-pass and low-pass filtering, as the concerned ${\displaystyle \overrightarrow{Q}}$-coordinates can be omitted from the calculation, hence reducing $N_{{\displaystyle \overrightarrow{Q}}}$ as well. A limitation to this practice is however that, in order to perform an extraction of the phase shift map from the WDD-retrieved transmission function, it should have a defined zero-frequency component, i.e. the mean of the complex numbers $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$ should not be equal to zero. For this reason, only the inherent high-pass filtering, and not both the high- and low-pass ones, may be used for the SFPA-based WDD calculation.

The employment of an orthonormal Fourier transform may moreover be useful in ensuring appropriate numerical normalization, as the number of pixels in the reconstruction grid and the scan grid are likely to differ. This proper choice of convention is important to conserve consistent results among slightly different scan grids among recordings.

Finally, the operations involving the dimensions ${{\displaystyle \overrightarrow{q}}}_d$ and ${\displaystyle \overrightarrow{R}}$, as in equations (14), (19), (21) and (24), are performed in the camera space. As such, ${\displaystyle \overrightarrow{R}}$ represents a spatially limited kernel, similarly to e.g. reference [70], in which numerical artefacts are prevented by a simple interpolation or zero-padding step. In that manner, the process can also precisely account for elliptical distortions observed in the far-field pattern [180] and prevent them from affecting the result, based on an initial calibration of the ${{\displaystyle \overrightarrow{q}}}_d$ dimension [166].

3 Atomic-resolution imaging of MoS₂

3.1 Conventional focused-probe conditions

3.1.1 Simulation and processing parameters

In order to test the dose-efficiency of the iCoM, SBI-D and WDD methods in the conventional focused-probe, high-resolution, condition of electron ptychography, an MR-STEM simulation was performed based on a monolayer MoS₂ specimen, for which the POA can reasonably be considered fulfilled. Diffraction patterns were calculated in a scan grid of 64 by 64 points, covering an area of 2 nm by 2 nm, hence with an interval of about 32 pm. Illumination conditions were chosen as representative of the capacities of a modern aberration-corrected microscope such as e.g. a Titan Themis 60–300 (Thermo Fisher Scientific). Specifically, the acceleration voltage U and the semi-convergence angle α were assigned values of 60 kV and 30 mrad respectively which, for reference, leads to a Rayleigh criterion δr_Rayleigh = 0.61/q_A of about 99 pm. As illustrated in Figure 3, the optical conditions described above lead to an area overlap ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}\approx 79.0\%$ between two scan points neighboring each other along a scan direction, and 62.3% along the diagonal. ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ is defined for higher distances $\delta {{\displaystyle \overrightarrow{r}}}_s$ as well. This highlights that a degree of redundancy remains beyond immediate neighbors, which is exploited by the ptychographic process as well. More details on the calculation of the area overlap and how it differs from the conventional approach [135] can be found in Appendix A.

Continuing, the propagation of the electron wavefunction through the specimen was modeled based on the multislice approximation [181–183] and the atomic potentials were calculated using parameterized hydrogen orbitals as described in reference [184]. The specimen potential was sliced below the atomic plane level and pixelated such that a maximum scattering vector of up to twice the range actually used could be included. Thermal motion within the lattice was accounted for by repeating the calculation for a total of 64 configurations of random lateral atomic shifts, and averaging the resulting distributions $I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)$. The random shift vectors were determined using the frozen phonon approximation [185,186] based on the Einstein model, i.e. assuming non-correlated atomic vibrations [187]. For simplicity, and also because this work aims at reproducing results obtainable with a Timepix3 chip [32,47] at a low acceleraton voltage U, thus in a condition where multiple counting is unlikely to occur [41,48,188], the simulation did not include an explicit MTF. As such, the values taken by the $M\left({{\displaystyle \overrightarrow{r}}}_d\right)$ function, included in practice in the SBI-D and WDD calculations, only encompassed the role of the finite pixel size of the simulated camera, as implied by the kernel size.

The results of the iCoM, SBI-D and WDD processes, implemented using the SFPA approach described in Section 2.5, are depicted in Figure 4. Specifically, the measurements of the projected potential, expressed in V⋅nm, are shown alongside the square roots of the corresponding Fourier transform amplitudes, for visualization of Fourier weightings along the two-dimensional ${\displaystyle \overrightarrow{Q}}$ coordinates. The calculations were done for a variety of average numbers of electrons per pattern $N_{e^{-}}$ and consequent doses D given in e⁻/Ų. To better highlight the non-linear relation between dose and contrast, $N_{e^{-}}$ was given values of 2 ^l with l ∈ [2, 3, ..., 10]. Dose-limitation was ensured by repeated random pixel selection, with the number of repetitions being probabilistically determined across the scan window by Poisson statistics. As such, the propagation of Poisson noise from detector space to the reconstruction window can be straightforwardly reproduced, while conserving a realistic sparsity in the scattering frames [105]. This approach is described in more details in Appendix B, alongside its wider interests.

Furthermore, for each dose-limited case, the generated sparse diffraction patterns were individually normalized by their sum, pre-treatment. This strategy was adopted for all reconstructions presented in this work and was chosen following the suggestion of reference [149]. This is equivalent to varying the normalization of the wavefunction, scan point-wise, to match the number of counts in each corresponding pattern. Importantly, taking this normalization choice into account will be required for any theoretical estimation of measurement precision in future work, as it leads, in effect, to a change in the variance of single patterns. This solution differs from the usual quantitative STEM approach [189,190], which would have consisted in uniformly normalizing by $N_{e^{-}}$. Finally, the $N_{e^{-}}=+\infty$ case corresponds to the direct use of the simulated $I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)$, where the intensity is implicitly normalized. The corresponding result can thus be understood as representing the experimental situation where the best achievable dose-dependent precision is reached, and hence where the noise level is negligible.

As a side-note, in the case of the WDD process, the projected electrostatic potential is obtained through a prior extraction of the phase shift map. Given that the [− π ; + π] range was not exceeded, no discontinuities were observed and thus no unwrapping was necessary.

3.1.2 Noise level in the micrographs

For the three methods, atomic patterns are already visible from $N_{e^{-}}=8$, hence with a dose below D = 81.92e⁻/Ų. Moreover, frequencies belonging to the specimen lattice are observed even in the Fourier transforms of results obtained given $N_{e^{-}}=4$. This first remark is particularly interesting for future applications of electron ptychography to beam-sensitive objects, as it empirically shows what is the true requirement in terms of dose, given a perfectly stable and coherent imaging system. As $N_{e^{-}}$ increases, the noise level in the images lessens and specimen frequencies become more dominant compared to the noise background. Such a dose-dependent precision in ptychographic computational imaging has been investigated empirically in the literature [76,88,91,191], as is done here as well, and its lowest achievable value can in principle be predicted by parameter estimation theory [192–197], in particular using the Cramér-Rao lower bound (CRLB) [198].

In this publication, the true frequency-dependent CRLB is not provided since, unless some simplifications such as the WPOA [195,196] are introduced, its formulation remains specific to the specimen [193,194]. The establishment of a general ${\displaystyle \overrightarrow{Q}}$-dependent metric, which would be dependent on the complete set of experimental parameters, is thus left for future work. Beyond that, the approximation made in reference [197], provided below, leads to a single number CRLB_RS representing the minimum standard deviation among distinct measurements, as induced by the propagation of Poisson noise [150], upon retrieving the phase shift map $\sigma \mu \left({\displaystyle \overrightarrow{r}}\right)$ in real-space. While it was derived in ideal illumination conditions which are not met here, e.g. the total illumination is not strictly restricted to the scanned area, this metric remains useful to establish a fundamental understanding of the concept.

$$\text{CRL}{\text{B}}_{RS}=\sqrt{\frac{N_{{\displaystyle \overrightarrow{Q}}}}{2N_s{N}_{e^{-}}}}\ge \sqrt{\frac{2\pi {q_A}^2-\frac{1}{2S}}{D}}.$$(30)

$N_s{N}_{e^{-}}$ represents the total number of probing electrons, while N_s is the total number of scan positions used. S is the surface covered by the scan window, necessarily smaller than the surface S_rec as explained in Section 2.5. The number $N_{{\displaystyle \overrightarrow{Q}}}$ of reconstructed frequencies was otherwise described by equation (29), and can here be understood as the number of useful pixels in the reconstruction. A “–1” term is added to account for the unmeasurable DC component. Given that the term 0.5/S is likely to be negligible compared to 2πq_A², CRLB_RS shows rather clearly that, in order to achieve a certain goal in measurement precision, the dose D has to be adapted to the aperture radius q_A, and thus implicitly to the spatial resolution [199].

3.1.3 Fourier ring correlations

In order to pursue the analysis further, Fourier ring correlations FRC ^m (k) [200,201], shown in Figure 5, were calculated from the projected potential results through

$$\begin{array}{ll}\hfill & \text{FR}{\text{C}}^m\left(k\right)=\frac{{\text{X}}^m\left(k\right)}{{\text{Y}}^m\left(k\right)}\hfill \\ \hfill & {\text{X}}^m\left(k\right)={\displaystyle {\displaystyle \sum }_{\parallel {\displaystyle \overrightarrow{Q}}\parallel \in R_k}}{{\displaystyle \tilde{\mu }}}_m\left({\displaystyle \overrightarrow{Q}}\right){{{\displaystyle \tilde{\mu }}}_{+\infty }}^{*}\left({\displaystyle \overrightarrow{Q}}\right)\hfill \\ \hfill & {\text{Y}}^m\left(k\right)=\sqrt{{\displaystyle {\displaystyle \sum }_{\parallel {\displaystyle \overrightarrow{Q}}\parallel \in R_k}}\left|{{\displaystyle \tilde{\mu }}}_m\left({\displaystyle \overrightarrow{Q}}\right){|}^2{\displaystyle {\displaystyle \sum }_{\parallel {\displaystyle \overrightarrow{Q}}\text{'}\parallel \in R_k}}\right|{{\displaystyle \tilde{\mu }}}_{+\infty }\left({\displaystyle \overrightarrow{Q}}\text{'}\right){|}^2}.\hfill \end{array}$$(31)

k is a spatial frequency modulus and R_k = [k − δk ; k] is the corresponding annular domain, with δk a case-dependent precision. ${{\displaystyle \tilde{\mu }}}_m\left({\displaystyle \overrightarrow{Q}}\right)$ is the Fourier transform of the measured projected potential, for the specific $N_{e^{-}}=m$ case. In this context, the calculated FRC can be interpreted as a frequency-wise measurement of the dose-dependent precision of each method, and thus provides a straightforward dose-efficiency metric. The closer FRC ^m (k) is to 1, for a given spatial frequency modulus k, the closer the corresponding R_k range of the signal is to reaching the best achievable precision. In this example, and for better visibility, the FRC curves are provided in a reduced selection of $N_{e^{-}}=4,16,64,256,1024$, for iCoM in Figure 5a, SBI-D in Figure 5b and WDD in Figure 5c.

A few observations can immediately be done from the calculated FRC. For all imaging modes, three peaks are observed, close to 0.5,1.0 and 1.5 times q_A, as well as an emerging fourth one. Those correspond to the hexagonal pattern of spatial frequencies belonging to the specimen, mirroring scattering orders of the atomic lattice, as observed in the Fourier transforms of Figure 4 as well. At the level of the peaks, perfect precision is reached at a much lower $N_{e^{-}}$ than in the rest of the k-axis. From a naive standpoint, this already tends to show that frequencies ${\displaystyle \overrightarrow{Q}}$ actually carrying information on the specimen are reconstructed much more efficiently than those containing no information, and which then end up reaching ${\displaystyle \tilde{\mu }}\left({\displaystyle \overrightarrow{Q}}\right)=0$ at infinite dose. This is expected in a situation where the spectrum of the illuminated object is sparse. In particular, ptychographically processed electrons end up contributing only to the recovered projected potential, i.e. its associated ${\displaystyle \overrightarrow{Q}}$-coordinates. This remains true as long as the illumination characteristics are known and no artefactual features are introduced, e.g. from an inaccurate interaction model.

In this context, if one considers the overall calculation as an additive inclusion of single counts' contributions to the measurement, with no question of normalization, the signal-to-noise ratio at ${\displaystyle \overrightarrow{Q}}$-coordinates belonging to the specimen is expected to directly improve for each dose increment, while the noise level at other frequencies remains the same. The reduction of the background noise is then due to the normalization, and thus occurs at a lower rate than the retrieval of specimen information in itself. Continuing, upon comparing the three frequency peaks mentioned above, it is also noticeable that FRC^m (k) ≈ 1 occurs with more difficulty as k increases, i.e. higher values of ${\displaystyle \overrightarrow{Q}}$ appear more dose-expensive at first sight. The practical reason for it is that the surface covered by an R_k ring increases with k, which thus implies a larger proportion of background noise compared to specimen frequencies. Similarly, the lower values of k lead to a less visually stable value of FRC, due to the low number of actual pixels in the corresponding R_k.

3.1.4 Contrast transfer capacities

Arguably the most important information to draw from Figure 5 is that, among the three investigated methods, the overall FRC profiles are rather similar. This is of particular interest, as it shows that, for a given frequency component, the dose-efficiencies of iCoM, SBI and WDD are more-or-less the same, in that they reach the best achievable result with comparable dose requirements. As such, what differentiates those imaging modes with regards to the measurement precision in real-space, and in particular to the noise background formed as a function of spatial frequency [149], is the existence of the CTF ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ for iCoM and ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ for SBI, displayed in Figures 2 and 1. Whereas those CTF lead to visually different micrographs, as is directly noticeable in Figure 4, the attenuation of frequency components also contribute to noise filtering. Consequently, ${\mu }^{SBI}\left({\displaystyle \overrightarrow{r}}\right)$ appears slightly, but noticeably, less noisy than ${\mu }^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$, e.g. for $N_{e^{-}}\le 64$. This is specifically related to the existence of high-frequency noise, as observed in the Fourier transforms of the WDD results, which is otherwise eliminated by the deconvolutive SBI process.

Beyond that, the underlying difference in the real-space measurement, between the two analytical ptychography methods, consists in an exaggerated dark halo around atomic sites, clearly visible at higher doses and only present in the SBI-D results. This feature is associated to ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ [152], which is thus shown to not intrinsically apply to the WDD result. As was explained in Section 2.3, this then constitutes a clear indication that the specimen is not a weak scatterer, i.e. ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ cannot be considered to constitute the PCTF of the experiment in general. Specifically, the WDD process only assumes the more general POA, here remaining reasonable, and should not be generally expected to show the specific frequency transfer met in the case of a weak phase object. The formulation of the SBI method, on the other hand, is still based on this assumption and will thus have ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ as a CTF in any case. Note that the violation of the WPOA is here further confirmed by the range of values covered by the WDD phase shift map, which is above 1.0 rad, as well as by the ground truth of the projected potential, for reference depicted in Figure 6 in a limited real-space window.

More fundamentally, the nonfulfillment of the WPOA means that, in the ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ distribution, the general sideband-like geometry arising from a weakly scattering specimen is not met in practice. As a consequence, this is not just noise, e.g. in the triple overlap areas, that is removed by the SBI process, but also potentially useful information on the specimen itself. This can be verified in Figure 4 as well, where the Fourier components of ${\mu }^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$, at high $\parallel {\displaystyle \overrightarrow{Q}}\parallel$ values, are visibly higher than those of ${\mu }^{SBI}\left({\displaystyle \overrightarrow{r}}\right)$. This thus leads, in addition to the absence of the artificial features mentioned above, to a slightly better resolution in the WDD measurement.

In parallel, the iCoM imaging mode is affected by the OTF ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$, whether the PCTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ of a weak phase object is applicable or not. As such, the higher weighting of low frequencies, with the rest of the spectrum being then attenuated, increases its susceptibility to long-range artefacts [162,169], as explained in Section 2.4 Consequently, low-frequency noise remains dominant up to e.g. $N_{e^{-}}\le 256$. This is verified by the Fourier transforms as well.

3.1.5 Phase shift value range

Continuing with the projected potential maps displayed in Figure 4, one last remark remains to be made. In high dose conditions, the range of values obtained with WDD is about twice larger than it is with iCoM and SBI-D. This is not due, for instance, to a normalization issue, as the WDD result is obtained by extracting the angle of the initially retrieved transmission function. Furthermore, as is shown in Section 3, the opposite situation can be met as well and, as proven in Section 2.2, the aberration function plays a role too. This thus points out the mismatch in value range as being a fundamental feature of the imaging method rather than a numerical issue. In that respect, it is also worth noting that such a mismatch was observed previously in the literature [106,108,169] as well. Finally, the effect is likely amplified by the higher resolution of the WDD reconstructions, itself due to the better transfer of high frequencies, leading to stronger atomic peaks in the image.

3.2 Overfocused illumination conditions

3.2.1 Simulation and processing parameters

To test analytical ptychography methods in overfocused illumination conditions, a second simulation was performed in the same conditions as described in Section 2.1, though with an added defocus of 40 nm, and only 16 by 16 scan positions leading to an interval of about 133 pm in the scan grid. The simulation and reconstruction windows were enlarged to avoid artefactual probe self-interference. Noteworthily, the reconstruction of frequencies exceeding the maximum that is allowed, in principle, by the scan interval [103,113] is enabled by the SFPA solution described in Section 2.5.

The illumination condition leads to ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}\approx 93.8\%$ between neighboring scan points, as shown in Figure 7. Whereas this somewhat high area overlap was found to be fully sufficient for the reconstruction, another attempt with only 8 by 8 points in the same region, which would have permitted 85.7% of overlap, was found to be insufficient. Empirically, this need for a significant value of ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ is expected from the literature [117,118,146,147] available on the use of defocused probes.

Continuing, because the correction of aberrations is not possible in the conventional framework of iCoM imaging, this method is not used in this subsection. As the number of distinct acquisitions is reduced here, the selection of $N_{e^{-}}$ values is adapted as well to include numbers 2 ^l with l ∈ [6, 7, ..., 14]. Noteworthily, the increment in the number of electrons per pattern is consistent with the higher complexity of those patterns, as they then constitute shadow images of the specimen [179].

Given the use of a large illumination including internal features, as shown in the inset of Figure 7, formula (33), used in the previous section to establish the dose, does not hold anymore. This is because a significant part of the incident intensity on the specimen surface ends up probing the area outside the scan window, hence the dose D serves to recover information from an inhomogeneously sampled surface, which is larger than S. This is demonstrated in Figures 8a and 8b, where the real and imaginary parts of the reconstructed transmission function $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$, for the $N_{e^{-}}=+\infty$ case, are shown in an extended field of view, though still contained in the normal reconstruction window. As can be seen directly, supplementary specimen information is obtained outside of the actual scanned area. A slight inhomogeneity may furthermore appear within the central scanned surface itself, with the regions close to the corners receiving less intensity overall. While it is not very striking, this noticeably occurs in the present case, as can be observed in Figure 8a, where a cross-like pattern is visible in the real part of the retrieved transmission function.

3.2.2 Comparison to the focused-probe case

Projected potential measurements by WDD and SBI-D, expressed in V⋅nm, are displayed in Figure 10, for the specified values of $N_{e^{-}}$, alongside the corresponding square roots of Fourier transform amplitude. Similarly to the conventional focused-probe case, specimen frequencies are already detected at $N_{e^{-}}=64$, although a clear observation of real-space features in the scan window is, arguably, only possible for $N_{e^{-}}\ge 128$. One further qualitative observation can be made on the noise in the micrographs, which seems more persistent than in those shown in Section 2.1. This is expected, since the use of a delocalized illumination implies an equivalent spread of retrievable information per recording, as explained in the previous paragraph.

Another important difference between the overfocused and the conventional cases is a slight loss of resolution, e.g. consistent with comparisons made in reference [147]. As the true aberration function of the illumination was included in the process, and since the number of single recordings and camera pixelisation were high enough, this cannot be attributed to an insufficiency of available scattering information [145,202] or a processing error. In particular, increasing the number of scan positions to 32 by 32 did not improve the resolution, hence showing no further need in overlap ratio. Another explanation can be found in the inherent information content of the acquired MR-STEM dataset, as determined by the CRLB [195,196]. In other words, different illumination conditions, including probe focus [196], may possess specific capacities to transfer specimen frequencies to the acquired data, hence leading to a supplementary ${\displaystyle \overrightarrow{Q}}$-dependent weighting in the result. In that context, ptychographic reconstructions, with their practical resolution limits, can be expected to remain probe-specific, even when this probe is known or refined in-process [69,70,107,148].

3.2.3 Fourier ring correlations

Moreover, Fourier ring correlations FRC^m(k) were calculated for the overfocused probe case and are displayed in Figures 10a for SBI-D and 10b for WDD. They show essentially the same features as were observed in Figure 5, in particular with four peaks at coordinates k corresponding to specimen frequencies. A difference is however found in the apparently lower dose-efficiency, which can be attributed to the wider illuminated area, as explained above. Noteworthily, this effect is likely amplified by the larger amount of pixels per frequency ranges R_k, which is due to the greater size of the reconstruction window, leading to a more important weighting of frequency coordinates without crystal lattice information. This is particularly visible when comparing the first and second peaks to the third and fourth ones.

3.2.4 Phase shift value range

Going back to the micrographs themselves, a few more remarks can be made on the ranges of value covered by the projected potential maps. First, in the case of the WDD result, a similar, though slightly smaller, range is obtained as in the focused-probe case. In that context, the reduction can be related to the loss of resolution, and thus to less strongly peaked atomic sites. The average potential in the scan window is also greater. This should nevertheless serve to highlight that the retrieved DC component, i.e. the mean phase shift in the reconstruction window, is arbitrary and only depends, numerically, on the size of the reconstruction window and on the sampling of specimen features, including beyond the scanned area.

The SBI-D result, on the other hand, shows a drastically lower range of values in comparison to Figure 4, nearly four orders of magnitude down. Interestingly, whereas ${\mu }^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$ is unaffected by this problem, as mentioned above, it is not the case for the transmission function itself. To understand this, an important difference between the two analytical ptychography methods should be highlighted again, which is that WDD measures the projected potential in an indirect manner, i.e. by extracting the angle of the initially retrieved $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$, post-use of equation (14). As such, it is determined by the ratio between its real and imaginary parts, irrespective of the amplitude. A comparison of the $T^{WDD}\left({\displaystyle \overrightarrow{r}}\right)$ map obtained in the conventional focused-probe case, as shown in Figures 8c and 8d, to the one retrieved in the overfocused case, in Figures 8a and 8b, is sufficient to confirm the role of this indirect measurement process in avoiding a similar defocus-induced value range issue as met for the SBI-D calculation. In particular, in the overfocused case, the amplitude of the measured transmission function is found to possess values more than an order of magnitude smaller than the focused-probe reconstructions, similarly to the SBI-D case though not as strongly, while the ratio of real and imaginary parts remains roughly the same.

As a supplementary note here, the slice-wise transmission functions actually used for the simulation, one set for each frozen phonon configuration, are all phase objects in the strict sense, i.e. with a constant unitary amplitude. This reflects the absence of absorption effects for the interacting electrons, which is assumed to be fully elastic. The measurement by ptychography, on the other hand, should in general not be expected to fulfill this condition, as mentioned in Section 2.1.

3.3 Role of the numerical aperture

3.3.1 Simulation and processing parameters

In analytical ptychography, the range of accessible frequencies is, outside of super-resolution [74,98,143–146], determined strictly by q_A = sin(α)/λ. In this subsection, an interest is thus taken in how the numerical aperture sin(α) affects the dose requirement of the reconstruction. Consequently, two supplementary simulations were performed given α = 15 mrad, in a scan grid of 32 by 32 points, and α = 60 mrad, with 128 by 128 scan points. For reference, the resulting unaberrated probes possess Rayleigh criterions δr_Rayleigh of 198 and 49 pm, respectively. Other than that, simulation parameters were identical as those described in Section 2.1.

Under those illumination conditions, the relation $\sqrt{N_s}\propto q_A$, with N_s the total number of scan positions, is fulfilled, which leads to approximately the same ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ values in all tested focused-probe cases, including in Section 2.1. The accessible frequency range is however twice smaller in the 15 mrad case, and twice larger in the 60 mrad one. Moreover, for both cases, the same selection of $N_{e^{-}}$ values was used as for the conventional focused-probe case, hence leading to comparable count sparsity in the exploited CBED patterns. The resulting doses nevertheless differ due to the change in the number of scan points.

Here, it should furthermore be noted that, for most instruments, using a semi-convergence angle of 60 mrad is either not technically possible or leads to an excessive loss of coherence due to chromatic aberration [203]. In this publication, the use of such large numerical aperture should thus be regarded as relevant for theoretical verification rather than an immediate experimental horizon, although some work has already been performed in that direction within the last few years [204–207].

3.3.2 Noise level in the micrographs

The results of applying the iCoM, SBI-D and WDD methods to the 15 mrad simulation are displayed in Figure 11, and those of the 60 mrad simulation in Figure 12. Owing to the different values of α, the indicated q_A differs among the two cases. An immediate consequence of the reduced frequency surface, for the 15 mrad case, is a resolution insufficient to clearly separate two neighboring atomic sites. As such, in Figure 11, the crystal lattice is visible only thanks to its hexagonal structure, i.e. the cavity in the middle of an given hexagon can be resolved, and no more than half of the first order of lattice-induced frequencies is transferred. In comparison, for the 30 mrad simulation presented in Sections 2.1 and 2.2, two orders, i.e. four hexagonal patterns of frequency peaks, were visible in the Fourier transform. Under α = 60 mrad, the resolution is significantly improved, and up to 11 hexagons can be seen.

Continuing, for both new values of α, the specimen frequencies are visible in the Fourier transform, even with excessive noise in real-space, already from the lowest doses introduced. Beyond this, as $N_{e^{-}}$ increases, the observed level of noise and relative strength of specimen frequencies evolve in a rather similar manner among the three focused-probe cases tested in this section, including Figure 4, with the reconstruction being nearly noiseless at $N_{e^{-}}=1024$. The same method-dependent frequency transfer capacities are furthermore observed in each case, in particular with an important presence of low-frequency artefacts in the iCoM result and slightly more persistent high-frequency noise in the WDD micrograph than for SBI, as explained previously.

3.3.3 Fourier ring correlations

Those first qualitative remarks are confirmed by the calculated Fourier ring correlations, shown in Figure 13 for the 15 mrad case and in Figure 14 for the 60 mrad one. The general behaviour described in Section 2.1 is observed for the two newly introduced numerical apertures too. In particular, the ${\displaystyle \overrightarrow{Q}}$-dependent dose-efficiency is higher for ${\displaystyle \overrightarrow{Q}}$-coordinates that are rich in specimen information and lower for the others, where noise reduction is only due to the normalization. Moreover, as is particularly visible in Figure 14, the manner in which FRC^m (k) is calculated for each method leads to an artificially lower dose-efficiency for higher k, owing to the more important weighting of noisy ${\displaystyle \overrightarrow{Q}}$-coordinates in the corresponding R_k range. It should finally be noted that the micrographs shown in Figures 11 and 12  confirm findings from Section 2.1 relating to the value range of the retrieved phase shift maps. Specifically, the WDD range is about twice as high as the SBI-D and iCoM ones. Beyond that, the resolution, in leading to more or less pronounced atomic peaks, contributes as well to this effect.

3.3.4 Dose requirement and reconstructed frequencies

At a fundamental level, the results presented in this subsection confirm, as was noted by equation (30), that the overall dose requirement of a ptychographic reconstruction, to obtain a specific precision determined by CRLB_RS, is proportional to q_A² and more generally to the surface covered by the reconstructed two-dimensional frequency space. This is expected, as a larger frequency surface implies a larger number of pixels to which Poisson noise [150] is propagated from the detector plane.

From a naive standpoint, as long as a well-focused probe is employed and that $\sqrt{N_s}\propto q_A$ is fulfilled, thus permitting the conservation of the same area overlap, this notion also implies that the relation between the average number of electrons per pattern $N_{e^{-}}$ and the noise level is not fundamentally dependent on the numerical aperture. In other words, it can be expected that, irrespective of the resolution, reconstruction can be performed with very low $N_{e^{-}}$, and thus with sparse CBED patterns [105], the dose being then fixed by the number of scan points. As such, count sparsity in itself is not a limitation for the reconstruction of the electrostatic potential by analytical ptychography. Noteworthily, if the other common strategy is adopted, consisting in recording in the defocused geometry [117,118], $N_{e^{-}}$ will need to be increased to match the information content of the focused-probe data. A revised area overlap is then also necessary, as shown in Section 2.2.

4 Imaging of apoferritin particles under high- and low-resolution conditions

4.1 Contrast predictions above vacuum

4.1.1 Difficulties in imaging light matter

In Section 2, interest was taken in the atomically resolved measurement of the projected potential in a monolayer 2D material which, though it often requires an acceleration voltage U below e.g. 80 kV to avoid excessive knock-on displacement of atoms [208,209], remains an experimentally realistic endeavour. On the other hand, the critical dose [25] of many beam-sensitive specimens, e.g. biological matter, is in practice too low to permit high-resolution imaging, unless done through the combination of a large number of images from identical objects, i.e. a single-particle analysis (SPA) [210,211] procedure. In particular, the dose requirement is proportional to the surface of reconstructed frequencies, as was exemplified in Section 2.3. Consequently, electron ptychography performed on viruses and proteins [1,61–64] has so far been focused on retrieving relatively small ranges of frequency components and limited resolutions.

It should also be noted that the amount of electrons needed depends on the imaged specimen itself. Specifically, it depends on the encountered atom types and their scattering cross-section [212,213], which determines the general amount of observed specimen-induced features in the scattering patterns. In other words, the heavier the imaged material is, the stronger the contrast ends up being in the retrieved phase shift map. This factor, as well as the high probability of radiolysis [27,214,215] leading to the dose limitations mentioned above, make biological specimens particularly difficult to investigate in STEM.

4.1.2 Simulation and processing parameters

In order to explore this topic further and, like in the MoS₂ case, empirically determine the dose requirement for the imaging of a macromolecule in the absence of further issues, e.g. the MTF of the camera, scan imperfections or an amorphous ice embedding, new simulations were performed based on an apoferritin particle in vacuum. The chosen acceleration voltage U was 300 keV. Two distinct semi-convergence angles α were furthermore tested, specifically 1.5 and 6.0 mrad, to verify the previously observed trends on the role of the numerical aperture. In both cases, a field of view of 15 nm by 15 nm, with the specimen in the center, was employed. This field of view was filled by 64² scan positions in the 1.5 mrad simulation, and by 256² positions in the 6.0 mrad one.

Importantly, those illumination conditions both permit an area overlap ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ slightly above 82%, when comparing immediately neighboring scan points, and lead to Rayleigh criteria δr_Rayleigh of about 801 and 200 pm, respectively. Moreover, following a suggestion made e.g. in reference [151,216], the probe focus was placed in the middle of the vertical distance covered by the specimen. Given the large depths of focus δz_DOF = λ/(2sin(α/2)²) [121,122] of about 1750 nm, for α = 1.5 mrad, and 109 nm, for 6.0 mrad, this is nevertheless not expected to be critical here. In particular, and also because such a light material is not expected to lead to e.g. strong channeling effects [123], the wave amplitude should remain sufficiently invariant throughout the propagation axis so that the POA can be considered fulfilled in any case.

Like in the previous section, the multislice method [181–183] was used to represent the elastic propagation of the electron wavefunction through matter, while the atomic potentials were parameterized according to reference [184]. As before, the specimen potential was pixelated in the two-dimensional plane such that a maximum scattering vector of up to twice the range actually used could be simulated. The object was furthermore sliced in the manner described in reference [110]. To perform the calculations in a reasonable time, owing to the large size of the simulation window, thermal vibrations were accounted for by multiplying the scattering amplitudes with an isotropic Debye–Waller factor, i.e. the wave was considered to interact coherently with a time-average of the atoms in motion. In general, this approximation may lead to errors for high enough scattering angles [133], e.g. above 40 to 50 mrad, but was not considered to be problematic here, as the extent of ${{\displaystyle \overrightarrow{q}}}_d$-space available was below this limit.

The results of applying reconstructions on the 1.5 mrad simulation are presented in Figure 15, and those of the 6.0 mrad one in Figure 16. Dose-limitation was imposed based on the procedure described in Appendix B, while employing average numbers of electron per pattern $N_{e^{-}}$ equal to 2 ^l , with l ∈ [2, 3, ..., 10]. The resulting doses are indicated in the figures. For all cases, both the measured projected potential map, expressed in V⋅nm, and the square root of its Fourier transform amplitude are displayed. For all cases, the single CBED patterns were normalized to their sum pre-treatment and the SFPA solution, described in Section 2.5 was employed for practical implementation.

4.1.3 Noise level and contrast transfer

As a first remark, specimen frequencies are not as obviously observable in the Fourier transforms as in the MoS₂ case. In particular, there are no lattice-induced frequency peaks with a width dependent on the size of the scan window to be observed, but rather a complex specimen pattern corresponding to this specific projection of the potential. The overall shape of the particle is also easy to notice in real-space, e.g. from $N_{e^{-}}=16$. This is especially true in the iCoM result, where the higher weighting of low-frequencies, high-frequency information having then been reduced, permits an easy detection of the edges [173,172]. This is nevertheless accompanied by a prevalence of low-frequency noise, as explained in the previous section. For the other imaging modes, the two distinct numerical apertures used permit the visibility of a varying degree of details in the inner structure of the specimen.

Furthermore, the CTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$, for SBI, and ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$, for iCoM, lead to clear differences between the different micrographs, as those two methods highlight specific information in the projected potential map. In contrast, in the case of an atomically resolved crystal where specimen frequencies are sparse, as mentioned above, those effects are not as striking. Upon comparing the two analytical ptychography approaches, it can thus be noticed that the presence of ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ leads to an exaggeration of the intermediary frequencies, e.g. close to $\left|\right|{\displaystyle \overrightarrow{Q}}\left|\right|=q_A$, as is especially visible in the high-dose micrographs.

Relating to arguments given in the MoS₂ case, as well as in Section 2.3, such a clear difference between the SBI-D and WDD results indicates that the imaged specimen cannot be strictly defined as a weak scatterer, i.e. ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ does not intrinsically apply while, under the WPOA, it should represent the information content of the scattering data itself and thus occur in all ptychographic imaging modes. This is furthermore confirmed by the value ranges of the WDD phase shift maps themselves, which are found above 0.4 rad for α = 1.5 mrad and 0.7 rad for α = 6.0 mrad, at infinite dose. Such a finding is of particular interest here, as it demonstrates that considering biological specimens as weak phase objects, even with a low value of σ following equation (2), may not be correct in the general case.

Outside of those aspects, the $N_{e^{-}}$-dependent measurement precision is found to be comparable among the two q_A cases, as explained in Section 2.3. This thus implies, following the necessary adaptation of N_s for the conservation of the area overlap, a proportionality between the required dose and the frequency surface being reconstructed, which here extends to q_A² directly as was illustrated by equation (30). The count sparsity of the CBED patterns used is furthermore not a limitation for the reconstruction itself, as expected from previous results. While $N_{e^{-}}$ increases, internal features of the particle become better resolved, thus here providing a direct empirical verification of the dose requirement of specifically targeted structural information.

4.1.4 Fourier ring correlations

Those observations are confirmed by the calculated FRC^m(k), as shown in Figure 17, for the 1.5 mrad case, and in Figure 18, for the 6.0 mrad one. Like in Section 2, no striking differences of dose-efficiency are observed among the three imaging modes used in this work, which again is reflective of the comparable needed dose to reach the best achievable precision for a given frequency component. This furthermore illustrate the preponderant role of the CTF in making high- and low-frequency noise more persistent in the WDD and iCoM micrographs, respectively.

As was highlighted previously as well, the FRC shows much more difficulty in reaching high dose-dependent precision for higher frequencies k in general, which can be related to the higher number of pixels in the corresponding R_k range. Moreover, the absence of frequency peaks, owing to the difference of structure in the Fourier transform, observed in the previous paragraph, can be noted here as well and leads to a more homogeneous variation of the FRC, though a specific frequency response is still visible. Relating to Figures 15 and 16, it is noteworthy that this specific frequency response, upon comparing the dose-limited cases to the infinite dose reconstruction, can arguably be noticed from e.g. $N_{e^{-}}=64$.

4.1.5 Phase shift value range

Finally, going back to the micrographs themselves, it should be noted that the different imaging modes, much like for the MoS₂ simulations, do not lead to the same general value ranges, though this time the effect depends on the semi-convergence angle α as well. This confirms the hypothesis made in Section 2.1 that disagreements in the range of values covered by the retrieved phase shift map, outside of the role of the resolution and of the reconstructed frequency surface, are related to the methods themselves.

4.2 Influence of protective amorphous ice

4.2.1 Simulation and processing parameters

While Section 4.1 was sufficient to provide a general prospect for the imaging of biological objects via analytical ptychography, it left out an important practical aspect of such experiments. When performing TEM imaging on this type of specimens, it is common to first embed them in a relatively thick layer of amorphous ice, in order to permit stability in a vacuum environment [214]. This experimental protocol is known for leading to a so-called structural noise [217] effect in the result. In particular, as the micrograph should represent a vertical projection of the illuminated object [218], the frequency distribution of the ice directly adds up to the image spectrum, potentially creating difficulties of interpretation.

In order to investigate this effect further in the case of analytical ptychography, the simulations described in Section 4.1 were repeated with a specimen consisting of the same apoferritin particle, though this time embedded within a representative amorphous ice layer. The ice and particle ensemble, having a total thickness of about 50 nm, was relaxed via molecular dynamics, as described in reference [110], before performing the actual multislice calculation. The focus point of the probe was placed in the middle of the object, and other simulation parameters were chosen identically to the in-vacuum simulation cases.

4.2.2 Noise level and contrast transfer

Results of applying the imaging methods to the new simulations are displayed in Figure 19, for α = 1.5 mrad, and in Figure 20, for 6.0 mrad. In general, the remarks made in Section 4.1 can be transferred to this second situation as well. In particular, the visualization of the overall specimen structure in real-space is possible from about $N_{e^{-}}=16$ and the CTF of iCoM and SBI lead to very different final images among the methods. The main difference is the presence of the projected potential of the ice, superposed to the contribution of the apoferritin particle and thus leading to the deterministic noise-like effect [217] mentioned in the last paragraph.

An important difference in its overall influence should furthermore be noted between the two semi-convergence angles employed here. Specifically, while the 1.5 mrad case does not show a very striking loss of contrast due to the ice, as can be observed by comparing it to Figure 15, the 6.0 mrad one is affected much more strongly. In particular, in the SBI-D and WDD images, the particle is nearly not visible at all anymore, in contrast to Figure 16 where its inner structure was well-resolved even at relatively low doses.

As can be directly noticed in the Fourier transforms, the frequency spectrum of the amorphous ice, which possesses a ring-like shape owing to its amorphous structure [218,219], is added to the specimen frequencies, thus obstructing them in the resulting image. It is then clear that, while the value of 1.5 mrad is sufficiently small to mostly cut off the affected ${\displaystyle \overrightarrow{Q}}$-coordinates, then found beyond the 2q_A limit, it is not so in the 6.0 mrad case. In this context, due to its large thickness compared to the apoferritin itself, the amorphous ice furthermore ends up dominating the projected potential measurement, and thus preventing a direct interpretation of the micrograph. Noteworthily, the iCoM method, in Figure 20 is the least affected of the three, which is related to its CTF attenuating higher frequency components.

In general, those findings show that, for the imaging of a single ice-embedded biological object, not only does the frequency distribution of the amorphous ice has to be known beforehand [218], but the numerical aperture may need to be adapted as well in order to obtain an interpretable micrograph. That is, unless further post-processing is employed like in the context of SPA [210,211]. It is moreover clear that, if the WPOA does not strictly apply to a protein particle standing in vacuum, it will be the case as well for its ice-embedded version.

4.2.3 Fourier ring correlations

For completeness of the arguments, FRC calculations were performed for both cases and are displayed in Figures 21 and 22, respectively for 1.5 and 6.0 mrad. As expected, the results displayed in Figure 21 do not show obvious differences from those in Figure 17, while the FRC profiles in Figure 22 have completely changed from their in-vacuum version, as displayed in Figure 18. Specifically, at the coordinates k where the ice-induced ring of frequencies is highest, the fine structure of specimen-related information has been largely replaced by a near-homogeneous response.

5 Discussion

5.1 Comparison of frequency transfer capacities and role of the interaction model

5.1.1 Generalizable contrast transfer function for a weak scatterer

The three STEM-based phase retrieval techniques used in this work, encompassing the iCoM approach and the two existing analytical ptychography methods, can be distinguished in how well they recover specimen information at reconstructible spatial frequencies ${\displaystyle \overrightarrow{Q}}$. In particular, they may be attributed contrast transfer functions, denoting a ${\displaystyle \overrightarrow{Q}}$-dependent attenuation of signal-to-noise ratio. In principle, such frequency-wise reductions of the object spectrum can be solved by deconvolving the result with the predicted point-spread function. This is however difficult in the low-dose case, as the concerned ${\displaystyle \overrightarrow{Q}}$-coordinates may then have been brought below the noise level, hence resulting in significant noise amplification upon deconvolution. As a result, the dose-efficiency is then expected to be much worse for frequency components where the CTF has a low value.

A first step towards determining the CTF is the derivation of the intrinsic phase contrast transfer function ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ [104], depicted in Figure 1, occurring when the illuminated specimen is a weak phase object and being then reflective of a sideband-like geometry in the acquired scattering data. Owing to this geometry, the SBI method [103,104] constitutes an optimized approach for the reconstruction of the phase shift map, where parts of the data containing only noise are excluded as explained in Section 2.3.

More generally, when the WPOA is fulfilled, the PCTF ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ is applicable to all methods investigated in this work, which implies the equality of the WDD and SBI results, given a high enough dose. The known form of the PCTF, and the resulting ${\displaystyle \overrightarrow{Q}}$-dependent noise level in SBI [149], furthermore make it possible to establish a noise normalization strategy [106,149], rendering the noise level homogeneous across the full frequency spectrum of the retrieved object.

In the particular case of iCoM imaging, and as explained in Section 2.4, a supplementary frequency weighting is imposed, following the optical transfer function ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ [114,168] shown in Figure 2. In contrast to the PCTF arising in the case of a weak scatterer, this OTF is due to the much simpler measurement method based on the prior calculation of the average momentum transfer at each scan position, and does not represent the information content of the experiment. As such, under the WPOA, both ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ and ${\displaystyle \tilde{\gamma }}\left({\displaystyle \overrightarrow{Q}}\right)$ can be expected to apply to the iCoM result.

5.1.2 Specimen-dependent frequency transfer

Continuing, from the results presented in this work, it should nevertheless be clear that the WPOA is inappropriate in the general case. Specifically, it was shown, for both the monolayer MoS₂ and the apoferritin model-objects, that not only did the SBI and WDD results differed in a clear manner, but also that the retrieved ranges of phase shift exceeded those of weak scatterers. What this then implies is the inapplicability of the derived ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ as a general PCTF in focused-probe ptychography. In this context, it is noteworthy that, due to its process still being based on the assumption of a weak phase object, the SBI method may then remove useful specimen information from the available scattering data rather than just noise, and still imposes ${\displaystyle \tilde{\zeta }}\left({\displaystyle \overrightarrow{Q}}\right)$ as a CTF [106,151].

This is not the case of the WDD method, which only assumes the more general POA, and consists in a complete deconvolution of the four-dimensional Wigner distribution $\Gamma \left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ of the illumination from the scattering data. As such, and in contrast to iCoM and SBI, its process may be expected to lead to no supplementary ${\displaystyle \overrightarrow{Q}}$-wise reduction in signal-to-noise ratio, at least outside of the frequency-dependent CRLB [195,196], then representing the fundamental information content of the data.

At a deeper level, what this means is that, while the WPOA allowed for the extraction of a specimen-independent frequency transfer capacity in ptychography, the POA alone does not permit this simplification. This then implies that the empirical CTF of the WDD method should always be expected to depend on the illuminated specimen and its scattering power [193,194]. In general, more work on the CRLB as a theoretical precision metric will be needed in the future, as it can be used to derive inhomogeneities along ${\displaystyle \overrightarrow{Q}}$-space, for different specific cases. Noteworthily, this encompasses the role of the aberration function which, even when included in the process, still affects the reconstruction [196], as exemplified in Section 2.2.

5.1.3 Further limitations

Continuing on the WDD method, the absence of a truly process-induced frequency attenuation effect may remain true only as long as the assumed interaction model, i.e. the fully coherent POA, is fulfilled. In this condition, equation (12) is correct and may be used as an accurate basis for the treatment of the scattering data. A first practical limit is the partial coherence of the illumination, which imposes an envelope effect in $\Gamma \left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ [112].

Furthermore, in the case where the specimen is too thick to be accurately described as a phase object, but a ptychographic reconstruction is still performed on the basis of a single transmission function, artificial features may then be introduced in the result [108,151,220]. In this context, it is also noteworthy that, even with a thicker specimen, using an optimally focused illumination has been shown to partly alleviate the artificial features mentioned above [151,216], as is also well-known in the case of CoM and DPC [163–167]. In iterative ptychography, another increasingly popular solution is the inclusion of a multislice propagation within the process [221–224], i.e. the use of a more accurate interaction model.

5.1.4 Use of dark field electrons and super-resolution

One more advantage of the WDD method, in comparison to SBI, is its ability to exploit the dark field electrons, which are otherwise neglected under the WPOA, as shown by equations (17) and (20).

Importantly, in focused-probe ptychography, the use of the scattering vectors above q_A is also necessary to achieve super-resolution [74,98,143–146], i.e. the ability to enhance resolution in real-space by accessing spatial frequencies that extend beyond the diffraction limit. The spectrum of the object is then completed outside of the conventional 2q_A range, by exploiting the relation of the intensity scattered outside the primary beam with those initially missing frequencies. In practice, for analytical ptychography, this is done through the stepping out approach [98,113].

Because of its ability for super-resolution, WDD has the potential to access much larger frequency ranges than SBI and iCoM, though this requires a significant amount of electrons to be present in the dark field, which makes this approach very expensive in terms of dose [147]. As the present work focuses on the imaging of beam-sensitive specimens, that topic was left out of it.

5.2 Other aspects of the reconstruction strategy for low-dose imaging

5.2.1 Dose requirement

Overarchingly, this publication verifies that the dose requirement of ptychography is proportional to the frequency surface to be reconstructed. In the case where no super-resolution [74,98,143–146] is sought, this proportionality extends directly to q_A².

As illustrated by equation (30), this furthermore implies that the total number of detected electrons, needed to reach a certain accuracy in the real-space measurement, should be expected to be proportional to the quantity of reconstructed frequency pixels $N_{{\displaystyle \overrightarrow{Q}}}$. What this then means is that, in general, to achieve a pre-defined signal-to-noise ratio in the micrograph, the numerical aperture has to be adapted to the critical dose of the imaged specimen [25,199].

Furthermore, while super-resolution may be of interest for many other applications, it should realistically not be relied on for the low-dose imaging of beam-sensitive objects, as it is based on exploiting the least intense scattering vectors across the far-field. Consequently, enhancing the resolution when imaging such specimens should rather be done by enlarging the numerical aperture itself, as this then represents the least dose-expensive option.

5.2.2 Normalization

An appropriate normalization choice is also important in analytical ptychography, especially in the case where the acquired data is sparse [105], and thus where large changes in the variance of single patterns occur across the scan window. For this purpose, the strategy proposed in reference [149], consisting in dividing the acquisitions by their individual sums pre-treatment, was adopted in this work.

Noteworthily, this pattern-wise approach is equivalent to performing the calculation while adapting the normalization of the electron wavefunction itself to the number of counts in each corresponding pattern. Further investigations on the normalization strategy may otherwise be relevant in the future, which will also need to be correctly accounted for in any estimation of the theoretical measurement precision, e.g. using the CRLB [198].

5.2.3 Use of a Wiener filter

Finally, in the case of the WDD and SBI-D processes, the use of a Wiener filter [139] as a deconvolution method implies the introduction of a parameter ϵ to avoid divisions by zero, as included in equations (13) and (19). Whereas, at infinite dose, this number may be considered as a simple numerical precision term, it in practice needs to be adapted to the noise level in the distributions ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ and $J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)$ to avoid its amplification in the final result. This is however at the cost of accuracy for the reconstruction itself, i.e. the range of phase shift values and the transfer of higher frequencies are affected, as noted in more details e.g. in reference [225].

In this work, a single ϵ = 10⁻⁶ was consistently used in all deconvolutions, which was found sufficient to avoid noise amplification, including in the condition of highest count sparsity, or an unwanted modification of the value range in the retrieved phase shift maps, as verified in the infinite dose case.

High values of ϵ, up to 10⁻³, were tested as well in the WDD reconstruction, with no further reductions in the noise level, but leading to excessive modifications of the result in high-dose conditions. In the case of SBI, comparisons were performed between the deconvolutive form, using the Wiener parameter of 10⁻⁶, and the summative form. This was done for the few lowest doses considered and did not show a better signal-to-noise ratio in the SBI-S result, hence confirming the stability of the deconvolutive process in this case. There too, higher values did not permit an improved noise level compared to the 10⁻⁶ case. For both the WDD and SBI-D processes, ϵ = 10⁻⁷ furthermore led to a slight, but clear, amplification of the noise.

Importantly, the stability of both the WDD and SBI-D calculations, even with very sparse scattering data and given no adaptation of the parameter to the dose, can be related in large part to the scan position-wise normalization strategy chosen here. In particular, it ensures that each treated CBED pattern has a total value of one, hence leading to a reduction in their individual variances [149]. As such, similar value ranges are consistently found in the amplitudes of the distributions ${{\displaystyle \tilde{J}}}_{{\displaystyle \overrightarrow{Q}}}\left({{\displaystyle \overrightarrow{q}}}_d\right)$ and $J_{{\displaystyle \overrightarrow{Q}}}\left({\displaystyle \overrightarrow{R}}\right)$. Other practical choices included the use of orthonormal Fourier transforms. In this context, the precise selection of the Wiener parameter ϵ becomes less critical for the reconstruction, which then permits more reproducible performances for the analytical ptychography procedures.

As a final note on this topic, reaching a correspondence between the ranges of phase shift covered by the SBI and the WDD results, which was an argument used in reference [225], should not be an objective in choosing a value for the Wiener parameter ϵ. As was extensively discussed in Section 4.1 and noted e.g. in reference [108], the two methods are generally expected to lead to different results, since their fundamental assumptions on the specimen differ. On the other hand, once ϵ has been elevated sufficiently above zero to make sure that numerical divergence and noise amplification are avoided even with very sparse data, there is little justification in further increasing it. In fact, the higher this parameter becomes, the farther away the numerical operation goes from an actual inversion, and the less representative of the real interaction the deconvolved Wigner distribution ${\rm Y}\left({\displaystyle \overrightarrow{Q}};{\displaystyle \overrightarrow{R}}\right)$ becomes. As such, the change in values for the retrieved WDD phase shift, induced by an increment of ϵ, should be seen as artefactual rather than a possible validation of the SSB/SBI result.

6 Conclusion

Analytical ptychography methods present several advantages for the imaging of beam-sensitive materials. Specifically, they are direct, fast and relatively easy to implement. Their requirements in terms of computer memory can furthermore be reduced to allow efficient parallelization and GPU implementation, for instance through the scan-frequency partitioning algorithm introduced in this publication. They also do not require a specific choice in reconstruction parameters for particular cases, such as e.g. a coupling of loss and regularization functions, the batch size or an update strength, which could otherwise be needed to obtain a satisfying output. Finally, they permit the treatment of sparse scattering data [105] with no risk of numerical divergence.

Overall, those advantages make analytical ptychography especially relevant for the low-dose investigation of beam-sensitive objects, where the same measurement often has to be repeated multiple times to reach an accurate result, thus creating a need for streamlined acquisition and reconstruction procedures, even encompassing live processing [100,101]. In that manner, the direct form of analytical ptychography would also facilitate the inclusion in a more complex experimental protocol, such as e.g. three-dimensional structure retrieval based on a single-particle analysis [62,210,211]. Moreover, in this context, potential reproducibility issues could be prevented, given the uniqueness of the processing compared to the high variability among iterative algorithms and their parameter sets.

Acknowledgments

The authors acknowledge financial support from the Horizon 2020 research and innovation programme, as well as from the European Research Council within the Horizon Europe innovation funding programme.

Funding

H.L.L.R. and J.V. acknowledge funding from the Horizon 2020 research and innovation programme (European Union), under grant agreement No 101017720 (FET-Proactive EBEAM). M.L.L. and K.M.-C. acknowledge funding from the European Research Council within the Horizon Europe innovation funding programme under Grant Agreement 101118656 (ERC Synergy project 4D-BioSTEM).

Conflicts of interest

The authors declare no conflicts of interest.

Data availability statement

No experimental data was required in this publication.

Author contribution statement

H.L.L.R.: conceptualization, methodology, software, simulation, data treatment, writing. M.L.L.: simulation, draft review. K.M.-C.: draft review, supervision. J.V.: conceptualization, draft review, supervision.

References

Appendix A: redundancy condition and illuminated area overlap

For a successful ptychographic reconstruction, given a well-focused and unaberrated probe, the overlap between successively illuminated regions is usually required to be 70–80%. In particular, this condition leads to the necessary degree of redundancy in the four-dimensional STEM dataset $I_{{{\displaystyle \overrightarrow{r}}}_s}^{det}\left({{\displaystyle \overrightarrow{q}}}_d\right)$, i.e. specific object locations are probed as part of multiple recordings, thus creating common recognizable features among neighboring scan positions and making the correlative measurement of $\mu \left({\displaystyle \overrightarrow{r}}\right)$ possible.

The illumination overlap ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ between two scan positions distant from a vectorial distance $\delta {{\displaystyle \overrightarrow{r}}}_s$ can be determined numerically using a normalized autocorrelation metric, e.g. given by

$${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}=\frac{{{\displaystyle \sum }}_{{{\displaystyle \overrightarrow{r}}}_0}\left|P\left({{\displaystyle \overrightarrow{r}}}_0-\delta {{\displaystyle \overrightarrow{r}}}_s\right){|}^2\right|P\left({{\displaystyle \overrightarrow{r}}}_0\right){|}^2}{{{\displaystyle \sum }}_{{{\displaystyle \overrightarrow{r}}}_0\text{'}}|P\left({{\displaystyle \overrightarrow{r}}}_0\text{'}\right){|}^4}.$$(A.1)

In comparison to simpler approaches [135], the ratio described by equation (A.1) has the advantage of being calculated in two dimensions rather than just one, which can be expected to make it more accurate. ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ is furthermore defined for an arbitrary aberration function or aperture shape, which would be useful e.g. for future work involving phase plates [116,153–155]. In the case of a strongly overfocused probe, it encompasses the influence of near-field propagation on the precise incident intensity distribution in real-space, rather than simply assuming it to be an homogeneous disk. Examples of calculation are provided in Section 2.

Continuing, as long as the probe remains well-focused and if the area overlap ${\beta }_{\delta {{\displaystyle \overrightarrow{r}}}_s}$ is sufficient, it can be assumed that the scanned area is homogeneously illuminated and that the region outside does not receive any electrons [25]. The dose is then simply estimated by

$$D=\frac{N_s{N}_{e^{-}}}{S},$$(A.2)

where N_s s the total number of scan points and S is the surface of the square scan window. Importantly, equation (A.2) assumes a perfect detection probability, such that all electrons sent to the specimen, and thus contributing to the dose D, end up being measured. In practice, this is not necessarily the case, as higher energy thresholds for electron detection may be imposed when the acceleration voltage U s above e.g. 200 kV [226]. This is then done to prevent multiple counting [48,188], which tends to lower the effective DQE of the camera. This effect nevertheless cannot be represented in the presented dose-limitation process, except by correcting the assumed dose value post-calculation. It may also be that the maximum collection angle of the camera is too low to include every strongly scattered electrons, though this is easily prevented by a correct choice of camera length.

Appendix B: inclusion of dose-limitation in simulated scattering patterns

In order to include dose-limitation in the simulations of scattering data while reproducing the single electron sensitivity, and thus the resulting count sparsity, of a hybrid-pixels DED [226], this publication proposes the following approach. First, the user defines an average number $N_{e^{-}}$ of electrons sent on the specimen at a given scan position, i.e. a ${{\displaystyle \overrightarrow{r}}}_s$-wise expectancy of the incident intensity. Each CBED pattern is then attributed a random number of counts $n\left({{\displaystyle \overrightarrow{r}}}_s\right)$ following the Poisson probability

$$p\left[n\left({{\displaystyle \overrightarrow{r}}}_s\right)|N_{e^{-}}\right]=\frac{{N_{e^{-}}}^{n\left({{\displaystyle \overrightarrow{r}}}_s\right)}{e}^{-N_{e^{-}}}}{n\left({{\displaystyle \overrightarrow{r}}}_s\right)!}.$$(B.1)

A dose-limited intensity $I_{{{\displaystyle \overrightarrow{r}}}_s}^{N_{e^{-}}}\left({\displaystyle \overrightarrow{q}}\right)$ is then obtained, for each scan point ${{\displaystyle \overrightarrow{r}}}_s$. This is done through the random selection of a single pixel, with probability weighted by the underlying pre-calculated $I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)$ and repeated $n\left({{\displaystyle \overrightarrow{r}}}_s\right)$ times, thus generating a new count at every step.

In that manner, the simulation results can be made to encompass Poisson statistics [150] in the amount of counts per scan points, while providing a faithful representation of the detection process involved in devices such as e.g. the Medipix3 [30], Timepix3 [32] and Timepix4 [36] chips, which is itself represented by a multinomial distribution [149]. This then constitutes an alternative to the more conventional approach, which would simply consist in adding noise over the simulated CBED patterns, hence with no direct representation of sparsity other than rounding pixel values to the closest integer.

The newly introduced dose-limitation procedure possesses an additional advantage, in that it provides an opportunity to make predictions on multiple counting [40,42,48,188] in hybrid-pixels DED [226], which is normally due to single electrons depositing an amount of energy above the detection threshold in more than one location. As each wavefront collapse on the camera is here represented individually, it becomes in principle possible to model the stochastic travel among pixels, e.g. through a Monte-Carlo calculation informed on the varying velocity of the incident electron [227] and encompassing a choice of threshold energy. This would in turn lead to a more realistic representation of the resulting information spread effect than a direct convolution of $I_{{{\displaystyle \overrightarrow{r}}}_s}\left({\displaystyle \overrightarrow{q}}\right)$ with a known isotropic $M\left({{\displaystyle \overrightarrow{r}}}_d\right)$, either post-pixel selection or pre-noise supplementation. In particular, multiple counting manifests as non-isotropic clusters unique to each incident electrons [41,50,228,229], whose sizes and shapes depend on the acceleration voltage and which maintain a constant value of 1 among activated pixels. Hence, for a ptychographic calculation based on a collection of sparse diffraction patterns such as the ones generated in this work, it can be expected that those subtleties become important. This topic will thus be critical for future work on low-dose ptychography making use of the Timepix3 [32] or the Timepix4 [36] chips.

All Figures