Background

The technique of core-level photoelectron diffraction (PD) has been applied to surface structure determination on the atomic scale for more than 20 years. This local diffraction technique probes short-range order around the photoemitter. A broad variety of surfaces have been successfully studied, including metals, semiconductors, oxides, systems exhibiting surface core-level shifts, adsorbed atoms and molecules, epitaxial overlayers, and atoms at buried interfaces. A number of research groups have performed photoelectron diffraction experiments to study surface and interface structures, using both laboratory X-ray and synchrotron radiation sources, and doing the measurements in both scanned-angle and scanned-energy modes. Several reviews of this field have appeared in recent years. A further element which has recently been added to this technique is the application of holographic inversion methods to photoelectron diffraction data so as to more directly yield atomic structures in three dimensions. Thus, the accurate theoretical modeling of photoelectron diffraction is crucial to the development and use of this technique in its various forms, and we here discuss one method for accomplishing this.

In core-level photoemission, a photon illuminates an emitting atom at or near the surface or in the solid and excites an electron from an atomic core-level, ejecting the electron to a detector far away from the surface. The resulting photoelectron wave components can occur both by direct propagation to the far-field detector and via a number of scatterings from atoms in the neighborhood of the emitter. The quantum interference (diffraction) between the different pathways depends sensitively on the relative atomic positions. The intensity associated with this process, when measured as a function of photon energy and/or emission angles, thus includes information about the atomic structure around the emitter. The energy- and angle-resolved photoemission intensity I(k,q, f) at the detector can be written in general as

(1)

where k is the final electron wave vector (related to the final kinetic energy), q and f are the polar and azimuthal angles of photoelectron emission, respectively, f₀ is the wave-component representing travel along a path directly to the detector without being scattered by another atom, and f_sj is the wave-component representing travel via paths involving single or multiple scattering by one or more atoms, respectively. The multiple-scattering order is defined as the number of scattering atoms in a path: j = 1 = single, 2 = double, etc. Higher-order scattering processes are less important because of the roughly 1/r falloff of the outgoing distorted spherical wave represented by f₀, damping due to inelastic processes and loss of diffraction modulation due to thermal vibrations (Debye-Waller effects). The explicit sum in Eq. (1) includes all possible scattering paths inside the solid.

To briefly review the history of such theoretical modeling of PD effects, the first qualitative analyses of higher-energy data made use of Kikuchi-band theory. An early quantitative theoretical explanation of these diffraction data was provided by Liebsch in 1974 and then improved by him and others with multiple-scattering effects in 1976; it was based on LEED theory. In the next few years, Pendry, as well as Li, Lubinsky and Tong, put forward similar theories involving the more complex nature of both the initial and final states. These theories require the assumption of full translational symmetry parallel to the surface for the system under investigation.

A single-scattering cluster-based theory based on prior work in extended x-ray absorption fine structure (EXAFS) and preliminary work on Auger electron diffraction (AED) was first applied to scanned-angle photoelectron diffraction by Kono et al; this was based on the plane-wave approximation. The next major advance in the cluster approach was by Barton and Shirley, who included both spherical-wave corrections and multiple scattering effects. The cluster-based theories are inherently more suitable for photoelectron diffraction modeling in view of the point-source nature of the problem, the spherical outgoing waves involved, the short inelastic attenuation lengths, and the sensitivity to short-range rather than long-range order. In particular, one key advantage of PD as a structure probe is in not requiring long-range translational order, and so being able to model such short-range order structures conveniently without having to resort to some sort of repeated structure with fictitious long-range order is crucial. Cluster-based methods are ideal for this. It has also been found that the maximum cluster size required for accurate simulations is about 100 atoms.

The multiple scattering expansion for spherically-symmetric scatterers is conveniently formulated in terms of diagonal plane-wave scattering t-matrices with elements t_l = sin(d_l)exp(i d_l) and matrix elements of the free-particle propagator G_L,L'(r) = <L,R|G|L',R'> in an angular momentum L = (l,m) and site R basis. Here r = k(R-R') denotes the dimensionless vector between two sites, k = |k| is the wave number. G_LL'(r) is defined by the following integral involving spherical Bessel functions j_l(kr) and spherical harmonics Y_L():

(2)

where is a unit vector along k, r and r' are arbitrary displacements. An (n-1)-atom path (Fig. 1), including the emitter atom at R₀, n-1 scatterers at R₁, R₂, ..., R_n-1, and the detector at R_n, is represented by the following total propagator in this exact multiple-scattering expansion:

