## research papers

## X-ray absorption spectroscopy: state-of-the-art analysis

^{a}Laboratori Nazionali di Frascati, INFN, CP13, 00044 Frascati, Italy, and ^{b}Universita' dell'Aquila, Via Vetoio, loc. Coppito II, 67100 L'Aquila, Italy^{*}Correspondence e-mail: natoli@lnf.infn.it

State-of-the-art techniques for analysing X-ray absorption spectra are reviewed, with an eye to biological applications. Recent attempts to perform full spectral fitting of the XANES energy region and beyond for the purpose of structural analysis have met with encouraging success. The present paper analyses the theoretical motivations behind this success and indicates routes for future improvements. The theoretical background is not entirely new, although the point of view is, and some sections and appendices present material that the authors believe has never been published before. The aim of this paper is to provide a theoretical analysis that is as self-contained as possible.

Keywords: X-ray absorption spectra; XANES; biological applications.

### 1. Introduction

In the past 20 years, *via* subtraction of an approximate atomic background can be performed with a certain confidence and reliability only in the region, whereas the low-energy part cannot be adequately background subtracted.

There exist several analysis packages based on spherical-wave multiple-scattering theory (MST) and the complex optical potential of the Hedin–Lundqvist type (Hedin & Lundqvist, 1969, 1971). These packages are able to reproduce satisfactorily the signal [for a review see Rehr & Albers (2000)]. The same codes even offer the possibility of fitting the entire spectrum; however, it does not seem that this potential has been seriously exploited as yet. Some of the factors that have probably deterred the various practitioners of MST from extending the fitting procedure used in the region to the edge region, *via* the calculation of the total are well known: the inadequacy of the muffin-tin approximation to the potential at low photoelectron energies, the lack of a satisfactory description of the screening and relaxation processes following the sudden creation of the core hole, and the need to include electronic correlations and, in particular, two-electron excitations (since they change the slope of the background absorption). Therefore, the analysis of the near-edge region has remained at a semi-quantitative level.

Nevertheless, quantitative analysis of the X-ray near-edge structure (XANES) in order to obtain structural and electronic information can be very relevant in many fields of scientific application, like extra-dilute systems, surface spectra, real-time measurements of dynamic systems, trace-element analysis, the local investigation of materials under extreme conditions and especially biological systems (enzymes), where the low S/N ratio and the weak scattering power of the light elements constituting the organic material limit the *k*-range of the available experimental data. In all these instances, the region of the spectrum cannot be adequately exploited, since the usable data are very often below 200 eV from the absorption edge.

Very recently, a method for performing full spectral fitting of the XANES energy region and beyond (up to ∼200 eV above the edge) for the purpose of structural analysis has been proposed and applied to the *K*-edge spectra of a number of transition metal compounds, both organic and inorganic (Benfatto & Della Longa, 2001; Della Longa *et al.*, 2001). The encouraging success of this attempt requires a reconsideration of the theory with a view to understanding why the method works and highlighting areas of possible improvement.

The present formulation reflects our long experience of the performance of multiple-scattering theory with a complex optical potential and might not duly acknowledge other groups' contributions to the subject. We apologize if this is the case.

### 2. Theoretical background

In this section we shall present the derivation of the photoemission *k*^{2} along the direction and illustrate the reduction of this many-body problem to an effective one-particle problem with complex energy-dependent optical potential. This process of reduction will help us to understand the validity of the necessary approximations to the optical potential and will give us guidance in choosing among various approximation schemes. The photoabsorption is nothing more than the integration of the photoemission over all the emission angles and all the final channels (elastic plus inelastic) with the same final energy. The reason we treat both cases together is duplex. Firstly, the mathematical formalism is the same, as is the reduction process to a one-particle problem; secondly, and more importantly, on purely physical grounds we can think of photoabsorption as a kind of photoemission that has the same electron source (the photoabsorber), in which the detector, instead of being outside the measured system, coincides with the source. The validity of this assumption will be apparent from the mathematical formalism in the following. By treating photoabsorption and photoemission together, we can judge the sensitivity to structural details of a particular potential in absorption by looking at its performance in photoelectron diffraction. In photoabsorption, because of the obvious impossibility of controlling the measured variables (except total energy), we are obliged to sum over all final states at a given energy.

#### 2.1. Reduction of the many-body problem to an effective one-particle problem: the photoemission case

The photoemission *k*^{2} along the direction can be written in the dipole approximation as

where is the many-body final scattering state, normalized to one state per energy interval unit, for the *N*-electron system with one electron of momentum travelling to infinity and is its ground state; these states have respective energies *E*_{k}^{N} and *E*_{g}^{N}; is the incoming photon energy and is its polarization. Energy conservation imposes that . According to Breit & Bethe (1954), in order to satisfy the correct boundary conditions for the ejected photoelectron (no electron in a continuum state in the remote past), we must take the time-reversed scattering state by application of the time-reversal operator . Throughout this paper we shall use of length and Rydberg units of energy. The quantity is the fine-structure constant, equal to .

The photoabsorption *E*_{k}^{N} = *k*^{2} + *E*_{n}^{N - 1}, where *E*_{n}^{N - 1} is the energy of the remaining *N*-1 electrons of the system. By definition, the elastic channel is the channel for which *E*_{n}^{N - 1} = *E*_{g}^{N - 1}, the ground state of the (*N*−1)-electron system. The states with energy *E*_{n}^{N - 1} can have one or more electrons in the continuum, as in the case of double (shake-off channels), or one excited electron in a valence that remains in the system (shake-up channels).

Since we sum over all final states (*f*) we can eliminate the time-reversal operator and write the photoabsorption in the usual way as

In the case of photoemission from a deep core state of angular momentum *L*_{0} = (*l*_{0}, *m*_{0}), we assume that, to a good approximation,

where *A* is the usual anti-symmetrization operator [*A* = , with *A*^{2} = *A*], and are Slater determinants describing the configurations present in the ground state of the system. Normalization imposes if , and for simplicity we shall omit spin variables, since we are going to deal with non-magnetic systems.

In a similar way, we can write without loss of generality

where the functions , ignoring exchange effects, can be thought of as describing the excited photoelectron while the states are eigenstates of the Hamiltonian *H*^{N - 1} describing the remaining (*N*−1)-electron system with eigenvalues :

The tilde stands as a reminder that in the expansion (4) the relaxed states around the core hole are dominant. If needed, they can in turn be expanded in terms of Slater determinants that describe the intervening configurations in the final state. Borrowing the term from many-particle scattering theory, we can call the states final-state channels. Here and henceforth, the lower index *f* in the final state can be replaced by whenever we deal specifically with the scattering state .

The wave function is an eigenstate of the total Hamiltonian *H*^{N} with eigenvalue , *i.e.*

Moreover,

where is the interaction potential of the excited photoelectron with the rest of the system.

By inserting (4) into (6), projecting onto the states and using (5), we obtain for the amplitude functions the set of coupled equations

where

*I*_{c} being the for the core state and Δ*E*_{α} the left behind in the (*N* − 1)-particle system. The non-local interchannel potentials are the matrix elements between states and of the interaction potential and include local terms coming from the Coulomb interaction as well as non-local exchange terms originating from the exchange interaction. The set of equations (8) is supplemented with the boundary conditions that are related to the behaviour of the photoelectron at infinity and to the state of the (*N* − 1)-electron system according to the partition of the total energy between the photoelectron and the system. Each different partition has a different set of boundary conditions, which lead to a different solution of (8). For example, if we are interested in a particular photoemission channel with leaving behind the energy in the system, in the limit we should impose the scattering conditions

where we have made explicit the dependence of on as an argument rather than an upper index. Hence, as usual, δ_{αβ} is the Kronecker symbol, and *f*_{α} is the scattering amplitude. The factor normalizes the photoelectronic plane wave at the detector to one state per Rydberg. In the case of double electron transitions with two electrons in the continuum, (10) can easily be changed accordingly. For more details we refer the reader to Natoli *et al.* (1990), where a formal solution of (8) is given in the framework of the multichannel multiple-scattering theory. For the sake of the argument developed here, we need only the expression for the photoemission which is obtained from (1), (3) and (4) as

