One-Step Theory View on Photoelectron Diffraction: Application to Graphene

Diffraction of photoelectrons emitted from the core 1s and valence band of monolayer and bilayer graphene is studied within the one-step theory of photoemission. The energy-dependent angular distribution of the photoelectrons is compared to the simulated electron reflection pattern of a low-energy electron diffraction experiment in the kinetic energy range up to about 55 eV, and the implications for the structure determination are discussed. Constant energy contours due to scattering resonances are well visible in photoelectron diffraction, and their experimental shape is well reproduced. The example of the bilayer graphene is used to reveal the effect of the scattering by the subsurface layer. The photoemission and LEED patterns are shown to contain essentially the same information about the long-range order. The diffraction patterns of C 1s and valence band photoelectrons bear similar anisotropy and are equally suitable for diffraction analysis.


Introduction
Photoelectron diffraction (PED) [1][2][3] has long been established as a powerful tool for the determination of the surface crystal structure along with the low-energy electron diffraction (LEED) [4,5] and very low-energy (VLEED) I(V) technique [6]. A conclusive interpretation of the diffraction patterns and a reliable determination of the atomic geometry depend on an adequate theoretical modeling of multiple scattering. Powerful computational methods have been developed that treat the multiple scattering in real space and are indispensable for the study of the short-range order of surfaces and adsorbates [7][8][9]. Recent progress in the fabrication of two-dimensional (2D) materials (such as graphene or h-BN) draws attention to methods of determining the geometry of long-range-periodic 2D crystals. The recently developed angle-resolved reflected-electron spectroscopy (ARRES) technique [10,11] has allowed angular scanning of electron reflection at very low energies and revealed rich structure of the R(E, k || ) diffraction patterns carrying information about the surface geometry and the electronic structure. The angular distribution of photoemitted electrons I(E, k || ) presents an alternative diffraction pattern. Currently, the experimental focus has been on ultra-high kinetic energies [3], where the theoretical analysis allows certain simplifications [12]. On the other hand, low energies have some advantages, such as higher sensitivity to the local environment and lesser effect of atomic vibrations. In the present work, the low-energy regime is theoretically explored with the aim to interpret the electron reflection and photoemission maps on the same footing and understand the relation between PED and LEED.
This calls for theoretical methods capable of treating periodic systems with minimum adjustable parameters and providing a straightforward connection to the underlying band structure. The basis for such approaches is the one-step photoemission theory, where the photoemission intensity equals the probability of the transition to the timereversed LEED state Φ * k || [13][14][15][16][17]. The most common implementations of the one-step theory employ the multiple-scattering formalism within the layer-KKR [17,18] or real-space cluster approach [19]. The present study employs a one-step theory realized through the variational embedding method [20], where the scattering properties of a 2D slab are represented by the band structure of a supercell containing the slab. In some contexts, the term PED implies photoemission from atomic-like core states, so the outgoing wave can be calculated as the scattering of a spherical wave emanating from the atom by the surrounding atoms. In this case, one does not need to construct Φ * k || for each detection angle, which is exploited in the PED theories based on Refs. [7,8].
Nevertheless, the application of the one-step theory is instructive because it uses the same function Φ k || to describe scattering by the crystal potential in the LEED and in the PED experiment and is not limited to atomic-like states. The aim of the present work is to compare the diffraction patterns (energy-dependent angular distributions) observed in VLEED and in PED. Monolayer and bilayer graphene are used as examples. The dependence of the PED diffraction patterns on the initial state will be analyzed and the effect of backscattering will be considered.

