Application of iterative phase-retrieval algorithms to ARPES orbital tomography

Electronic wave functions of planar molecules can be reconstructed via inverse Fourier transform of angle-resolved photoelectron spectroscopy (ARPES) data, provided the phase of the electron wave in the detector plane is known. Since the recorded intensity is proportional to the absolute square of the Fourier transform of the initial state wave function, information about the phase distribution is lost in the measurement. It was shown that the phase can be retrieved in some cases by iterative algorithms using a priori information about the object such as its size and symmetry. We suggest a more generalized and robust approach for the reconstruction of molecular orbitals based on state-of-the-art phase-retrieval algorithms currently used in coherent diffraction imaging. We draw an analogy between the phase problem in molecular orbital imaging by ARPES and of that in optical coherent diffraction imaging by performing an optical analogue experiment on micrometer-sized structures. We successfully reconstruct amplitude and phase of both the micrometer-sized objects and a molecular orbital from the optical and photoelectron far-field intensity distributions, respectively, without any prior information about the shape of the objects.


Introduction
Organic semiconductors play a key role in modern devices such as organic light-emitting diodes and photovoltaic cells [1,2]. More recently, organic molecules have been used as catalysts in photolytic water splitting, a promising route towards production of hydrogen as renewable energy source [3]. Tailoring the physical properties of molecular optoelectronic devices [4,5,6] crucially depends on a deep understanding of the charge transfer mechanisms at metal-organic interfaces. The time-resolved spatial visualization of such processes would hence be highly desirable.
The frontier orbitals, i. e. the highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals, largely determine the chemical reactivity and electronic properties of molecular systems. Detailed information about the electronic structure of molecular systems can be inferred from angle-resolved photoelectron spectroscopy (ARPES) of well-ordered molecular layers on single-crystalline substrates [7,8,9,10,11,12]. The photoemission intensity I(k f , E kin ) is derived from Fermi's golden rule as where ψ i and ψ f denote initial and final state wave functions with corresponding momentum components k i and k f parallel to the surface, respectively. The delta functions in the second line comprising photon energyhω, sample work function Φ and reciprocal lattice vector G ensure energy and momentum conservation in the photoemission process. The transition matrix element is given in the dipole approximation, where p and A denote the momentum operator and the vector potential of the exciting light. The photocurrent I(k f , E kin ) is obtained by summation over all transitions from occupied initial states ψ i to the final state ψ f characterised by the kinetic energy E kin and the parallel component of the final state momentum k f of the photoelectron. The photoemission final state ψ f can be approximated by a plane wave ∝ e ik f r provided the following conditions are fulfilled [7,8,13]: (i) photoelectrons are emitted from π-orbitals of large planar molecules, for which all the contributing orbitals are of the same p z character; (ii) the molecules consist of mainly light atoms (H, C, N, O) and final state scattering effects can thus be neglected. Under these assumptions, the measured ARPES intensity becomes proportional to the squared modulus of the Fourier transform of the initial state wave function ψ i weakly modulated by a slowly varying angle-dependent envelope function [7,8]: The recorded intensity pattern, however, does not contain any information about the phase of the complex-valued electron wave distribution in the detector plane, which inhibits the direct reconstruction of the molecular wave function via computation of an inverse Fourier transform. In certain cases, phase information can be inferred from the parity of the wave function [7] or from dichroism measurements [11] and be imposed onto the measured data. However, the reconstruction of the molecular wave functions in such a way is not applicable to the most general type of problems when the phase distribution cannot be deduced from symmetry considerations. This issue was addressed by Lüftner et al. [10] by suggesting an iterative phase retrieval procedure similar to the Fienup algorithm [14]. In the suggested procedure, one iterates back and forth between real and reciprocal spaces by computing Fourier transforms and satisfying the constraints in both domains. In real space, the wave function is confined to a rectangular box which roughly corresponds to the van der Waals size of the molecule and thus represents the support of the object. The absolute value of the wave function is reduced to 10% outside this confinement box at each iteration step. In reciprocal space, the computed value of the amplitude is replaced by the measured one and the phase is kept.
In this work, we suggest that the phase problem in ARPES-based molecular orbital imaging can be solved in a more robust manner by utilizing the analogy to the phase problem in coherent diffraction imaging (CDI) [15]. Both in CDI and orbital imaging, the far field pattern in the detector plane is proportional to the squared modulus of the Fourier transform of the object distribution. Provided the far-field intensity pattern is measured at the oversampling condition [16,17], both the amplitude and the phase of the object can be reconstructed from the experimentally available modulus of its Fourier transform using the phase retrieval algorithms [14] as it is done in CDI [15]. Therefore, we suggest to directly apply state-of-the-art phase-retrieval algorithms, currently used in CDI, for the reconstruction of molecular orbitals. These algorithms were specifically optimized for objects described by a complex-valued transmission function [18], which makes them ideal for the reconstruction of electron wave functions. Moreover, recent advancements in CDI allowed for the solution of the phase problem without need for the precise knowledge of the shape of the object, which is instead found in the course of the reconstruction using the shrinkwrap algorithm [19]. To facilitate a better understanding of the CDI phase-retrieval algorithms in view of their applicability to reconstruction of molecular wave functions, we designed an optical analogue experiment and performed CDI on micrometer-sized structures produced by means of photolithography. Available CDI phase-retrieval algorithms [14,18,19] were employed for the reconstruction of the micrometer-sized object. Eventually, the same algorithms were applied to a set of ARPES data and the lowest unoccupied molecular orbital (LUMO) of pentacene was reconstructed.

