A minimalist approach to 3D photoemission orbital tomography: algorithms and data requirements

Photoemission orbital tomography provides direct access from laboratory measurements to the real-space molecular orbitals of well-ordered organic semiconductor layers. Specifically, the application of phase retrieval algorithms to photon-energy- and angle-resolved photoemission data enables the direct reconstruction of full 3D molecular orbitals without the need for simulations using density functional theory or the like. However, until now this procedure has remained challenging due to the need for densely-sampled, well-calibrated 3D photoemission patterns. Here, we present an iterative projection algorithm that completely eliminates this challenge: for the benchmark case of the pentacene frontier orbitals, we demonstrate the reconstruction of the full orbital based on a dataset containing only four simulated photoemission momentum measurements. We discuss the algorithm performance, sampling requirements with respect to the photon energy, optimal measurement strategies, and the accuracy of orbital images that can be achieved.


Introduction
Image reconstruction based upon sparse Fourier space measurements is a problem that recurs in a wide range of imaging modalities, with well known examples being X-ray computed tomography [1,2] and magnetic resonance imaging [3].In fields such as coherent diffractive imaging [4,5] and photoemission orbital tomography [6], the challenge of image reconstruction is further increased due to the inherently incomplete measurement data: whereas image reconstruction requires the complex-valued field to be known, only the amplitude can be measured.For the latter two, iterative image reconstruction methods nowadays enable highresolution imaging with only limited prior information [7,8], but these methods still require strongly oversampled data to achieve this.In many cases, particularly for photoemission orbital tomography which concerns itself with the quantum-mechanical electronic orbitals of molecules, the object of interest is intrinsically sparse, and this can be efficiently exploited using the ideas of compressive sensing surveyed in [9].
Photoemission orbital tomography (POT) is well established as a probe of the electronic structure of thin molecular films.Until recently, the conventional approach of this method to achieve information on the molecular orbitals would be the comparison of angle-resolved photoemission spectroscopy (ARPES) data to simulations generated using density functional theory (DFT) [10][11][12].Iterative phase retrieval, however, allows one to reconstruct the orbitals directly from the measurements [6,8], and this is independent of DFT or similar calculations (see Fig. 1).While spatially resolved access to the electron density of molecular orbitals can also be achieved, e.g., by scanning tunneling microscopy [13,14], the unique nature of POT provides several advantages to probe the molecular orbitals at high resolution.Most notably, together with gas-phase molecular orbital tomography [15,16], POT using photon-energy-dependent ARPES measurements is uniquely suited to image the full threedimensional (3D) molecular orbitals, and POT remains so far the only technique that can do so for larger, nanometer-scale molecules [17,18].Such 3D-POT is highly promising, as the combination with time-resolved photoemission orbital tomography [19][20][21][22] will access spatiotemporal properties of the dynamic electronic wavefunction that are so far inaccessible by experimental means.Additionally, 3D-POT promises a crucial insight at hybrid interfaces, where strong modifications of the electronic structure are known to occur [23,24].
However, despite the rapid development of powerful multi-dimensional photoelectron detection schemes [25][26][27], their application for 3D-POT has remained limited for a number of reasons: First and foremost, the application of phase retrieval algorithms requires a highresolution 3D momentum distribution, implying the time-costly measurement of an ARPES dataset with a large number of photon energies, where, furthermore, the photon flux must be accurately calibrated.Moreover, the plane-wave model of photoemission, upon which phase retrieval methods are based, is only an approximation with respect to the overall photon-energy dependence of the photoelectron intensity [28,29].In this article, we deal with these challenges using a fundamentally new approach for 3D-POT: by implementing ideas from sparse optimization [30], we perform an image reconstruction on a pixelized object that is defined by the combination of several qualitative constraints in addition to sparsity constraints.
In Section 2, we review the fundamentals of three-dimensional photoemission orbital tomography, and identify the challenges that must be solved to facilitate an efficient access to the 3D molecular orbital.We translate this theoretical description in Section 3 to an image reconstruction problem, and explain how this nonconvex problem can be solved by a cyclic projection algorithm with several qualitative constraints.The orbital reconstruction from simulated data for an isolated pentacene molecule is demonstrated in Section 4. After establishing the convergence behavior of cyclic projections in this setting, we focus on two questions: first, how many photon energies must be measured for accurate orbital reconstruction, and second, we investigate how experimental noise such as counting statistics and inaccurately estimated light intensity calibration affect the reconstruction.

