Extracting Dynamic Information from EXAFS: Simultaneous Analysis of Multiple Temperature-Dependent Data

A new approach to the extraction of dynamic information from extended X-ray absorption fine-structure (EXAFS) spectra has been developed. With this method, a complete set of temperature-dependent spectra are fit simultaneously to one of a variety of pair-distribution functions. Distributions are calculated in r-space using the appropriate absorber–scatterer pair potential. The temperature-dependent EXAFS spectra are calculated by summing k-space models over a range of distances and angles weighted according to the relative contribution of each geometry to the distribution. This approach allows refinement of data using a full multiple-scattering analysis with only modest computational time.

1. Introduction

Over the last twenty years, extended X-ray absorption fine structure (EXAFS) has come to be recognized as one of the premier tools for determining local atomic structure. For most applications, the information of interest is the static structure. In this case, the EXAFS can be defined as

where χ(k) is the fractional modulation in the absorption coefficient above the edge, k is the photoelectron wavevector, F(k) is the theoretical backscattering amplitude, r is the scatterer distance, λ(k) is the photoelectron mean free path, φ(k) is the phase shift encountered by the photoelectron on passing through the potentials of the absorber and scatterer, and g(r) is the radial distribution function (Teo, 1986).

The EXAFS equation can be simplified by assuming small disorder in the system (Beni & Platzman, 1976). This allows the pair-distribution function to be approximated as a Gaussian distribution, giving the more familiar

where N_s is the number of scatterers, S₀² is an inelastic loss term, A_s(k) is the backscattering amplitude and is a disorder term also known as the Debye–Waller factor, and the sum is taken over every shell of scatterers.

In (2), the Debye–Waller factor gives a measure of the disorder in the absorber–scatterer interaction. Disorder can arise either from static disorder (i.e. from a range of bond lengths within a single shell of scatterers) or from the dynamics of the system (i.e. the motion of the absorber–scatterer pair). From temperature-dependent EXAFS measurements it is possible, at least in principle, to extract information about the dynamics of a system.

In cases where the disorder (static, dynamic or both) is large, analysis using the traditional EXAFS equation (2) can lead to errors in bond length, coordination number and as a result of breakdown of the small disorder assumption (Hayes et al., 1978; Eisenberger & Brown, 1979; Crozier & Seary, 1980; Balerna & Mobilio, 1986). Large disorder can occur in a variety of situations but is particularly likely to be important in temperature-dependent measurements where the temperature is high compared with the characteristic energy of an absorber–scatterer interaction. This frequently results in broad anharmonic pair-distribution functions that cannot be described by (2). As part of a study of dinuclear metal clusters, we have developed a computational approach to extract dynamic information from EXAFS in cases where anharmonicity is important.

2. Previous analyses

2.1. Temperature dependence

Many approaches have been used to describe temperature-dependent effects in EXAFS. For infinite lattice systems, the temperature-dependent mean-square disorder is often calculated from the phonon distribution. Two of the models commonly used to describe the distribution of phonon modes are the Einstein model (one characteristic phonon frequency, ω_E) and the correlated Debye model (a distribution of frequencies from 0 to ω_D) (Beni & Platzman, 1976; Sevillano et al., 1979). These approaches require the somewhat arbitrary assumption of a particular phonon mode distribution and can only account for Gaussian distribution functions. These are most useful for near-perfect crystal lattices and cannot model systems with severe disorder.

A more involved approach for analyzing temperature dependence in EXAFS is to include contributions from all of the vibrational normal modes that cause relative motion of the absorber and scatterer (Boland & Baldeschwieler, 1984). In this method, the disorder of the system is described in terms of amplitude vectors and the mean-square amplitude of vibration, and the sum of atomic motion is taken over all normal modes. The equation used in this derivation includes single scattering pathways as in the familiar EXAFS equation (2), but can be expanded to include double scattering pathways and higher-order scattering processes as necessary, although the treatment becomes increasingly complex. This approach readily describes the scattering processes of simple systems and distinguishes the effects of bond angle and temperature on multiple scattering and EXAFS amplitude. In principle, this could be extended to include anharmonic contributions to one or more of the normal modes. However, only harmonic distributions have been analyzed.

2.2. Anharmonic distribution functions