which is valid if we orthogonalize the excited channels to all the one-particle states belonging to the configurations that are present in the ground state . Here we have introduced the overlap integrals of the `passive' electron and indicated by the time reversal of (in practice, the complex conjugate, if spin is neglected). Spin degeneracy has been taken into account by including an extra factor of 2.

The set of equations in (8) contains the complete description of all the outcomes of the photoemission process, be they of intrinsic (*i.e.* consequent to the relaxation of the system around the core hole) or extrinsic (excitations created by the photoelectron on its way out of the system) origin. The complete solution of (8) is out of the question; however, one can analyse the consequences of some solutions in particular cases. Since we are mainly interested in structural analysis, both in photoemission and in photoabsorption, we need only consider the completely relaxed or elastic channel (*i.e.* that for which ), since this is the only channel capable of carrying structural information. For simplicity we shall attribute to this channel here and in the following the index . This channel carries most of the weight since, as a typical value, –0.9 (Mustre de Leon *et al.*, 1991).

We can then think of solving the set of coupled Schrödinger equations (8) by eliminating all unwanted channels in favour of the elastic one. The result is a single equation for the channel function with an effective complex energy-dependent non-local optical potential of the kind

where we have isolated the local Coulomb part of the potential (*V*_{c}) and indicated the energy dependence coming from the eliminated channels by the argument in . Once this equation is solved, we can write each in terms of through a relation of the type , which involves complicated inversions of the operators in (8). Again, a formal scheme of solution is given by Natoli *et al.* (1990). We can therefore write (11) as

so that everything is expressed in terms of . Notice that in the summation over α the most important term is , since by definition .

Our task is then to solve (12) with the asymptotic boundary conditions given by (10). This can be achieved in the framework of MST by transforming the integro-differential equation (12) into a Lippman–Schwinger equation with non-local potential, following the method illustrated by Natoli *et al.* (1986, 1990). In these papers the derivation of the MS equations is given for local potentials, but it is clear that the same derivation is valid for non-local potentials as well. The essence of the method rests on the partition of the space into Voronoi polyhedra (equivalent to Wigner–Seitz cells for periodic systems) such that the diameters of the polyhedra are always smaller than the nearest distance between their centres. In each such polyhedra, a local solution of (12) is obtained, which behaves like near the origin, where *j*_{l} (*k**r*) is the usual spherical Bessel function of order *l* and are the spherical harmonics of type (we shall use a real basis throughout). These functions are used to expand locally in the *i*th cell the overall solution of class *C*^{1} (continuous together with first derivatives) in the whole space, satisfying the boundary condition (10). Thus we can write locally,

provided the amplitudes satisfy the compatibility equations

where is the exciting amplitude originating from the plane wave in (10). A derivation is provided in Appendix *A*. Here, and are surface integrals over the surface *S*_{i} of the *i*th cell,

and are the KKR structure factors, which have the well known form

where , *h*_{l}^{+} is the Hankel function with outgoing spherical-wave behaviour, is the vector connecting the origins of the two cells centred at and , and connects site *i* with the origin of the coordinates *o*, which is assumed to coincide with the photoabsorber. The quantities are known as Gaunt coefficients.

In order to be able to use the physical language of MST, we introduce the quantities , and we define the suitably normalized local basis solutions , so that . The coefficients are easily seen to satisfy the MS equations

where we have introduced the quantity . In the case of cells of spherical shape it can be shown (Natoli *et al.*, 1990) that is the of the non-spherical potential located inside the *i*th cell and is a total scattering amplitude impinging on that cell in response to a plane-wave excitation of momentum . Indeed, (19) are the self-consistent equations for these amplitudes. We shall assume in the following that we are dealing with spherical cells, in order to make the discussion more physical. In the general case, (15) should be used.

For convenience, we define the cell *T*-matrix and the matrix of the free-spherical-wave propagator. Then (19) are easily solved for the amplitudes to give

where we have introduced the inverse *τ* of the MS matrix; , which is known as the full-scattering path operator and gives the total amplitude of propagation from site *i* to site *j*, starting with angular momentum *L* and arriving with angular momentum .

We now have all the ingredients that are needed to calculate the photoemission

Introducing the energy-dependent matrix element for the creation of the photoelectron,we easily find that our quantity of interest takes the simple form

Remembering (20), we then see that the photoemission current along direction is the modulus squared of the sum of all of the possible composite amplitudes that are obtained as products of an amplitude *M*_{L0 L}^{c} for exciting a core electron at site *o* (the origin), times the amplitude of propagation from this site to any other site (cell) *j* with final angular momentum , times the amplitude for emission along , times the phase factor that takes into account the phase relation of the electronic wave between sites *o* and *j*. The photoemission current is therefore the result of a complicated interference process that requires the reliability of the optical potential.

Expression (22) has constituted the basis for the interpretation of photoelectron diffraction data at low photoelectron kinetic energies with some success (Gunnella *et al.*, 1998). Computer programs in the MT approximation are also available (Gunnella *et al.*, 2000). The emitted current depends on three variables: the photon energy and the two polar angles of the direction . Therefore, as anticipated at the beginning of this section, by fixing the energy we can test the performance of a particular model optical potential in reproducing polar and azimuthal diffraction scans.

#### 2.2. The photoabsorption cross section

Turning now to absorption, we need to integrate (11) over and sum over all the final channels that have the same total energy . To this purpose it is expedient to make explicit the energy-conservation condition and to label the final state with *f*. By defining the Green's function matrix with outgoing wave boundary conditions,

the photoabsorption

is easily seen to bewhich is the same as the expression that we would have obtained if we had started from (1), (3) and (4). In Appendix *C*-4 of Natoli *et al.* (1990) it is shown that this matrix satisfies the set of coupled equations,

where we have written *E* for . Following the same steps as in the photoemission case, the elimination of all channels in favour of the relaxed one () leads to the following expression

where the non-local operator is the same as before and *G*_{00} obeys an equation corresponding to (12) with the same optical potential:

The solution to this equation within MST, for in cell *i* and in cell *j*, can be written as (Faulkner & Stocks, 1980)

remembering that the functions do not carry the normalization to one state per Rydberg. The second term on the right-hand side is the singular part of Green's function, is the lesser (greater) of and is the solution of (25) inside cell *i* that is irregular at the origin and matches smoothly to at the boundary. Making the reasonable assumption that the range of the functions is of the order of the atomic dimensions (this assumption is obviously true for *A*_{0}) and using (24) and (26) in (23), we finally obtain for the photoabsorption cross section

This expression arises because of the localization of the core initial state. *M*_{L0 L}^{c} is given by equation (21) and we define an `atomic' absorption given by

where

Equation (27) is the expression that we wanted to arrive at. It is valid under quite general conditions, as is apparent from our derivation, and shows the natural partition into an `atomic' contribution and an MS contribution. By removing all cells of the cluster except that containing the photoabsorber, reduces to so that the MS contribution is zero. Even though in a multiatom system we cannot define precisely an atomic entity, we shall continue to use the term `atom' for the central cell, because it is after all a reasonable approximation to the isolated atom. The fact that the second site index is equal to the first means that the electron source and the detector coincide, as already anticipated. Since the escape direction of the photoelectron has been integrated out, the only variable left is the energy and we can only test energy-dependent diffraction patterns. Therefore, these patterns might be distorted (compared with photoelectron diffraction patterns) by the energy dependence of the matrix elements that are necessary to create and detect the photoelectron at the absorption site.

However, we can hope to recover the structural information that we are interested in if we can neglect electronic correlations in the final state and if only one configuration (Slater determinant) is predominant in (3). Both of these conditions seem to be true, to a good approximation, for the *K*-edge spectra of metal ions in materials of biological interest. Indeed, the final photoelectronic states reached in these instances possess *p* symmetry around the photoabsorber and therefore are sufficiently delocalized not to suffer correlation effects with the electrons of the system. Moreover, very often the spin-restricted (unrestricted) Hartree–Fock approximation for the initial state provides a reasonably good description of the ground state for closed (open) shell configurations. Only in special cases do we have examples of configuration mixing in the ground state, as in the case of the Cu^{2+} ion in La_{2}CuO_{4} (Wu *et al.*, 1996) or valence fluctuating compounds. In such instances, a better strategy would be to solve the multichannel equations for the configurations that enter in the ground state and to eliminate all the other channels by expressing them as a function of the configurations of interest.

The predominance of one particular configuration implies that the spectral shape of the excited channels is rather featureless and similar to the ground state, since most of them consist of particle–hole or plasmon-like excitations, which do not drastically change the ground-state potential. Therefore, to a reasonable approximation, we can put . We can also write , because is the preponderant term of the series [by definition ]. Thus, we obtain from (24)

This equation, together with (23), tells us that the effect of the eliminated channels results in a shape function that modulates the originating from the primary channel

The expression for is easily obtained by observing that, in the same approximation, (29) gives

whereas (21) becomes, dropping ,

so that insertion of (33) into (27) and (28) provides the desired result, simply by replacing the transition matrix elements *M* with ,

Even with these approximations, the first-principle calculation of is not an easy task. The most noticeable effect is provided by double electron transitions (that are included in our scheme), since they show up as a change of slope in the `atomic' background component of the *b*), which occurs at definite energies that depend on the photoabsorber. An analytical or numerical modelling of these effects would be highly desirable, although for our purposes it is sufficient to know that such `kinks' can be factorized into this shape function.

We are left with the problem of determining the nature of the optical potential in (25). This potential contains the effect of both the intrinsic channels (excitations induced in the system by the sudden creation of the core hole) and the extrinsic channels (excitations created by the photoelectron on its way out of the system). Model systems to describe both types of processes, and their interference patterns have been studied by some authors (Fujikawa, 1999; Campbell *et al.*, 2002, and references therein). However, a practical scheme for realistic calculations has not yet been devised. The only optical potentials currently in use are those based on the Dyson self-energy of the photoelectron propagation in the system. This type of self-energy clearly accounts only for the extrinsic losses. Depending on the systems, various reasonable approximations have been devised.

It is well known, for example, that for metals we can obtain very good agreement with the observed absorption spectra by using a one-particle approach with an *X*-α potential and convoluting the calculated spectrum with a Lorentzian broadening function, which has an energy-dependent width that is related to the mean-free path of the photoelectron in the system by the relation (Müller & Wilkins, 1982; Müller *et al.*, 1984)

where is the full width at half-maximum.

In the framework of the above multichannel approach, this finding can be rationalized by observing that in a metal the completely relaxed channel together with the plasmon excitation channels (whether intrinsic or extrinsic) almost completely exhaust the sum rule

