Relativistic Corrections to Phase Shift Calculation in the GNXAS Package

: Modern XAFS (X-ray Absorption Fine Structure) data-analysis is based on accurate multiple-scattering (MS) calculations of the X-ray absorption cross-section. In this paper, we present the inclusion and test of relativistic corrections for the multiple-scattering calculations within the GnXAS suite of programs, which is relevant to the treatment of the XAFS signals when atoms with high atomic number are contained into the system. We present a suitable strategy for introducing relativistic corrections without altering the basic structure of the programs. In particular, this is realized by modifying only the Phagen program calculating the atomic absorption cross sections and scattering t -matrices for the selected cluster. The modiﬁcation incorporates a pseudo-Schrödinger Equation (SE) replacing the Dirac relativistic form. The phase-shift calculations have been put to a test in two known molecular and crystalline cases: molecular bromine Br 2 and crystalline Pb. Calculations in an extended energy range have been shown to be very close to the non-relativistic case for Br 2 (Br K-edge) while corrections have been found to exceed 25% for amplitude and phases of the XAFS multiple-scattering signals (Pb L 3 -edge). Beneﬁts in the structural reﬁnement using relativistic corrections are discussed for crystalline Pb at room temperature.


Introduction
The GnXAS package [1,2] is an advanced software for XAFS (X-ray Absorption Fine Structure) data analysis. The suite of programs implements state-of-the-art methods for the computation of the XAFS signal in the framework of the multiple-scattering theory (MST), using a complex effective optical potential for the photoelectron moving in a cluster of atoms modeling the system under study. The XAFS signal so obtained constitutes the input for a rigorous fitting procedure to the raw experimental data in order to derive structural and dynamical information on the system [3,4]. The underlying theory, basic methodology and practical applications have been widely discussed elsewhere. The GnXAS package is currently available for scientific research and numerous successful applications of the GnXAS package have been published so far by different research groups.
In this approach for XAFS data-analysis, the effective optical potential is assumed to be a functional of only the electron density of the selected cluster of atoms, which is obtained by superposition of self-consistent relativistic atomic charge densities in the muffin-tin (MT) approximation [5]. For these atomic potentials, one solves the local non-relativistic Schrödinger Equation (SE) to calculate the atomic t-matrices (t l = e i δ l sin δ l , where δ l are the phase-shift functions) describing the scattering process for each atom and relevant l.
These quantities, together with the spherical wave free propagators for the photo-electron, constitute the ingredients to calculate the wanted XAFS signal.
One of the present limitations of the package is the use of the non-relativistic SE (and consequently t-matrices) which may prevent accurate XAFS data-analysis when atoms with 'high' atomic number Z are involved in the system under study. By 'high' we mean atoms for which relativistic corrections are important (non negligible) for the definition of physical properties, such as, for example, scattering power and transition probabilities. The introduction of relativistic effects in XAFS multiple-scattering calculations has been described in several papers (see Refs. [6,7] for example) but few examples of applications and detailed studies of the effect of the various approximations were given so far.
The purpose of this paper is to present a suitable strategy for eliminating this limitation without altering too much the present non-relativistic structure of the GnXAS suite of programs, putting to a test the multiple-scattering (MS) simulations in selected molecular and crystalline systems. In particular, the introduction of relativistic corrections is realized by modifying the program (Phagen, part of the GnXAS suite) that generates atomic absorption cross sections and scattering t-matrices for the selected cluster. This goal is made possible by deriving a pseudo-SE from the Dirac Equations that takes into account all of the relevant relativistic corrections. Model XAFS relativistic and non-relativistic calculations are performed for the Br K-edge in molecular bromine (Br 2 ) and for the Pb L 3 -edge of crystalline lead (at room temperature). The latter is a system for which relativistic corrections are expected to be higher. The importance of the relativistic corrections and the impact on structural refinement for those systems can be thus evaluated by direct comparison with experimental data.
The paper is organized as follows: (i) derivation of the equations used by the modified Phagen program; (ii) presentation of phase-shift calculations using non-relativistic and relativistic calculations for the selected systems; (iii) comparison and discussion of results.