Optical CDI of a microstructure
The microstructures for the optical CDI experiments were patterned in a 105 nm-thick Cr film deposited on a 1.7 mm-thick fused silica substrate, thus providing transparent objects in a non-transparent medium. The individual microstructures had an identical shape but different sizes and were separated from one another by several millimeters to avoid interference between the neighbouring objects. The size of the microstructures was selected in such way that the ratio between microstructure length (e.g., 15 µm) and employed laser wavelength (0.532 µm) was comparable to the ratio between length of pentacene molecule (≈ 1.5 nm) and de Broglie wavelength of the electrons (≈ 0.17 nm) at the used photon energy (50 eV). The experimental setup for optical CDI is shown in Fig. 1. The laser beam profile had a Gaussian distribution as shown in the inset. For CDI experiments, the laser beam is usually spatially filtered and then expanded using two lenses, which ensures that the intensity profile in the object plane is constant [20]. In our experiment, we employed the laser beam without expansion because the light intensity variations on the length-scale of the microstructure were negligible. The microstructure was illuminated from the side of the Cr film. The far-field distribution of the scattered wave was imaged onto a semitransparent screen and the diffraction patterns were recorded with a 10-bit CCD camera (Hamamatsu C4742-95) placed behind the screen as shown in Fig. 1. In order to increase the dynamic range, we recorded several diffraction patterns at different exposures by using a rotatable neutral density filter with optical densities ranging from 0 to 4.0. The recorded images were then combined into one high-dynamic-range (HDR) image by a procedure proposed by Debevec [21].

ARPES of pentacene/Ag(110)
A well-ordered sub-monolayer of pentacene molecules adsorbed on Ag(110) served as model system for orbital tomography. Pentacene ARPES data has been acquired during a beamtime of A. Schöll and coworkers (University of Würzburg) at the NanoESCA beamline at Elettra synchrotron (Trieste, Italy) and has been provided to us for validation of our phase retrieval algorithm [22]. The crystal was prepared according to standard procedures [23] and pentacene molecules [7] were deposited from a homebuilt Knudsen cell [11]. ARPES constant binding energy (CBE) momentum maps of the pentacene LUMO were recorded with the p-polarized light at a photon energy of 50 eV using the photoemission electron microscope (PEEM) [24,25]. The setup of the PEEM and the experimental geometry are shown in Fig. 2(a) and Fig. 2(b), respectively. The microscope was operated in the momentum mode and allowed for detection of electrons with the acceptance angle of α = ±90 0 corresponding to slightly less than ±3Å −1 at 50 eV photon energy without any sample rotation. The CBE map was integrated over a 200 meV energy window, which is of the order of the electron analyzer resolution and of the full-width at half-maximum of the pentacene LUMO at the binding energy of 0.1 eV.