which holds because of the completeness of the intermediate relaxed states . Indeed, the intensity of the double-electron excitation channels is of the order of 10^{−2}–10^{−3} times that of the main relaxed channel (Filipponi, 1995*b*). Therefore, an optical potential given by is able to give a satisfactory picture of the absorption process for metals. The only discrepancy with experiment is that the calculated absorption maxima fall short of the observed maxima because of the energy independence of the *X*-α exchange.

A better approximation is provided by the Hedin–Lundqvist (H–L) potential (Hedin & Lundqvist, 1969, 1971), owing to its energy-dependent exchange and its imaginary part that is able to reproduce rather accurately the observed in metals (Penn, 1987). The H–L potential is the self-energy (based on the GW approximation and calculated at the local density of the system under study) of an electron that is propagating in a homogeneous interacting electron gas. Although initially devised to describe exchange and correlation corrections to the Coulomb potential due to the valence charge, Lee & Beni (1977) have extended its validity in the atomic-core region as well.

By neglecting the effect of the intrinsic processes, we can approximate as the Dyson self-energy of the photoelectron in the final state. This approximation is consistent with the physical picture of the photoabsorption process, in which we add an electron to the ground state of the (*Z* + 1)-equivalent atom. From this point of view, describes the propagation amplitude of the excited photoelectron from point to point . is the probability amplitude that the added electron remains in the state in which it was added to the system. The imaginary part gives the total probability amplitude for scattering out of this initial state. In this scheme the Dyson self-energy acts as a complex optical potential that describes the reduction of the wave-function amplitude of the elastic channel due to the transitions to all the other channels. The localization of the initial core state has the consequence that the optical paths of the photoelectron in the final state begin and end at the photoabsorbing site. We expect that neglecting the effects of the intrinsic processes on is a reasonably good approximation, since their main effect is already incorporated in the shape function .

We can therefore interpret the H–L potential as an effective optical potential that controls the propagation and damping of the excited photoelectron everywhere in the system. In this approach, this potential can be viewed as a local-density approximation to the self-energy of the photoelectron in real systems. Nowadays it is the potential that is most widely used in the calculation of the absorption and photoelectron diffraction cross sections of many systems, ranging from metals and semiconductors to ionic and covalent systems with varied success. In Appendix *B* we derive in some detail analytical expressions for the H–L potential, both for the sake of completeness and in order to make contact with other kinds of potentials, like the *X*-α and the Dirac–Hara potentials.

Finally, Fujikawa *et al.* (2000) have improved on the H–L potential by restricting its validity to the valence charge, as originally devised. They relied on the GW approximation for the photoelectron self-energy, in the solid, where , *V* is the bare Coulomb interaction and is the dielectric response function of the system, and they split both the polarization propagator *P* = *P* ^{vo } + *P* ^{c} and the one-electron Green's function *G* = *G* ^{v} + *G* ^{c} into core and valence parts. Since the core polarization was assumed to be much smaller than the valence polarization, Fujikawa *et al.* (2000) obtain an expansion in powers of *P* ^{c} for . Here *G* ^{v} *W* ^{v} is the self-energy for the valence electrons, which, when calculated *via* the plasmon-pole approximation for the dielectric function, is equivalent to the H–L potential; *V* _{ex}^{c} = *G* ^{c} *V* is the bare Hartree–Fock exchange potential and *G* ^{v} *W* ^{v} *P* ^{c} *W* ^{v} is the screened polarization potential for the ion cores. Details of this latter potential are given in Appendix *C*. Preliminary calculations for photoelectrons with greater than about 100 eV show that this non-local potential gives the same scattering amplitude as the total H–L potential at large scattering angles and provides a better description of small-angle scattering (Fujikawa *et al.*, 2000). However, more work is needed to establish its performance at low photoelectron energies.

Once the optical potential has been specified, we can proceed to the calculation of the key ingredients in (34), namely the cell functions and in the absorbing sphere, which are needed to calculate the atomic absorption and the transition-matrix elements that create the photoelectron, and the scattering-path operator . The inversion of the MS matrix (*T*_{c}^{ - 1} + *G*) becomes time consuming at energies greater than about 50–150 eV, depending on the number of cells in the cluster. Fortunately, in most cases, above ∼50 eV we can invert the MS matrix by series expansion, *i.e.*

so that the total absorption is seen to be made up of an atomic smooth contribution plus an infinite series of oscillatory signals (Benfatto *et al.*, 1986; Tyson *et al.*, 1992). This observation has been the basis for the development of the packages mentioned in §1, based on the empirical extraction of a structural signal from the experimental data and the comparison with a theoretical signal.

A major ingredient in the extraction procedure is the definition of a background atomic

such that the structural signal is defined asBy identifying with (ignoring the usual proportionality factor , where *N* is the Avogadro number) we are led by (34) to the identification

In this formula, the shape function seems to have dropped out from the ratio. However, the function's disappearance is a consequence of the identification of with , which is only approximate because the empirical definition of does not take into account the exact form of . The nearer is to , the less dependent becomes on the shape function. This is the reason why in many instances of

analysis seems to be constant and very near to 1. Also note that the atomic absorption does not factorize from the structural signal . Indeed, only if the optical potential is real can we show thatso that the structural signal is proportional to the atomic absorption. However, because of the presence of inelastic processes, the potential is complex, and the different behaviour of and should be taken into account in a refined treatment of the experimental data.

Therefore, in order to extract structural and electronic information on the system under study, we should not rely on an empirical definition of , which becomes increasingly difficult to define in the XANES part of the *et al.*, 1998; Gunnella *et al.*, 1998). Regarding the matrix elements for the creation of the photoelectron and its energy dependence, past and present experiences indicate that the smoothing action of the complex part of the potential helps us in reproducing the correct behaviour. Remaining discrepancies can be easily absorbed into a redefinition of this complex part. Even the effect of the shape function on the under the assumptions described above, can be mimicked *via* an additional damping. These statements will be substantiated by the applications of the method to particular systems in a following paper (Benfatto *et al.*, 2003).

#### 2.3. The mean free path

Expression (19) tells us that the structural part of the photoabsorption is proportional to the imaginary part of the product of the amplitude *M* for emitting the photoelectron, times the full scattering amplitude of propagation from the photoemitter site and back, excluding the contribution of the photoemitter itself, times another amplitude *M* for detecting the photoelectron. Since the optical potential is complex, we expect that the coherent part of the propagation (the part that only carries the structural information) will be attenuated by a damping factor related to the of the photoelectron in the system, whereas incoherent terms will appear that describe the effects of scattering out of the coherent channel due to the presence of inelastic processes. These incoherent terms are described by the complex part of the potential.

In this section we shall show explicitly that the mathematical structure of the theory contains the intuitive picture of the photoelectron damping in the expected way, by making contact with existing expressions for the

. At the same time, we shall find out how the incoherent terms due to the presence of inelastic processes find their place in the theory and, as a by-product, how we can mimic their effects when we ignore the details of their manifestation.For the sake of illustration, we shall assume that the MS series for in (35) converges and that the optical potential is local and of the muffin-tin form. The general *n*th term of the series, dropping one factor *T*_{c}^{o} that factorizes into the amplitude for emitting and detecting the photoelectron, has the form

where is the atomic *i* and angular momentum *l* in terms of the corresponding phase shift ; are the spherical wave propagators that are introduced in Appendix *A*. Without loss of generality, for simplicity we shall consider in (37) the case *n* = 3 with three sites *o*, *i* and *j*, so that *j* = *k*, *L*_{n} = *L*_{2}.

Since the potential is complex, the generic atomic phase shift is also complex, so we can write . Therefore,

which reduces to the known expression with real in the limit . In an electron–atom scattering process with complex potential (Landau & Lifshitz, 1966), it is known that, within a factor *k*^{-2 }, Im *k**t* represents the total scattering (elastic plus inelastic), whereas gives the elastic (without energy loss for the impinging electron). Based on the definition (25), the relation

holds, so that is the inelastic *kt*, since the propagator *G* in (37) is proportional to *k*, the photoelectron momentum.

From (38), it is clear that the coherent signal is obtained by choosing, for all the *kt* factors appearing in (37), the term . All the other terms containing at least one factor describe inelastic processes with loss of coherence for the photoelectronic wave. Indeed, we have for sufficiently small , which is usually the case.

Concerning the coherent signal of interest, we see that it has the same form and value as the structural signal obtained for a real potential, except for a damping factor coming from each *kt* and each propagator *G*. To proceed further, we need to find an expression for this latter factor. Again, for illustrative purposes, we shall work in the so-called plane-wave approximation (PWA), but it can be shown that all the following considerations are valid even when spherical-wave corrections for the photoelectron propagating wave are taken into account. In the PWA we can write (Natoli *et al.*, 1986)

where *R*_{ij} is the distance between sites *i* and *j* and is the photoelectron momentum in the interstitial region, in which the potential *V*_{I} is constant and can be complex.

Calling *R*_{i} the sphere radius of the potential at site *i*, it is easily seen that in the case of the three sites *o*, *i* and *j* the damping is represented by an exponential factor with exponent

We now need an expression for the phase shift , which in the WKB approximation is given by (Hara, 1967)

where the potential *V*_{i} (*r*) inside the sphere *R*_{i} is assumed to be complex. Notice that in keeping with the plane-wave approximation for the propagator *G* we have dropped the centrifugal term in the WKB expression, so that in this case the phase shifts are independent of the angular momentum *l*. However, as already emphasized, the argument remains valid even when spherical-wave corrections are taken into account.