Derivation of the Pseudo-Schrödinger Equation
In the present form, Phagen numerically solves the non-relativistic SE by the Numerov method for the continuum states in the effective optical potential V e f f (r, ), a sum of the electrostatic (Coulomb) potential V c (r) and the energy-dependent exchange-correlation potential of the Hedin-Lundqvist type [8]. For short, we shall write in the following V e f f (r, ) as V(r).
The Dirac Equation (DE) is the starting point for a fully relativistic calculation of the atomic t-matrices. In atomic units for lengths and Rydberg units for energies, indicating by α = e 2 /(hc) = 1/137.036 the fine structure constant, the DE for an electron with total energy (including rest mass) E in presence of the central potential V(r), specified above, writes Here, Ψ is the Dirac bi-spinor where g κ (r) and f κ (r) are the upper and lower radial function components, and are the velocity and mass term operators, σ denoting the Pauli matrices. In the bi-spinor Equation (2), χ µ κ is the spin-orbital labeled only by the orbital projection µ and the relativistic quantum number κ. l µ − ν|J µ is the usual Clebsh-Gordan coefficient. Since the spin-orbital is an eigenstate of J 2 , L 2 and S 2 , one can easily derive that Therefore, the relativistic quantum number κ determines both J = |κ| − 1/2 and l according to κ = l for J = l − 1/2 and κ = −(l + 1) for J = l + 1/2. We refer the interested reader to [9] for the steps leading from the bi-spinor DE to the following radial equations where now = E − 2/α 2 is the electron energy less its rest mass. This is positive for continuum scattering states. Following Ref. [10], we use the second equation of (6) to replace v κ (r) into the first equation, obtaining By performing the derivative on the r.h.s. and introducing the quantity we find which is the sought pseudo-SE. The H m (r) term is known as the mass relativistic correction, the H D (r) one as the Darwin correction and H so (r) is the spin-orbit potential. Indeed, according to Equation (5), 1 + κ acting on a spin-orbital state is equivalent to −2 s · l. To make contact with the non-relativistic SE, we observe that κ(κ + 1) = l(l + 1) ∀κ. The associated indicial equation is obtained by inserting the assumed low r behavior of u κ (r) ∼ r ρ [a 0 + a 1 r + ...] into Equation (10) and equating to zero the coefficient of the lowest power of r, (r ρ−2 ), taking into account that lim r→0 B(r) ∼ 2r/(Zα 2 ). We have assumed, as usual, that V(r) ∼ −2Z/r at the origin, where Z is the atomic number. One finds which gives ρ = [κ 2 − (Z α) 2 ] 1/2 (taking only the positive root for regular behavior at the origin). This is obviously the same low r behavior of the DE. The solution of Equation (10) provides the upper component of the radial DE. The lower one can be obtained by inserting this solution into the second equation in (6). We shall neglect the lower component, so that only dipole transition matrix elements will be affected, whereas the upper component is sufficient to provide the exact relativistic atomic t-matrix.
The Equation (10) is numerically solved by the Numerov procedure after elimination of the first order derivative (as detailed in Ref. [10]), by applying the method of Gaussian elimination [11]. We also use the more flexible linear-log mesh [12] ρ =αr + β ln r (12) instead of the usual Herman-Skillman (HS) mesh, since the relativistic corrections take effect at much smaller values of r than the initial point of the HS mesh. In Equation (12), α and β are constants to be conveniently determined, as illustrated below. We should, therefore, perform a change of variable from r to ρ before applying the Numerov procedure, which requires constant spacing of the integration variable. A constant mesh size of ρ can be constructed in the interval ρ 0 ≤ ρ ≤ ρ N as follows. The initial value of ρ 0 is chosen according to the empirical formula which gives a minimum value of r already of the order of 10 −6 for Z = 1, taking β = 1, as usually done in atomic physics. The final ρ N is determined by so thatα is given byα where the value of r N is taken to be slightly greater than the MT radius of the atom under consideration. A convenient value of N is given by N = 100 r N , so that Equation (14) determines the mesh size h. The value of r = r(ρ) corresponding to a given value of ρ can be readily found by application of the Newton technique. Ref. [13] describes in great detail the Numerov method of integration together with Gaussian elimination, using a linear-log mesh, in the case of a non-relativistic SE in three dimensions. The method applied here is an extension of that method to the relativistic pseudo-SE of Equation (10) (and all the related approximations as illustrated below) in one dimension, after elimination of the first order derivative. In order to find the atomic t-matrices in the relativistic framework we need the free radial solutions of the DE for positive energies, both regular and irregular. They are obtained from Equations (6) (via Equation (10)) by setting V(r) to zero and are given in Ref. [9].
For a MT potential, as in the non-relativistic framework, the t-matrix is obtained by matching smoothly at the MT radius R s the solution of the DE inside the MT sphere with the free solution outside the sphere where C = α/[2 (1 + 1 4 α 2 )], S κ = κ/|κ| and p 2 = (1 + α 2 4 ) ∼ is the free electron energy. Therefore, indicating by l κ the l-value associated to κ so that l = l κ − S κ , we find As a consequence, t κ is determined solely by the knowledge of the upper component g κ (r). As a check, we have verified that the second expression in terms of the lower component f κ (r) provides the same t-matrix. Note that the smooth matching conditions renormalize the pair of functions g κ (r) and f κ (r) into new ones g n κ (r) and f n κ (r) that are continuous with the free solutions at R s , where the superscript n stands for "normalized".
In the course of a fitting procedure, it is useful sometimes to establish the relevance of the various potential terms in the effective Equation (8). To this purpose, Phagen generates three kinds of t pa κ -matrices, according to the kind of potential approximations (pa) retained. In Equation (10), the potential V(r) ≡ V e f f (r, ) appears in four terms. Retaining, from left to right, only the first, one recovers the usual non-relativistic SE (pa = nr), integrated numerically on a linear-log mesh as for the relativistic calculations. According to Equation (17), the corresponding radial solution generates the non-relativistic t nr κ -matrices. Since the program keeps the same quantities generated with the HS mesh, the matching comparison is a good check that there are no unexpected surprises.
By adding to the non-relativistic potential V(r) the second H m (r) and the third H D (r) term we obtain the scalar relativistic Equation (pa = sr) that provides the t sr κ -matrices in the same approximation.
Finally, when the potential includes also the spin-orbit term H so (r) we recover the fully relativistic potential, including spin-orbit effects (pa = so). In this case, the value of the quantum number κ determines both J and l. Correspondingly, according to the value of the quantum number J = |κ| − 1/2, the program generates two sets of t-matrices, t sojm for J = l − 1/2 (κ positive) and t sojp for J = l + 1/2 (κ negative).
The wave-functions used in the generation of the t-matrices in the various approximations serve also to calculate the corresponding atomic absorption. These are calculated from the knowledge of the atomic Green Function (GF) as σ abs (ω) = −4 π αh ω φ c (r)|(ε · r) * | G + at (r, r ; )| ε · r |φ c (r ) dr dr (18) wherehω is the energy of the incoming photons with complex polarization ε, |φ c (r ) is the initial core state and the scalar product ε · r is between spherical components of the respective complex vectors. Restricting oneself to non-magnetic materials, it is intended that one should perform the sum over the magnetic quantum number projections of the initial core and final states, and all the components of the dipole (l = 1) transition operator C lm (r) = 4π/(2l + 1) Y lm (r). The expression of the upper component of the relativistic GF, in analogy with the non-relativistic one, is where, as usual, r < (r > ) is the lesser (greater) of r, r and h n κ (r) is the irregular upper solution of the DE matching smoothly to the Hankel function h + l (r) of the first kind at R s . The double integral in Equation (18) breaks up into a radial part and an angular part. The radial part has the standard expression [14] 2 where the radial functions g n κ and h n κ take the appropriate form according to the approximation retained.
The angular part varies according to whether in Equation (8) one retains or not the spin-orbit potential. Without spin-orbit interaction, one can work in the l, σ representation, so that, according to the Wigner-Eckart theorem [15] and in the standard notations, Therefore, taking into account the spin and orbital degeneracy of the initial state 2(2l c + 1), transitions to a final state l f = l c + 1 bear a weight 2(l c + 1), whereas transitions to l f = l c − 1 bear a factor 2l c . This result is valid both for the non-relativistic and scalarrelativistic cases. Since they are already coded in, there is no need to change anything in case of scalar-relativistic corrections.