Algorithms
Prior to reconstruction of the pentacene LUMO, we tested the performance of the algorithms on the optical CDI data set, taking advantage of the high dynamic range of these data. We employed a combination of the phase-constrained [18] hybrid inputoutput [14] (PC-HIO) and error reduction [14] (ER) algorithms. The usage of both algorithms in an alternating scheme has been shown to eliminate stagnation problems and to provide faster convergence [14,18,26]. The support of the object was found using the shrinkwrap algorithm [19]. Following the conventional procedure of this algorithm, the initial estimate of the object support was obtained from the autocorrelation of the object by computing the inverse Fourier transform of the experimental diffraction pattern I(X, Y ), convolving it with a Gaussian function (width σ = 5 pixels) and applying a threshold at 10% of its maximum. The pixel values below the threshold were zeroed. The reconstruction began with 40 iterations of the PC-HIO algorithm followed by 2 iterations of the ER algorithm. We found that this number of iterations is sufficient to yield a resonable estimate of the object shape and thus to perform the first update of the object support by using the shrinkwrap procedure [19] described in detail below. The scheme of the iterative phase retrieval procedure is shown in Fig. 3, which included the following steps: (i) In the first iteration k = 1, the experimental amplitude |F (X, Y )| = I(X, Y ) was combined with a random phase and the inverse Fourier transform supplied an initial input object distribution g k (x, y), where (X, Y ) and (x, y) denote the coordinates in the detector and object planes, respectively. We assume the most general case of a complex-valued object distribution and keep both its real and imaginary parts.
(ii) By computing the Fourier transform of g k (x, y), we obtain the complex-valued  (iii) By replacing the calculated amplitude |G k (X, Y )| with the experimental amplitude |F (X, Y )|, while keeping the calculated phase distribution, we obtain an updated complex-valued field distribution in the detector plane G ′ k (X, Y ). (iv) Inverse Fourier transform of G ′ k (X, Y ) provides the output object distribution g ′ k (x, y). (v) In the PC-HIO algorithm [14,18], the input object for the next iteration g k+1 (x, y) is obtained as where β = 0.9 is a feedback parameter and γ corresponds to a set of points which comply with the object domain constraints (belong to the support region and have their phases within an expected range). In the ER algorithm [14], the object distribution g k+1 (x, y) is calculated as where γ fulfills the same criteria as in the PC-HIO algorithm.
The output object distribution g ′ k (x, y) obtained in the last iteration of the ER cycle was used to update the object support. This was done by convolving g ′ k (x, y) with a Gaussian function and setting a threshold at 12% of its maximum, as it is typically done in the shrinkwrap algorithm [19]. The width of the Gaussian was initially set to 2.5 pixels. Upon the first update of the support, the algorithm continued with alternating cycles of 20 iterations of the PC-HIO algorithm followed by 2 iterations of the ER algorithm [14,18,26]. The end of each cycle was finalized by computing a new distribution of the object support. The threshold value and the Gaussian width were chosen empirically so that no part of the reconstructed pattern was truncated, but instead the support converged smoothly towards the shape of the object. The latter requirement was ensured by reducing the width of the Gaussian at every support update by 1% as it is conventionally done in the shrinkwrap algorithm [19]. The quality of the reconstructions was estimated by computing the mismatch between the iterated and the experimental amplitudes [14,27,28]: where |F (X, Y )| is the experimental amplitude, |G it (X, Y )| is the iteratively obtained amplitude.

Oversampling requirements
The solution of the phase problem requires the fulfillment of the oversampling condition [16,17]. Given an N × N pixel sampled amplitude |F (X, Y )| = | N −1 X,Y =0 f (x, y)e −2πi(xX+yY )/N | in reciprocal space, we obtain a set of N 2 equations, which have to be solved in order to find both the amplitude and phase of f (x, y). Miao et al. [16] defined the oversampling ratio as where N total is the total number of pixels and N unknown is the number of pixels with unknown values. The set of equations is solved by dense sampling of the diffraction pattern so that the object distribution is surrounded by a zero-padded region with σ > 2 [16]. In each dimension of a 2D data set, we can define a linear oversampling ratio where N is the linear number of pixels, ∆r is the size of the pixel in the object domain and a is the largest extent of the object. The oversampling requirement then corresponds to Ø > √ 2 [16]. Fig. 4 shows the results of the reconstruction of the micrometer-sized structures. In optical CDI, we employed micrometer-sized structures of 30 × 12 µm 2 (sample 1) and

