Inverse Multislice Ptychography by Layer-wise Optimisation and Sparse Matrix Decomposition

We propose algorithms based on an optimisation method for inverse multislice ptychography in, e.g. electron microscopy. The multislice method is widely used to model the interaction between relativistic electrons and thick specimens. Since only the intensity of diffraction patterns can be recorded, the challenge in applying inverse multislice ptychography is to uniquely reconstruct the electrostatic potential in each slice up to some ambiguities. In this conceptual study, we show that a unique separation of atomic layers for simulated data is possible when considering a low acceleration voltage. We also introduce an adaptation for estimating the illuminating probe. For the sake of practical application, we finally present slice reconstructions using experimental 4D scanning transmission electron microscopy (STEM) data.


Introduction
One of the fundamental challenges in electron microscopy is dealing with phase retrieval from the intensity of diffraction patterns. The reason for this problem is that current detectors used in electron microscopy are unable to record phase information, which is necessary for example to improve the image resolution, to understand the interaction of electrons and atoms within a material, and in particular to recover the electrostatic potential of the specimen.
Several methodologies and approaches have been developed for solving the phase problem. One of the most prevalent techniques is ptychography. Instead of only exploiting the intensity of a single diffraction pattern, ptychography takes advantage of a large set of subsequently recorded diffraction patterns stemming from multiple, partly overlapping illuminations of the object. Here, the illumination, given by the electron wave incident, is sometimes also called as the probe. In essence, the acquisition of diffracted intensities from adjacent scanning positions provides additional information enabling to solve the phase problem [1][2][3]. Following on from these original approaches, various contributions to the phase retrieval from a single diffraction pattern have led to the introduction of new algorithms, i.e. enhanced methods adapted for ptychographic reconstructions. For instance, adoptions from alternating projectionbased algorithms like classical Gerchberg Saxton (GS) [4] and Fienup Hybrid Input-Output (HIO) [5] are referred to under the term of Ptychographic Iterative Engine (PIE) algorithms [6][7][8]. Another approach for solving the phase problem by direct inversion has been proposed in refs. [9,10]. It utilises the property of the ambiguity function, sometimes also called Wigner function, which naturally appears by reformulating the equation for deriving the intensity in terms of the probe and object transfer functions.
In addition to the aforementioned methods, further approaches for modelling ptychography as an optimisation problem have been developed over the last few years. As the phase retrieval problem is generally non-convex, there is no certainty that the global optimum can be attained. However, several contributions [11][12][13] manage to achieve the convergence to a local optimum. The crucial fundamental assumption for most studies is the single multiplicative approximation used for modelling the interaction between the electron beam and a thin specimen. However, this assumption does not necessarily hold when investigating thick specimens due to strong dynamical electron scattering effects [14,15]. For this purpose, one should take into account the theory of multiple scattering and propagation when solving the phase problem for thick specimens, e.g. via the multislice approach [14,16], Bloch waves [17,18], i.e. scattering matrix-based formulations [19].
Several attempts have been made to adapt the phase retrieval model for thick specimens by incorporating the scattering matrix, as discussed in [19][20][21]. In [19], the authors developed an iterative projection algorithm called N − phaser for estimating the scattering matrix. The key idea stems from the specific eigenvalue structure of the scattering matrix. It can therefore be used for estimating the object transfer function from a thick specimen while eliminating the unwanted scattering artefacts in the recorded diffraction patterns. Another approach has been proposed in [21], where the authors used optimisation methods, e.g. Alternating Direction Method of Multipliers (ADMM) and block coordinate descent, in order to estimate both the scattering matrix as well as the probe. A similar approach to estimate object and probe is presented in [22], where the reconstruction is done iteratively by a modified Gauss-Newton method.