(3)

where r_i = k(R_i+1-R_i), and L₀ and L_n denote fixed initial and final angular momenta.

Fig. 1 A photoelectron excited from the core-level of an emitter atom at R0 propagates to the detector at Rn via n-1 scatterers R1, R2, ..., and Rn-1. Each straight segment linking consecutive atoms represents a Green's function propagator GLL'(k(Ri+1-Ri)). The emitter can become a scatterer after emission.

Several computer programs have been written to successfully calculate photoelectron diffraction intensities based on this exact cluster formalism, for example, by Chassé and co-workers^, and by García de Abajo. However, without special optimization (as we consider below) such exact programs in general require considerable computing times for large clusters, due to the large matrix dimensions involved and the large number of paths that are explicitly summed over. That is, the matrix G_LL' has dimensions (l_max+1)² x (l_max+1)², where l_max is the largest angular momentum value needed to accurately represent scattering; it can be roughly evaluated using l_max ~ k_maxR_mt, where R_mt is the muffin-tin radius of the scattering potential. For photoelectron diffraction in the energy range 50-1500 eV, this l_max lies between 5 and 30, respectively, of course depending on the electron energy. The computing time can be considerably reduced by an iterative summing over the multiple-scattering paths, as proposed by García de Abajo, which leads to proportionality to the number of multiple scatterings, rather than the much steeper behavior as in the traditional summing method. More precisely, the total calculation time in traditional summing is T µ T₂ ^nmax/2, while the iterative summing has a total time of T µ n_maxT₂, where T₂ is the time for a double scattering calculation, and n_max is the maximum multiple-scattering order included. A convergent calculation is found to requires n_max up to about 6, or even higher in cases with long straight chains for strongly forward scattering atoms. We use these two exact formalisms here as references with which to compare calculations performed with the Rehr-Albers expansion.

Several methods have been proposed previously for approximating the exact scattering formalism so as to reduce computation times in cluster calculations.

• Barton and Shirley first derived a representation of the exact formalism with their Taylor-series magnetic-quantum-number expansion (MQNE). This reduces the matrix dimension that needs to be used, since only the first 3 Taylor orders are usually required.
• Fritzsche, Rennert and Chassé provided a reduced angular momentum expansion (RAME) approximation. The incoming spherical waves are approximated by a limited set of spherical harmonics with low angular momenta (quantum numbers (l,m) = (0,0), (1,-1), (1,0) and (1,1)). Improvement has been achieved by taking into account two more components of higher angular momenta (l,m) = (2,0) and (3,0). A generalized scattering amplitude of each scattering process contains curved wave corrections and a rotation operator which rotates the z-axis between the wave vectors before and after the scattering.
• As another approach to reduce computation times, Rehr and Albers developed a theory, denoted R-A in this paper, that is based on a separable representation of the exact formalism; it was developed primarily to calculate NEXAFS (near-edge x-ray absorption fine structure) and EXAFS (extended x-ray absorption edge fine structure). Curved-wave multiple-scattering contributions can in this approach be calculated with an efficient formalism similar to that based on the plane-wave approximation, but with scattering amplitudes replaced by distance-dependent scattering matrices that act like effective scattering factors. A useful property of this method is that the size of the scattering matrices can be systematically increased as needed for sufficient accuracy, ultimately recovering the exact result. Thereby, the successive orders of scattering can be built up and convergence achieved in a convenient and efficient way. The first application of this theory to photoelectron diffraction was in a computer program developed by Kaduwela, Friedman, and Fadley, and this was based on the second order of the R-A approximation, i.e. scattering matrices of dimension (6x6).
• In a more recent development based on the R-A approach, Wu, Chen and Shirley developed a faster program that benefits from a novel reverse-summation method over paths and an iterative summing technique, but is again limited to second order in R-A. Chen and Van Hove have further developed this method to enable the inclusion up to 4th order in R-A (up to (15x15) matrices). By an efficient path-cut process, this program runs even faster. The resulting computer package, called MSCD, is also very portable, running on a variety of sequential and parallel computers.

By applying object oriented programming within C++, the MSCD code provides a better user interface (particularly important for convenient use by experimentalists), general data formats, and fully functional subsidiary utilities. The simulation program is also designed to be very portable, with very simple changes from one computer type to another. The sequential version of the code has been implemented on supercomputers (under Unix), on Sun workstations (Unix), as well as on PC and Macintosh desktop computers. The parallel version has been implemented on a massively-parallel supercomputer (Cray T3E) and on a cluster of multiprocessor systems (COMPS, a type of network of workstations) under MPI.