Optical CDI: Reconstruction of the micrometer-sized structures
14.8 × 6 µm 2 (sample 2). The scanning electron microscope (SEM) images of samples 1 and 2 are shown in Fig. 4(a,b) next to the experimental diffraction patterns (Fig. 4(c,d)). The size of the diffraction patterns sampled with 1000 × 1000 pixel was 40 × 40 cm 2 in each case, thus giving the size of the pixel in the detector plane ∆p = 400 µm. The size of the pixel in the object plane ∆r can be related to the distance z = 22.5 cm from the object to the detector plane and to the employed laser wavelength λ = 532 nm [29]. The linear oversampling ratio defined by Eq. 7 can be rewritten as [16]: For the samples 1 and 2 with the lengths a 1 = 30 µm and a 2 = 15 µm, the linear oversampling ratios fulfilled the oversampling condition and were Ø 1 ≈ 10 and Ø 2 ≈ 20.2, respectively.  Prior to application of the phase retrieval algorithms, the experimental diffraction patterns were pre-processed: First, each of the recorded 1000 × 1000 pixel images was centered. Centering of the experimental diffraction pattern was shown to have a strong effect on the quality of the reconstruction in CDI [30]. The noise of the CCD camera (average count rate of 50 counts) was subtracted from each pixel and the images were truncated to 500×500 pixels around their centers because of the low signal-to-noise ratio at the peripheral parts. The central part of each diffraction pattern was dominated by an intense laser spot due to the partial transparency of the chromium film to the laser beam. Pixel values exceeding the thresholds of 1.5 · 10 5 counts (sample 1) and 6 · 10 3 counts (sample 2) were defined as missing and their values were updated in the course of the reconstruction by using the corresponding pixel values of the calculated amplitudes in the detector plane [19,31]. In each case, the square root of the resulting diffraction pattern was fed into the algorithm. We found that 10 alternating cycles of the PC-HIO and ER algorithms, each followed by an update of the support, were enough to achieve a stable reconstruction. Further increase in the number of the reconstruction cycles was not necessary since it did not improve the quality of the reconstructed object distribution. At the end of 10 cycles, each reconstruction was stabilized by 100 iterations of the ER algorithm [31]. In total, we performed 1000 independent reconstructions by employing a random phase distribution for each reconstruction run. Eventually, the 50 reconstructions with the smallest error E as defined by Eq. 5 were selected and averaged [31] and are shown in Fig. 4(e-h). The reconstructed amplitudes correctly reproduce the shape and dimension of the microstructures. Furthermore, as it was expected for a purely transmitting object illuminated by a Gaussian beam with an almost planar wavefront at the object site, the phase distributions turned out to be almost constant. The lower quality of the reconstructed amplitude of sample 2 ( Fig. 4 (g)) can be attributed to the low signal-to-noise ratio in the respective diffraction pattern.