Computational Methodology and Approximations
The LEED wave function Φ k || (r) is a scattering solution for a plane wave incident from vacuum with the energy E and surface parallel wave vector k || . It is a Bloch function with the crystal momentum k || , and inside the graphene slab it satisfies the Schrödinger equation with the HamiltonianĤ = −∆ + V(r). In the present calculation, the imaginary potential responsible for the inelastic scattering is not included. The crystal potential V(r) is the self-consistent all-electron Kohn-Sham potential obtained within the local density approximation of the density functional theory using the numerical procedure described in Ref. [21].
The scattering calculation setup for the graphene bilayer is shown in Figure 1. (The parameters for the monolayer are presented in Ref. [22].) The plane wave normalized to unit flux is incident from the right vacuum half-space z > z R , and the current in the left vacuum half-space, z < z L , is the transmitted current T(E). In the scattering region [z L , z R ] the wave function is a linear combination of the eigenfunctionsĤξ n = n ξ n of an auxiliary three-dimensional z-periodic crystal (lattice constant c = 39.6 a.u.), which contains the graphene slab as part of the unit cell [20]. The solution of the scattering problem is sought as a linear combination of the basis functions ξ n that satisfies the Schrödinger equation for the energy E in the scattering region z ∈ [z L , z R ] and at the boundaries z L and z R matches the function and derivative of the plane-wave representations in the respective half-spaces (where it satisfies the Schrödinger equation by construction). The problem is solved with the variational embedding method in the augmented plane waves (APW) representation [20]. The basis set in the scattering region comprises the ξ n functions with energies up to max ∼130 eV above the Fermi energy, which for the bilayer graphene amounts to around 400 ξ n functions.
To calculate the momentum matrix elements, a Laue representation of the LEED state is constructed by a straightforward expansion of the all-electron wave function in terms of 11,999 plane waves (with wave vectors smaller than 9.8 a.u. −1 ). It comprises 19 surface reciprocal lattice vectors G || : The photocurrent I(E, k || ) is calculated as the probability of transition from the eigenstate of the graphene slab | i k || to the time-reversed LEED state | f k || : where r | f k || = Φ * −k || (r) andô is the dipole operator,ô = −i∇ · e. The light polarization e in the present calculation is chosen to be along the surface normal. In order to draw the connection with Ref. [23], let us write the Bloch state r | i k || of an isolated initial-state band λ as a lattice sum of Wannier functions: ψ λk || (r) = ∑ R exp(ik || R)η λ (r − R). (The final state cannot be written in this form because Φ k || is not k || periodic.) Then it follows immediately from Equation (2) that the transition matrix element is Φ −k || (r)ô η λ (r) dr, which is physically the same as Equation (5) of Ref. [23]. Note that this expression is not a result of a Brillouin zone averaging. For the emission from core states, the Wannier function is clearly close to a linear combination of the atomic-like states, while for the valence band η λ (r) may have a complicated shape and be rather extended.

Structure of the LEED Wave Functions
Two examples of scattering solutions for the AB-stacked bilayer graphene for the normal incidence (k || = 0) are presented in Figure 1a,b for E = 22 and 14 eV, respectively, as the density profiles The energy E = 14 eV lies at the lowest transmission minimum, while E = 22 eV is close to a transmission resonance T(E) = 1 [10,11,24], see the inset in Figure 1. With increasing the number of layers, the interval around the T(E) minimum becomes a k || -projected gap, while the number of the T(E) resonances grows, and in graphite they merge to form a conducting band, see Ref. [25]. At the T(E) minimum, the wave function is seen to rapidly decay into the interior of the bilayer, while at the (almost) full transmission the density distribution is almost symmetric relative to the z = 0 plane. The symmetric ρ(z) at T = 1 is the consequence of the fact that in the absence of reflection, the scattering states for the wave incident from the right and from the left are related by the time reversal.
The red-shaded areas in Figure 1a,b show the contribution of G || = 0 harmonics (surface Umklapps) to the Φ k || function, see Equation (1). Somewhat surprisingly, this contribution is much smaller at the high reflection (which implies strong scattering) than at the transmission resonance. This means that while the scattering resonances can be qualitatively reproduced in a 1D model crystal, the 1D wave function would be rather unrealistic, see also the analysis in Ref. [26].