Related work
Implementations of inverse multislice ptychography have been applied for instance in [23][24][25]. The key idea in these studies bears a strong resemblance to extending the established algorithms, such as the extended PIE. As three-dimensional specimens were investigated, these algorithms were named 3PIE. The forward model deals with the propagation of the specimen entrance wave (probe) and calculates the complex wave function for the observed specimen at a specific thickness. The backward model constructs an estimation of the entrance wave by applying an inverse Fourier transform to the product of the estimated phase and the intensity of diffraction patterns acquired by the measurement. However, mentioned works focus on the reconstruction of visible light and x-ray datasets. The same algorithm has been applied to reconstruct images from the LED microscope data in [26]. Another approach are gradient-based methods, where the gradients are calculated over the whole multislice model all at once. Examples can be found in [27][28][29][30].
In order to address the inverse multislice ptychography problem for electron microscopy data sets we present two different approaches, an adaptation from the Amplitude Flow method and a matrix decomposition, respectively. Amplitude Flow is a gradient-based method, which has been analysed for a randomised one-dimensional phase retrieval [11]. This analysis was later enhanced for arbitrary measurements and in particular for ptychography [31]. In the second approach, a matrix factorisation technique adopted from the field of optimisation and dictionary learning [32] is incorporated into the estimation of the matrix from intensity measurements. Apart from proposing different techniques to solve inverse multislice ptychography we have also reformulated the forward multislice model. This adjustment enables separation of the effect of the illuminating probe from interaction with the specimen, which in turn yields only the construction of a thick object transmission function in respect of a single matrix. Additionally, we outline a methodology to reconstruct each atomic plane by applying the proposed algorithm to synthetic data simulated for a low acceleration voltage. However, it should be noted that increasing the number of slices affects the decomposition performance.

Summary of Contributions
• The forward multislice model is reformulated to disentangle the effect of the probe from the object transfer function at any thickness. This allows modelling a thick object as a matrix consisting only of the product of phase gratings and Fresnel propagators.
• Two approaches for estimating the slices of a thick object are proposed, namely layer-wise optimisation and sparse matrix decomposition. In the first approach we cycle over slices applying the Amplitude Flow algorithm, only optimising with respect to a single slice. The matrix which represents the phase gratings and the Fresnel propagations were recovered by applying a second algorithm. Further factorisation of this matrix was then carried out in order to extract the slices of the object.
• Simulations of diffraction data of specimens (MoS 2 , SrTiO 3 and GaAs) with different crystal structures have been carried out. These simulations were performed for different energies of the incoming electrons, i.e. different wavelengths. They serve as the ground truth for determining the error of the reconstructions, which were carried out using the proposed algorithms. Since the depth resolution in an electron microscope is limited, we investigated the impact of the slice thickness and under what conditions a unique reconstruction is possible. We also adapted the algorithm to estimate the illuminating probe and present the probe reconstruction.
• To highlight the practical use of our proposed method we provide results for the first applications of the algorithms in respect of experimental data, notably in our reconstructions of the object transfer function of a MoS 2 specimen using a four-dimensional data set acquired by scanning transmission electron microscopy (STEM).

Notations
Vectors are written in bold small-cap letters x ∈ C L and matrices are written as a bold big-cap letter A ∈ C K×L for a complex field C and for a real field R. A matrix can also be written by indexing its elements The set of integers is written as [N ] := {1, 2, . . . , N } and calligraphic letters are used to define functions A : C → C. Specifically, we denote the discrete two-dimensional Fourier transform by F. For both matrices and vectors, the notation • is used to represent element-wise or Hadamard product. The A H is used to represent conjugate transpose. For a vector x ∈ C L , the p -norm is given by |x |. For a matrix X ∈ C K×L , the Frobenius norm is denoted by X F := K k=1 L =1 |X k | 2 and the spectral norm is given by The trace operator Tr (.) is the operator to sum all elements in the diagonal of square matrices. Indices are wrapped around, so that x −p = x N −p .
2 Problem Statement

Forward Multislice Model
The forward multislice model is based on the idea that a thick object can be approximated by multiple thin slices stacked on top of each other. For each slice of the specimen the interaction between the electron wave incident on this slice and the potential of the slice can be modelled by a multiplicative approximation in analogy to the standard model in ptychography. Furthermore, as the illumination progresses through the object, the exit wave of the previous slice propagates through potential-free space to the subsequent slice, where it acts as the new illumination for this slice, as schematically shown in Figure 1.
Considering an aberration-free probe, the two-dimensional probe P ∈ C N ×N can be described for different aperture sizes q max with entries given by where J 1 is a Bessel function of the first kind of order 1. The intensity of this function is called Airy disk. In general, this function can be derived analytically by applying a two-dimensional inverse Fourier transform to a circular aperture. The probe shifted to the scanning position (x s , y s ), s ∈ [S] is denoted by a matrix P s ∈ C N ×N with entries p s x,y = p x−xs,y−ys . In general aberrations exist and affect the probe formation. In this study the focus is on the aberration-free condition when generating the simulated data. For a more general treatment of this subject please refer to [15].
The interaction between the probe P s ∈ C N ×N at scanning point s ∈ [S] and the first slice X 1 ∈ C N ×N is given by the element-wise product and in turn produces an exit wave of slice 1 After passing through the first slice the propagation of the exit wave between the slices is modelled by the Fresnel transform V z which is given by The parameters q y , q x denote the discrete grid in the reciprocal space and hence represent spatial frequencies, ∆ m is the distance of the wave propagation, and θ x , θ y are the two-dimensional tilt angles. In this article, the illumination was set to be perpendicular to the object surface along a major crystallographic axis, i.e. tilt angles are zero. As the beam reaches the second slice, it is described by V z (E s 1 ) and the next exit wave is given by Consequently, the general representation of the m-th observed exit wave is written as , m = 1, Finally, the intensity of the Fraunhofer diffraction pattern that is recorded by a detector in the far field is given by In scanning transmission electron microscopy (STEM) the illuminating probe is rastered across the specimen. Therefore, the illumination is varied to yield a set of S diffraction pattern intensities collected throughout an experiment. This four-dimensional data set is then subjected to phase retrieval by multislice ptychography.