ARPES orbital tomography: Reconstruction of the pentacene LUMO
We then applied the same algorithm to the ARPES data. Fig. 5 shows the results of the reconstruction of the pentacene LUMO. The experimental CBE map is shown in Fig. 5(a). Given the resolution in reciprocal space of ∆k ≈ 0.01Å −1 and the length of the pentacene molecule a ≈ 15Å, the linear oversampling ratio in the ARPES experiment can be calculated using Eq. 7. Taking the relation ∆r∆k = 2π N between the pixel size in object space, ∆r, and reciprocal space, ∆k, into account, the linear oversampling ratio can be expressed as The linear oversampling ratio was Ø ≈ 42 and thus fulfilled the oversampling condition [16]. The experimental CBE map was pre-processed following similar steps as those applied to the reconstruction of the micrometer-sized objects: First, the image was centered and the quasi-constant noise of the CCD camera (average count rate of 50 counts) was subtracted from each pixel. To ensure a sufficient number of pixels allocated per unit length of the molecule, we zero-padded the experimental CBE map to 2000 × 2000 pixels around its center. The square root of the processed CBE map was fed into the algorithm with the same parameters as used for the reconstruction of the micrometer-sized structures. Varying these parameters did not lead to any substantial improvements in the quality of the reconstruction. In total, we performed 1000 reconstructions of the pentacene LUMO. About 56% of the reconstructed objects g(x, y) were reconstructed together with their conjugate g * (−x, −y) or twin images [28]. The identification of the twin images could be automated by a procedure proposed by Fienup [28], but here they were easily identified by visual inspection and discarded. From the remaining reconstructions, 50 with the smallest error E as defined by Eq. 5 were selected and averaged. The reconstructed amplitude and phase of the pentacene LUMO are shown in Fig. 5 (b-c) together with the overlayed carbon frame of the molecule for comparison. It should be noted that we did not perform a normalization of the ARPES intensity by the angle-dependent factor |A · p| 2 nor did we enforce any symmetry constraints in the course of the reconstruction onto the amplitude and phase shown in Fig. 5 (b-c). The object distribution was let to freely evolve until the stable solution was reached, which makes the utilized algorithm independent of any symmetry properties imposed onto the object under reconstruction. Furthermore, we note that the recorded CBE map shown in Fig. 5(a) contains features coming from the Ag(110) substrate (mostly at high momenta), but they do not seem to have a profound effect on the results of the reconstruction. By comparing our results with the literature, we find that the phase distribution weighed with the correspondent amplitude values as well as the shape of the orbital correctly reproduce the DFT calculations [10,32] as well as the data reconstructed by Lüftner et. al [10].
Finally, in order to assess the robustness of the algorithm in terms of the quality of reconstruction from the unsymmetrized CBE map, we made use of the symmetry properties of the pentacene LUMO amplitude and phase and symmetrized the CBE map shown in Fig. 5 (a) around its center. The symmetrical version is shown in Fig. 5 (d). In optics, the far field diffraction pattern is symmetric only in two cases: either due to the real-valued nature of an object or in case of an even complex-valued object distribution with an even amplitude and an even phase. In the latter case, the Fourier transform of the even complex-valued function is an even function as well and the far field intensity distribution is therefore symmetric. In the case of the complex-valued wave function of the pentacene LUMO, the symmetrization is justified purely due to the symmetry of the LUMO amplitude and the phase as it is known from the DFT calculations [10,32]. The symmetrized CBE map was pre-processed following the same procedure as described above and the results of the reconstruction are shown in Fig. 5(d-f). Qualitatively, the reconstructions from the unsymmetrized CBE map are as good as the reconstructions from the symmetrized data set, except for some minor differences in the shapes of the lobes due to the intrinsic asymmetry of the CBE in Fig. 5(a). This agreement further proves the robustness of the employed algorithm for the reconstruction of molecular orbitals with arbitrary symmetry properties.

Summary and Conclusion
In this work, we show that the state-of-the-art phase retrieval algorithms currently employed in CDI can be successfully used for the reconstruction of complex-valued wave functions of molecules adsorbed on single-crystalline substrates. We tested and applied these algorithms in an optical analogue experiment and then successfully applied them to the reconstruction of the LUMO of pentacene adsorbed on Ag(110). The advantage of using modern CDI algorithms and in particular the shrinkwrap algorithm for the reconstruction of molecular orbitals is that they do not require any a priori information about the shape of the object. Instead, they smoothly converge to the correct shape of the object in the course of the reconstruction. In case of molecular wave functions, this is highly important, since precise estimation of the object support is difficult and cannot be guaranteed in every case. This applies, for instance, if the orbital tomography technique aims at visualizing chemical reactions or following the dynamics of excited states, where effective electronic wave functions are unknown. The availability of a general and robust reconstruction algorithm is thus an important step for further advancement of orbital tomography.