In the path going from site *o*, the photoabsorbing site, to site *i* and *j* and back to *o*, each atom is traversed for a length 2*R*, if *R* is the radius of the corresponding sphere. Writing *R*_{path} = *R*_{oi} + *R*_{ij} + *R*_{jo} for the total length of the path, we see that the damping factor is given by , where

since in the interstitial region and, for example, inside sphere *j*, . Equation (41) is exactly the same as the expression that we would have obtained had we studied the propagation of an electronic wave in the potential *V*(*r*) by solving the problem in the WKB approximation. The is accordingly given by

which is consistent with the fact that is an attenuation factor for an amplitude of propagation. A further simplification is achieved if we take into account that everywhere in the system , *V*_{1} (*r*) and *V*_{2} (*r*) being the real and imaginary parts of the potential. Then, expanding the square root in (41) as

we obtain

remembering that *k* = *E*^{1/2}. Therefore, from (42), we finally get an expression in atomic units,

where , since *V*_{2} (*r*) is actually a self-energy. Equation (43) is the same as that given by Penn (1987), who has, in conventional units,

taking into account that *E*_{k} = *k*^{2} (Ryd). The only difference is that is defined as a volume average of the imaginary part of the H–L self-energy (see Appendix *B*) instead of the line average found in this paper.

From the above discussion, we see that what in practice determines the photoelectron damping is a line average of the various MS paths. If the system is homogeneous enough, this average will not depend on the path, so in practice we can ignore the position dependence of the self-energy and replace it by a function of only the energy of the photoelectron: . Since the *E*, neglecting the cut at the we can approximately take into account the effect of the absorptive part of the potential on the , calculated with the real part of the potential, by performing the following convolution

which holds if we are away from the edge, so that we can extend the integral to , and if does not vary too rapidly with *E*. Explicit calculations *via* both methods substantiate this relation.

The effect of the core-hole lifetime can be easily incorporated into the theory by adding to , where is the full width at half-maximum of the core hole and is related to its lifetime by . We assume here an , so that

of the core hole, which is true in most of the cases that we are interested in. To estimate the effect of this added damping on the we simply have to add to in (43)where is the inelastic

and is the corresponding to the finite lifetime of the core hole.In the H–L approximation for the self-energy at energies lower than the plasmon energy, the *B*. Damping originating from channels other than plasmon excitations (creation of electron–hole pairs) cannot be calculated by the H–L approximation. However, these damping terms have been evaluated (Hedin & Lundqvist, 1969) and found to be negligible compared with the core-hole lifetime and the experimental resolution. The damping terms can be introduced back into the equations using the formula given by Quinn & Ferrell (1958), as implemented in the *FEFF* code (Rehr & Albers, 2000). The photoelectron, on reaching the plasmon energy, interacts with the electron gas by creating plasmon excitations and suffering a reduction in the coherent wave-function amplitude. This fact would entail a rather sharp decrease in the photoelectron which should show up as a localized feature in photoabsorption spectra. In reality, the setting in of the plasmon damping is not so sharp because of the quantum interference between intrinsic and extrinsic losses at the plasmon edge (Fujikawa, 1999; Campbell *et al.*, 2002, and references therein). Therefore, we introduce into our fitting analysis an empirical model that takes into account both the smooth opening of the plasmon-loss channel and the effect of the electron–hole excitations.

#### 2.4. The muffin-tin approximation

As a result of the reduction operated in §2.2, Green's function obeys an effective one-particle Schrödinger equation, (25), which is better known as the Dyson equation. In (25) we have made explicit the appearance of the usual Coulomb or Hartree potential , which is given by

Here and *Z*_{k} indicate the position and the charge of the *k*th atomic nucleus and is the charge density of the system under study, which, in an independent-particle scheme or in a local-density approach, is given by

In general, the solution of the Schrödinger–Dyson equation in three dimensions does not pose any particular problems in the framework of the MST, as illustrated in Appendix *A*. However, the generation of the local solutions and the calculation of the surface integrals and for all site and angular momentum indices are rather time consuming, especially if we have to generate many theoretical signals to compare with experimental data in a fitting procedure. A great simplification is achieved if we spherically average both the Coulomb potential and the self-energy inside the atomic spheres. Both quantities depend on the electron density of the system under study, so, clearly, a self-consistent charge density would be highly desirable; however, this is often complicated to obtain in an amount of time that is appropriate to a fitting procedure. Therefore, usually the electron density is approximated by a superposition of spherically symmetric self-consistent atomic charge densities that are generated by currently available atomic programs. In order to make this approximation, we need to know how to expand a spherically symmetric function referred to one centre *j* around another centre *o* and take only the *L* = 0 component. This expansion is easily obtained based on the procedure proposed by Löwdin (1956) and later utilized by Mattheiss (1964). If is the atomic radial charge density around centre *j*, normalized so that , then the component of the charge density at distance *r* from the centre *o* is given by

where *R* is the distance between the two centres. As a result, the total overlapped spherically averaged charge density around centre *o* is given by

where *R*_{j} is the distance of neighbour *j* from centre *o* and is the charge density of the atom at this centre. Similarly, because of the linearity of the Poisson equation, it is easy to see that the spherically averaged potential around centre *o*, *V*_{c}^{o} (*r*), generated by this overlapped charge density is given by

where *V*^{j} (*r*) is the atomic potential generated by the radial charge density of the atom at site *j*:

Now, the spirit of the muffin-tin (MT) approximation consists of partitioning the molecular space (*i.e.* the space occupied by the atomic cluster under consideration) into three regions. Region (I) is made up of atomic spheres around the physical atoms, and the radii of these spheres are determined as described below. Region (III) lies outside an outer sphere circumscribing the cluster, and the interstitial region (II) lies between the outer and the atomic spheres. Empty spheres can be added to (II) to minimize the and to better account for the local variation of the true potential. This latter is spherically averaged into the spheres and approximated in regions (II) and (III) to a constant, which is set equal to the volume average of the potential in region (II). The true long-range variation of the potential in region (III) (Coulomb tail, in the case of isolated molecules) or a different constant (if the cluster is embedded in a different medium) can be easily incorporated into the scheme without any problems, although these corrections are usually neglected because of the damping effects of the imaginary part of the optical potential when the radius of the outer sphere is greater than the photoelectron The interstitial constant is easily calculated by expanding the potential around the centre of the outer sphere (and centre of the cluster). Assuming that *o* is such a centre and denoting as the volume of the interstitial region, we can use (48) to obtain

where *R*_{o} is the radius of the outer sphere. Similarly for the constant interstitial charge we find

Equations (47)–(51) define the Coulomb potential and the charge density throughout the molecular system in its muffin-tin form. *Via* the density, the exchange-correlation part of the electron self-energy can be calculated at each photoelectron energy, either in its local (H–L) or in the more elaborate non-local (GW) form.

In order to determine the MT radii of the atomic spheres, we can follow one of two prescriptions: one given by Norman (1974) and the other by Wille *et al.* (1986). In the Norman scheme, a Norman radius *R*_{i}^{Nor} is determined for each site *i* such that

which is roughly proportional to the ionic radius of the atom in the periodic table.

Then, given two nearest-neighbour sites *i* and *j*, the touching-sphere MT radii are given by the following formula,

which implies that their ratio equals the ratio of the corresponding Norman radii. Empirically, this prescription works in the case of covalent bonds, so that both spheres should be full. In the case of several neighbours at slightly different distances, we can take average interatomic distances in place of *R*_{ij}. The use of averages might entail some overlap between the various MT spheres, but that overlap is not considered to be a problem. Actually, in practical calculations, we can allow an overlap of about 10–15% of the MT radii to empirically simulate the and make up for the approximation of the potential in the interstitial region. In many cases this procedure seems to improve the agreement between theory and experiment, even though it is arbitrary on theoretical grounds. In the case of empty spheres or H atoms, *ad hoc* prescriptions can easily be devised.

On the other hand, Wille *et al.* (1986) suggest choosing the MT radii in such a way that the potential discontinuities at the boundaries of the MT spheres are minimized. This prescription seems to work for ionic compounds. In any case, even though the potential discontinuities can be reduced to an acceptable number, discontinuities remain between the boundaries of the atomic spheres and the interstitial region. These latter produce unphysical scattering, especially at photoelectron energies near the threshold, and can deform the shape of the calculated spectrum. Fortunately, by using the freedom of slightly varying the interstitial Coulomb potential and charge density and the overlap factor between MT spheres, we can nearly always reach a good agreement with experiments. However, the entire situation is highly unsatisfactory and calls for the elimination of the MT approximation.

In the generation of the charge density , care must be taken in many cases to reproduce the relaxation of the atomic electrons following the creation of the core hole and the screening action of the valence electrons. This relaxation is a dynamical process that, for deep core holes, is simulated in a static way by considering for the photoabsorbing atom the SCF charge density obtained by promoting the core electron to the first non-occupied valence orbital. This simulation seems to work for many systems (metals, semiconductors *etc.*) for which electronic correlation in the final state can be neglected, including also the biological compounds that we are interested in.

Given the MT potential and exploiting the ensuing simplifications, it is then easy to generate the local regular and irregular solutions of the Schrödinger–Dyson equation and Green's function, as detailed in Appendices *A* and *C*. In this way, polarized and unpolarized spectra can be calculated to fit to the experimental data. In the above procedure it would be a step forward to use self-consistent charge densities (and therefore potentials) in the MT form. Indeed, starting from Green's function, a self-consistent loop could be initiated by first generating a new charge density based on the formula