Theory of 3D photoemission orbital tomography
In image-reconstruction photoemission orbital tomography, we consider an ordered molecular layer in which the unit cell consists of a single molecule (e.g., a monolayer of pentacene on Ag(110)), and use the plane-wave model of photoemission to calculate the ARPES signal.In this case, for a given molecular orbital ψ the ARPES signal I for the photoelectron momentum ⃗ k can be expressed as Here ⃗ A is the vector potential of the ionizing radiation, ψ := Fψ, the Fourier transform of the molecular orbital, and | ψ( ⃗ k)| denotes the absolute value of ψ at the point ⃗ k in momentum space.The Dirac δ function stems from Fermi's golden rule and ensures energy conservation, taking into account the ionization potential E b of the molecular orbital, the final state kinetic energy E kin = ℏ 2 k 2 /2m e and the photon energy ℏω.
Note that E kin is proportional to the magnitude squared of the momentum ⃗ k; hence, a single ARPES measurement provides the amplitude of the orbital's momentum distribution along a hemispherical shell with a radius that depends directly on the photon energy ℏω.In order to probe the full 3D momentum distribution | ψ( ⃗ k)| of the molecular orbital, it is therefore necessary to measure ARPES patterns for multiple photon energies.This naturally leads to a question that is central to the work that we present here: How many photon energies are required to yield a reliable orbital reconstruction from ARPES momentum maps?Before it is possible to address this question, however, first a strategy must be established that enables the reconstruction of a 3D molecular orbital from a sufficiently large dataset.
A reconstruction strategy for 3D photoemission orbital tomography (3D-POT) must solve two problems: 1) Starting with an ARPES dataset that only contains discrete cuts through the 3D photoelectron momentum distribution, the complete 3D photoelectron momentum distribution must be determined, and 2) the phase distribution corresponding to the measured amplitudes must be resolved.Up until now, these steps have been performed sequentially: first, the 3D momentum distribution is completed using an interpolation algorithm [17].Subsequently, the full 3D momentum distribution can then be used as input of an iterative phase retrieval algorithm, which is the standard approach to determine the phase pattern in 2D orbital imaging as well.The interpolation step has significant implications: here, it is (implicitly) assumed that the momentum distribution is of a shape, or regularity, that can be interpolated and furthermore that it is properly sampled.Such Fourier-space interpolation, however, is prone to errors, particularly when also the intensity-only nature of the measurement is taken in account.Consider for example a zero crossing in the momentum distribution ψ.Here, intensity measurements at opposite sides will both yield Figure 1: Photoemission orbital tomography provides access to the three dimensional structure of molecular orbitals (top left) by combining photon-energy-and angle-resolved photoemission spectroscopy (ARPES) with image reconstruction methods.Here, the photoelectron distribution of an ordered molecular layer measured at different photon energies ℏω (top right) corresponds to the amplitude squared of hemispherical cuts through the three dimensional Fourier space (center).Based on the sparsely sampled Fourier space (bottom right) the initial molecular orbital can be reconstructed using a cyclic projection algorithm (bottom left).non-zero intensity, and interpolation of this data will lead to the false conclusion that the intensity in-between must also be non-zero.The size of the error introduced here decreases as the sampling density increases, meaning that an accurate reconstruction of 3D molecular orbitals through this sequential approach relies on large datasets in which many photon energies have been measured.
In order to reduce the number of ARPES measurements needed and, more importantly, to prevent errors arising from the interpolation procedure, it is thus clear that a completely different reconstruction strategy must be developed.To that end, we now consider the phase retrieval algorithms that are already commonly used in orbital reconstruction photoemission orbital tomography.Generally, these are mathematical optimization algorithms that aim to find a solution that best satisfies a set of constraints.In photoemission orbital tomography, the solution corresponds to the image of a molecular orbital, while the set of constraints consists of the ARPES measurements and the prior knowledge about the molecular orbital.These algorithms are commonly referred to as phase retrieval algorithms because they reconstruct the phase that is lost upon the measurement, but the recovered information can include intensity as well as phase.For example, phase retrieval algorithms in coherent diffractive imaging are already routinely used to also recover the intensity in parts of the diffraction pattern that cannot be measured [31,32].As such, a promising approach to reconstruct 3D molecular orbitals from a discrete set of ARPES measurements is to generalize the already-necessary phase retrieval algorithms to 3D reconstruction algorithms.In this article, we build upon the sparsity-driven phase retrieval algorithm for orbital imaging [8] to develop a 3D-POT reconstruction algorithm.In particular, using the compact and voxelsparse shape of the molecular orbital as a real-space constraint, we present an algorithm that efficiently and accurately interpolates the momentum distribution based upon ARPES measurements at just a few photon energies.
3 Basic approach

From data to phase retrieval problem
To lay the foundation of an image reconstruction method based on the description of the 3D-POT measurement in (1), we first define the mathematical image reconstruction problem that follows from this description.Given that the objective is to reconstruct a 3D image of the molecular orbital, we denote by D ⊂ R 3 the "object domain" of the 3D ARPES data.An example of the object (i.e., a molecular orbital) in the domain D is shown in Fig. 2(a)-(b).Note that throughout this paper, the presented slices of the orbital or reconstructed orbital in the objective domain are truncated for visualization purposes only.We denote by (x i , y j , z l ) the point in D corresponding to the index (i, j, l), for (i, j, l) ∈ I := {1, . . ., N x } × {1, . . ., N y } × {1, . . ., N z }, where N * denotes the number of voxels respectively in the x, y, and z directions.We will denote the total number of voxels |I| by N = N x × N y × N z .The number of voxels and their dimensions are indirectly determined by the measurement resolution and maximum measured photoelectron momentum, respectively.Let D denote the momentum space, i.e., the Fourier domain, where the ARPES data is measured.The momentum space is discretized in the same fashion as the object domain, though perhaps with different numbers of pixels/voxels.We denote by ⃗ k (i,j,l) the point in D corresponding to the index (i, j, l) ∈ I := {1, . . ., N x } × {1, . . ., N y } × {1, . . ., N z }, where N * denotes the number of voxels respectively in the k x , k y , and k z directions in momentum space.We will denote the total number of voxels | I| by N = N x × N y × N z .We take as the starting point of our algorithm a photon-energy-dependent set of momentum maps.The details on how these momentum maps are acquired differ from one experimental facility to the next; for example, the data structure as retrieved using a toroidal analyser [33] can be subtly different [34] to what would be measured using a time-of-flight momentum microscope [19-21, 25, 35].For simplicity and generality, we will therefore assume that we extract from the measurement a two-dimensional intensity distribution I ℏω for the coordinates (k i , k j ) with for (i, j) ∈ I := {1, . . ., N x } × {1, . . ., N y }.Furthermore, we will assume that the measured intensity distribution is already corrected for the polarization factor ⃗ A • ⃗ k (c.f.eq. 1).Each of these momentum maps is a measurement of the intensity on a hemisphere, which can be completed to a sphere by applying the symmetry expected from a real-valued object, i.e. the molecular orbital.The momentum maps are therefore in the data domain D with k z > 0 and with radius r = ∥ ⃗ k∥ depending on the photon energy ℏω used for photoemission; these data spheres are denoted S r ⊂ D. The value of r is determined by the experimental parameters and must be mapped to the voxel grid described above.We index the various measurement spheres with the symbol ρ so that the voxels corresponding to the measurement sphere with radius r ρ are indicated by S ρ .Since a sphere with radius r will run through voxels with an inherent thickness, the measurement sphere will inherit this thickness through the discretization.Complicating matters further is the fact that the measurement device is a planar array, i.e., that k z is not directly measured.All of these effects of discretization and projection onto a planar grid can be accounted for approximately by an appropriate scaling of the measurements.Accounting for this explicitly in the discussion below adds notational clutter that obscures the main ideas, so we will suppress these specifics in our presentation.Instead, interested readers may refer to the related data repository (ref.[36]), which includes the script that was used to cast ARPES measurements to D.
A 3D-POT data set containing momentum maps at m photon energies then provides measurement data on m spheres with radius r ρ for ρ = 1, 2, . . ., m.The measurement spheres taken together are denoted by is the corresponding collection of voxels.The equation ( 1) can be reformulated in the following form, as presented in, for example, [37,38]: where The problem we are faced with is to find ψ that is consistent with the relationship (3), i.e., find ψ when we know only the amplitude of its Fourier transform on finitely many spheres in momentum space.This is a type of phase problem where the data is sparse and restricted to spheres in the Fourier domain.The new challenge is to retrieve not only the phase of the measured data but also the phase and Fourier magnitude on the missing voxels I \ S.Only by accomplishing this can the complete molecular orbital ψ be calculated.