Reformulation of Multislice Ptychography
The measurement model in (2) can be further reformulated in order to separate the probe-and the object-related terms. This reformulation is based on the following property of the Hadamard product.
The notation vec : C N ×N → C N 2 is an operator that vectorises the matrix and diag : C N 2 → C N 2 ×N 2 constructs a diagonal matrix by placing the elements of the given vector on the main diagonal. Additionally, the second equality in (3) is valid from the commutative property of the Hadamard product. By using (3) the first exit wave for s-th position of the probe can be rewritten as vec (E s 1 ) = vec (P s • X 1 ) = diag (vec (X 1 )) vec (P s ) .
For convenience, the following notations are introduced and As a result, the first exit wave can be simplified as For the second exit wave the action of the Fresnel propagator is required. Again, by facilitating (3), it is given by where matrices F 2D , F −1 2D ∈ C N 2 ×N 2 are two-dimensional Fourier and inverse Fourier matrices, respectively. The term one-dimensional Fourier matrix represents the discrete implementation of the Fourier basis, i.e. when the Fourier basis is sampled and stored as a matrix. Along the same line, the two-dimensional Fourier matrix can be constructed by using the Kronecker product between two one-dimensional Fourier matrices.
Hence, the Fresnel propagator is a multiplication of the vectorised exit wave E s m with a matrix Substituting the obtained representation to the second exit wave results in Consequently, the M -th exit wave is given by and the resulting far-field vectorised intensity is By combining all vectorised intensities i s and probes p s as columns of the matrices respectively, can be simplified to There are at least two benefits of this reformulation. Firstly, the object transfer function for an arbitrary thickness M is now represented by the matrix A M ∈ C N 2 ×N 2 . Note that A M purely represents the properties of the thick specimen without being affected by the probe, in opposition to the model (2), where the probe's illumination is entangled with the slices. Secondly, the matrix A M decomposes into the product according to (4) and each slice O m can therefore be separated from other multipliers, i.e. Fresnel propagator G m , which will be convenient in the next section where the recovery of a thick object is discussed. With this reformulation the matrices now have an ambient dimension of N 2 × N 2 which increases the computational complexity for processing data of this form. However, because most of these matrices are diagonal matrices they may therefore allow a more efficient treatment and storage in comparison to, e.g. the Bloch wave method, in which an eigenvalue decomposition is directly performed on a scattering matrix of dimension N 2 × N 2 .

Methods and Algorithms
This section considers the recovery of a thick object and the probe from intensity measurements (5) in diffraction space. Firstly it is posed as a constrained optimisation problem, then two algorithms are proposed for solving the optimisation problem under the assumption that the probes are known. Finally concepts for incorporating the probe estimation into the suggested methods are provided.

Inverse Multislice Ptychography as an optimisation problem
One of the standard approaches for the recovery of the object from intensity measurements is the optimisation of the data fidelity, which is represented by the least squares problem The challenge in obtaining the minimiser of (6) is the non-convexity of the objective function, which results from the absolute value and product representation of the matrix A M as well as the multiplication with the probes P. In general, non-convex functions are known to require non-polynomial time to find the global optima [33]. One common method for tackling a non-convex minimisation is the alternating projections method [34,35]. It is based on minimising the objective function with respect to a single selected unknown at a time, while other unknowns remain fixed. This usually results in simpler intermediate subproblems which can be solved efficiently. Afterwards, the next unknown is chosen for optimisation and this process is continued until the minimisation with respect to any of the variables does not improve further. Whereas the alternating minimisation is well-understood for convex functions [34], it often acts as a heuristic for non-convex ones. Nevertheless, in applications such as ptychography [36], the estimates obtained by alternating minimisation are quite accurate, which motivates applying this technique to inverse multislice ptychography.