where *L* is any path in the complex energy plane containing the cut on the negative real axis due to the valence states and the poles due to the inner core states. By deforming the path so that it comprises two vertical lines at (the Fermi level) and (the energy value below any core state) and two horizontal lines at in the upper and lower planes with an appropriate , the contour integral could easily be calculated. The is easily determined by the condition

where *V* is the molecular volume (outside which the charge is negligible) and *N*_{el} is the total number of electrons. The new charge so obtained could serve to generate a new Coulomb and exchange-correlation potential (always of the MT form) that would generate a second Green's function, so as to repeat the loop until self-consistency was reached. This procedure would substantially improve the reliability of the fitting procedure at the cost of increasing the fitting time, although we could choose to repeat the loop only when really necessary. The self-consistent procedure has been implemented in the package *FEFF*8 [Ankudinov *et al.* (1998), to which we refer the reader for details] and seems to generate good final charge densities.

### 3. Future improvements and conclusions

The current status of techniques for calculating the inner-shell photoabsorption ), whereas XANES analysis has been treated in the reviews by Kizler (1993), Binsted & Hasnain (1996), Ankudinov (1999) and, in the special case of narrow-band highly correlated materials, de Groot (1996). In particular, Binsted & Hasnain (1996) have advocated the fitting of the entire spectrum rather than its components, and XANES, and have reached the conclusion that this is a viable procedure. To our knowledge, they were the first to indicate such a procedure. These authors analyse four representative model compounds and by a comparison of matrix-inversion and finite-path-sum methods they reach the conclusion that the latter method is more promising for fitting the edge region. Our preference is for the matrix-inversion method, because the MS series might not converge near the edge, in which case the maximum order of scattering becomes a further parameter without physical significance, which could lead to false minima in the *R* factor. Rather, we should attempt to reduce the impact of the many theoretical approximations present nowadays in the calculation of spectra. We summarize in the following the motivations behind this choice, having in mind not only biological applications but also the many X-ray spectroscopies that are used to obtain information on the many aspects of condensed matter physics.

Indeed, experience has shown that the sophistication of present-day spectroscopies has reached a point where a simple interpretation of the data is not sufficient to extract all the information that they conceal. Structural *et al.*, 1995), and in the domain of surface physics photoelectron diffraction is another example (Chen *et al.*, 1998). In such cases the possibility of calculating a good theoretical signal has been of paramount importance in establishing reliable tools for structural analysis, and the success of these techniques is due to the fact that the corresponding theories are asymptotic in the energy parameter. In fact, at high photoelectron energies simple approaches, like the description of the photoemission process in terms of one electron moving in an optical potential of the local density type and the muffin-tin approximation, are usually sufficient for a realistic description of systems of general type.

Unfortunately, while structural information dominates the high-energy spectral region, both structural and especially electronic information are present in the low-energy region. It would be extremely useful to have a reliable theory to help us extract the relevant information from the data, especially since the most interesting and promising spectroscopies for unraveling electronic structure (like X-ray magnetic and natural dichroism, resonant elastic and . Indeed, electron-correlation effects in the final state and the coupling of the excited photoelectron with the core hole concur to complicate the description of the near-edge region in many spectroscopies. Last but not least, the effects of atomic vibrations should be included in the theoretical simulations if we extend the fitting procedure beyond the near edge, since the damping effect of these vibrations interferes with the electronic damping due to the imaginary part of the optical potential.

both magnetic and non magnetic) give most of their electronic information near the excitation edge. In this energy region, the muffin-tin approximation is no longer adequate for describing the geometrical details of the potential, since the excited electron is quite sensitive to them. Moreover, complicated many-body processes intervene to screen the core hole and to interfere with the direct excitation channel, as described in §1Despite many efforts in the past to tackle these problems, they have not yet received a satisfactory solution. This is the reason why in many spectroscopies the near-edge region is not adequately exploited to extract the rich information, both structural and electronic, present in it. This latter information would be extremely valuable in the case of biological systems, concerning, for example, the spin and charge state of metal centres. Based on the above considerations, efforts should be concentrated in the following areas.

#### 3.1. Inclusion of vibrational effects in spectra

In the structural analysis of *t*-matrix elements that take into account in an approximate way the effect of atomic vibrations [see de Andres & King (2001) for details and references]. It would be very useful to incorporate such developments into the photoabsorption programs and continue the research in the field. Also, the availability of an experimentally determined vibrational spectrum would be highly desirable, since this would reduce the number of quantities to fit. Note that the experimental data can support only a limited number of fitted parameters (Stern, 1993; Filipponi, 1995*a*; Michalowicz & Vlaic, 1998). In this respect several attempts have been made to avoid the need to consider the Debye–Waller (DW) factors as parameters in the fitting procedure. These DW factors can be either calculated using semi-empirical and *ab* *initio* methods that consider a small cluster around the absorbing atom (Dimakis & Bunker, 1998; Poiarkova & Rehr, 1999) or derived from other types of experiments as in the work of Loeffen & Pettifer (1996). These authors performed an calculation for the zinc tetraimidazole molecular cluster, which included thermal four-body correlations established from inelastic neutron scattering. In this case full vibrational information has permitted a quantitative comparison between theory and experiment, without recourse to fitting. They found that differences between theory and experiment persist, setting the limits on the systematic errors still present in the theory. Presumably the MT approximation (not good for short atomic bonds) and the over-damping of the imaginary part of the H–L potential (not suitable for such low-*Z* covalent systems) are the causes of the remaining discrepancy.

#### 3.2. Elimination of the muffin-tin (MT) approximation in the solution of the Schrödinger–Dyson equation

This stage is preliminary to all spectroscopies, since a reliable solution of the one-electron Dyson equation moving in an effective optical potential is necessary for obtaining the related response functions. A method based on the finite-difference method (FDM) (and therefore free of the MT approximation) and a program for solving the associated equation are already available. This method has been developed by Joly (2001), and a collaboration with the Frascati group in the past few years has led to interesting applications, ranging from the electron population analysis of near-edge absorption spectra (Joly *et al.*, 1999) to the interpretation of resonant in systems showing orbital ordering (Benfatto *et al.*, 1999). However, such applications have also shown the method's practical limitations, which are due to the large memory and the long CPU time required by the computer code, especially for low-symmetry systems. In this scheme, fitting procedures to experimental data become prohibitive and many applications to interesting systems in biology, earth sciences, chemistry and physics are out of reach. An alternative method has been developed by Foulis *et al.* (1990), based on the non-MT (full potential) multiple-scattered-wave theory of Natoli *et al.* (1986). This method it is not completely free from approximations in treating the interstitial region. Nevertheless, this method, together with self-consistent-field electron densities, has demonstrated that quantitative accuracy for molecules and clusters of modest size can be achieved (Foulis *et al.*, 1995). Probably, a combination of the two methods, as illustrated in Appendix *A*, is the way to eliminate the drawbacks of both. Moreover, in view of applications to magnetic dichroism, spin–orbit interaction and spin-polarized potentials can be included and indeed have been incorporated by Joly into his code (to be published).

#### 3.3. Implementation of the multichannel MST (MMST) in the photoemission process

As already discussed, in the near-edge spectral region (either in absorption or in photoemission) intrinsic inelastic phenomena (shake-up, shake-off) very often modify the transition amplitude due to the *et al.* (1990). This theory is, in principle, `*ab initio*' and is able to incorporate into its description the local atomic multiplet aspect of the problem and the extended character of the excited photoelectron wave function. However, the theory has never been implemented in a computer code. If successful, this implementation should lead to a breakthrough in the analysis of near-edge spectroscopies concerning electronic properties.

#### 3.4. Core hole dynamic screening and the problem of the optical potential in the final state

As we have shown, the calculation of the total absorption *ad hoc* rules to construct such a potential. In particular, the most widely used rule is based on the work of Von Barth & Grossman (1982), which prescribes the use of a relaxed fully screened charge to calculate both the Coulomb and the exchange-correlation potential of the H–L type (the so-called *Z* + 1 screened approximation). However, in order to convince ourselves that the situation is more complicated, it suffices to consider some absorption spectra. The *Z* + 1 approximation may be necessary for the description of the *K*-edge absorption in diatomic molecules (N_{2}, O_{2}) in order to bind the resonance in the final state, since this state lies in the continuum in the ground state of the molecule. However, this approximation is not necessary for the *L*_{II}, *L*_{III} circular absorption edges in the ferromagnetic phase of nickel. In this case the dichroic effect would disappear if we were to use a self-consistent spin-polarized relaxed potential, as has been checked by direct numerical calculations. Clearly the optical final-state potential is the result of a dynamic screening process that depends on the system and on the type of final state reached. In particular, in the case of nickel, the self-screening of the core hole by the excited photoelectron and the Pauli principle hinder the electronic relaxation and screening by the other electrons of the system, so that the optical potential closely resembles the initial-state unrelaxed and unscreened potential. In general, in the case of electronic correlations in the final state, the true amount of screening charge should be determined based on a as discussed in §3.3 above. It is very important to have a theoretical model for the screening process since it affects the relative position of the pre-edge and rising-edge peaks in the absorption spectra.

In summary, even though decisive advances in the understanding of inner-shell X-ray absorption have been made over the past few decades, much remains to be done in order to obtain a versatile tool for extracting electronic and structural information from experimental spectra. However, we have also given theoretical arguments, supported by the applications made up to now (Benfatto *et al.*, 2003), that even with the present approximations structural analysis by full-spectrum fitting is possible. Electronic analysis, like spin and charge state, electron population analysis *etc.*, should probably await further improvements.

### APPENDIX A

### Derivation of the MS equations

We derive here the solution of the Dyson-like equation

subject to the asymptotic scattering boundary conditions

in the framework of multiple-scattering theory and in the non muffin-tin case.

The main ingredients for this derivation are the single and two-centre expansions of the free Green's function with the following outgoing wave boundary conditions:

(*a*) Around one centre located at the origin,

with and the same definitions as in §2.1.

(*b*) Around two centres located at and . Defining and , we have, provided ,

where

are the real-space KKR structure factors already introduced by (18) in §2.1.

(*c*) If the two centres are such that , where *o* is the origin of coordinates, we find

with

We assume here that in order to treat scattering states; however, the above relations still hold for , if we take the analytic continuation of the functions of energy appearing in the above equations and if we set the incoming plane wave to zero in (53). In this way, within the same formalism we can also treat bound states.

For simplicity we consider only short-range potentials (the extension to long-range potentials is, however, straightforward, as will be apparent from the following). There exists a volume that encloses all the atoms of our system, and this volume is sufficiently large that in the solution of (52) is of the form of (53). We partition this `molecular' volume in *N* non-overlapping cells with centres at (these cells might even be empty, *i.e.* not enclosing a physical atom). The cells completely fill (*i.e.* without interstices, so that ) in such a way that:

(*a*) There is a finite, but small, neighbourhood around the origin of each cell that lies completely within the cell. If there is an atom inside the cell, its nucleus should coincide with the origin.

(*b*) The shortest inter-cell vector, joining the origins of nearest-neighbour cells, is larger than any intra-cell vector.

We then start from the obvious equality involving surface integrals

in which the meanings of the symbols are obvious. This equality is valid for all , provided that is continuous with its first derivatives, as required by the solution of the Schrödinger-like equation.

By taking inside such that , we use (54) in the surface integral over *S*_{i}, (55) in the integral over *S*_{j} and (57) in the integral over *S*_{o}. We then find

so that, because of the angular completeness of the set , the expression inside the brackets {} should be zero for each *L*.

Now, inside each cell , we introduce basis functions , which are solutions of the Dyson-like equation (52) and behave at the origin like . These basis functions constitute a complete set so that the general solution can be locally expanded as . Likewise, in the outer domain , in order to impose the boundary conditions (53), we can take

where , owing to the well known decomposition of a plane wave. Inserting these expressions into (60) using the relation

and the identity

which is obtained from (58) by observing that , we derive (15) of §2.1,

where the definitions (16) and (17) give the surface integrals and .

In the case of long-range potentials (or if the short-range potential is substantially different from zero in the outer region) we have to supplement (62) with a similar equation obtained from (59) by taking to describe the effect of the scattering of the external potential. However, the solution inside the central cell is very often rather insensitive to the potential in the external region, provided the radius of the cluster is greater than the of the photoelectron at the energy considered. Clearly the above formalism can be easily extended to treat bound states.

The solution of (62) proceeds *via* the calculation of the surface integrals and and the structure factors . These latter are routinely calculated and, in principle, there should be no problem in calculating the surface integrals and . However, their practical calculation for polyhedral cells of general shape, though feasible, might be rather complicated. For this reason, it is preferable to use cells of spherical shape, inscribed in the polyhedra, especially since the quantity , which was introduced in §3.1, can then be interpreted as the scattering amplitude of the cell. The filling of space in this case leaves the problem of the interstitial region (IR), which can be minimized by increasing the number of empty spheres. However, in order to treat the irregular shape of the remaining IR we are obliged to approximate the potential by a constant (*e.g.* its volume average). By absorbing this constant into the definition of energy the interstitial potential becomes zero, so that the surface integrals and and (62) remain unchanged by application of Green's theorem. A further simplification, though not a necessary one, is achieved by replacing the potential inside each sphere by its angular average, so that we have to deal with spherically symmetric potentials. In this case the basis functions can be written as , where *R*_{l}^{i} (*r*_{i}) is the solution of the radial Schrödinger equation inside sphere *i*. Consequently the surface integrals and reduce to

and

where we have introduced the Wronskian of two functions *f*(*r*) and *g*(*r*), , calculated at the sphere radius *R*_{s}. The atomic scattering amplitude , also known as the atomic *t*-matrix, takes the form

which is the expression that would be found by solving the scattering problem for a spherical wave of angular momentum *l* impinging on a spherical potential truncated at radius *r* = *R*_{s}. From scattering theory we also have , where is the phase shift of the radial solution *R*_{l} (*r*) with respect to the free solution *j*_{l} (*r*) caused by the potential. If this latter is real, conservation of requires that (optical theorem), so that is real; otherwise , which implies that is complex and that the difference is related to the loss of due to the absorptive part of the potential, which becomes the source of the damping of the electronic wave (see §2.3).

### APPENDIX B

### The Hedin–Lundqvist potential

Much work has gone into approximating the exchange-correlation part of the optical potential in a way suitable for numerical applications. Hedin & Lundqvist (1969, 1971), by incorporating the Sham–Kohn (Sham & Kohn, 1966) density-functional formalism for excited states within the single-plasmon-pole approximation of the electron gas dielectric function, have produced a useful local-density approximation to :

Here is the self-energy of an electron in a homogenous interacting electron gas with momentum , energy and density , the local density of the actual physical system.

Since , neglecting the small exchange and correlation correction, we can write for the local exchange and correlation potential (following Lee & Beni, 1977)

The local momentum is defined as

where *k*^{2} is the photoelectron measured from the for an extended system or from the first for finite systems (Hara, 1967) (atoms and molecules) and is the local Usually, we can omit to perform the self-consistent procedure implicit in (66) for the determination of , except perhaps near the Fermi energy.

We follow Lee & Beni (1977) in extending the validity of (64) to the atomic core region, although the H–L approach was initially devised to describe exchange and correlation corrections to the Hartree potential due to the valence charge.

To calculate , we use equations (25.1), (25.14) and (25.15) of Hedin & Lundqvist (1969), which are based on the so-called GW approximation:

Equation (67) corresponds to the self-energy of the test electron interacting with the charge fluctuations of the medium. Here is the propagator for the test electron at energy measured from the ; is an infinitesimal positive quantity that is necessary for imposing the correct boundary conditions so that the integral converges when .

The

of the electron gas is given by , whereis the Fermi momentum, which corresponds to the constant density of the homogeneous electron gas [eventually to be taken equal to the local density of the inhomogeneous system]. The mean inter-particle distance *r*_{s} is given by

so that

Moreover, the function is the frequency-dependent dielectric function of the electron gas in the plasmon-pole approximation (Lee & Beni, 1977):

where

is the plasmon energy and

is the plasmon-pole dispersion relation. We have taken the coefficient 4/3 following Filipponi *et al.* (1988) rather than Lee & Beni (1977).

By performing the ω′ integral of (67), we obtain

and hence

where we should take the principal part of the integral in (68); is the Fermi distribution function. To proceed further we introduce the dimensionless quantities

so that

with

Calculating the above integrals at (*i.e.* *E*_{r} = *X*_{k}^{2}) we obtain for

where

and

Note that we have been considering excitations above the the usual static and energy-dependent exchange obtained by Dirac (1930) and discussed by Hara (1967):

so that . We recognize in the first term of (70)This comes from the constant part of the dielectric function in the first integral of (68). When multiplied by , where is a parameter, we obtain the usual form of the *X*- Slater static exchange potential (Slater, 1979; Schwarz, 1972),

Usually a parametrized form is also assumed for *V*_{ex}^{DH},

in such a way that *V*_{ex}^{DH} reduces to the *X*- form for *X*_{k} = 1.

As already noted by Chou *et al.* (1985),

for excitation energies much larger than the *V*_{ex}^{DF} becomes proportional to the density and decreases inversely with energy, instead of being proportional to *k*_{F} (or ) and independent of energy, as for the *X*- exchange *V*_{ex}^{S}. This behaviour is to be expected since, with increasing energy, the Pauli exclusion principle is less and less active, owing to the fact that at high energies the wave function of the test electron is already nearly orthogonal to that of the electrons in the system. Several authors have found that this energy dependence is crucial for a good calculation of phase shifts (Chou *et al.*, 1985, 1987, and references therein).

The next contribution in (70) still comes from the first integral in (68) but contains the effect of the frequency dependence of the dielectric function. Therefore, this contribution is a kind of screened exchange. The contribution is positive and roughly constant in the interval and then decays rapidly like

since

(This asymptotic relation actually holds already for ).

Finally, the last contribution in (70) comes from the last term in (68), which in turn arises from the poles of the inverse dielectric function. Therefore, this contribution is associated with a Coulomb correlation hole and stems from the fact that the test electron tends to keep away electrons of both spins because of the Coulomb repulsion. The contribution has a negative sign, like the exchange, but cannot be given an analytic form in terms of elementary functions. However, we can calculate the behaviour of the term at high energies,

so that the net contribution to *V*_{exc} is

This behaviour shows that the effect of the Coulomb correlation hole is more persistent at high energy than the D–H exchange. This effect causes the difference in the phase shifts that are calculated with the two potentials.