Results and Discussion
The k || distribution of the specular reflectivity R for the monolayer graphene is compared in Figure 2 to the photoemission intensity from the 1s carbon core band and from the 7 eV-wide lowermost valence band (VB) (see, e.g., Figure 1 in Ref. [27]) for three finalstate energies: E = 20, 29, and 54 eV. The VB disperses from −19.3 eV at Γ to −12.4 eV at K, so the intensity distribution I VB (k || ) implies a scan over the initial state energies. Note that to facilitate the comparison, constant-E maps are shown for VB PED rather than constant-hω ones.
The three maps in each row manifest similar gross features, which reflect the structure of the Φ k || states, although the details may be rather different. For example, for E = 20 eV, the central hexagon of side 0.7 Å −1 has very similar shape in R(k || ) and I VB (k || ), but it looks very differently in I 1s (k || ). This hexagon originates from a sharp scattering resonance that is specific to graphene: it was theoretically predicted in Ref. [27] and first experimentally observed in Ref. [28]. The present results agree well with the measurements in Ref. [28], regarding both the orientation and dispersion of the hexagon: around 20 eV, its size is about 7% smaller in the present theory than in the experiment (according to the low-transmission signature of the resonance), which corresponds to an upward shift by about 1 eV of the measured unoccupied band structure relative to the calculated one. (This is a typical value of the self-energy correction for the Kohn-Sham states at such energies.) The concave-sided larger hexagon formed by arc-shaped narrow stripes is present in all three E = 20 eV spectra, as well as in the measurements of Ref. [28]. At E = 29 eV, the small hexagon transforms into a star with rays in the ΓK direction, which again is similar in the R and I VB maps but appears different in the I 1s map. On the other hand, the elongated petals pointing in the ΓM direction bear the same shape in all three maps. Apart from this, the VB PED map manifests a hexagonal shape that coincides with the 2D Brillouin zone (BZ) of graphene. The well-visible BZ contour is due to the highly dispersive initial states, whereby the initial state rapidly changes character in crossing the BZ boundary. This feature persists over an energy interval from 20 to 39 eV, see the movie film1-MLRCV.mov provided in the Supplementary Material. Clearly, it provides the most direct information about the lattice constant. Finally, the E = 54 eV maps present an opposite example, where the reflectivity and 1s emission produce very similar patterns, and the VB one is different.
The question of interpreting the PED patterns obtained from an extended initial-state band has been first raised in Ref. [29], where a strong similarity in the anisotropy of the core and VB X-ray emission from Al(100) was observed. A strong PED effect was observed in the ultraviolet emission from Cu(100) and Cu(111) [23], with a significanthω dependence. An instructive aspect of the present calculation is that no band-averaging is involved as the extended initial-state band is strictly two-dimensional and energetically isolated. It can be concluded that the three diffraction experiments manifest similar anisotropy, see Figures 2 and 3. For a reliable determination of structural parameters, it may be useful to compare diffraction maps of different incident wave sources in order to reveal common features. Electron diffraction from the monolayer graphene. Constant energy k || distribution of the specular reflectivity R, photocurrent from the C 1s core states I 1s , and from the highly dispersive lowermost band I VB are shown in the 1st, 2nd, and 3rd columns, respectively. Three energies are presented: 20, 29, and 54 eV in the 1st, 2nd, and 3rd rows, respectively. The ΓM line is placed horizontally. The value at the top of the color scale bar indicates the map maximum, which presents the fraction of the specularly reflected current in the R maps, and arbitrary units in the I maps (the same units in all I maps). The full movie for E = 5 to 59 eV can be found at film1-MLRCV.mov.
The electron diffraction maps from the AB bilayer are presented in Figure 3 for the same three energies. The bilayer bears C 3v symmetry, and thus the PED maps naturally do not exhibit vertical-mirror symmetry axis. At the same time, the R(k || ) maps show C 6v symmetry, which is the consequence of 3D inversion symmetry of the bilayer. Indeed, one can prove (see Appendix A) that the transmission and reflection coefficients, t and r, for the incident plane waves with opposite k || are related as r − = r + and t − = −t * Thus, the reflectivity R = |r| 2 is the same for k || and −k || , whereas the wave functions are generally different. Similar to the monolayer graphene, the map at 20 eV is dominated by the scattering-resonance hexagon, which is to be expected since this resonance originates from coupling of the surface-perpendicular and in-plane motion within one layer [27]. Generally, the large-scale features in the monolayer maps have their counterparts in the bilayer maps, but at smaller k || the differences are significant, which points to a more important role of the interlayer scattering at smaller incidence (emission) angles.  Let us now consider the role of backscattering in the bilayer graphene. One should keep in mind that the time-reversed LEED state is not the true final state of the photoemission process (as seen from the presence of the incoming wave ψ * r traveling from the detector in Figure A1b), so it does not provide detailed information about the multiple scattering of the photoexcited wave. However, for the C 1s states, we can make use of the fact that the coherence of the 1s orbitals at different sites is physically unimportant and consider separately the emission from the front layer and from the back layer, see Figure 4. In terms of multiple scattering, the electrons photoemitted from the front layer experience extra backscattering, while those emitted from the back layer are additionally scattered by the front layer. For all energies considered, the intensity distribution of the front-layer emission is rather close to the true full bilayer pattern, see Figure 4, while for the back-layer emission this depends on energy (see the full movies at film3-ABBFT.mov for AB and film4-AABFT.mov for AA stacking). In the example shown in Figure 4, the I 1s (k || ) patterns from the back and front layers are surprisingly similar. The total photoyield from the front layer is everywhere significantly larger than from the back layer, see Figure 5. This result is in accord with an intuitive single-scattering picture: the photoelectrons emitted from the front layer in the direction away from the detector are reflected from the back layer to add to those traveling towards the detector without backscattering, and the photoelectrons emitted from the back layer are partially reflected from the front layer and do not reach the detector. Indeed, the photoyield from the stand-alone graphene monolayer is larger almost everywhere than from the back layer and smaller than from the front layer for both stackings, see Figure 5.  Diffraction patterns are seen to depend substantially on the number of layers and on the stacking order, concerning not only the higher symmetry of the I 1s (k || ) distribution for AA stacking but often the overall shape of both the I 1s (k || ) and R(k || ) patterns (as in the example of Figure 4). However, at certain energies there appear thin arc-like stripes that are present both in I 1s (k || ) and R(k || ) maps in all the systems, for example, in Figures 2 and 3 one can see the lines of enhanced (dark lines at 20 eV) or reduced (light lines at 54 eV) reflectivity. The light lines form the transmission slits typical of 2D crystals. They are best seen in the R(k || , E) plots in the upper row of Figure 6. For both bilayers, two transmission slits due to interlayer scattering appear around Γ at low energies: the large concave-up arc (well-known experimentally [10,11,24]) disperses from 6 to 12 eV and the small concave-down arc from 21 to 23 eV. Both transmission resonances correspond to an enhanced photoemission intensity, see lower row in Figure 6. Apart from this, in both directions, R(k || , E) manifests almost straight high-transparency lines with negative dispersion: from 55 eV at 1.2 Å −1 (1.4 Å −1 ) to 31 eV (43 eV) at 2 Å −1 along ΓK (ΓM). This feature is present both in the bilayers and in the monolayer, which points to its in-layer-scattering origin. In the I 1s (k || , E) maps, this corresponds to a reduced intensity, in contrast to the features of interlayer origin. Interestingly, in the I 1s (k || , E) maps it is seen more distinctly than in the R(k || , E) ones.
We have seen that the difference in stacking results in quite a different structure of the PED and reflectivity maps. However, they rapidly vary with energy, which implies certain difficulties in using them as structure fingerprints. At the same time, in the present example there exist gross features that clearly discriminate between AB and AA stacking: Note the rapid growth in reflectivity above 50 eV around the normal incidence for AB stacking, which does not occur for the AA stacking. It is accompanied by an even stronger increase in C 1s emission from the AB bilayer that is not seen in the AA one. (Interestingly, for k || = 0, both R(E) and I 1s (E) spectra are virtually identical for the two stackings below 30 eV.) To summarize, the scattering of the incoming (LEED) and outgoing (PED) electrons share many common features; however, the respective structures may bear the same or inverse (maxima instead of minima) character. Figure 6. Energy-momentum distribution of reflectivity (upper row) and C 1s photoemission (lower row) for graphene monolayer and bilayers of AB and AA stacking, 1st, 2nd, and 3rd columns, respectively. Shown are ΓK (negative k || ) and ΓM (positive k || ) directions.