Phase-Shift Calculations for Br 2 and Pb
We have explicitly calculated the t-matrices using different approximations for selected well-known systems and in this paper we present the results for molecular Br 2 and crystalline Pb, for which average interatomic distances are known and previous EXAFS experiments and data-analysis were already carried out (see Refs. [16,17] for Br 2 and [18] for Pb). In those systems, the first-neighbour interatomic distances were rather precisely defined at room temperature using non-relativistic models and compared with previous experiments (R BrBr = 2.289 ± 0.001 Å [17] and R PbPb = 3.496 ± 0.008 Å [18]). We refer to the original publications for details about experimental data, modelling of the X-ray absorption background including double-electron excitation channels and full (non-relativistic) structural analysis.
Structural models for those systems can be easily introduced in GnXAS and calculations of the t-matrices are performed using the present new version of Phagen, containing the relevant relativistic corrections. For Br 2 and Pb phase-shift calculations have been carried out using a complex Hedin-Lundqvist energy-dependent optical potential. In Br 2 , we have considered a typical interatomic Br − Br distance R BrBr = 2.286 Å and R MT = 1.143 Å for the muffin-tin radius (Br 1s excited atom). For Pb (2p 3/2 excited atom), we have considered a face-centered-cubic structure (lattice spacing a = 4.941 Å at 300 K) consistent with the present accepted value around room temperature (see Ref. [18] and references therein). In this case the muffin-tin radius for the excited Pb atom was R MT = 1.222 Å corresponding to about 80 electrons for the charge density integral. In both cases, atomic phase-shifts have been calculated for the photoabsorbing (excited) and neutral atoms, following standard procedures used for muffin-tin MS calculations.
The Phagen program has been used to produce various non-relativistic and relativistic outputs including the scalar relativistic and the full-Dirac calculation. Following the outlines of the previous section, these are denoted, respectively, as 'nr', 'sr' and 'so'. In particular, the fully-relativistic t-matrices and atomic absorptions of the 'so' case are written into two different outputs, classified as 'sojp' and 'sojm' following the notation introduced in the preceding Section 2.
The results of these phase-shift calculations for Br 2 and Pb are reported in Figure 1 for the lower angular momenta l = 0, 1, 2 as a function of the photoelectron kinetic energy (results shown are only for the photoabsorbing Br and Pb atoms, since no substantial differences have been observed for neutral Br and Pb). As expected, minor corrections are found for the Br case, while substantial changes are observed in Pb. Those corrections are especially important at low kinetic energies and low angular momenta as it is shown by differences between non-relativistic and scalar-relativistic phase shifts. Differences up to about 0.5 rads are found in Pb for l = 0 in the entire energy range, although relativistic corrections tend to be less important for increasing l. This trend is more evident looking at the difference curves |δ (nr) l − δ (sr) l | shown in the lower panels of Figure 1. The energy dependence of the difference curves is rather smooth, showing only weak oscillations, but the increasing importance of relativistic corrections for low angular momenta is clear. This is in line with what is known to take place for electrons bound within central-field potentials [19].