One of the first papers to recognize the importance of anharmonicity was a study of Zn EXAFS near the melting point of Zn (Eisenberger & Brown, 1979). An anomalous decrease in the Zn—Zn distance was observed as the temperature was increased to the melting point, in contrast with the bond length increase that was expected based on bulk thermal expansion. This apparent decrease in Zn—Zn distance was attributed to the loss of low-k information (k < 3 Å⁻¹, where near-edge effects prevented EXAFS analysis). The low-k data contains much of the phase information needed to describe the anharmonic Zn—Zn distribution function. Near the melting point of Zn, the Zn—Zn distribution function is strongly skewed to high r, leading to incorrect fitting results when modelled with a symmetric distribution function. Using a simple correction term for the scattering phase, it was possible to model the data and extract limited information about the asymmetry of the distribution.

An alternative approach for non-Gaussian distributions is to analyze the data in r-space (Hayes & Boyce, 1980; Hayes, 1984). Analysis in r-space allows the deconvolution of scattering shells radially without having to re-Fourier transform the data, thereby further introducing distortions due to the transform windows. Analysis in r-space also reduces the number of transforms necessary to compare the data with model parameters. Most importantly, however, is the freedom to utilize any type of atomic potential to describe the absorber–scatterer interaction.

It is possible, in principle, to use anharmonic distribution functions to model EXAFS. The Morse potential,

where D is the depth of the potential well, a = (8π²cμωx/h)^1/2, μ is the reduced mass and ωx is the anharmonicity constant (Morse, 1929), is often used in vibrational spectroscopy to model anharmonic distributions (Wilson et al., 1980). This can be expanded in a power series to give (Miyanaga & Fujikawa, 1994)

This expression has been used in combination with displacement vectors to obtain anharmonic corrections to harmonic displacements in terms of a, D and T for mon­atomic and diatomic chains.

Many model-independent approaches have been used in the analysis of anharmonic distributions (Babanov et al., 1981; Bunker, 1983; Gurman & McGreevy, 1990; Yang & Bunker, 1996), with use of cumulants being the most common. The cumulant expansion (Bunker, 1983) of the Fourier transform of a pair-distribution function is

where λ is the photoelectron mean free path and n = 0, 1, 2 etc. is the number of the cumulant. The even cumulants, C_2n, depend only on the amplitude of the EXAFS while the odd cumulants, C_2n+1, depend only on the phase. For a Gaussian distribution, all cumulants of n > 2 are zero. Cumulants can be used in combination with the ratio method, in which the high-temperature data set of interest is compared with data measured at low temperature (Bunker, 1983). A plot of the natural log of the amplitude ratio versus k² is a straight line for a harmonic distribution but will curve if higher-order cumulants, which arise from anharmonicity, are important. Similarly, a plot of the difference in phase between two data sets versus k can be used to determine the importance of odd cumulants. The cumulants derived in this way can be used to reconstruct the radial distribution function using (5).

One of the difficulties with the cumulant expansion is that it often diverges at high k, and thus the cumulants are only accurate at low k. One solution to this is the `splice' method (Stern et al., 1992) in which the cumulants are used to recreate the missing low-k data but are not used to generate distribution functions directly. The low-k data are instead `spliced' with the measured data at higher k (k > 3–5 Å⁻¹), in principle giving a more accurate EXAFS function over the complete k range. This gives the low-k data that are crucial for analysis of anharmonic distributions without a loss of accuracy at high k. The data can then be used to calculate the radial distribution function.

Cumulant analysis has been widely used in the study of disordered systems. The method, however, has several limitations. Cumulants do not provide a good description of systems with large noise (Stern et al., 1992), and accurate extraction of phase and amplitude information can be overshadowed by uncertainties in the analysis procedure (Monti, 1996). As the disorder increases, the number of cumulants needed to describe the data increases and the usable data range for fitting decreases. Most importantly, cumulant analysis is impossible to use in systems with mixed coordination shells and systems displaying multiple-scattering effects (Dalba & Fornasini, 1997).

All of the approaches to anharmonicity that have been described thus far give correction terms that can be plotted versus temperature to extract information about the temperature dependence of the asymmetry. This is useful for describing anharmonicity, but fails to provide a complete description of the temperature dependence. A comprehensive treatment of both temperature dependence and anharmonicity has been given by Mustre de Leon et al. (1992). In this work, a double-well potential,

where a and b define the depth of the two potential wells and r_eq is defined by requiring the distribution to be continuous, was used to describe the motion of an O atom relative to two Cu absorbers, where r_eq1 and r_eq2 are the centres of the two possible O positions. Empirical EXAFS parameters for A_s(k) and φ_as(k) were used in (1) to model the data. A quantum-mechanically derived radial distribution function, g(r), was fit to the data with the constants a and b and the distances r_eq1 and r_eq2 being refined to give the best fit to the EXAFS. This method can be readily extended to an arbitrary potential. More importantly for the present purposes, this method can accommodate simultaneous analysis of multiple temperature data sets by refining g(r,T) rather than g(r), although this has not previously been reported.