Model
The phase retrieval problem has a rich history and many different approaches to its solution.The algorithm we use here is the most prevalent and successful for phase retrieval and is based on a nonconvex feasibility model [38].As explained above, the 3D-POT data records the momentum distribution of a specific molecular orbital on set of spheres.The first step in the feasibility model approach is conceptual: we view the data as specifying a set of orbitals ψ that could have produced the data; i.e. any ψ consistent with the data is acceptable.We won't belabor the transition in thinking about ψ and ψ as functions taking values on D (respectively D) to their discretized versions as vectors in C N and C N .This allows us to write ψ (i,j,l) instead of ψ( ⃗ k (i,j,l) ); likewise ψ (i,j,l) corresponds to ψ(x i , y j , z l ).The set of points ψ satisfying (3) in this discretized form is thus defined by where S is defined by ( 2).This of course leads to too many possibilities, so we narrow the options down by adding other constraints or prior knowledge about the orbitals.The first constraint concerns the experimental resolution in real space.To ensure that the full measurement data can be included in the analysis, the domain D is chosen to be large enough to fully contain the data for the highest photon energy, which lies on the shell S r with the maximum radius r.Consequently, the corners of D lie outside the measured area in momentum space.As the momentum distribution in the corners is not constrained by the measurement, and an overestimation of the intensity can lead to high-frequency noise on the reconstructed orbital, we set the intensity outside of the shell S r to zero.This reflects simultaneously the experimental limitations as well as bandwidth limitations on the types of orbitals which can be successfully recovered.More will be said about this below when we impose symmetry constraints on the orbital, but imposing such a bandwidth constraint is equivalent to applying a low-pass filter in momentum space at the determined experimental resolution.Let B r be a ball in momentum space D that contains all measured spheres S j (j = 1, . . ., m) in the Fourier domain.The set of orbitals satisfying the bandwidth constraint is given by The following constraints concern the properties of the molecular orbital in the object domain D. Let s > 0. In the object domain D, we utilize a sparse-real constraint which is defined as follows: where, Im(ψ) is the imaginary part of the complex vector ψ, and ∥ • ∥ 0 denotes the counting function that counts all non-zero voxels.We employ the sparse-real constraint since the molecular orbital that we are recovering consists of a small number of positive and negativevalued lobes (e.g., Fig. 2a).Sparsity constraints were first employed in phase retrieval by Marchesini in [39], and have previously been applied successfully in orbital imaging [8].The sparsity constraint imposed by the restriction on the value of the counting function allows us to impose a movable support constraint on the orbitals; this is very convenient since we don't know the molecular orbital initially, and the counting function allows the support to evolve during the reconstruction process.This constraint therefore also eliminates the need for a shrink-wrap procedure to fine-tune an a-priori unknown object support [5,40].
It is not uncommon that the symmetry of a given molecular orbital is known in advance.Based on the properties of the orbital shown in Fig. 2a and Fig. 2b, we introduce additional assumptions, namely symmetry and anti-symmetry constraints.The set of points that are symmetric about the X-axis, anti-symmetric about the Y-axis, and anti-symmetric about the Z-axis is denoted as SYM and is expressed as follows: Imposing a symmetry constraint in the object domain D is made more challenging by the fact that the magnitude of the Fourier transform is invariant under spatial shifts in the object domain.Reconstructions that are consistent with the measured data, and even the sparse-real constraints can appear anywhere in the object domain D (see Fig. 3e).While we might see symmetries in the reconstruction with our eyes, to impose these as a constraint we need to align the object with an axis of symmetry (see Fig. 3d).We accomplish this by simply applying a conventional, fixed support constraint (as opposed to the sparsity support constraint above) centered at the origin that is symmetric in each of the x, y, and z directions separately.We emphasize that the size of this support is chosen small enough to ensure proper alignment of the object, but also large enough that it does not affect the shape of the ultimately reconstructed object.This is shown in Fig. 3b.Given a symmetric binary mask Ω ∈ {0, 1} N , the support constraint set is denoted SUPP and given by We refer to the mask Ω as the support region.By using the support constraint, we can not only avoid spatial translation of the recovered object, but also speed up the recovery process by reducing the number of indices that need to be estimated.
Having defined the qualitative constraint sets above, we can now formulate the feasibility model as finding a point that lies in all these sets: This is an inconsistent problem because the qualitative constraints are idealized and the intersection above is empty.From a physical perspective, the most significant impediment is due to the fact that the true molecular orbital in Fourier space is actually non-zero outside the ball B r .From a numerical perspective, different sparsity parameters s in (7) or support regions Ω in (9) may also cause inconsistency, as they introduce hard edges which are suppressed by the low-pass filter LF.Inconsistent feasibility has been studied in [41] where the interpretation of "solutions" to (10) is explained in terms of difference vectors or gaps of smallest magnitude between successive sets [41, Lemma 3.2]; this is depicted in the case of two sets in Fig. 4. Generally, it is found that the size of the gap correlates with the distance to the 'true' solution, such that the reconstruction with the smallest gap provides the best approximation.