Results and Discussion
The atomic t-matrices calculated by Phagen using non-relativistic and relativistic approximations have been used to compute the irreducible two-body γ (2) and three-body γ (3) (and η (3) ) MS signals relevant for the final XAFS data-analysis carried out by GnXAS. We refer to the original publications [1,2] for details about this procedure. In Figure 2, we report the results for the calculations of the single-scattering signals χ 2 (first MS contribution to the two-body term γ (2) ) related to Br 2 and Pb, as a function of the photoelectron wave-vector k. Those χ 2 signals are related to the first-neighbor shell in the two systems, corresponding to typical interatomic distances R BrBr = 2.286 Å and R PbPb = 3.494 Å.
Looking at Figure 2, upper panels, it is possible to observe that the χ 2 signal in Br 2 is visibly affected by relativistic corrections only at moderate wave vector values k < 4 Å −1 , while substantial modifications are found in Pb at least up to 8 Å −1 . The χ 2 signals can be expressed as the product of amplitude and phase functions χ(k) = A(k) sin(φ(k)) (see for example Refs. [1,2]), both generated by GnXAS and shown in the center panels of Figure 2 (as usual, phase functions are shown removing the linear distance factor: Ψ(k) = φ(k) − 2kR). Those two functions are related to the backscattering amplitude and to the phase shifts experienced by the photoelectron and are obviously affected by the approximations used in the calculation. While differences are very limited in Br 2 , we observe some important effects in Pb. In this case, the amplitudes shows fluctuations and minima close to zero in the typical region used for XAFS data-analysis (k > 4 Å −1 ). The k regions for which amplitudes are near zero (highlighted in figure) or shows minima, are also characterized by changes in phase and may appear as typical beats in the XAFS χ(k) signal. Figure 2. Results of calculations of the single-scattering χ 2 (k) signals in Br 2 molecules (left panels) and in solid Pb (right-hand panels, first shell) as a function of the photoelectron wavevector k. From top to bottom, the total single-scattering signals χ(k) = A(k) sin(φ(k)), amplitudes A(k), and total phase functions Ψ(k) = φ(k) − 2kR (where R is the first-neighbour distance), are reported. Nonrelativistic (nr) and scalar relativistic (sr) corrections are indicated, respectively, by blue and orange (continuous) lines; spin-orbit relativistic (sojp and sojm) are denoted, respectively, by green and red (dashed) lines. The k regions for which amplitudes are near zero (highlighted) or show minima, are also characterized by changes in phase and may appear as typical beats in the XAFS χ(k) signal. The lower panels display the difference of 'nr' and 'sr' reported as normalized residuals, A (sr) −A (nr) . Deviations well below 5% are found for Br 2 in the typical XAFS range k > 4 Å −1 , while substantial corrections are needed for Pb also in this high-energy regime. Average deviations were found to be 6% and 25% for Br 2 and Pb, respectively. Corrections may be critical near amplitude minima and in the near-edge region (low k values), reaching deviations more than 30%.
We have performed a full structural refinement of the Pb L 3 -edge experimental signal [18] of face-centered-cubic (fcc) crystalline Pb at room temperature (299 K), using the GnXAS data-analysis method and the new Phagen program including relativistic corrections. The individual irreducible multiple-scattering signals related to a typical fcc structure [1,2] were calculated up to the fourth coordination shell. The most important contribution is of course that related to the first neighbors two-body γ (2) 1 signal, while MS effects are limited in Pb due to a combination of effects: (i) large atomic vibrational amplitudes at room temperature for this low-melting point metal; (ii) large interatomic distances; (iii) short lifetime of the excited state (effect of the ∼3 eV HWHM core-hole width and electron mean-free-path). However, weak MS high-frequency signals can be introduced, related to first-neighbor equilateral configurations γ An important aspect of the XAFS data-analysis is the modeling of the atomic background that for the Pb L 3 -edge is particularly complex, being affected by the presence of multielectron channels involving 4 f electrons [20]. In GnXAS, XAFS data analysis is carried out including realistic models for the absorption background that are refined with the structural signal (fitheo refinement program). The absorption background used in the present structural refinement includes a smooth polynomial spline and the contribution of double-electron excitation 2p4 f channels.
The improvement in XAFS refinements obtained using relativistic approximations has been clearly assessed by our tests and allowed us to obtain a very good agreement with the experimental data, using a realistic background model of Pb L 3 data and structural values compatible with the known crystalline structure. In fact, the quality of the refinement measured by the residual value R (see Ref. [2]) was found to improve by a factor of 2, a statistically significant value. The agreement of the scalar relativistic approximation with experimental data has been always found to be better than that of non-relativistic calculations, for any tested background model. The better quality of the relativistic XAFS calculations can be appreciated looking at Figure 3, where we compare total best-fit multiple-scattering (MS total, scalar relativistic) and experimental signals. The agreement with experimental data is particularly improved using relativistic approximations in a wide wave-vector region k < 6 Å −1 . Similar good results were found using the MS signals generated by spin-orbit relativistic (sojp, see above) phase-shifts. We have also verified that the best-fit structural parameters found in this study are compatible with previous results. In solid Pb at room temperature, the first neighbor distance distribution can be reproduced using a Γ-like function (see Ref. [18] and refs. therein) which depends on four parameters: N, the number of neighbours (coordination number, 12 in this case); R, the average bond distance; σ 2 , the distance variance (the so-called Debye-Waller-like disorder factor); β, a dimensionless asymmetry parameter which tends to zero at low temperatures (limit of a Gaussian function). For scalar relativistic MS calculations, the best-fit structural parameters resulted to be R = 3.498 Å, bond variance σ 2 = 0.032 Å 2 , β = 0.38, in line with previous results [18].
4 shells corresponding to first-neighbors configurations at 90, 120 and 180 degrees. The improvement in XAFS refinement obtained using scalar relativistic phase-shift calculations can be appreciated looking at the comparison between the total multiple-scattering (MS total) and experimental (expt. data, blue) signals and the residuals (red color).