We have developed an approach to the analysis of dynamic information from EXAFS that builds on that reported by Mustre de Leon et al. (1992). We use k-space analysis to facilitate calculation of multiple-scattering effects and to avoid introducing artifacts into the data during Fourier transformation. We have adopted a model-dependent approach in order to allow refinement of all data measured at different temperatures simultaneously. In order to speed up computation, all of the necessary k-space multiple-scattering calculations are performed at one time and this library of spectra is used to construct the EXAFS during the curve fitting. A brief account of this program has been given (Daly & Penner-Hahn, 1997). In the following we provide a more detailed description and evaluation of the approach.

3. Program method

The EXAFS simulations are based on FEFF 6.01 calculations for the amplitude and phase (Rehr & Albers, 1990). In contrast with traditional simulations, which typically use a single FEFF calculation to define φ_as(k,r) and A_s(k,r) and then refine r, N_s and in using (2), we calculate the FEFF simulations for a range of r. For systems where multiple scattering is important, the bridging angle, θ, can also be varied. These are weighted according to an r-space distribution function to give k-space EXAFS. This amounts to replacing the integral in (1) with a summation,

where N_mod is the number of models with spacings of δr = r_i + 1 − r_i. This type of integration by summation calculation of EXAFS is the same as that used in the methods by Mustre de Leon et al. (1992) and originally used by Benfatto et al. (1989) in multiple-scattering calculations. The r-space distribution function, g(r), is refined so as to optimize the k-space simulation relative to the measured EXAFS. A flow chart summarizing the computation is shown in Fig. 1.

Figure 1 Flow chart of the program method.

There are four principal advantages to this approach. First, this method provides the advantages of r-space refinement without the need for Fourier transformation during the fitting. Second, this approach gives a slight increase in accuracy by handling the 1/r² weighting appropriately. It has previously been noted (Brouder, 1988) that in conventional analyses the 1/r² term in (2) is held constant at . For narrow distributions, the use of has no effect on the fit. However, for broad distributions, it leads to a slight over-weighting of the high-r portion of the distribution. This can result in small (typically 0.010–0.015 Å) errors in distance and 2–5% errors in amplitude. A third advantage of our approach is that r-space refinement allows us to model several temperature-dependent data sets simultaneously, optimizing only a single pair potential function. Any type of radial distribution function or, equivalently, any pair potential function can be used to describe the atomic interactions. In particular, asymmetric distributions are straightforward to handle. Finally, multiple-scattering refinements require far fewer computational resources in the present refinement scheme. Scattering pathways with greater than two legs are highly dependent on the geometry of the system. In traditional refinements (either k- or r-space) it is necessary either to assume that the multiple scattering does not change, and thus that a single amplitude and phase function can be used throughout the refinement, or to repeat the multiple-scattering simulation at each step of the refinement. The latter is a computationally intensive process. In contrast, the present approach requires only one time calculation of the complete multiple scattering for a relatively small number of geometries (typically 20–30).

The fit is optimized using a modified Levenberg–Marquardt non-linear least-squares algorithm (Garbon et al., 1979) with a standard goodness-of-fit parameter

where N is the number of data points. In simultaneous refinements of multiple temperature-dependent data sets, it is the high-temperature data sets that contain the most information concerning the distortion of the radial distribution function. Unfortunately, these data sets have the lowest EXAFS amplitude and thus contribute the least to the goodness of fit, equation (8). In this case we have found it advantageous to use the weighted goodness-of-fit parameter F_w in (9), in which each data set contributes equally to the optimization, and the overall goodness-of-fit parameter is summed over each of the M data sets,

The weight, W_j, for the jth data set is taken as the reciprocal of either maximum EXAFS amplitude or the integrated EXAFS amplitude for the data set,

Both methods of weighting were effective at increasing the importance of the high-temperature data sets. This weighting is necessary for accurate modelling. In some systems, asymmetry only becomes evident at temperatures above 100–150 K where EXAFS amplitudes are often reduced by 50% or more relative to the low-temperature data. Without weighting, these systems are best fit with symmetric distributions, since the more symmetric low-temperature data sets contribute most to the goodness-of-fit parameter in (8).

4. Radial distribution functions

A variety of distribution functions have been implemented with this fitting approach. All of the distributions that we have used are in the form of effective pair potentials for a diatomic system. This is, of course, only an approximation of the motion that is possible for a polyatomic system. However, this is expected to be a good approximation for molecular systems in which the temperature-dependent motion is dominated by a single low-frequency normal mode.