Method
Algorithms based on metric projections are the most successful methods for finding points where the gap between the sets is smallest [38].The Cyclic Projection algorithm (CP) is one of the simplest, and more effective approaches.Since the focus of this work is not on the algorithms, we will limit ourselves to CP; if the reconstructions are successful for cyclic projections, they will also be successful for more sophisticated algorithms.The algorithm CP is given by: Cyclic Projections: Given an initial point ψ (0) ∈ C N , for n = 1, 2, . . ., do Here ψ (n) denotes the reconstructed orbital at the n-th iteration and P C denotes the metric projector onto a set C. The formulas for the individual projectors are given below.We use character "∈" instead of "=" since the projection of a point onto the sets M and SR can be multi-valued.
Luke et al. [41,Theorem 3.2] proved guarantees on local linear convergence of CP to fixed points, i.e. points x * with T CP x * = x * , for non-convex problems and global linear convergence for convex problem under regularity assumptions of the sets.In short, for our non-convex problem, as long as the starting point is close enough to a fixed point of the algorithm, then the iterates of the algorithm converge linearly to a fixed point.
A projection of ψ onto the set M with infinite precision arithmetic is computed by [37, Theorem 4.2] where F −1 is the discrete inverse Fourier transform and i 2 = −1.As argued in [37, Corollary 4.3], this formulation is numerically unstable, so instead we approximate the projection by where v = (v (i,j,l) ) (i,j,l)∈ I is defined by The first case of (13) states that if the triple index (i, j, l) corresponds to voxels that belong to the measured spheres, then the updated magnitude of the Fourier signal must be equal to the measured information.The parameter ϵ is introduced to prevent division by zero.
The choice of ϵ depends on the machine precision.For our examples ψ is a complex vector of length of about 10 6 elements, so with double precision arithmetic, we shouldn't expect better than 10 −9 accuracy for vector-vector multiplications.To be on the safe side we choose ϵ = 10 −8 to ensure numerical stability.For comparison, the measurements b (i,j,l) represent photoelectron counts and are on the order of 10 3 or more.The second case in (13) means that if a voxel does not belong to the set of measured spheres, its updated Fourier value at that point is simply the Fourier value of the previous iteration's reconstruction, taken at that index.This implies that we allow the unmeasured signal to be free as we have no information about it.
The projectors onto the sets LF and SUPP, satisfying the bandwidth constraints in the Fourier domain and the support constraints in the object domain, respectively, are just well known projections onto supports of vectors in the respective domains.Indeed, recall that the ball B r with radius r represents the measurement bandwidth of the image in momentum space.The projector P LF is given by The projection ψ SUPP of a point ψ ∈ C N onto the set SUPP (i.e., ψ SUPP ∈ P SUPP ψ) is given by The projection of a point ψ onto the sparse-real constraint set, SR, with some fixed sparsity parameter s is defined by first projecting onto the set of real-valued vectors, and then projecting onto the sparsity constraint by keeping the s components of Re(ψ) with the greatest absolute values, and then setting the other components to zero.It is important to observe this order of operations since first applying the sparsity constraint, then the realvalued constraint does not yield the nearest point in the intersection of these two sets to the original point.This projector is multi-valued whenever there are ties for the sth largest component.
A projection ψ SYM of ψ onto SYM (i.e., ψ SYM ∈ P SYM ψ) is calculated via the following steps: (i) Take u as the average of the ψ and its reflection along the X-axis: (ii) Take v as the average of u and opposite sign of its reflection along the Y-axis: (iii) Take ψ SYM as the average of v and the opposite sign of its reflection along the Z-axis: A straight-forward argument shows that this is indeed a metric projector; we leave this to the reader.