Conclusions
In this paper, we have described the inclusion of relativistic corrections in the multiplescattering calculations used for XAFS data-analysis within the GnXAS suite of programs. This task has been accomplished by realizing a new code for the program generating the scattering t-matrices (Phagen), without altering the basic structure of the programs. The modification incorporates a pseudo-Schrödinger equation containing several terms including scalar relativistic and fully relativistic spin-orbit potentials. The new scheme for phase-shift calculations has been put to a test in two known molecular and crystalline cases: molecular bromine Br 2 and crystalline Pb. Non-relativistic, scalar relativistic and fully relativistic scattering t-matrices (and phase shifts) have been calculated for the two cases showing the importance of taking into account relativistic effects especially for Pb. Phaseshift corrections are extended to the whole photoelectron energy range and have found to be more important for low angular momenta. Therefore, relativistic corrections may be less important for high photoelectron energies where high angular momenta components are needed to describe the final state. Relativistic corrections were found to be much more pronounced in Pb than in Br (up to about 0.5 rads in δ 0 (E) for l = 0 in Pb) as expected due to the higher atomic number.
Multiple-scattering calculations have been then carried out for the Br K-edge and Pb L 3 -edge in molecular bromine and solid lead showing that clear modifications are related to relativistic effects especially at lower photoelectron energies (or wave-vector). In particular, large deviations in amplitudes and phases are observed in a wide vector range k < 8 Å −1 in Pb, while corrections are limited to k < 4 Å −1 in Br. Amplitude deviations are found to be 25% and 6% for solid Pb and Br 2 , respectively. The new multiple-scattering signals generated using relativistic approximations for the scattering t-matrices have been used for direct structural refinements of the XAFS experimental signal of crystalline Pb. We have observed significative improvements in the agreement between the experimental and calculated XAFS signals using relativistic phase shifts, quantified by a decrease in the residual value up to a factor of 2 with a Pb-Pb distance distribution compatible with previous works. In view of these results, we can conclude that the new scheme implemented in GnXAS for relativistic MS calculations can be safely used for XAFS structural refinements of systems containing atoms for which relativistic corrections are important.