Phase retrieval via Amplitude Flow
As observed throughout this section, the minimisation of the objective function in (6) with respect to a single unknown, either O m , A M or P, leads to the phase retrieval problem. It concerns the recovery of an unknown vector z ∈ C L from the measurements of the form with the measurement matrix Q ∈ C K×L . One popular approach for the reconstruction of z is the Amplitude Flow algorithm [11,31]. It applies the gradient descent in order to minimise the least squares objective The generalised Wirtinger gradient of the function A is given by where each element k ∈ [K] in the fraction Qz |Qz| k is set to 0 whenever (Qz) k = 0. Then, starting from a position z 0 , the t-th iteration is obtained via the gradient step with learning rate µ = Q −2 given by the squared inverse of the spectral norm of the matrix Q. This approach with the chosen learning rate µ was adopted such that the convergence of the algorithm to the critical point of the objective function (7) can be guaranteed [31].

Layer-wise Optimisation
For the layer-wise optimisation, the alternating minimisation was adopted so as to optimise (6) with respect to a single slice O m at a time, while keeping all other slices unchanged. Note that the objective in the optimisation problem (6) can be understood as a sum of errors for each scanning point s, Using initial guesses of the object transfer functions O 0 1 , . . . , O 0 M , the alternating minimisation technique is employed in order to optimise with respect to a single slice O , ∈ [M ] by solving where the supporting prefix and suffix matrices are given by Once an estimate for the -th slice is produced, the algorithm continues with the + 1-th slice. After the L-th slice is estimated, the estimation process is repeated from the first slice until a desired stopping criterion is reached. By applying (3), the intensity measurement for a single probe p s , s ∈ [S] can be rearranged according to The optimisation problem (10) is equivalent to the phase retrieval problem with the measurement matrix and the measurements given by respectively. The problem in (13) can be solved by running the gradient descent method as discussed in Section 3.2 for a fixed number of iterations.
Overall, it grants us the Algorithm 1 summarised below.

Algorithm 1 Layer-wise Estimation
• Number of iterations T and number of gradient steps K. Compute prefix and suffix matricesR t andŜ t via (12).

Sparse Matrix Decomposition
Another approach for solving the optimisation problem (6) is separating it into two subproblems. At first, the matrix A M ∈ C N 2 ×S , which represents the complex object transfer function of a thick specimen is estimated from the measurements. In the second step, this estimated matrix A M is to be decomposed in order to determine the slices O 0 m , m ∈ [M ], which is sparse in the sense that only diagonal elements are non-zero. Both steps are then reiterated. The detailed procedure is described below.
In the first step, the object transfer function is estimated by solving the optimisation problem In case the matrix A is vectorised, (15) can be treated as a phase retrieval problem of the form (7), which gives the gradient step with O 0 m , m ∈ [M ] being the intialisations for each slice. Once matrix A is estimated usingÂ, its sparse decomposition [32] is achieved by solving the following problem where the objective is reformulated in terms of the prefix and suffix (12) matrices. Thereby, proximal gradient descent methods [37] can be applied, which grants an update of the form with the projection operator P acting onto the space of diagonally normalised matrices given by and the learning rate µ = 1 c where c ≥ (λ t R t S t ) 2 , as discussed in [32].
For the minimisation with respect to λ, estimates of the M -th slice are combined according tõ and λ is updated by minimising it according to the one-parameter least squares problem Therefore, the update for λ t+1 is given by which concludes the second step of our method. These two steps are reiterated by using the new initialisation for the first step The summary of the procedure is given in Algorithm 2. Estimate matrixÂ with initial A t by solving (15).

5:
Compute prefix matrix : R t

6:
Compute suffix matrix : S t

Probe Reconstruction
After estimating the object i.e. the phase gratings of a specimen the optimisation method can be adapted related to Amplitude Flow in (7) in order to estimate the centered probe p c ∈ C N 2 , by utilising the intensity of diffraction patterns at the same position, i.e. i c ∈ R N 2 , Additionally, the gradient update for the t-th iteration is similar to (8), where there is In this case, the learning rate µ is calculated by using the spectral norm of the estimated matrix F 2DÂ , i.e. F 2DÂ .