What to monitor?
To decide when to stop the algorithms, we define: where ∥ψ∥ := (i,j,l)∈I ψ 2 (i,j,l) for all ψ ∈ C N .The iterations are terminated once the distance between successive iterates computed by ( 16) reaches a tolerance, for example TOL = 10 −14 .The theory [41,Example 3.6] establishes that algorithm (11) converges locally linearly for this problem; how long one needs to wait until a region of local linear convergence is found, however, is undetermined.The simplest solution is therefore to wait until the expected local convergence behavior is observed.Figure 5a) and b) show that the range of possible iteration counts before local linear convergence is observed remains relatively stable under changes in other model parameters, like the number of spheres, or the sparsity parameter.
Iterates of the cyclic projection algorithm (11), and in particular final iterate, always belong to the set SYM but not necessarily to the other sets.However, we can always look at the intermediate points from the sets SR, SUPP, LF or M simply by storing these intermediate projections.Specifically, if ψ is the retrieved orbital obtained after completing the iterative process, the "shadow" of this point on the set M is calculated by P M ψ; the corresponding shadow on the set LF is computed by P LF P M ψ, and so forth.Depending on the application, one might prefer one "shadow" over other possible shadows for the reconstructed orbital.In this paper, we opt for reconstructions on the SYM set because they do not retain artifacts due to, e.g., measurement noise.Some practitioners take averages of reconstructions on all sets; we don't recommend this because the mathematical interpretation could be lost, but in the end, this decision is application dependent.
For the simulated 3D-POT analysis we denote the truth, which is the Kohn-Sham molecular orbital that is used for the data simulation, by ψ * .As mentioned earlier, problem (10) is inconsistent, and thus, the point that solves (10) in some sense is not ψ * , but rather a point close to it; this "best numerical solution" is denoted ψ * .There are different ways to define this point.We explain our approach to this issue in more detail in the next section.Regardless of how one defines ψ * , we have the following errors: and where ψ (n) denotes the preferred shadow of the reconstructed orbital at the n-th iteration.
The best numerical solution ψ * is constructed so that it belongs to the same set.The error is the minimum of either the sum or the difference of the reconstruction and the reference "truth" to account for the unavoidable global phase ambiguity, e.g., the change in sign, shown in Fig. 3(d).The formula (17a) is the error between the reconstruction and actual orbital.This metric helps us to estimate reliability of the model.On the other hand, formula (17b) tells us how far the algorithm is from the best numerical truth to the model problem.Thus, the model error (17b) could be small, while the physical error (17a) is large; in this case, it may be necessary to consider adding additional constraints or modifying the model.It is easy to see that 0 ≤ E (n) , E (n) ≤ 1 for all n ∈ N.Moreover, when the model fits perfectly with the truth, i.e, ψ * = ψ * , we get The accuracy of the reconstructed orbital is also influenced by the initial starting point, given the non-convex nature of the problem.While in simulation it is trivial to compare the reconstructed orbital and the truth, this is not possible in experiment.Thus, there is a challenge in selecting the best reconstruction among the multiple possibilities.To address this challenge, we monitor the sum of the gaps between the sets at any given iterate: where gap When the algorithm converges, the gap remains constant but nonzero.For a sufficient number of trials, the reconstruction with the smallest gap corresponds to the point where the different constraints are nearest to each other.If the model of the molecular orbital is reliable, then it is likely that this point provides the best approximation of the true solution (see Fig .4).However, a smaller gap does not necessarily have to imply a smaller error.In our numerical experiments, we will investigate to what extent the gap can be used as a measure for the reconstruction accuracy.

Numerical results
The implementation of algorithms described in this section can be downloaded from the link [36].
Starting with the highest occupied molecular orbital of an isolated pentacene molecule as calculated by density functional theory (Fig. 2a), we generated 3D-POT data for up to 56 photoelectron kinetic energies in the range up to 110 eV such that the spheres are equidistantly spaced in momentum.The individual momentum maps span a momentum range of 11.7 × 11.7 Å-2 with N x = N y = 120 pixels (although our simulated measurement data only covers the central part of the image up to a = 5.5 Å-1 radius).From these momentum maps, we calculate the measured intensities I( ⃗ k (i ′ ,j ′ ,l ′ ) ) in the data domain D, which has dimensions N x = N y = N z = 120.Each voxel in D therefore has a volume 0.098×0.098×0.098Å-3 .Since the orbital is real-valued, we set the intensities I(− ⃗ k) = I( ⃗ k).To define the symmetric support region Ω in the object domain D, we chose a shape of 28 × 12 × 10 voxels (= 15 × 6.4 × 5.4 Å3 ) located at the center of the image.The radius of the low-pass filter B r is set to r = 56 (= 5.5 Å-1 ), as shown in Fig. 3a for the 2D case; as claimed, the bandwidth contains all the spheres for which we simulated measurements to perform the reconstructions.We start with ideal data, generated using a fixed total number of 10 12 detected photoelectrons with constant, well-calibrated vector potential A (i.e., the probe light intensity).In Fig. 5, we investigate the convergence behavior under various scenarios.Frame a) shows 100 globally random instances (i.e., the starting points are chosen randomly over the entire domain space) of the algorithm behavior with a fixed sparsity value s = 3360 and fixed data including all m = 56 measured spheres.All instances exhibit local linear convergence, as expected.The variation of the slopes of the graphs in this frame indicates considerable heterogeneity of local minima.In our experiments, we observed no correlation between the rate of convergence of the algorithm and the quality of the local minima.It is interesting to observe in frame b) the convergence plots of 1000 locally random instances using data m = 56 and s = 3360, i.e., instances whose starting points are randomly chosen in a small neighborhood of a fixed point.Here, we specifically select the neighborhood around the starting point of the best numerical reconstruction, ψ * .Most of the graphs have the same slope in the final iterations, indicating a constant of linear convergence concentrated around the value ≈ 0.38, and its variance is only 6.12 × 10 −5 .This indicates that the geometry of the problem around this local minimum is not wildly varying; a mathematical statement about why this should be the case is beyond the scope of this work.As argued in [42], under the assumption of Q-linear convergence, the rate of convergence r of each instance can be estimated empirically from the observed local numerical behavior.Here, we observe that the mean of r is approximately 0.28 with a small variance of about 1.85 × 10 −5 .Moreover, the average error E over 1000 reconstructions is 3.28 × 10 −12 with a variance of about 4.66 × 10 −33 , suggesting that the fixed point with the smallest gap is an isolated point, which further bolsters the hypothesis that numerical convergence is Q-linear.
In order to benchmark the reconstruction depending on various parameters such as the sparsity parameter, the number of photon energies and the intensity calibration, we now define the best numerical solution ψ * , which may be referred to as the ideal reconstructed orbital for the inconsistent problem (10).This is constructed with the fixed point of algorithm (11) obtained using data from 56 spheres of ideal data with the setting s = 3360 in Fig. 5a); note that we constructed ψ * so that it belongs to the set SYM, and from 100 globally random initializations it was chosen as the solution that achieved the smallest gap.The errors that we report are the distances between points computed by the algorithm in the set SYM.The distance between the truth and the best numerical solution, i.e., 1 2 min{ ψ * ∥ψ * ∥ ± ψ * ∥ ψ * ∥ }, is 7.398%.
In the following subsection (4.2), we present exemplary numerical results that demonstrate the successful reconstruction of the orbital using as little as m = 4 photon energies.In subsection 4.3, we then conduct a more comprehensive investigation of the effects of the sparsity parameter s, the number of spheres m, the starting points, and calibration errors (defined below), exploring their impact on the reconstructions.

