Chapter 6.11. GNXAS. I. Phase shifts and signal calculations

The GNXAS package has been under development since 1990 and its theoretical background has been thoroughly described in two main articles (Filipponi et al., 1995; Filipponi & Di Cicco, 1995). The suite of programs implements state-of-the-art approaches for computation of the XAFS signal in the framework of the multiple-scattering theory in a complex effective potential for the photoelectron within the muffin-tin approximation. A flow chart of the input/output information through the main and ancillary programs of the suite is shown in Fig. 1.

Figure 1 A flow chart of the main input/output information (dashed rectangles) through the PHAGEN and GNXAS codes is shown, together with the roles of the auxiliary CRYMOL and GNPEAK programs that are useful for preparing the required input files.

A specific feature of this code is that the signals are directly associated with a list of specific relevant pair, triplet, … atomic configurations involving the photoabsorbing atomic species in the model structure. The related multiple-scattering expansions or the cumulative irreducible n-body signal can be directly computed using a continued fraction algorithm (Filipponi, 1991). The output signals are tabulated in a suitable energy mesh with separate amplitude and phase information, together with the derivative with respect to all geometrical parameters. This approach allows efficient energy interpolation, and accounts for possible displacements from the reference structures and thermal or structural disorder in the subsequent comparison with experimental data, which is performed using the FITHEO program, without the need to recompute the signals. The reliability of the underlying approximations for the XAFS range has been validated by numerous successful applications of the method to solid, liquid and molecular structures. Applications to the XANES range are possible, especially for systems where the muffin-tin approximation is still a reasonable model.

The aim of the GNXAS package (Filipponi et al., 1995; Filipponi & Di Cicco, 1995, 2000) is to generate an accurate multiple-scattering simulation of the EXAFS signal χ(k) related to a model cluster of atoms relevant to the structure under investigation that can be used for structural refinement (fitting) of the raw experimental data. The experimental EXAFS signal χ(k) is defined as the relative oscillation of the absorption cross section σ(E) with respect to a smooth total atomic cross section normalized to the atomic cross section of the edge under consideration σ₀(E), i.e. χ(k) = , where k = [2m(E − E_e)/ℏ]^1/2 is the modulus of the photoelectron wavevector and E_e is the threshold energy of the edge.

The procedure can be divided into the following steps: (i) construction of the model potential for the final-state photoelectron, (ii) calculation of the atomic phase shifts, (iii) signal calculation and (iv) configurational averaging. The program PHAGEN takes care of steps (i) and (ii), while the model cluster can be generated using the CRYMOL program.

The absorption cross section for a cluster of atoms as a function of the incoming photon energy ℏω is given in terms of the full Green's function (GF) of the system bywhere α is the fine-structure constant, = is the wavefunction of the excited core electron with angular momentum L_c ≡ l_cm_c and spin σ, E = ℏω − E_e is the photoelectron energy and G(r, r′; E) is the solution of the equationHere, V_opt(r) is the (complex) potential obtained from the reduction of the photoabsorption many-body problem to an effective one-electron problem moving in an effective optical potential. In nonmagnetic materials V_opt(r) is spin-independent, so that in equation (1) we can drop the spin summation and multiply the cross section by 2.

In multiple-scattering theory (MST) one divides the space into non-overlapping domains Ω_j (cells) and introduces a partition of the potential V_opt(r) = , where v_j(r) coincides with V_opt(r) within the cell Ω_j and is zero outside. The Schrödinger equation (SE) is then solved in each cell and these local solutions are assembled together into the desired global solution.

One of the advantages of MST is that one can write an explicit form of the GF asin terms of the scattering-path operator and the local regular and irregular solutions Φ_Lⁱ(r_i; E) and Ψ_L(r_i; E) of the SE at energy E inside the cell located at i. Ψ_L(r_i; E) smoothly matches the Hankel function = at the boundary of the cell (k = E^1/2 in atomic units), whereas is regular at the origin. Here, are spherical harmonics, r_i = r − R_i is referred to the origin R_i of cell i and r_< and are the lesser and greater of r and r′, respectively.