Simulation and experimental details
In this section, information regarding the simulated dataset of a specimen used as the ground truth, including its type and crystal structure, is given. Furthermore, the microscope and a description of the experimental conditions used for obtaining actual experimental diffraction data are provided.

Simulated data sets
Intensities of simulated diffraction patterns from Molybdenum disulfide (MoS 2 ), Strontium titanate (SrTiO 3 ), and Gallium arsenide (GaAs) specimens with elevated thicknesses were generated by using a forward multislice algorithm. In Figure 2, their 3D structural representations as well as 2D projections along [0 0 1] of the unit cells are shown. The structural as well as simulation parameters are given in Table 1. In addition, the parameter unit cell presents the most simple repeated lattice point in the crystal. The collection of several unit cells is called a supercell. At last, the semi-convergence angle represents the semi-angle that appears in a cone shape when a convergent electron beam illuminates a specimen.

Experimental dataset
Besides simulated datasets, numerical evaluations of experimental datasets were also performed. From a bulk crystal of 2H-MoS 2 , sheets were exfoliated by using a poly-dimethylsiloxane elastomeric film supported on a glass slide and transferred onto a holey silicon nitride membrane for the use in transmission electron microscopy (TEM). Experimental data of MoS 2 was acquired using a probe corrected Hitachi HF5000 field emission microscope in STEM mode and with an acceleration voltage of 200 keV as well as a beam current of about 7.4 pA. Intensities of diffraction patterns were recorded by using a Medipix3 Merlin4EM camera with 256 × 256 pixels. The distance between neighbouring scan points was 26.5 pm in x, i.e. horizontal, and y, i.e. vertical, scanning directions. In addition, the acquisition time per diffraction pattern was 0.5 ms and data was acquired using a dynamic range of 6 bit. The positionaveraged convergent beam electron diffraction pattern (PACBED) is depicted in Figure 3, where the intensity of all diffraction patterns from 128 × 128 scanning points is averaged.
(a) (b) Figure 3: (a) PACBED of an experimental data set of MoS 2 acquired using a Hitachi HF5000 microscope in STEM mode with a Medipix3 Merlin camera, (b) Amplitude of the probe initialisation, or the socalled Airy disk, which is generated by taking the absolute value of the two-dimensional inverse Fourier transform of the circular aperture generated from the PACBED.

Numerical Results
Several numerical evaluations that measure performance of the proposed algorithms are presented in this section. The object O 0 m and the probe are initialised by using an identity matrix and an Airy disk previously outlined in Section 2. Initially the error metric used for measuring the quality of a reconstruction will be defined.

Error metrics for Evaluation of the Algorithm
The error metric used for evaluating the reconstruction of each slice is calculated as the mean square error with O m andÔ m being the ground truth and the estimated object at the m-th slice, respectively. The objective of the presented optimisation problem is to minimise the error between measured and estimated intensities of diffraction patterns. Accordingly, it is necessary to introduce an additional error metric withÂ M being the total estimated object transfer function at slice M . In (21) the error metric is referred to as the relative reconstruction error whereas, the error metric in (22) is referred to as the relative measurement error. Two settings, the reconstruction of an arbitrary thickness and the reconstruction for decomposing into the different atomic planes, are evaluated. Both differ by the conducted Fresnel propagation distance. Using the latter a reconstruction of the phase grating for each slice with the actual Fresnel propagation distance was attempted, i.e. the thickness resulting from the crystal structure of the specimen.