Example: 4 photon energies
Figure 6 shows the results of typical reconstructions using m = 4 spheres (radii: 1.9, 3.2, 3.9 and 4.6 Å-1 , photoelectron kinetic energies 13.1, 39.6, 58.2, and 80.3 eV, respectively) of ideal data and sparsity s = 3360 at two different starting points, in comparison to the ideal reconstructed orbital.Here, one point, denoted ψ 1 , achieves an error of E = 0.25, while another point, denoted ψ 2 , achieves E = 0.04.These reconstructed orbitals are selected from a pool of 100 experiments with varying starting points.The reconstruction ψ 2 is the result with the smallest gap among the 100 random initializations; ψ 1 is just a randomly selected trial.
Frame (a) of Fig. 6 shows that the algorithm is converging at a locally linear rate on termination for both instances, as expected.In frames (f) and (k), it is clear that the quality of the reconstructed orbital is affected by the choice of starting points, as measured by the gap and the error, respectively.The fact that the gaps and errors appear unchanged from iteration 500 onward does not necessarily indicate that the solutions have stopped changing, as indicated in frame a).Our objective is to find a point whose corresponding difference vectors are as small as possible; by this reasoning, the better solution should be ψ 2 (orange) as it corresponds to the reconstruction with the smallest gap.Indeed, visual inspection of the constant height (z) maps in frames (b),(g), (l),the constant k z maps in frames (c), (h), (m) and the isosurfaces in frames (d), (i), (n), (e), (j), (o) confirms that ψ 2 is a better approximation of the molecular orbital.Further experimental results regarding the relationship between the gap and error can be found in the next Subsection.

Reconstruction accuracy in experiment
In subsection 4.2 we have seen that, although good orbital reconstructions can be achieved using the CP algorithm for 3D-POT, also but not all random starting points yield the same reconstruction accuracy.In order to investigate how the most optimal results can be achieved, it is first useful to characterize the impact of the sparsity parameter s on the reconstruction result since the sparsity parameter is one of the few control parameters that the user can easily adjust.In this context, also the size (whether larger or smaller) of the support or low-pass filter LF may impact the quality of the reconstruction.In particular, without a low-pass constraint at all, the algorithm tends towards orbitals with un-physically strong high-frequency components.However, a detailed discussion of the effects of the lowpass and support constraints is beyond the scope of the present paper.

The sparsity parameter s
In Fig. 7a), we present the dependence of errors E and E on the sparsity parameter s for a dataset consisting of m = 56 input spheres of ideal data.Each s value underwent testing through a total of 100 globally random trials [i.e., we use the same reconstructions as in Fig. 5a)].We then calculated the averages of E and E across these starting points.We also identify the minimum values of E and E among the 100 trials for each value of s.Here, the minimum achieved error provides an indication of the reconstruction quality that can be achieved, while the mean and spread show that a significant number of trials is necessary to achieve the best result.
Figure 7a) shows a trend that the number of required iterations increases as the number of nonzero voxels s increases, the error decreases as s increases.This holds for both the mean and the minimum errors of E and E. The data shows that a wide spread of reconstruction errors is to be expected for all s.Therefore, it is necessary to perform reconstructions with multiple globally random instantiations and select the best one.Having established that larger values for s generally lead to better reconstruction results, with diminishing returns for s > 2500, we choose to set the sparsity parameter to s = 3360 for the remainder of our  Figure 7b) displays the average number of iterations until convergence over 100 trials, corresponding to the data presented in Fig. 7a) for each sparsity level s.A smaller value of s results in a faster reconstruction process and a smaller spread.The average number of iterations remains the same when s is large, but there is a larger relative variation.