The numerical calculation of (70) does not pose particular problems, although it is a little time consuming. A good numerical approximation, valid in the entire energy range of the excited photoelectron, has been given by Mustre de Leon *et al.* (1991). Both options are incorporated into our computer code.

It is found that the imaginary part in (69) is amenable to quadrature. In fact, for and the integral becomes

where

is the maximum real root of the equation

Therefore, is non-zero only if , *i.e.* for

In other words there is no damping if the , page 93). The function is defined to be 1 for and zero outside. The interval [*X*_{1}, *X*_{2}] is any interval, lying inside [0,*X*_{M}], where or

For , this cubic polynomial has three real roots, of which only one, *X*_{1}, lies in the interval [0,*X*_{M}] and is such that in [*X*_{1}, *X*_{M}]. Therefore, the integral becomes

where

The root *X*_{1} is given by

where

We could also exploit the facts that *X*_{M} is a solution of the equation *w*(*X*) = (*X*_{k}^{2} - 1)/2 and *X*_{1} is solution of *w*(*X*) = *X*(*X*_{k} - *X*/2) to eliminate *w*(*X*) in (73).

Equation (73) is the same as the expressions that were obtained by Penn (1987), apart from a volume average and a correction for some misprints, and by Mustre de Leon *et al.* (1991). In the spirit of the local-density approximation (66), the expressions (70) and (71) for the real and imaginary parts of become functions of the position through the same dependence on the effective density of the system under study.

### APPENDIX C

### Non-local core polarization potential

In order to calculate the screened polarization potential *G*^{v} *W*^{v} *P*^{c} *W*^{v} we need the explicit expression for the full RPA (random phase approximation) polarization propagator (Hedin & Lindquist, 1969),

Here includes both space and spin variables, and are the one-electron excitation energies of the system. The sum over *k* runs over unoccupied electron states, while *l* runs over the occupied core and valence electron states. By splitting the summation over *l* into core and valence contributions, *P* can be written as a sum of core and valence parts,

where v_{o} means occupied valence states. Similarly, we can split the summation over *k* in the expression for the one-electron Green's function,

to obtain

In this case, the index c runs again over the core states while v runs over both the occupied and the unoccupied valence states.

The question to consider is the effect of screening in the core-polarization term *G* ^{v} *W* ^{v} *P* ^{c} *W* ^{v}. For a free atom, the screening is small and we have to a good approximation *G* ^{v} *V**P* ^{c} *V*, where *V* is the bare Coulomb potential. For long distances from the ion core, the core-polarization term further reduces to the well known local potential , where is the dipole polarizability and *e* is the (Hedin, 1965*a*,*b*). In a solid we expect, from simple physical considerations, that for long distances we should have a statically screened polarization potential, *G* ^{v} *W* ^{v} (0)*P* ^{c} *W* ^{v} (0). Since an *r*^{ - 4} potential is already very weak, the additional screening should make it negligible outside the Wigner–Seitz cell of the ion under consideration. Inside the Wigner–Seitz cell, on the other hand, we do not expect much screening to take place because of the cost in that is needed to localize the screening charge.

The symbol *G* ^{v} *V**P* ^{c} *V* stands for a convolution in energy space (Fujikawa & Hedin, 1989), which can be performed analytically, giving

where and . The more tightly bound the core level *l*, the smaller its contribution to , because the overlap with the unoccupied function *k* is smaller. Thus the outermost core level will give the dominant contributions.

We replace by a constant , the average , 1977*a*,*b*). We define a function ,

where the last equality follows by closure and and are the one-electron density matrices for all electrons and for core electrons, respectively. With we then have

With the use of closure we avoid the summation over the unoccupied states. Still, the density matrix contains a sum over the occupied extended states. We will here take a simplified approach and represent the sum over Bloch functions in each atomic cell *i* by a sum over localized functions . This approach is well motivated for rare-gas solids but a more serious approximation for, say, metals. In future, a more accurate treatment should be considered.

From (80) and the expression for Green's function (26), we see that the range of is apparently not limited to atomic dimensions. However, there are damping effects present in Green's function, and the energy-independent term defined by (79) decreases fairly rapidly with the distance [see (101) below]. Therefore, in (26) we can neglect, to a good approximation, terms coming from spheres other than the sphere under consideration (labelled *i*). Since, for simplicity, we shall spherically average the optical potential inside each sphere, we can then write (see Appendix *A* and §2.2)

with

where is the regular solution of the radial Schrödinger equation with the optical potential matching smoothly to at the sphere radius *R*_{s}; is the irregular solution matching smoothly to , , and *t*_{l}^{i} is the atomic *t*-matrix for the atom at site *i*.

Therefore, in order to obtain an expression for the atomic optical potential at this site, we need to calculate the energy-independent term . In each atomic region *i* we average the core charge density to obtain a spherically symmetric function

Though both and are in the same atomic region, and are not necessarily in the same atomic region. From the first term in square brackets in (79) we have an intraatomic contribution to :

where both and belong to atomic region *i* and is also in region *i*. is

where

and *R*_{i} is the radius of the atomic region *i* (muffin-tin radius).

In addition to this intraatomic contribution to the first term of (79) we have the interatomic term, where and are still in region *i*, whereas is in region *j* . This interatomic term is given by