Reconstruction of Arbitrary Slice Thickness
The intensities of diffraction patterns of GaAs, MoS 2 and SrTiO 3 specimens were generated by using the forward multislice method for a thickness of 20 nm. The simulation parameters are given in Table 1. In this case, the reconstructions of 1 nm and 4 nm Fresnel propagation distance in the inversion process are evaluated. It should be noted that, depending on the thickness, the reconstructions will accumulate all the atom positions of each phase grating into one slice. Thereby, all atom positions are projected onto one image. This setting is necessary to evaluate experimental data as the Fresnel propagation distance of the specimen cannot be exactly determined at the atomic scale. Thus, one can heuristically approximate the correct reconstruction by using an initial, comparably large slice thickness and evaluating the atomic positions.
In Figure 4 the phase reconstructions of five slices with a Fresnel propagation distance of 4 nm are shown for the layer-wise optimisation as well as the sparse matrix decomposition. The sparse matrix decomposition yielded a higher range of the phase reconstruction and in comparison outperformed the layer-wise optimisation in terms of the accumulated atom positions of the MoS 2 and SrTiO 3 specimens. Due to the large slice thickness the reconstructed phases differed significantly from the ground truth GT, in particular when using the sparse matrix method, although the average structure of each slice was correctly reconstructed. It was observed that a 20 nm thick specimen leads to very strong multiple scattering effects.
Numerical evaluations with 1 nm slice thickness, i.e. Fresnel propagation distance, have been performed as depicted in Figure 5. Similar to the 4 nm case, except for MoS 2 where the phase reconstruction appears unstable after the second slice, the sparse matrix decomposition generally performed better than the layer-wise optimisation. This confirmed that adding another constraint to impose the solution as a diagonal matrix can significantly improve the reconstruction. The reconstruction performed best for the GaAs and SrTiO 3 cases, this could be explained due to their slice thicknesses being closer to an integer multiple of the lattice parameter in comparison to MoS 2 . This caused a slight beating effect in dependence of the thickness. In these evaluations knowledge of the correct lattice parameter from Table 1 was intentionally not used as a prior to verify the outcome for the realistic case where the structure of the investigated material is also unknown.

Atomic Plane Decomposition
The reconstruction of each atomic plane, i.e. direct recontruction of each slice in the inverse process at the Fresnel propagation distance similar to the forward multislice model described in Table 1, was also examined. This approach proved quite challenging due to ambiguities concerning the atom positions along the direction of the electron beam which occurred in the reconstruction. The relative measurement error is depicted in Figure 6.
In terms of the error between the ground truth and the reconstructed intensity of the diffraction patterns it was observed that the sparse matrix decomposition converged faster than seen in the layerwise optimisation. The reconstructions of each slice for both algorithms are presented in Figure 7. While attempting to uniquely decompose each slice in the atomic plane, at the original Fresnel propagation  therefore is simply the product of each atomic plane. Consequently it was difficult to uniquely decompose each slice. Since the intensity of the diffraction patterns and the depth resolution depend highly on the convergence angle of the electron beam, they also affect the reconstruction. Conceptual approaches for improving the phase reconstruction in order to resolve each atomic plane are provided below.

Increasing the Fresnel propagation distance
A huge advantage of conducting a conceptual study is the ability to set the Fresnel propagation distance in the forward multislice model to any desired value within the simulation in order to improve the convergence of the sparse matrix decomposition, as presented in Figure 8. The relative reconstruction error, as in (21), was evaluated for each slice. Compared to an arbitrary thickness that also determines the Fresnel propagation distance in the inversion process, the Fresnel propagation distance for both the forward and the inverse processes were deliberately increased.
Reconstructions with both layer-wise optimisation and sparse matrix decomposition for the first three atomic slices of GaAs, MoS 2 and SrTiO 3 specimens are presented in Figure 9. The reconstructions resulting from the layer-wise optimisation still suffer from ambiguities even after the Fresnel propagation distance was increased. In contrast to the layer-wise optimisation, sparse matrix decomposition could uniquely reconstruct each slice after increasing the Fresnel propagation distance to 2 nm.
The investigation of the capability of slice-wise reconstructions for crystals with larger thicknesses as a function of the Fresnel propagation distance was also conducted. This distance was artificially increased to 1 nm and 2 nm, respectively. From the results in Figure 10 it is apparent that increasing the specimen thickness also affects the performance of the sparse matrix decomposition, since some ambiguities appear despite a Fresnel propagation distance of 2 nm. Apart from ambiguities that occurred in the phase   Figure 9: Slice phase reconstruction, in radian, of Fresnel propagation distance 2 nm, using an acceleration voltage of 200 keV, by layer-wise optimisation and sparse matrix decomposition, observed at 100 iterations, for the following specimens: (a) GaAs, (b) MoS 2 and (c) SrTiO 3 . GT is the ground truth, L is the layer-wise optimisation, and S is the sparse matrix decomposition. retrieval problem, the performance of both algorithms generally depends on a trade-off between the Fresnel propagation distance and the number of slices to be reconstructed.
Increasing the Fresnel propagation distance would correspond to artificially increasing the lattice parameter in electron beam direction. This is only possible in a simulation study. However, the conceptual insight is that the phase fronts of electron waves with 200 keV energy would only change significantly after propagating 2 nm for the algorithms to separate the Fresnel propagation from the interaction with the Coulomb potential of the slices.