Number of spheres
As the number of photon energies and thereby measured spheres increases, the algorithm is provided with more information and needs to fill fewer voxels.In Fig. 8a), we run a total of 500 reconstructions using ideal data with a fixed sparsity parameter s = 3360 for five various settings of m: 4, 7, 13, 26, 56.The m = 4 data includes the same spheres as described in subsection 4.2, the m = 7 data covers radii of 1.2 to 5.2 Å-1 (5.2, 13.1, 24.6, 39.6, 58.2, 80.3, and 106.1 eV kinetic energy), and the m = 13, 26, 56 data sets cover radii of 0.5 to 5.4 Å-1 .Again, each value of m was tested with 100 globally random instances.
Similarly to Fig. 7a), we measure the mean values of E and E, along with their spreads and respective minimum values, and find that the mean errors decrease as the number of measurements increases.Moreover, the decrease in spreads indicates that the errors exhibit less variation when more data is used.Also as expected, the error measured against the best numerical solution, E, goes to zero as the number of shells increases to 56.
Rather surprisingly, whereas the mean E and E increases significantly as the number of spheres is reduced, the minimum error increases much less.For all numbers of spheres, the minima of E and E stay below 0.1.This observation leads to the conclusion that a dense input is not necessary; using 13, 7 or even 4 shells can still yield a close approximation of the orbital as long as a sufficient number of trials is performed.Finally, we note that the simulated data for m ≥ 13 also includes photoelectron kinetic energies below 5 eV which are  often inaccessible in experiment, e.g., due to strong final-state effects or the presence of a strong inner potential [43].To verify that the inclusion of these shells in the reconstruction does not affect our conclusions, we have also performed reconstructions where we exclude shells for which the photoelectron kinetic energy is less than 14 eV.From this analysis, we find that the exclusion of these shells leads to an increase of the best feasible error of less than 10 −3 from E = 0.07398 to E = 0.07415.We thus conclude that the the lower kinetic energy shells do not affect the reconstruction strongly.In all of these numerical experiments, we chose the photon energies such that the radii of the corresponding spheres are evenly spaced in momentum (for example, see in Fig. 2e) and f) for the case of 7 spheres).It is possible that the energy spacing can be adapted to the specific molecular orbital of interest to provide better reconstruction results, however this falls outside of the scope of the present study.
In subsection 4.2, we have observed an example where examining the gap can offer insights for selecting a reconstruction when both the ground truth and the ideal reconstruction are unknown.Using the full dataset of 500 numerical reconstructions with different m, we can investigate more quantitatively to what extent the gap can be used as an indication for the error E (or E).In Fig. 8b), we depict the dependence of the gap and E for each random instance and input data m.First, we note a clear trend in the data, where a larger error E generally corresponds to a larger gap.Surprisingly, this trend is highly similar between datasets with different m.We emphasize that this match between different m is not expected to be general, but in this case it further supports that there is a direct correlation between the error E and the gap.However, we also observe that the correlation between the gap and E weakens as the amount of data decreases; this is particularly true for m = 4.This is because the problem becomes more ill-defined and, therefore, more non-convex.However, even for m = 4 the smallest gaps do result in the smallest E as claimed.

Calibration error
In experiment, measuring 3D-POT data requires to tune the photon energy to a range of values and to measure the photoemission rate normalized to the amplitude of the incident light field.Consequently, a precise calibration of the -often fluctuating -extreme ultraviolet (EUV) light intensity if necessary.If fluctuations in the EUV power are not properly accounted for, calibration errors may arise which scale the measured intensity on individually measured spheres with respect to the true momentum distribution of the orbital.Furthermore, it is well known that the plane-wave model of photoemission (Eq. 1) does not always provide a good prediction of the photon energy dependence of the photoelectron spectrum [28,44], which can lead to similar errors in the measured momentum distribution.To assess how sensitive the algorithm is to such model miss-specification, we now add EUV intensity calibration errors to the data.
For each m, we generated 10 samples {L ρ } m ρ=1 yielding 50 datasets for each value of m across the respective levels of α.We add 1 ideal dataset (with α = 0) for each m.This results in a total of 255 datasets.
We study the behavior of the algorithm, for fixed algorithm parameters, under the realizations of calibration error as described above.
The results displayed in Fig. 9 show the recovery performance, indicated by the error with the best numerical solution E defined by (17b), over 25500 reconstructions with a fixed s = 3360 (i.e., 100 globally random initializations per simulated dataset).The average number of iterations required for convergence for m = 4, 7, 13, 26, 56 was about 2409, 485, 238, 155, 147 iterations, respectively.Fig. 9a) provides an overview of the average of E depending on m and α.As expected, we observe a decrease in mean error both as the data becomes more accurate (i.e. as α decreases) and as the data becomes richer (i.e. as m increases).Note that for m = 4, the dependence of the mean E on α is relatively weak.This surprising result has a relatively simple explanation: at m = 4, the model already allows for so many false local minima, that the added complexity due to model miss-specification does not worsen the results significantly.
However, these results also show that the mean error is not enough to characterize the algorithm performance.In the following, we define that a given reconstruction is successful if it has an error E ≤ 0.1 (i.e., when the error is below 10%).At this error level, the reconstructed isosurfaces all match well to the true molecular orbital.Strikingly, in Fig. 9c), we find very similar success rates for α ≤ 0.3 as for ideal data with α = 0.The success rates are slightly reduced, but a series of 100 globally random initializations is still   likely to retrieve several successful reconstructions, which can then be selected by of the minimum gap.For larger values of α, the data becomes significantly more inaccurate, resulting in a marked decrease in the percentage of success.Nevertheless, even in these cases a few successful reconstructions are found.From a numerical perspective, this shows that systematic model errors like calibration error do not eliminate physically meaningful reconstructions (interpreted as global minima), but such errors do introduce more local minima that are less meaningful.Furthermore, the miss-specification of the model affects the correlation between the gap and the error, as shown in Fig. 9b.For larger miss-specification (here, α = 0.6 with m ≤ 13), the gap can no longer be used to select the best reconstructions, however we emphasize that this problem only arises for very large miss-specification level which are not realistic in experiment.We therefore conclude that the algorithm is robust in the presence of the calibration error.Optimal results of course require the best possible intensity calibration, but significant errors up to 30% in the intensity calibration are not detrimental.