Comprehensive test calculations agree well with more exact (but much slower) calculations performed by other programs. A new version is currently being generalized to take into account photon and electron-spin polarization effects such as circular dichroism (the latter is obtained by contrasting measurements performed with right- vs. left-circular polarized light), magnetism (including magnetic circular dichroism), and relativistic effects such as spin flip.

In the next sections, we will briefly introduce both the exact and R-A formalisms, discuss the methods that have been used to increase the calculational efficiency of R-A, present a number of results calculated with the MSCD program based on R-A, and compare these results with exact cluster calculations to quantitatively evaluate the performance and limitations of the R-A approximation as carried out to different orders. We will discuss a few more necessary details about the exact cluster-based formalism, and briefly review the R-A method and its properties. The methodology and convergence of the implementation of R-A used in the MCSD will be presented, as well as the overall reliability of the MSCD approach for structural determinations.

Table of some photoelectron diffraction simulation programs (as of 8/1/98)

The programming of the MSCD package began on December 9, 1996, using the ANSI C++ object oriented programming (OOP) technique and message passing interfacing (MPI) parallel communication.

• Version 1.37, June 10, 1998
  Make the source code compatible to Macintosh CodeWarrier Compiler.
  Fixed a bug in net search mode.
• Version 1.36, June 1, 1998
  Impoved emission angle search feature.
• Version 1.35, April 21, 1998
  Number of emitters in different logical layers is considered. Total number of data points can be up to 30000.
• Version 1.34, February 22, 1998
  Fixed a bug in calculation of large cluster.
• Version 1.32, January 16, 1997
  Fixed a bug in calculation of reliability factors for two hologram data files. Calculate Fadley and Van Hove group's reliability factors in program caldif
• Version 1.31, August 31, 1997
  Expand the Rehr-Albers approximation order from 3 (10x10 matrix) to 4 (15x15 matrix). To choose 4th order of Rehr-Albers approximation, set raorder=4 in the input file.
• Version 1.30, August 30, 1997
  Expand the Rehr-Albers approximation order from 2 (6x6 matrix) to 3 (10x10 matrix). To choose 3rd order of R-A approximation, set raorder=3 in the input file.
• Version 1.24, June 30, 1997
  1) Add usermac.cpp file. The program works on Macintosh computer. 2) Add two functions in userutil.cpp, skipendline and getendline, making the program package support PC, Unix and Macintosh ascii text format data file. 3) Add usercomp.cpp file. The program works on COMPS network of workstation in parallel mode. 4) Add cpu time analysis function for each processor in parallel mode. The processor times are divided into four categories: computation, sending, receiving and idle. To activate this feature, set display mode to 10. The total processor time and its distribution are shown in output file mscdlist.txt.
• Version 1.23, June 12, 1997
  The maxmum number lmax of angular momentum is expanded again to 60. To make it possible, the makecurve function in rotamat.cpp has been changed, using double precission floating point numbers.
• Version 1.22, June 10, 1997
  The maximum energy and number of angular momentum are expanded. The maximum kmax set to 25.0, corresponding to energy 2.38 keV. The maximum lmax set to 30.
• Version 1.21, June 3, 1997
  Removed a bug in fitting mode and a possible bug in test-running mode.
• Version 1.20, May 28, 1997
  First parallel version, working on Cray T3E in MPP parallel mode.
• Version 1.10, April 1, 1997
  Introduce the rotation transfer method to save computation time for circular polarized light and magnetic circular dichroism.
• Version 1.00, March 1, 1997
  First version, working on PC, Sun workstation, Cray J90 supercomputer, and Cray T3E supercomputer in sequential mode, covering all the features of the SCAT Version 3.65 package.

We thank Dr. H. Wu, Dr. Z. Hussain and D.A. Shirley for their collaboration in developing the SCAT package. The MSCD package uses the same method and similar data format, and this documentation takes some contents from the SCAT package user guide.

We thank our colleagues in D.A. Shirley's group at Lawrence Berkeley National Laboratory, E. Moler, S. Kellar, W.A. Huff, X. Zhou, and B. Peterson, for their detailed discussions and many test calculations. Special thanks also go to C.S. Fadley, J.J. Rehr and A. Kaduwela for many helpful discussions regarding photoelectron diffraction theory and its computation. R.X. Ynzunza, F.J. García de Abajo and A. Chassé ran many time consuming test calculations and comparisons with other codes, which helped us to debug and improve the software.