4.1. Individual data set analysis

In order to validate our program and to determine the necessary r-space grid for acceptable refinement, we performed initial tests using a symmetric Gaussian distribution function. The Gaussian width, σ², is equivalent to the Debye–Waller factor, , in a traditional k-space calculation, equation (2). A simple, albeit crude, approach to asymmetry is to allow the width of the Gaussian to vary independently on each side of the equilibrium value

where r_eq is the equilibrium bond distance, is the width on the low-r side and is the width on the high-r side. More strongly skewed asymmetric distributions as described by Teo (1986) can also be incorporated.

4.2. Temperature-dependent analysis

Given an absorber–scatterer pair potential function, it is possible to calculate the pair-distribution function, g(r), at every temperature. If the spacing of energy levels is small in comparison with k_bT, a simple Boltzmann factor can be used to convert from potential energy, V(r), to the probability of finding an arbitrary absorber–scatterer distance, g(r), using

This is an appropriate distribution within the classical approximation but is expected to be a poor description of real systems, especially at low temperatures. It is necessary to consider the quantization of the energy levels in order to accurately compare data sets measured at different temperatures. In this case, a harmonic oscillator and a Morse potential function are among the easiest distributions to model.

For the harmonic oscillator (Noggle, 1989), the wave functions are given in terms of the Hermite polynomials, and the energy is E = (v + 1/2)ω, where v = 0, 1, 2, . . . are the vibrational levels and ω = (2πc)^−1/2(k/μ)^1/2 is the vibrational constant for each mode. The populations of each vibrational level can be calculated from the Boltzmann factor

The Boltzmann weighted sum of the squares of the individual wave functions gives the temperature-dependent g(r). The normal mode parameters are calculated using Gaussian 94 (Frisch et al., 1995).

The Morse potential, equation (3), gives a similar expression, but incorporates asymmetry in the distribution. Solving the Schröedinger equation for (3) gives wave functions based on the associated Laguerre polynomials (Morse, 1929). The energy of the vibrational levels is given by E =−D + ω(v + 1/2) − ωx(v + 1/2)^1/2, where ωx = a²h/8π²cμ is the anharmonicity constant and ω = (4Dωx)^1/2.

It is instructive to compare the temperature-dependent analyses for a one-shell traditional optimization, where only two or three parameters are varied, to the fit of a series of seven temperature data sets where a minimum of 14 parameters would be required (r and at each temperature). The present approach provides a much better description of the potential distribution function by relating all temperatures to a single set of parameters.

5. Results

The program is relatively memory intensive. For example, if 50 FEFF models are used to define the EXAFS at 200 points for each of five shells of scatterers, a total of ∼1 Mbyte of storage is required. The benefit of large memory usage is that the optimization is simplified to a straightforward summation with shorter refinement times. For single-shell non-multiple-scattering fits, the current approach ranges from 30–60% faster than conventional refinements (see Table 1). Much larger improvements would be observed for full multiple-scattering refinements.

Table 1
Comparison of computation times

	Time per iteration (s)
Type of fit	Traditional	This paper	% Faster
One shell	0.153	0.110	28
Three shell	0.410	0.140	66
One shell, 8 data sets	8 × 0.153 = 2.05	0.825	60
Two shell, multiple scattering	0.138 + 2.19 min	0.195	99.9

For each shell, FEFF models are calculated for a range of distances. The number of model calculations, N_mod, and the r-space grid, δr, required to accurately represent the distribution is dependent on the breadth, ΔR, of the distribution function to be modelled. Fig. 2 is a generalized distribution function showing the definition of these parameters. We simulated EXAFS scattering using a simple Gaussian damping factor as in equation (2) and modelled these data using a Gaussian distribution as in equation (7). As seen in Fig. 3(a), the root-mean-square percent error drops as ΔR increases relative to σ and is less than 0.1% for ΔR ≥ 4σ. In Fig. 3(b), the error remains less than 1% as long as δr ≤ σ. In the absence of static disorder (Teo, 1986), distribution widths as small as σ ≃ 0.033 Å could be found for a shell of scatterers strongly bonded to the central atom (estimated for a metal–ligand stretch with μ = 50 a.m.u. and  = 300 cm⁻¹). In practice, the σ values that we have found have ranged from 0.05 to 0.12 Å. To provide a margin of safety in simulations of experimental data, we have kept δr < 0.5σ_min and ΔR > 4σ_max, giving models spaced every 0.025 Å and ±0.5 Å from r_eq.