Conclusions
The present proof-of-concept calculation demonstrates the prospect of employing the PED technique at low energies in application to 2D crystals and interpreting energydependent intensity maps within the one-step photoemission theory. In particular, in the example of graphene monolayer and bilayers, the LEED and PED patterns are shown to contain essentially the same information about the long-range order, the main difference being that for the AB-stacked bilayer, the LEED pattern bears higher symmetry than the PED pattern (owing to the combination of time-reversal and inversion symmetry). The calculations confirm the important role of in-layer scattering in graphene multilayers: many characteristic features are already present in the graphene monolayer. The layer-resolved C 1s PED patterns are found to be sensitive to the location of the initial state, but this depends on the final-state energy. On average, the photoemission from the back layer is weaker than from the front layer, which is the result of solely elastic scattering. The PED patterns from the atomic C 1s and extended VB states differ in detail but have similar anisotropy and common large-scale features and are equally suitable for diffraction analysis.
Funding: This research was funded by the Spanish Ministry of Science, Innovation and Universities (Project No. PID2019-105488GB-I00).

Data Availability Statement: Not applicable.
Acknowledgments: The author is grateful to Aleš Paták and Roman Kuzian for the careful reading of the manuscript and valuable comments.

Conflicts of Interest:
The author declares no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Consider the incident ψ i , reflected ψ r , and transmitted wave ψ t in Figure A1a, ψ r (r) = r exp i(k || r || − k ⊥ z), (A2) ψ t (r) = t exp i(k || r || + k ⊥ z).  Figure A1b shows the time-reversed wave function and Figure A1c the scattering of the wave with opposite momentum −k || (the electron source is rotated by π around the surface normal). Finally, Figure A1d is the spatial inversion of Figure A1a:ψ(r) = ψ(−r). We are interested in the reflection and transmission coefficients for the case "c": ψ π r (r) = r π exp i(−k || r || − k ⊥ z), ψ π t (r) = t π exp i(−k || r || + k ⊥ z).
Case "d" can be represented as a linear combination of cases "b" and "c":ψ = bψ * + cψ π . By equating the coefficients of each of the four waves in both half-spaces, we can write: incident from the left half-space, (A6) br * exp i(−k || r || + k ⊥ z) = − c exp i(−k || r || + k ⊥ z) incident from the right half-space, (A7) outgoing in the right half-space, (A8) ct π exp i(−k || r || + k ⊥ z) = r exp i(−k || r || + k ⊥ z) outgoing in the left half-space.
This gives r π = r and t π = −t * r r * , which explains the C 6v symmetry of the reflectivity patterns in Figure 3.