Number of electrons
Next to the calibration errors that were discussed above, an important source of model inconsistency is intensity noise in the measured momentum maps.Precisely how the momentum maps are acquired differs per experimental setup; however a common factor to all ARPES experiments is electron counting statistics.Here, we compare the effect of reducing signal-to-noise ratio of the data on the quality and reliability of the reconstructed orbitals by simulating datasets with varying total photoelectron counts.Other possible sources of experimental noise, such as background noise or the aforementioned calibration error, are neglected.Each of the datasets in the previous numerical experiments were generated with 10 12 photoelectrons, which is large enough that the data can be assumed to be noise-free.For this experiment, we now generate a dataset including m = 7 photon energies using instead a total of 10 5 detected photoelectrons (see Fig .10 for an exemplary ARPES simulation).
At this signal strength, the measurement data is significantly affected by Poissonian noise, and a significant fraction of the voxels that would otherwise indicate a low intensity now are zero: where in the noise-free case 6.4% of the voxels in the data domain are non-zero, this reduces to 1.6% at 10 5 counts.Nevertheless, the orbital reconstruction is not significantly impacted: As before, we fixed s = 3360, conducted reconstructions with 100 globally random starting points, and found that in this case a success rate 7% (i.e. with %( E ≤ 0.1)).With 10 6 incident electrons, where 3.1% of the voxels is nonzero, the percentage success increases slightly to 8%.We emphasize that these total photoelectron counts are very easily achievable with all of the currently used photoemission orbital tomography detectors.Our results therefore indicate that efforts towards achieving higher data quality should be aimed at reducing other experimental noise factors, such as the determination and subtraction of background signals.

Conclusions
In this paper, we have addressed the challenge of reconstructing full 3D molecular orbitals from sparsely sampled 3D photoemission orbital tomography data.We propose a feasibility model and utilize the cyclic projections method to reconstruct the orbital.The model incorporates prior constraints such as low-pass filtering, support, voxel sparsity, and symmetry properties.Although this leads to an inconsistent feasibility model, we show that the global minimum closely approximates the true solution, which can be found with reasonably high probability by initializing the algorithm at several randomly chosen points and selecting the reconstruction with the lowest gap-value.
Remarkably, we find that the iterative orbital reconstruction method can deal with very sparse data and is robust in the presence of crucial experimental noise factors such as the light intensity calibration and the total signal strength.In particular, a 3D-POT dataset covering only four photon energies in the EUV range already allows for accurate orbital reconstruction with over 8% of the reconstructions achieving an error E within 90% of the best possible value.This indicates that the orbitals to be reconstructed are sparse in some representation (e.g.spherical harmonics).Future research will investigate how to exploit such sparse representations.
Our results are especially promising in the context of laboratory-based, table-top photoemission orbital tomography, where high-harmonic generation EUV light sources provide access to only a limited number of photon energies.The reduced measurement requirements also provide other opportunities, including potentially the extension to time-resolved 3D-POT [19][20][21][22].A limitation of the method remains the relatively large variation in the results, as indicated by the standard deviation of errors in Fig. 7.This reflects a rather large variety of fixed points of the algorithm.Future work will investigate the performance of the this algorithm on experimental data, as well as exploring alternative numerical methods for solving the feasibility problem whose fixed point sets are verifiably smaller.

Figure 3 :
Figure 3: Visualization of (a) the ball B r bounding the low pass region in the Fourier domain.The updated Fourier modulus on the dark background will be set to zero, see (6); (b) the support region which limits the updated value only inside the light box and is applied in real space; (c) sparse real constraint removes values as low as zero, allowing us to obtain a tighter support; (d) shows a change in sign and (e) is an example of a translation, which makes it difficult to use symmetric constraints.(f) shows the true molecular orbital.

Figure 5 :
Figure 5: Algorithm performance of locally and globally random initialization.(a) Iterate differences of 100 globally random initializations (b) Iterate differences of 1000 locally random initializations.

Figure 6 :
Figure 6: Comparison of retrieval results ψ 1 , ψ 2 at two different starting points with the best numerical solution ψ * using ideal data on m = 4 spheres.The sparsity is set to s = 3360.The reconstructions have been normalized by dividing them by their norms for the purpose of visualization.

Figure 7 :
Figure 7: Dependence of (a) errors E and E on the sparsity parameter s, and (b) the average number of iterations needed to achieve convergence using ideal data, setting m = 56.The data points indicate the mean error and minimum achieved error (out of 100 trials), while the shaded regions indicate the standard deviation: in frame a) the green shaded region is the standard deviation of E, the pink-shaded region is the standard deviation of E and the brown shading indicates the overlap.How these results were obtained is described in subsection 4.3.

Figure 8 :
Figure 8: Correlation between a) number of spheres m and errors E, Ẽ and b) gap and error Ẽ using a varying number (m) of spheres, setting s = 3360, 100 globally random initialization for each m, for well-calibrated data.

Figure 9 :
Figure 9: a) The average value of E as a function of calibration error α; b) The correlation of gap and E with α = 0.6; c) The dependence of success percentage ( E ≤ 0.1) on m and α.

Figure 10 :
Figure 10: Exemplary simulated ARPES momentum maps from the m = 7 data for different number of photoelectron counts.