Figure 2 Description of the parameters used in the calculation for a Gaussian distribution with width σ centred at req. Model calculations are performed for a grid of Nmod, represented by circles, at R values from req − ΔR to req + ΔR with a spacing of δr.

Figure 3 Comparison of EXAFS simulation with this method using a Gaussian function. (a) Logarithm of the root-mean-square (r.m.s.) percent error versus the width of the distribution in units of σ for a model spacing interval of 0.25σ. (b) Logarithm of the root-mean-square (r.m.s.) percent error versus the spacing of the models in units of σ for a distribution width of 4σ.

As an illustration of the benefits of this analysis protocol, preliminary results for a rhodium chloride dimer are described below. Complete results for this system will be published later (Daly et al., 1998). The system consists of two Rh(I) cations bridged by two chloride anions forming a four-membered ring. Each Rh atom is also bonded to two ethylene molecules in a pseudo square-planar environment. The system exhibits significant temperature dependence in the EXAFS scattering, especially in the Rh⋯Rh interaction (Barnes et al., 1995).

Independent fits to each of the different data sets were sufficient to model the low-temperature data but did not work well at high temperatures. The fits to the 4 K and 300 K filtered Rh⋯Rh interaction using a symmetric Gaussian function are shown in Fig. 4. Slightly better fits could be obtained for the 300 K data by using the ad hoc asymmetric Gaussian. However, even these fits were not particularly good. Worse, the widths of the two half Gaussians did not vary monotonically with increasing temperature, thus illustrating the difficulty of this approach for modelling temperature-dependent data.

Figure 4 The filtered Rh⋯Rh interaction (full line) for Rh2Cl2(ethylene)4 and the symmetric Gaussian fit (dashed line) at (a) 4 K and (b) 300 K.

When simultaneously fitting all temperatures, the Morse potential is significantly better than the harmonic oscillator potential as shown in Fig. 5. The need for a weighting factor, equation (9), is evident in Fig. 5 where the amplitude of the EXAFS at 300 K is an order of magnitude smaller than that at 4 K. The Morse potential gives a fit to the 4 K data that is comparable with, and perhaps slightly worse than, the corresponding symmetric Gaussian fit [compare Figs. 4(a) and 5(a)]. The slightly worse Morse potential fit reflects the fact that this Morse potential fit has been optimized against seven different data sets simultaneously. The Morse potential has a slightly larger amplitude than the experimental data, but continues to fit the phase and amplitude up to k = 17 Å⁻¹ where the symmetric Gaussian (Fig. 4) and the harmonic oscillator (Fig. 5) both fail. Fitting statistics show that the Morse potential, which models all seven data sets at once, describes the motion of the system better than individual fits using symmetric Gaussians. These results highlight two of the main advantages of the fitting protocol: asymmetric distributions can be used to accurately describe molecular systems and multiple temperature data sets can be analyzed simultaneously to allow the extraction of dynamical information.

Figure 5 The filtered Rh⋯Rh interaction (full line) for Rh2Cl2(ethylene)4, the weighted harmonic oscillator fit (dashed line) and the weighted Morse potential fit (dotted line) at (a) 4 K and (b) 300 K.

We are currently exploring the ability of the temperature-dependent EXAFS analyses to distinguish between different bridging geometries (Daly, 1998). Possible applications include the Mn site in the oxygen-evolving complex of Photosystem II. EXAFS data collected at 10 K show a ∼3.3 Å Mn⋯Mn interaction (Yachandra et al., 1996; Penner-Hahn, 1998) that was not seen with data collected at 150–190 K (Cole et al., 1987) suggesting that the 3.3 Å interaction may show pronounced temperature dependence. Temperature-dependent studies of this system may be helpful in elucidating the nature of the Mn—Mn bridging ligands.

6. Conclusions

An EXAFS analysis method that permits rapid complete analysis of asymmetric and temperature-dependent systems has been described. The protocol has been implemented using a Morse potential, but could be readily generalized to other pair-distribution functions. The calculation handles multiple scattering more effectively and optimizes in a shorter period of time than conventional analyses.

Some other systems of interest for future applications include dinuclear iron–sulfur centres, such as those found in feredoxin where the Fe—S—Fe bridges exhibit temperature dependence, the dinuclear Fe site of hemerythrin, which shows different Fe⋯Fe distances at different temperatures (Riggs-Gelasco et al., 1995), and the Mn cluster in the photosynthetic oxygen-evolving complex, where the long-range Mn interactions show temperature-dependent behaviour.