Low electron energy
The Fresnel propagator on eq. (1) contains the product of the wavelength and the propagation distance as the governing parameters. Therefore, a realistic reconstruction with atomic layer sensitivity needs to use larger wavelengths if the Fresnel distance is reduced to atomic spacings in the range of 0.1 nm. Furthermore, the same potential in a given specimen leads to a larger phase change of low-energy electrons as compared to high energies due to the effect of the so-called interaction constant on the phase grating. The relation between electron acceleration voltage U and wavelength λ is given by λ = hc √ e 2 U 2 +2eU mc 2 , with c being the speed of light, h the Planck constant, m the electron mass, and e the elementary charge. For 200 keV electrons, the wavelength is approximately 2.5 pm, whereas it increases to 4.18 pm for electrons with 80 keV energy. The interaction constant increases by about 40%. Note that both acceleration voltages, 200 and 80 kV, are common settings in STEM such that atomic resolution can be obtained readily in aberration corrected machines. Figure 11 shows a reconstruction of the atomic planes using an acceleration voltage of 80 keV. Despite the low-signal artefacts related to the location of the atoms in the different slices, the exact location of the atoms in each slice can now be correctly determined. As discussed earlier, it should be noted that a unique decomposition highly depends on the number of slices to be reconstructed, as shown in Figure 12. It can be seen that having more slices affects the reconstruction since the ambiguities related to the atom positions were still present even after observing the reconstruction at 10000 iterations.

Probe Reconstruction
The reconstruction of the probe after estimating the matrix at the M -th slice, i.e. A M by using the optimisation problem in (19), is presented below. As discussed in the forward multislice model, the ground truth up to a global phase factor, which is in general an undefined quantity.

Experimental Data
Numerical evaluations on diffraction intensity measurements acquired experimentally from a MoS 2 specimen are presented below. The thickness of the experimental specimen was determined by comparison with simulated position-averaged diffraction patterns to be approximately 35 nm. Further experimental details are described in Section 4. Applying the same algorithms to experimental instead of simulated data is a crucial aspect to demonstrate the practical usefulness of the methods. However, reconstructing based on real data is also a critical point, because experiments are affected by additional parameters that are difficult or even impossible to include in the algorithmic setups above. For example, the recording is inherently containing Poissonian counting noise, the camera has a modulation transfer function which leads to a blurring of diffraction space features, and the projection system of the microscope can cause geometrical distortions of the diffraction patterns. Furthermore, the scan positions of the STEM probe usually deviate slightly from the ideal regular raster due to instabilities of the scan engine. Due to this a detailed analysis of the performance of layer-wise optimisation and sparse matrix decomposition algorithms in dependence of the experimental conditions are set aside for a future task. Instead this preliminary evaluation focuses on demonstrating the principal applicability by targeting the qualitative reconstruction of the MoS 2 structure using a relatively small number of five slices, similar to the example in Figure 4. In particular, this was necessary due to computational efficiency and the much higher dimensionality of the experimental data as compared to the simulations, i.e. the large number of probe positions and camera pixels.

Slice reconstruction
Consequently, the phases of the individual slice reconstructions in Figure 14 are not expected to quantitatively represent the actual phase gratings on the one hand. On the other hand, they are supposed to resemble the atomic structure of the specimen in the respective slices, taking a large portion of the dynamical scattering into account. Indeed, the atomic structure is consistently visible in all slices, opposite to single-slice models for which evaluations at thicknesses of tens of nanometers are by far out of range. The dynamic range of the phase is comparably low, most probably because more slices would be needed to disentangle the slice potentials and Fresnel propagation between the slices completely. The reconstruction is observed after 50 iterations. Figure 14 shows the reconstructions for both sparse matrix decomposition and layer-wise optimisation, in which both algorithms are able to reconstruct the atom positions of MoS 2 . In addition the object transfer function matrix A was generated by considering the products of the Fresnel propagation matrices and all reconstructed slices. Furthermore, the two-dimensional projection of all atoms can also be easily produced, as presented in Figure 15, where all atom positions of MoS 2 are present in the projection. A coloured overlay of the Molybdenum (Mo) and Sulfur (S) atoms were added to the figure in order to better visualise the reconstructed atomic arrangement.