where in terms of (Gaunt's integral)

We note that depends on *m*. To avoid this difficulty, we use the spherically averaged value of . Finally, we have the averaged expression of *F*_{L} in (89) in terms of the Clebsch–Gordan coefficient , which is denoted *F*_{l}^{0}:

The second term in square brackets in (79) is more difficult to calculate in general. We note that

where and we assume that the core orbital is localized on site *i*. We take both and to be in the same region *i*. For the occupied states, we will not use Bloch functions but instead take a simplified approach and use localized functions . This approach is well motivated for rare-gas solids but a more serious approximation for, say, metals.

By spherically averaging at each site *i*, we can obtain a simple representation for (94) as

where *m* and *n* stand for one-electron atomic states in atom *i* that have the same angular quantum number *L* = (*l*,*m*), *e.g.* 2*s*, 3*s* and so on. The quantity *d*_{m}^{i} ( = *d*_{mm}^{i}) is the spherically averaged electron density of the *m*th atomic function at site *i*, and *d*_{mn}^{i} is the cross charge from the *m*th and *n*th atomic functions on site *i* written in terms of the radial parts of the atomic wave functions, *d*_{mn}^{i} = *R*_{m}^{i} (*r*)*R*_{n}^{i} (*r*)^{*}. *R*_{m}^{i} refers to a localized function, which for metals has a fractional occupation number. When we use this simple approximation in (79), the *intraatomic* contribution to *A* can be written as

with

We see that *J*_{i} has only a spherically symmetric contribution.

The *interatomic* contribution *J*_{i}^{j} is not difficult to evaluate:

where, by taking *r*_{i} = *r*, *J*_{m}^{ij} (*r*) is given as

Here *n*_{m}^{j} is the number of electrons on site *j* and . Cross terms such as *d*_{mn}^{j} give no contribution to *J*_{m}^{ij} because of the orthogonality between the *m*th and *n*th shell functions. Therefore, can be written as

where *A* is a sum of one- and two-centre terms,

We find a good convergence for the two-centre sum, the second term on the right-hand side in (101), when we include the surrounding atoms up to the third shell for the systems considered here.

The optical potential can be given by (105) after the spherical averaging of the potential over and in the same atomic region:

where is expressed in terms of , and Clebsch–Gordan coefficients,

As previously demonstrated, it is enough to include only *A*_{0} and *A*_{1} in the expansion (106) (Fujikawa *et al.*, 1993). We thus obtain an explicit expression for ,

where is defined by

Each *g*_{l} includes the radial solution for the potential to be determined. Since depends on *g*_{l} and , in principle we have to solve coupled self-consistent equations.

### References

Andres, P. L. de & King, D. A. (2001). *Comput. Phys. Commun.* **138**, 281–301. Web of Science CrossRef Google Scholar

Ankudinov, A. L. (1999). *J. Synchrotron Rad.* **6**, 236–238. Web of Science CrossRef CAS IUCr Journals Google Scholar

Ankudinov, A. L., Ravel, B., Rehr, J. J. & Conradson, S. D. (1998). *Phys. Rev. B*, **58**, 7565–7576. Web of Science CrossRef CAS Google Scholar

Benfatto, M. & Della Longa, S. (2001). *J. Synchrotron Rad.* **8**, 1087–1094. Web of Science CrossRef CAS IUCr Journals Google Scholar

Benfatto, M., Della Longa, S. & Natoli, C. R. (2003). *J. Synchrotron Rad.* **10**, 51–57. CrossRef CAS IUCr Journals Google Scholar

Benfatto, M., Joly, Y. & Natoli, C. R. (1999). *Phys. Rev. Lett.* **83**, 636–639. Web of Science CrossRef CAS Google Scholar

Benfatto, M., Natoli, C. R., Bianconi, A., Garcia, J., Marcelli, A., Fanfoni, M. & Davoli, I. (1986). *Phys. Rev. B*, **34**, 5774–5781. CrossRef CAS Web of Science Google Scholar

Binsted, N. & Hasnain, S. S. (1996). *J. Synchrotron Rad.* **3**, 185–196. CrossRef CAS Web of Science IUCr Journals Google Scholar

Breit, G. & Bethe, H. (1954). *Phys. Rev.* **93**, 888–890. CrossRef Web of Science Google Scholar

Byron, F. W. & Joachain, J. C. (1974). *Phys. Rev. A*, **9**, 2559–2568. CrossRef CAS Web of Science Google Scholar

Byron, F. W. & Joachain, J. C. (1977*a*). *Phys. Rev. A*, **15**, 128–146. CrossRef CAS Web of Science Google Scholar

Byron, F. W. & Joachain, J. C. (1977*b*). *J. Phys. B*, **10**, 207–226. CrossRef CAS Web of Science Google Scholar

Campbell, L., Hedin, L., Rehr, J. J. & Bardyszewski, W. (2002). *Phys. Rev. B*, **65**, 064107-1–064107-13. Google Scholar

Chen, Y., Garcia de Abajo, F. J., Chassé, A., Ynzunza, R. X., Kaduwela, A. P., Van Hove, M. & Fadley, C. S. (1998). *Phys. Rev. B*, **58**, 13121–13131. Web of Science CrossRef CAS Google Scholar

Chou, S. H., Kutzler, F. W., Ellis, D. E., Shenoy, G. K., Morrison, T. I. & Montano, P. A. (1985). *Phys. Rev. B*, **31**, 1069–1076. CrossRef CAS Web of Science Google Scholar

Chou, S. H., Rehr, J. J., Stern, E. A. & Davidson, E. R. (1987). *Phys. Rev. B*, **35**, 2604–2614. CrossRef CAS Web of Science Google Scholar

Della Longa, S., Arcovito, A., Girasole, M., Hazeman, J. L. & Benfatto, M. (2001). *Phys. Rev. Lett.* **87**, 155501-1–155501-4. Google Scholar

Dirac, P. A. M. (1930). *Proc. Cambridge Philos. Soc.* **26**, 376–385. CrossRef CAS Google Scholar

Dimakis, N. & Bunker, G. (1998). *Phys. Rev. B*, **58**, 2467–2475. Web of Science CrossRef CAS Google Scholar

Faulkner, J. S. & Stocks, G. M. (1980). *Phys. Rev. B*, **21**, 3222–3241. CrossRef CAS Web of Science Google Scholar

Filipponi, A. (1995*a*). *J. Phys. Condens. Matter*, **7**, 9343–9356. CrossRef CAS Web of Science Google Scholar

Filipponi, A. (1995*b*). *Physica B*, **208**/**209**, 29–32. CrossRef Web of Science Google Scholar

Filipponi, A., Bernieri, E. & Mobilio, S. (1988). *Phys. Rev. B*, **38**, 3298–3304. CrossRef CAS Web of Science Google Scholar

Filipponi, A., DiCicco, A. & Natoli, C. R. (1995). *Phys. Rev. B*, **52**, 15122–15134. CrossRef CAS Web of Science Google Scholar

Foulis, D. L., Pettifer, R. F., Natoli, C. R. & Benfatto, M. (1990). *Phys. Rev. A*, **41**, 6922–6927. CrossRef CAS PubMed Web of Science Google Scholar

Foulis, D. L., Pettifer, R. F. & Sherwood, P. (1995). *Europhys. Lett.* **29**, 647–652. CrossRef CAS Web of Science Google Scholar

Fujikawa, T. (1999). *J. Phys. Soc. Jpn*, **68**, 2444–2456. Web of Science CrossRef CAS Google Scholar

Fujikawa, T., Hatada, K. & Hedin, L. (2000). *Phys. Rev. B*, **62**, 5387–5398. Web of Science CrossRef CAS Google Scholar

Fujikawa, T. & Hedin, L. (1989). *Phys. Rev. B*, **40**, 11507–11518. CrossRef Web of Science Google Scholar

Fujikawa, T., Saito, A. & Hedin, L. (1993). *Jpn. J. Appl. Phys.* **32**, 18–22. CrossRef CAS Web of Science Google Scholar

Groot, F. de (1996). *X-ray and Inner-Shell Processes: Seventeenth International Conference*, Hamburg, Germany, edited by R. L. Jonson, H. Schmidt-Böcking & B. F. Sonntag, p. 497. New York: AIP. Google Scholar

Gunnella, R., Bullock, E. L., Patthey, L., Natoli, C. R., Abukawa, T., Kono, S. & Johansson, L. S. (1998). *Phys. Rev. B*, **57**, 14739–14748. Web of Science CrossRef CAS Google Scholar

Gunnella, R., Solal, F., Sébilleau, D. & Natoli, C. R. (2000). *Comput. Phys. Commun.* **132**, 251–266. Web of Science CrossRef CAS Google Scholar

Hara, S. (1967). *J. Phys. Soc. Jpn*, **22**, 710–718. CrossRef CAS Web of Science Google Scholar

Hedin, L. (1965*b*). *Arkiv Fysik*, **30**, 231–258. CAS Google Scholar

Hedin, L. (1965*a*). *Phys. Rev.* **139**, A796–A823. CrossRef Web of Science Google Scholar

Hedin, L. & Lundqvist, B. I. (1971). *J. Phys. C*, **4**, 2064–2083. CrossRef Web of Science Google Scholar

Hedin, L. & Lundqvist, S. (1969). *Solid State Physics*, edited by F. Seitz, D. Turnbull & H. Ehrenreich, Vol. 23, pp. 1–181. New York: Academic Press. Google Scholar

Joly, Y. (2001). *Phys. Rev. B*, **63**, 125120-1–125120-10. Google Scholar

Joly, Y., Cabaret, D., Renevier, H. & Natoli, C. R. (1999). *Phys. Rev. Lett.* **82**, 2398–2401. Web of Science CrossRef CAS Google Scholar

Kizler, P. (1993). *Phys. Lett. A*, **172**, 66–76. CrossRef Web of Science Google Scholar

Landau, L. D. & Lifshitz, E. M. (1966). *Mécanique Quantique. Théorie non Relativiste*, p. 632 and *ff.* Moscow: Éditions MIR. Google Scholar

Lee, P. A. & Beni, G. (1977). *Phys. Rev. B*, **15**, 2862–2883. CrossRef CAS Web of Science Google Scholar

Loeffen, P. W. & Pettifer, R. F. (1996). *Phys. Rev. Lett.* **76**, 636–639. CrossRef PubMed CAS Web of Science Google Scholar

Löwdin, P. (1956). *Adv. Phys.* **5**, 3–172. Google Scholar

Mattheiss, L. F. (1964). *Phys. Rev.* **133**, A1399–A1403. CrossRef Google Scholar

Michalowicz, A. & Vlaic, G. (1998). *J. Synchrotron Rad.* **5**, 1317–1320. Web of Science CrossRef CAS IUCr Journals Google Scholar

Müller, J. E., Jepsen, O. & Wilkins, J. W. (1984). *Solid State Commun.* **42**, 4331–4343. Google Scholar

Müller, J. E. & Wilkins, J. W. (1982). *Phys. Rev. B*, **29**, 365–367. Google Scholar

Mustre de Leon, J., Rehr, J. J., Zabinsky, S. I. & Albers, R. C. (1991). *Phys. Rev. B*, **44**, 4146–4156. CrossRef CAS Google Scholar

Natoli, C. R., Benfatto M., Brouder, C., Ruiz Lopez, M. F. & Foulis, D. L. (1990). *Phys. Rev. B*, **42**, 1944–1968. Google Scholar

Natoli, C. R., Benfatto, M. & Doniach, S. (1986). *Phys. Rev. A*, **34**, 4682–4694. CrossRef PubMed Web of Science Google Scholar

Norman, J. G. (1974). *Mol. Phys.* **81**, 1191–1198. Google Scholar

Penn, D. R. (1987). *Phys. Rev. B*, **35**, 482–486. CrossRef CAS Web of Science Google Scholar

Poiarkova, A. V. & Rehr, J. J. (1999). *Phys. Rev. B*, **59**, 948–957. Web of Science CrossRef CAS Google Scholar

Quinn, J. J. & Ferrell, R. A. (1958). *Phys. Rev.* **112**, 812–827. CrossRef Web of Science Google Scholar

Rehr, J. J. & Albers, R. C. (2000). *Rev. Mod. Phys.* **72**, 621–654. Web of Science CrossRef CAS Google Scholar

Schwarz, K. (1972). *Phys. Rev. B*, **5**, 2466–2468. CrossRef Web of Science Google Scholar

Sham, L. J. & Kohn, W. (1966). *Phys. Rev.* **145**, 561–567. CrossRef CAS Web of Science Google Scholar

Slater, J. C. (1979). *The Self-Consistent Field for Molecules and Solids; Quantum Theory of Molecules and Solids.* New York: McGraw-Hill. Google Scholar

Stern, E. A. (1993). *Phys. Rev B*, **48**, 9825–9827. CrossRef CAS Web of Science Google Scholar

Tyson, T. A., Hodgson, K. O., Natoli, C. R. & Benfatto, M. (1992). *Phys. Rev. B*, **46**, 5997–6019. CrossRef CAS Web of Science Google Scholar

Von Barth, U. & Grossman, G. (1982). *Phys. Rev. B*, **25**, 5150–5179. CrossRef CAS Web of Science Google Scholar

Wille, L. T., Durham, P. J. & Sterne, P. A. (1986). *J. Phys.* (*Paris*), **47**, C8-43–45. CrossRef Google Scholar

Wu, Z., Benfatto, M. & Natoli, C. R. (1996). *Phys. Rev. B*, **54**, 13409–13412. CrossRef CAS Web of Science Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.