Moreover, introducing the matrix T = , which is diagonal in the site indices, and the site off-diagonal matrix G = , we haveThe series expansion in equation (4) is absolutely convergent if the spectral radius (maximum eigenvalue of TG) ρ(TG) < 1, otherwise it is only conditionally (or asymptotically) convergent. Usually, for small clusters of (∼10–15) atoms there is absolute convergence. In all other cases it is not difficult to find the order n that gives meaningful results (Sébilleau & Natoli, 2009).

The quantity describes the scattering amplitude due to the potential v_i(r) of cell i for an impinging spherical wave of angular momentum L into a spherical wave of angular momentum L′, whereas gives the amplitude of free propagation from site i to site j for spherical waves. Its expression is given in terms of 3-j symbols bywhere R_ij = R_i − R_j is the vector joining site j to site i. For further details, refer to Natoli et al. (2024) and references therein.

With the aim of coping in a simple way with the wide variety of systems that are usually encountered in practical applications, one makes the ansatz that the optical potential is a function of the local density ρ(r) of the system under consideration. Experience has shown that a potential made up of the sum of the Hartree Coulomb potential and the complex Hedin–Lundqvist (HL) exchange and correlation energy (Natoli et al., 2003), satisfying this requirement, constitutes a very good approximation to this kind of `universal' optical potential with satisfactory results.

Moreover, partitioning the space into non-overlapping domains in an automatic way is not straightforward for arbitrary atomic locations. This difficulty motivated the introduction of the muffin-tin (MT) approximation, whereby one draws a sphere around each atomic position, inside which the potential v_i(r_i) is spherically averaged, whereas it is approximated by a suitably chosen constant in the interstitial region between the spheres.

In this approximation, the atomic T matrix is diagonal in the angular momentum indices and is independent of m, so that = = . In turn, = and Ψ_L(r_i; E) = , where R_l(r_i) and R_l⁺(r_i) are regular and irregular solutions of the radial SE with angular momentum l. Then,where the functions are calculated at the sphere radius . Here, j_l(kr_i) is the Bessel function of order l and f′(r_i) = df(r_i)/dr_i. Notice that if the potential is complex, then the phase shifts are also complex.

Within this approach, the total density ρ(r) of the cluster of atoms is constructed by superimposing the self-consistent charge densities of the neighbouring atoms on the atomic charge density present in each sphere, and then spherically averaging the anisotropic charge thus obtained. Moreover, since the Poisson equation is linear in charge density, we can overlap atomic Hartree potentials in the same way to obtain the total Coulomb potential of the cluster. This procedure makes it possible to calculate an average charge density and an average potential in the interstitial region. Finally, from the knowledge of the total density at each point in space, one can add the complex HL exchange and correlation energy to the Hartree potential to obtain the total optical potential. Once this latter is known in each sphere, one can integrate the radial SE by the Numerov procedure to obtain R_l(r_i) and and calculate via equation (6). As a result of the MT approximation, these latter are the same for all atoms that are equivalent under the symmetry operation of the cluster. Further details of this whole procedure are given in Natoli et al. (2003).

By inserting equation (3) into equation (1), we obtain the absorption cross section as the sum of an atomic absorption originating from the photoabsorber located at site 0 [the second term in equation (3)] and a contribution due to scattering from the environment, owing to the fact that σ_sct(ω) becomes zero in the absence of a neighbourhood of the photoabsorber. From the expansion in equation (4), taking into account that , we obtainwhereand the sum over is over the final angular momenta selected by the dipole selection rule.

If we average over the polarizations of the incoming photon, then equation (8) becomeswhere now is the radial integral. For s core states (1s, 2s, 3s, …) only the second term (l_c + 1) survives, which usually is much stronger than the (l_c − 1) term.

The separation of the total cross section into a scattering part σ_sct(ω) and a core absorption σ₀(ω) corresponds exactly to what is performed in the definition of the experimental structural signal χ(k), so that the theoretical signal to compare with χ(k) is σ_sct(ω)/σ₀(ω).

Writing down the MS series for s core states, we obtainwhere = σ₀/8παℏω and specifically Ξ⁰ⁱ⁰ = G_0,it_iG_i,0t₀ for , Ξ^0ij0 = G_0,jt_jG_j,it_iG_i,0t₀ for and Ξ^0ijk0 = G_0,kt_kG_k,jt_jG_j,it_iG_i,0t₀ for . In this notation, it is understood that the internal angular momentum indices have been saturated.

Notice that in the case of complex potentials the atomic absorption σ₀(ω) does not factor out from σ_sct(ω) as it does in the case of a real potential. PHAGEN takes this peculiarity of the complex potentials into account.

The analysis of MS paths based on equation (8), although useful, is not very efficient and versatile. One can, however, reorder the MS series in such a way as to obtain a faster convergence rate and to relate each order n to an n-body distribution function. Indeed, on the basis of equation (7) it is possible to calculate cross sections σ(0, i, j, …) and the related XAFS signals σ(0, i, j, …)/σ₀ − 1 associated with a photoabsorber 0 surrounded by an arbitrary number of atoms i, j, …, and this approach can be used to define and calculate the irreducible n-body contributions to the XAFS signal (Filipponi et al., 1995). The pair contribution due to atom i surrounding the photoabsorber is given bywhich is an oscillating function of k with a leading phase equal to 2R_0i. It can be calculated directly using the full-matrix inversions of equation (4) for two atoms, or expanded into the corresponding multiple-scattering series, as indicated explicitly in equation (13).

The triplet contribution associated with the simultaneous presence of atoms i, j around the photoabsorber is given bywhich is the XAFS signal due to atoms i, j minus the irreducible XAFS signals associated with the independent presence of atoms i and j. γ⁽³⁾(0, i, j) is also an oscillating function of k with a leading phase R_0i + R_ij + R_j0 determined by the shortest multiple-scattering path involving all sites. The multiple-scattering expansion for γ⁽³⁾(0, i, j) contains many terms and its direct calculation is computationally convenient if its amplitude and phase are well behaved.

The four-body contribution associated with the simultaneous presence of atoms i, j, k around the photoabsorber is given bywhich is an oscillating function of k with a leading phase determined by the shortest multiple-scattering path involving all sites. The multiple-scattering expansions for the γ⁽ⁿ⁾(0, i, j, …) signals become progressively cumbersome.

The total XAFS signal for a generic fixed atomic configuration surrounding the photoabsorber can therefore be expanded in a series of irreducible contributions of increasing order,that are expected to become progressively weaker and contribute with progressively higher frequencies to the spectrum. This series has improved convergence properties with respect to the multiple-scattering series since it combines in each term an infinite number of terms in the latter. The irreducible γ⁽ⁿ⁾(0, i, j, …) signals can be computed by combining a finite number of total cross section calculations for finite-size clusters according to equations (13), (14) and (15), and this task can be conveniently and efficiently accomplished using the continued fraction algorithm (Filipponi, 1991). An example of signal decomposition into the main irreducible two-, three- and four-body signals is illustrated for the face-centred cubic (f.c.c.) structure of solid nickel in Fig. 2.

Figure 2 XAFS signal decomposition for the f.c.c. structure of solid nickel. The signals associated with the first shell, and the first four triplet configurations, involving first-neighbour bonds at angles of 60°, 90°, 120° and 180°, are reported. The latter three triangular configurations are associated with two inequivalent signals according to the photoabsorber atom position (red in the insets) included in the corresponding γ(3) signal and the pair contribution associated with a longer distance in the triangle corresponding to the second, third and fourth shells, respectively. The intensity of the signal is greater than for the focusing effect. The last signal at the bottom is associated with the first four-atom tetrahedron configurations. No account for the thermal average is included.

In particular cases each γ⁽ⁿ⁾(0, i, j, …) signal can be approximated/computed using the leading terms of the corresponding multiple-scattering series. For γ⁽²⁾ such a term corresponds to the well known single-scattering contribution given by Schaich (1984):This is an exact curved wave expression which is at the basis of any modern single-scattering EXAFS calculation. Efficient equations for higher multiple-scattering paths were first introduced by Gurman et al. (1984) and other authors (Brouder et al., 1989).