Supported data and manuscript "Layered multiple scattering approach to Hard X-ray photoelectron diffraction: theory and application"

npj | computationalmaterialsArticle Published in partnership with the Shanghai Institute of Ceramics of the Chinese Academy of Sciences https://doi.org/10.1038/s41524-025-01653-y Layered multiple scattering approach to Hard X-ray photoelectron diffraction: theory and application Check for updates Trung-Phuc Vo1,2,OlenaTkach 3,4, Sylvain Tricot5, Didier Sébilleau5, Jürgen Braun6, Aki Pulkkinen1, Aimo Winkelmann7, Olena Fedchenko3, Yaryna Lytvynenko3,8, Dmitry Vasilyev3, Hans-Joachim Elmers3, Gerd Schönhense3& Ján Minár1 Photoelectron diffraction (PED) is a powerful technique for resolving surface structures with subangstrom precision. At high photon energies, angle-resolved photoemission spectroscopy (ARPES) reveals PED effects, often challenged by small cross-sections, momentum transfer, and phonon scattering. X-ray PED (XPD) is not only an advantageous approach but also exhibits unexpected effects. We present a PED implementation for the spin-polarized relativistic Korringa-Kohn-Rostoker (SPRKKR) package to disentangle them, employing multiple scattering theory and a one-step photoemission model. Unlike conventional real-space approaches, our method uses a k-space formulation via the layer-KKR method, offering efficient and accurate calculations across a wide energy range (20-8000 eV) without angular momentum or cluster size convergence issues. Additionally, thealloyanalogymodelenables simulations offinite-temperatureXPD andeffects insoft/ hard X-ray ARPES. Applications include modeling circular dichroism in angular distributions (CDAD) in core-level photoemission of Si(100) 2p and Ge(100) 3p, excited by 6000 eV photons with circular polarization. Asaninfluentialtechniqueforprobingthe atomicstructurein thevicinityof a given emitter, photoelectron diffraction (PED) is widely used for various investigation purposes: crystal structures, bonding geometries of atoms, and the local environment of impurity or dopant atoms inside surfaces1–4.Itis considered analogous to angle-resolved photoemission spectroscopy (ARPES) at the fundamental level, specifically regarding the angular distribution of photoelectrons emitted from a crystal surface. However, the underlying physics and the investigative objectives of the two approaches are different. The angular distribution of the emitted electrons represents the momentum of the initial states in ARPES, while in PED it reveals the interference of photoelectron waves from final states. The observed change in PED intensities depends on theconstructive or destructive interference of allthecoherentlyemitted electrons.Thismeansthatdetails relatedtoatomic arrangements, bonding orientations, and distances, and electronic and chemical properties are analyzed. Depending on the utilized photon energies with respect to photoelectron kinetic energies, this tool can be termed either ultraviolet-PED (UPD) or X-ray-PED (XPD). In summary, photoemission can provide information about core-level electrons or delocalized valence band electrons. Accordingly, there are core-level PED (CL-PED) and valence-band PED (VB-PED). The main difference is that the angular distribution in VB-PED is caused by two interference phenomena: from primary waves emitted at different atomic sites and from the scattering of corresponding secondary waves at neighboring sites5.Athighphoton energies, XPD effects6–12 obscure ARPES data, along with other complications such as poor cross-sections, considerable momentum transfer from photons causing a pronounced modulation of ARPES patterns, and substantial phonon scattering. Given the significant contributions of PED (especially in the circular dichroism in the angular distribution (CDAD) interpretation and disentangling of XPD effects), the development of a suitable theoretical model of scattering is of paramount importance.Powerful computational approaches that address single and multiple scattering in real space and reciprocal space 1New Technologies-Research Centre, University of West Bohemia, Pilsen, Czech Republic. 2Institute of Physics, Czech Academy of Sciences, Praha 6, Czech Republic. 3Johannes Gutenberg-Universität, Institut für Physik, Mainz, Germany. 4Sumy State University, Sumy, Ukraine. 5Univ Rennes, CNRS, IPR (Institut de Physique de Rennes) - UMR 6251, Rennes, France. 6Department Chemie, Ludwig-Maximilians-Universität München, München, Germany. 7AGH University of Krakow, Academic Centre for Materials and Nanotechnology, Kraków, Poland. 8Institute of Magnetism of the NAS of Ukraine and MES of Ukraine, Kyiv, Ukraine. e-mail: [email protected] npj Computational Materials | (2025) 11:159 1 1234567890():,; 1234567890():,; havebeendeveloped earlier andwill besummarizedand reviewed inthe section 4.1. However, these packages share some limitations. First, most of them perform calculations in real space13–15,whichmightnotbe convenient for high kinetic energies where the effective cluster sizes need to be very large. Second, some models have to borrow the atomic potential from other sources. For instance, in the MSCD code16,itis taken from the muffin-tin potential in the database of Moruzzi et al.17, which neglects relativistic effects. For low kinetic energy, MsSpec18 utilizes the potential file from Munich SPRKKR19,20. In the case of the dynamical high-energy electron diffraction approach to XPD21,22,the real and imaginary parts of the crystal potential (assumed to have 3D bulk periodicity) are computed based on the Fourier coefficients of atomic potentials available in various parameterizations23.Third,some approachesaredesignedtocopewith aspecific regimeofkineticenergy, not flexibly covering the whole range from UV to hard X-ray. For instance, EDAC24 is designed to describe diffraction at lower kinetic energies in comparison to the method proposed by Winkelmann et al.21,22. Fourth, although hard X-ray calculations (e.g., a chromium K α source)25,26 attempt to simulate XPD patterns, not all of them are availableforCDADpatternsforcorelevelsinthisregime,tothebestof our knowledge. Indeed, the TMSP code15 developed by Matsushita and colleagues can work with circularly polarized light27 and interpret the CDAD in terms of the rotational shift of forward focusing peaks around theincident-light axis28,29. However, no simulated CDAD for hard X-ray core levels has been published. Fifth, Kikuchi diffraction30 is welldocumented in scanning and transmission electron microscopy (SEM and TEM) with energy sources up to 100 keV, where it manifests as electron backscatter diffraction patterns31. Nevertheless, it is challenging to observe this phenomenon in photoemission spectroscopy, and few simulations are able to reconstruct it32–34. Lastly, to accurately represent scattering by cluster approaches, there is a high need for large values for the maximum angular momentum lmax (the number of scattering phase shifts). In the case of 10 keV electrons scattering from an atomic potential with a radius of 1 Åradius, this parameter must be increased to 10022. In our earlier effort33,majorfine features of W 3d 5/2 at hν=6 keV were relatively well reproduced by calculations. Based on this work and previous research35, which highlighted a more than 50% contrast in the spin-orbit doublets and the discriminant diffractogram of sub-levels in the case of W(110), this study describes our method utilized in the above-mentioned works in details and resolves these combined obstacles for Si(100) and Ge(100) at hν=6 keV. Our focus will be on studying fine Kikuchi lines and CDAD for core-level emission XPD (CL-XPD), which is more widely utilized compared to its counterpart, valence-band XPD (VB-XPD)36–39. Results Convergence tests This study makes use of the layer Korringa-Kohn-Rostoker (LKKR) scheme in the spin-polarized relativistic Korringa-Kohn-Rostoker (SPRKKR) package, which is based on multiplescatteringtheory.The theory posits that a system consists of atomic layers, with the scattering characteristics of each layer being computed using a partial-wave basis set. Here, Green’sfunctions are expressed in terms of spherical harmonics and radial functions, allowing the problem to be separated into angular and radial parts. This is where the maximum angular momentum lmax comes into play. The layers are then connectedina plane-wavebasis tocreateasolid. The inter-layerscatteringis handled by expanding Green’s functions into the infinite set of reciprocal lattice vector ~ Ghkl.Atthisstage,itistheroleplayedbythenumberof~ Ghkl. Nevertheless, calculations are carried out by the finite values of these two basis sets. It is not straightforward to predict their proper values as the situationreliesontheexperimentalgeometry, photonandfinal-stateenergy, and atom types. As a result, the convergent process must proceed with trial tests. In the beginning, systematic calculations are conducted with varying numbers of ~ Ghkl vectors and final state partial waves (orbital angular momenta lmax). Hundreds of separate simulations using various parameter sets lead to the exploration of significant elements influencing the fine structure of CDAD. As seen in Fig. 1a, increasing the number of ~ Ghkl brings more valleys, which become more sharp, besides changing the intensity. The calculations within 49-89 ~ Ghkl are in poor agreement with the rest. It starts to converge from 137 ~ Ghkl due to the common tendency. More Kikuchi lines are introduced as a function of the number of ~ Ghkl in Fig. 2. On the other hand, under lmax variation, the position of intensity peaks and valleys relatively remain the same with respect to θangles in Fig. 1b. Taking a closer look, the shape of these peaks and valleys as well as their vicinity are smoothened, leading to smearing out diffraction patterns in Fig. 3. Convergence tests are done up to lmax =11andwe observe that qualitatively the position of the photoemission peaks as well as their positions do not change for lmax > 3, except the proximity of ϕ=4°. In an attempt to balance between accuracy and efficiency, lmax has to be truncated at 4 based on the basis of computational time and memory limitations. Fig. 1 | The SPRKKR convergence test for Si 2p 3/2 .aThe intensity computed as a function of emission angles with different numbers of~ Ghkl (from 49 to 197) with lmax ¼4. bThe intensity computed as a function of emission angles with different numbers of lmax values (from 3 to 11) with 101 ~ Ghkl. https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 2 Fig. 2 | Calculated total-intensity patterns as a function of ~ Ghkl .a–cTotal intensity for Si2p 3/2 .d–fTotalintensityfor Ge 3p 3/2 . The different numbers of~ Ghkl arementioned in the panel. Emission at the final-state energy E Final = 5880 eV and photon energy hν= 6000 eV with lmax ¼4. Fig. 3 | Calculated total-intensity patterns as a function of lmax.a–cTotal intensity for Si 2p 3/2 .d–fTotal intensity for Ge 3p 3/2 .The different numbers oflmax arementioned in the panel. Emission at the final-state energy E Final = 5880 eV and photon energy hν= 6000 eV with 37 ~ Ghkl. https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 3 Figure2showsa seriesoftotal-intensitycalculationsI TOT =I RCP +I LCP (which possess fingerprints of Kikuchi bands as mentioned in the work of Fedchenko and co-workers32,40,41) for the Si 2p 3/2 (first row) and Ge 3p 3/2 (secondrow)core levelsathν=6000eV.Onecrucialimpressionisobviously visible when enhancing the number of G ! hkl such that the diffraction patterns get increasingly complicated since each family of lattice planes creates anew“Umklapp channel" in the diffractogram. The fine structure of the pattern grows more sophisticated in the series 25, 97, and 185 G ! hkl-vectors, where Kikuchi patterns vary on a very tiny k-scale. The Kikuchi diffraction in PED is distinct from conventional electron diffraction in that the emitter atom inside the material serves as the source of the diffracted electron wave. The scattered wave is diffracted at the crystal’s lattice planes before the scattered electrons reach the detector on the vacuum side. The Kikuchi diffraction mechanism is distinguished by incoherent yet localized electron sources within a crystal. PED enables the differentiation of chemically distinct emitter sites in a crystal due to the element-specific binding energies of the photoelectrons, which is unique to the one in SEM and TEM. The photoelectron intensities are observed as a fine structure in their angular distributions when the measurements were conducted with high enough angular resolution32. In the experiment, these filigree structures have not been detected, probably due to the restricted k-resolution. Nevertheless, there is another issue to blame for this. By reducing Bragg angles, the distances of the electron trajectories in real space rise (see Fig. 3in ref. 33). The scattering cross section is big, especially for high-Z materials such as tungsten, and the Kikuchi bands get more and more suppressed when indices increase. Conversely, there are a finite number of contributing ~ Ghkl in the experiment. Methodical computations of Si 2p 3/2 (a–c) and Ge 3p 3/2 (d–f) in accordance with the diverse values of lmax are performed in terms of total intensity patterns I TOT inFig. 3.The results are obtained from hν=6000and 37~ Ghkl-vectors. In general, the main Kikuchi bands show some reduction in intensity and spread out more both horizontally and vertically. As lmax values change, the diffraction patterns from the background become more spread out. Furthermore, the contrast gradually increases, particularly for the overall intensity, making the black Kikuchi bands and bright center stand out more against the background. lmax affects the intensity of the background diffract as typically seen in thehigh-to-low magnitude variation from lmax ¼3tolmax ¼7 (Fig. 3d–f). This influence also happens for Si but not with displayed values of lmax. In cluster approaches, this convergence parameterisveryhigh(e.g.,100 22) and must be handled with care as the computation time required to calculate the scattered wave function is directly proportional to nN2ðlmax þ1Þ3,wherenrepresents the scattering order, and Nrepresents the number of atoms utilized in the cluster42. However, in ourmethod,wedo not havetosufferfromits high-value thanks to a mixed basis set of partial waves and plane waves. When using lmax ¼4, it is possible to achieve satisfactory agreement between simulations and experiments in subsection 2.3. Withtheaimoffurthercross-checkingthefeasiblelmax integer,inFig.4 we utilize a cluster approach based on multiple-scattering spherical-wave cluster basis sets from the MsSpec program18,43.Thephaseshiftshavebeen calculated by a sophisticated Hedin-Lundqvist exchange and correlation potential44–46 so as to describe the finite electron mean free path in the final stateviaimaginaryparts.Theincidentaswellas emissionanglesandphoton source are determined identically to experiments and the one-step model setup. We use unpolarized light and the non-relativistic mode. Being interested in high-energy PED, the deviations from bulk physics owned by surface atoms are sparse47. As a result, we take into account bulk-terminated surfaces. The radius of our spherical cluster is set to 27.15 Å(Fig. 4a). Here, we investigate the test of cross-section, an intensity-related quantity, by tuning lmax values (4, 16, and 24) in the azimuthal scan of Si 2p 3/2 at θ=5°. Apparently, there is a gap among 3 cross-section behaviors, such as magnitude overall and the shape of the curve (especially in the vicinity of ϕ= 30° and ϕ= 60°). The result from lmax ¼4 seems a bit far from convergence as compared to the other two. Even so, its peaks and valleys share good relative positions (namely, in the surrounding of ϕ= 30° and ϕ= 60°) in common with the rest. More examples are addressed when ϕis in between 50° and 60°. Final-state energy dependence The final-state energy E Final is defined as the kinetic energy of the photoelectrons inside the crystal and the quantity in charge of diffraction dynamics. The expected variation amongst calculated patterns at individual final-state energies [4690, 6690] eV is seen in Fig. 5for Si (top row) and Ge (bottom row), respectively. From the simulation point of view, when energy is raised, a pronounced system of Kikuchi bands dominates the diffraction patterns, forming a rich fine structure. Each band consists of a projection plane sandwiched by two lines and its width is determined by the respective reciprocal lattice vector. For consistency in comparison, the entire calculations are designed by the identical value of momentum range and resolution. Overall, there is an agreement concerning the width of observed bands and the position of lines. Fine structures possess a horizontal and a vertical mirror plane that is parallel to the [010] and [001] directions, respectively. Another general feature is that the distinct central zone indicated by a green dash square, can be noticed easily regardless of final state energies. Figure5a–dshowscomputeddiffractogramsofSi2p 3/2 core-levelatthe (100) plane with 161 ~ Ghkl and lmax ¼4. Via visual inspection, there are several noticeably rapid alterations besides the above-mentioned common Fig. 4 | The convergence test done by the MsSpec package. a A clustermodelofSi (100). Theredcircle depicts the emitter.bThecross-sectioncalculatedasa function oflmax values for Si 2p 3/2 . https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 4 behaviors. Firstly, the texture inside the central zone evolves by enhancing E Final .Fourvarioussymbolsarefoundinthesurroundingsofthecenterzone axis [100], running from left to right: bow-tie dartboard (a), shield (b), 4-pointed star cross (c), and butterfly (d). However, all geometrical patterns still follow 4-fold symmetry, i.e., the dartboard pattern is obtained by rotating a bow-tie by 90°. Another shape turn can be noticed through the distance decrease of a line pair (marked as a yellow dash with a label (1)). Startingfrom4690 eV,it isdefinite toseethespacebetweenthem. Then,this gap continuously reduces in between (b-c) and seems to become 0 at 6690 eV. Consequently, some bright spots (e.g., four magnificent bright ones labeled (2) in (b)) become hidden. One more typical instance is the deformation of the bright area (orange dotted curves labeled 3 in (a)) along diagonal directions at four corners. These phenomena are derived from the bound relation between the angular range and Bragg angles. In the context where forward scattering appears overriding, the greater the final-state energy is, the smaller the Bragg angle is. Thus, there are some shifts and rearrangementsoffineKikuchilinesmeanwhilemainfeatures(horizontaland vertical bands as well as their intersections) remain unchanged. Analogously, the diffractogram of Ge 3p 3/2 is analyzed in Fig. 5e–h. Because of the mismatch related to photoelectron wavelength (caused by different core-level binding energies), the patterns are markedly dissimilar from those in Fig. 5a–d. The number of ~ Ghkl (113 in this case) is also one of the essential causes of reshaping patterns (as discussed in Section 2.1). The major Kikuchi grid is recognizable and aligned with the one from Si. For consistency and convenience in qualitatively comparing with Si outcomes, the colorbar is fixed. Hence, plots are not in good contrast and it is a bit hard to identify diffraction which is contributed by thin Kikuchi lines. Despite that, a striking difference regarding global intensity values is captured when bright areas are larger in general. As the energy rises, pattern adjustment is forcibly determined by the diffraction angle. As an illustration, the angle between two narrow fine lines which are in yellow dash (e) declines and its vertex position moves up when switching from 4690 eV to 5190 eV. Continuouslywith theenhancement ofenergyto 6190eV, thebrightspot tagged by (4) in (f) is fading. In Kikuchi photoelectron diffraction, there exists an orientation of higher and lower intensity. This bright-dark distribution can be interpreted in terms of the reciprocity principle and the coupling probability of photoelectron from localized emitters with outgoing wave22. Experiments meet theories Figure 6displays a comparison of measured and calculated CDAD in pairs for the Si 2p 3/2 . The photoelectron diffraction pattern is simulated over a polar angle range of 0 ° ≤θ≤10° and over a full 360° azimuthal range (ϕ) with 180 points for both angles at thephoton energy hν= 6 keV. The pattern centers are indicated by +. Band edges are marked by full arrows on the lefthand side of the figure, and band centers are marked by dotted arrows. The CDAD signal (CDAD = I RCP -I LCP ) is a difference between the intensities of right and left circular-polarized light that emphasizes faint details in the patterns. When this difference is normalized, it leads to a so-called CDAD asymmetry A CDAD =(I RCP (k x ,k y )-I LCP (k x ,k y ))/(I RCP (k x ,k y )+I LCP (k x ,k y )) reflecting the symmetry behavior of CDAD. This normalized factor is a good choice for experiment-theory comparisons due to its autonomy of the spectrometer transmission function and the free-atom differential photoelectric cross-section. As anticipated, the A CDAD is antisymmetric referring to the horizontal mirror plane. This “up-down antisymmetry" is an explicit consequence of the atomic CDAD’s symmetry (see Fig. 1in ref. 33)andthe Kikuchi diffraction’s characteristics (see Fig. 2in ref. 33). The crystal-lattice mirror plane is positioned horizontally to maintain the antisymmetry. The agreement between observed and computed intensity (a-b), CDAD difference (c,d), and A CDAD (e,f) looks quite reasonable. Particularly, for the intensities I RCP +I LCP the principal Kikuchi bands (as labeled by black arrows) and the 4-pointed star cross clearly show up in Fig. 6a, b. Two elliptic shapes (marked as green) are also captured. Compared to the measured total intensity, the contrast is lower and there are more visible minor Kikuchi bands from calculations. The reason for these differences originates from the lmax effect (above-mentioned in Fig. 3a–c). The issue is solvable by increasing lmax but for the sake of clearly displaying fine structures in detail, the small converged value is opted. It is likely that the CDAD image (Fig. 6c) shares commonly observed patterns with the total intensity (Fig. 6a), such as a large diamond shape around the center and two ellipses (markedasgreen).Thesesimilaritiesarealso nicelycapturedbytheone-step Fig. 5 | Sequence of calculated total intensity as a function of final-state energies. a–dThe diffractogram of Si 2p 3/2 .e–hThe diffractogram of Ge 3p 3/2 . Calculations areperformed with161~ Ghkl forSi and113~ Ghkl forGeatlmax ¼4andhν=6keV. The final-state energy E Final is between 4.69 and 6.69 keV. All plots are made at the same color scale. Green dash squares indicate the central zone of studied structures. The zone axis is [100], directly out of the page. https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 5 photoemission model (Fig. 6d). The CDAD difference I RCP −I LCP depicts a higher contrast than the sum of intensities I RCP +I LCP . However, the 4-pointed star cross, surrounding the CDAD center, is practically hidden in the experimental image (c). From the computational point of view, this footprint is reproduced charmingly in Fig. 6d. Several band edges are present in agreement as well, in particular two crossing points (yellow arrows in (c-d)). Likewise, the latter emerges from the A CDAD patterns as the vertex of small blue and red triangles, identified by the black arrows in (e,f). There is no doubt that the resolution of A CDAD from measurements is not in good quality in comparison with the two above-mentioned ones. Overall, the consensus between the experimental and simulated A CDAD is far from perfect. But the mirror plane and some characteristics (namely, Kikuchi bands forming a central diamond) still emerge. Calculated spectra have sharp features that need to be broadened to match the momentum resolution ofthe experiment.To maketheoreticalcalculations morecomparable to experimental data, we apply convolutions with a Gaussian function in which the standard deviation σ= 2 (Fig. 6b*,d *,f*.ForFig.6b*,thePeronaMalik (anisotropic diffusion) filter is employed for I TOT to smooth the resultswhilepreservingthemain spectralfeatures.The“weight"parameterλ is set to 1 in this case. In general, convolution steps significantlyimprove the visual comparison. To further clarify the parallels between measuring and calculating, we conducted the same measurement as shown in Fig. 6abutwithhigher resolution in Fig. 7a.Next,weexaminethespecifics from a quadrant of A CDAD in the second row of Fig. 7. The superior Kikuchi bands are marked with dashed lines in the experimented intensity (a) and A CDAD (c). These line patterns are then mapped to the corresponding calculated ones (b, d). The Kikuchi-band edges and crossing points take place in a well-founded arrangement in the matching experiment (a) and simulation (b). Dark (addressed by a dashed arrow) and bright (pointed by a solid arrow) spots located in the corner of the diamond shape are observed in the relevant calculation (Fig. 7b). The computed and measured patterns of the total intensity are arranged in such a way that features can be reflected via a mirror plane. The rough red-blue texture of experimental A CDAD is analogous to the theoretical one with respect of not only color but also relative positions. Next, we turn our attention to the Ge(100) cases in Fig. 8which can be interpreted by the formalism like in Si(100) diffraction. In Fig. 8a, b, the intensity from two light helicities are summed for Ge 2p 3/2 . An inset of a higher-resolution measurement is added to the experimental pattern. The relevant area in the simulation is depicted via a blue dashed circle. Figure 8c, d shows the difference between these helicities while Fig. 8e, f belongs the CDAD asymmetry. The “four-fold" geometry marked with blue curves in the intensity difference indicates the matching. The grid made up of Kikuchi lines (dash dots) is superimposed on the diffractogram to conveniently clarify the parallels between measuring and calculating A CDAD .Notable features are easy to map between experimental and theoretical interrelationships: vertical (except CDAD patterns, as they are almost invisible) and horizontal Kikuchi bands and the pattern at four corners (indicated by green circles). More specifically, the color flipping of the pattern from upper to lower corners in Fig. 8c–f is triggered by the horizontal mirror plane which is a typical feature between CDAD and A CDAD . Besides, black and white spots of the intensity difference (blue and red ones in case of A CDAD ) formed by several visible band edges are sighted by green arrows. Although Fig. 6 | Comparison between measured and calculated patterns of Si 2p 3/2 at hν=6 keV. The top, middle, and bottom rows indicate the total intensity I TOT , intensity difference CDAD and CDAD asymmetry A CDAD . From left to right, there are measured (a,c,e), convoluted-calculated (b*,d*,f*) and purely calculated patterns (b,d,f). Computational results are performed at lmax ¼4 with 193 ~ Ghkl. To smooth the computed data, convolution is applied by the Gaussian filter with the standard deviation (σ) and the Perona-Malik filter (anisotropic diffusion) with the “weight" parameter (λ). https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 6 recognizing several systematic hallmarks, it seems we are far from reaching good consensus in comparison. From measured outcomes (especially in Fig.8c,thereisagenericdiamond-shapedpatternwhichiswidewithasharp rim. Nevertheless, it appears quite differentiable in simulation. Another distinguished discrepancy is the distorted version of the bow-shaped feature 1 (yellow dash in Fig. 8c, e). Also, we cannot ignore the fact that in CDAD andA CDAD fine diffraction attheupperandlower cornersisnot inharmony with respect to fine details. To enhance the comparability between theoretical calculations and experimental data, convolutions are applied with a Gaussian function characterized by a standard deviation of σas detailed in Fig. 8b*,d *,f*. Subsequently, we present the diffraction pattern of another core level, Ge 3p 3/2 ,asshowninFig.9. The similarities between measured and calculated diffractograms are easily found. For example, the Kikuchi lines, characterized by their very low intensity, are detected (denoted by black solid arrows in Fig. 9a, b). In addition, the bright fourfold cross (indicated by green dashed arrows) is reproduced though its contrast to the background is notashighasinexperiments.Theinner diamond shape surrounding the center of this cross also appears in the simulation. The nice agreement is further confirmed by considering the signal subtraction (Fig. 9c, d). The faint patterns are likely symmetrized through a vertical mirror plane separating two diffractograms. Obviously, the edges of the big diamond can be recognized via blue dashes. In pursuit of fleshing out the similarities between measured and calculated patterns, it is worth looking at their A CDAD in depth (Fig. 9e, f). The agreement is not ideal overall. However, the mirror plane and certain features (specifically the formation of the “redblue"hourglasses insideand outsidethe bigdiamond)stillbecome apparent. Taken together, the main findings indicate a greater result of Si2p 3/2 and Ge 3p 3/2 comparedtothe oneof Ge2p 3/2 .Despite owningthesame information in terms of crystal structure and period in the periodic table, it is understandable to experience the difference between their diffraction patterns due to dissimilar electron configurations and investigated core levels. To improve the consistency between theory and experiment, the calculated spectra are convolved with a Gaussian function of standard deviation σ=2, as illustrated in Fig. 9b*,d *,f *. This approach enhances the agreement between theory and experiment. It is distinguished between the instrument resolution and the width of the features in a specific k pattern. The former refers to the intrinsic performancelimitoftheexperimentalsetupandhasbeendeterminedtobe0.03 Å−1. The latter is the width of the features observed in a specific measurement or image, which is affected by several factors. This width can be much worsethan theinstrument resolutionwhen thesample qualityisnotgoodor the X-ray beam hits a “bad spot" on the sample. After cleaving, it is often complicated to find a good spot. In some cases, there may also be some misalignment, e.g., incorrect sample distance, or the photon beam is not optimally adjusted (e.g., too large a footprint, which also spoils resolution). Inmanycases, the statisticsarenot sufficient,so wehave to applya Gaussian blur. This is analogous to the width of a core-level signal in comparison with the resolution of an electron spectrometer. Fig. 8a and the inset of Fig. 9aare shown with better resolution and contrast because the corresponding measurements were done after we improved the setup and conducted the experiment in a smaller k-range. The rest of the figures come from the initial measurements. The instrument resolution is at least as good as the narrowest lines we observe. When lines appear broader, they may arise from several other factors, as discussed above. It is important to note that not the eye-catching broadest features reflect the instrument resolution, but the narrowest lines (which may be rather weak in intensity) in the pattern. While the overall agreement between theory and experiment remains qualitatively reasonable, we acknowledge that some quantitative discrepancies persist. From our perspective, several factors may contribute to Fig. 7 | Comparison between measured and calculated patterns of Si 2p 3/2 at hν=6 keV. aMeasured total intensity pattern I TOT by higher resolution, compared to Fig. 6a. bTotal intensity calculation by OSM. c,dCDAD asymmetry displayed by the quadrant of Fig. 6(e,f). Computational results are performed with 193 ~ Ghkl. https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 7 these differences. First, the choice of effective potential and exchangecorrelation functional plays a crucial role. We currently use the local density approximation(LDA)to describe the exchange-correlationpotential,which may not be optimal. Previous studies haveshown thatthe choiceof potential can strongly influence theoretical predictions. For example, McGovern et al.48 reported discrepancies arising from a suboptimal Te potential, which were later addressed by Wendin49 using the random phase approximation with exchange (RPAE)50, yielding more accurate atomic photoionization cross-sections for Te 4d. Additionally, while the influence of the exchangecorrelation functional is often minor in small systems51, its impact becomes more pronounced in complex materials52. The LDA is known to overbind, potentially leading to systematic deviations53,54. In our calculations, we use the atomic spheres approximation (ASA), which is similar to the muffin-tin approach but ensures full space-filling by allowing overlapping spheres and introducing empty spheres in the interstitial regions. However, ASA imposes shape-related constraints on the potential, and its accuracy can be improved by adopting the full-potential version of the SPR-KKR code55. Second, inelastic scattering effects may contribute to the observed differences. Inelastic scattering corrections to the elastic photocurrent (see, e.g., Borstel56 and Braun et al.57) are included via a parametrized, complex inner potential. While these effects are considered in the current work, they have not been systematically investigated. A more detailed study of their influence will be presented in a forthcoming publication58. Third, lattice vibrationsmustalso betakenintoaccount.Theexperimentalmeasurementswere conducted at 30 K, while our calculations were performed at 0 K. Although this temperature is low, vibrational damping effects are still relevant. We plan to analyze this temperature dependence more systematically using the Debye-Waller model and the alloy analogy approach within the coherent potential approximation (CPA) in a future publication aimed at reproducingthe SihardX-ray experiments34. Finally,crystalimperfectionsrepresent another possible source of discrepancy. Our theoretical simulations assume a perfect, defect-free crystal structure. However, real experimental samples maycontain hidden structural imperfectionssuch as local strain ordisorder, which are known to significantly affect photoelectron diffraction patterns3. These effects are not captured in our current model but are important to consider in explaining the remaining discrepancies. Discussion The combination of CDAD and XPD techniques provides comprehensive information on the crystal structure and chemical composition of the system being studied. Most of the time, it is rarely accessible to obtain all of this information readily. The challenges in interpreting emitted-electron spectra evolve from the intricacy of the photoemission itself. To extract it, one must make assumptions about a particular physical model and apply the relevant theoretical framework. Even so, a good qualitative analysis demands ab initio computations, and this process may include several levels of complexity. In this study, we concisely outline the interplay and importance of XPDand CDAD.Furthermore, severalstandardtheoreticalmodelsfor PED are reviewed with an emphasis on physical content and terminology rather than concentrating on mathematical formalism. In the context of high photon energies, we emphasized our fingerprint that can tackle the disadvantages of other codes. XPD diffraction computed in reciprocal space is Fig. 8 | Comparison between measured and calculated patterns of Ge 2p 3/2 at hν=6 keV. The top, middle, and bottom rows indicate the total intensity I TOT , intensity difference CDAD and CDAD asymmetry A CDAD . From left to right, there are measured (a,c,e), convoluted-calculated (b*,d*,f*) and purely calculated patterns (b,d,f). Computational results are performed with 177 ~ Ghkl.To smooth the computed data, a convolution is applied by the Gaussian filter with the standard deviation σ. https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 8 convenient to compare with momentum microscope measurements. The atomic potential is obtained from our own code, making the XPD results more controllable in a closed simulation loop. The observed outcomes were reproduced in the hard X-ray regime. The upcoming Ge work concerning the full-field photoelectron diffraction and circular dichroism texture will complete the energy range of applications.CDAD signals and asymmetries in the manner of CL-XPD were simulated at 6 keV. The diffraction patterns are significantly altered by the Kikuchi diffraction process. Furthermore, simulating Kikuchi patterns in spectroscopy encourages additional measurements and energizes the current trend in XPD59. The representation of scattering by our approach requires smaller values for the maximum angular momentum than others (e.g., cluster approaches), helping to avoid computationally intensive calculations. A series of convergence parameters (the number of ~ Ghkl and lmax ¼4) was tested. Next, the optimal values were selected for the calculations for Si and Ge core levels, yielding fairly nice harmony with experiments. Fine patterns are well-mapped from the total intensity and CDAD signal. The final-state effect was systematically investigated and had a strong influence on the pattern evolution. Thus, careful selection and control of the kinetic energy are essential for optimizing photoelectron diffraction studies. With this endeavor, we aim to bring another tool for XPD analysis and CDAD interpretation. Methods Brief review of PED simulation progress Basically, there exist three classifications for adequate PED simulation for solids as presented in Fig. 10: cluster-based models, layer-by-layer approaches, and the so-called “lattice-plane”methods. The first group(Fig.10a)is arightcall forshort-range probes (e.g.,PED or Auger electron diffraction) and broken crystal symmetry (for instance, by disorders). Because of the electron’s short inelastic scattering length, scattering can be limited to a cluster of individual atoms surrounded by the photo-emitter.Amodel canbeformulated by either single-scattering cluster (SSC)60–66 or multiple-scattering cluster (MSC)13,67,68 methods. For comparison in this subset, the most precise interpretation is MSC with spherical waves14. Specifically, the multiple-scattering effects have to be taken into account while investigating single-crystal because the XPD patterns originate from photoemission of ten (and more) surface layers. Inspired by the work of Kaduwela and co-workers14 where multiple scattering (MS) is derived from the Rehr-Albers (R-A) separable propagator approximation69 and applied for initial and final states, more CD measurements are nicely reproduced70,71. Furthermore, the experiment-theory match is improved by utilizing an accurate representation of the Green’s function propagator for evaluating MS expansion24,72. Thereafter, more indispensable modifications and implementations come into existence as seen in the following MSC packages: MSCD16,PAD273, MSPHD74,EDAC24 and MsSpec18. Recently, Rehr et al.75 has discussed the development of real-space Green’sfunction approach for PED based on the R-A formalism and boosted the importance of sharing theoretical software among groups. Workflow tools such as AiiDA76 and Corvus77 were proposed to enable advanced and efficient computation of PED without requiring significant changes to the original code. Besides, there is another cluster code called TMSP15,78, which relies on photoelectron holography (PEH), a rapidly growing and interconnected field of PED. In TMSP, multiple scattering is fully employed, but Fig. 9 | Comparison between measured and calculated patterns of Ge 3p 3/2 at hν=6 keV. The top, middle, and bottom rows indicate the total intensity I TOT , intensity difference CDAD and CDAD asymmetry A CDAD . From left to right, there are measured (a,c,e), convoluted-calculated (b*,d*,f*) and purely calculated patterns (b,d,f). Computational results are performed with 241 ~ Ghkl.To smooth the computed data, convolution is applied by the Gaussian filter with the standard deviation sigma (σ). https://doi.org/10.1038/s41524-025-01653-y Article npj Computational Materials | (2025) 11:159 9