Probe Reconstruction
A direct implementation of the Amplitude Flow as in (7) was adopted in order to reconstruct the illuminating probe after estimating the object transfer function matrix A. At this setting the focus was solely on the intensity of the diffraction patterns acquired at the center position of the illuminated area on the specimen. The resulting reconstructions of amplitude and phase of the probe are presented in . This shows that the data was taken with a well-focused probe as indicated by the sharp peak in the amplitude and a flat phase except for the noise. Note that the reconstruction of the probe is in general a robust check whether the algorithm and the parameters used for the reconstruction are suitable to separate illumination and specimen. In the present case, one can, therefore, conclude that the large slice thickness did not affect this, because no specimen details are visible in the reconstructed probe.

Discussion
Two algorithms based on optimisation methods for inverse multislice ptychography have been presented, namely sparse matrix decomposition and layer-wise optimisation, which are derived from the reformulation of the forward multislice model. The connections between reformulation of the multislice and other models to represent thick specimens are discussed and the possible direction of future research is outlined.
The numerical observations showed that the type of specimens and the number of slices impinges on the reconstruction performance of the algorithms. Theoretically it would be interesting to examine the fundamental limit of the algorithms with respect to the number of slices needed to disentangle interaction and Fresnel propagation sufficiently well. Moreover, a systematic study addressing the impact of the probe semi-convergence angle, the electron energy, aberrations of the electron-optical system and coherence effects could shed light on the robustness of the presented methodology in respect to the multitude of experimental parameters in real measurements. If one wants to work with low dose data, where Poisson noise is dominant, one has to modifiy the objective function. Poisson maximum likelihood has been used with success in case of single slice ptychography [41] and can be adapted to our layer-wise estimation. For the sparse matrix decomposition, one can even use a Poisson phase retrieval method [42] in the first step, estimating A, without changing the decomposition method at all.
Since the reformulation of the forward multislice model in this article yielded purely a matrix representing the transfer function of a thick specimen, it would be possible to relate such a matrix to a scattering matrix constructed from the Bloch wave method and observe the differences between both approaches. The latter requires an intensive computational effort of eigenvalue decomposition for huge scattering matrices.
Apart from the comparison between the proposed algorithms to the eigenvalue decomposition with the Bloch wave method, it should be possible to estimate the specimen thickness directly from the algorithms. One possibility could be to incorporate information from the high-angle intensity of diffraction patterns, or to start from a coarse slicing first with large slice thicknesses, and then increase the number of slices subsequently. In case a sufficiently high total thickness is assumed, empty slices should emerge, indicating that the specimen is actually compact along the electron beam direction. In general, a suitable regularisation should be developed and applied in future works. An interesting approach is a suitable sparsity model, as applied in [43] for the case of single slice ptychography.
Finally, the reformulation of the multislice scheme as a simple, though large, one-step matrix multiplication circumvents the successive forward and backward Fourier transform, which characterises the conventional multislice implementations that compute both the Fresnel propagation and its interaction with slice potentials in real space. In that respect, studying the capabilities and performance of the reformulation is not only relevant for solving inverse problems, but also interesting with respect to conventional forward simulations.

Conclusion and Summary
We proposed reformulation of the forward multislice method such that the transfer function of a thick specimen can be directly determined. In combination with the ptychographic approach we presented two optimisation models for solving the inverse multislice ptychography problem for both arbitrary thickness and atomic plane decomposition. In the first case, given the intensity of diffraction patterns, several atomic planes were jointly processed into a single reconstruction in order to show the total potential. In the atomic plane decomposition each slice was reconstructed at its atomic plane and given only the intensity of the diffraction patterns the results showed the unique atomic positions in each slice.
Although the resulting phase reconstructions by layer-wise optimisation still contained ambiguities, both algorithms could recover the locations of the atoms in the inversion process. Furthermore, this showed that for simulated data the sparse matrix decomposition could reconstruct the atom locations unambiguously for each atomic layer given the intensity of diffraction patterns with low acceleration voltage 80 keV. However, it could be observed that the reconstruction using both algorithms were dependent highly on the number of estimated slices. After the thickness of the specimen was increased in terms of the number of slices it resulted in ambiguities of the locations of the atoms. This indicated that the algorithm failed to accurately reconstruct each atomic layer in respect of the different slices.
We also supported our numerical observations with the reconstruction of the object given the intensity of diffraction patterns acquired from the experimental data set of MoS 2 . It was shown that both algorithms can reconstruct the phase of the specimen. Additionally, by using the reformulation of the forward multislice method, a matrix was constructed that represented the thick object transfer function, i.e. scattering matrix. This reformulation can be used to directly generate the two-dimensional projection of the atom arrangement.