Determining the radial pair distribution function from X-ray absorption spectra by use of the Landweber iteration method

The Landweber iteration approach is used to construct the radial pair distribution function (RPDF) from an X-ray absorption (EXAFS) spectrum. The physical motivation for the presented investigation is the possibility to also reconstruct asymmetric RPDFs from the EXAFS spectra. From the methodical point of view the shell fit analysis in the case of complicated spectra would be much more eased if the RPDF for the first shell(s) are computed precisely and independently. The RPDF, as a solution of the fundamental EXAFS integral equation, is examined for theoretical examples, and a detailed noise analysis is performed. As a real example the EXAFS spectrum of curium(III) hydrate is evaluated in a stable way without supplementary conditions by the proposed iteration, i.e. by a recursive application of the EXAFS kernel.

1. Introduction

The atomic pair distribution function describes the density of interatomic distances in matter. The radial pair distribution function (RPDF), which is independent of orientation, is of special practical importance for the analysis of EXAFS spectra. It is a major descriptor for the atomic structure of amorphous materials and liquids.

A RPDF has by definition the following properties (Babanov et al., 1981):

The first condition is the positivity constraint, which is valid for each probability distribution. The second condition states that, for large distances where the particles are uncorrelated, g(r) tends to 1.¹ Thirdly, for distances smaller than the sum of ionic radii of the two atoms (r_min), where the impenetrability of the particles becomes apparent, g(r) is equal to zero. Fourthly, g(r) is normalized to 1.

The conditions in (1) must be incorporated in the construction procedure of the function g(r), thus ensuring that the sought solution is a normalized RPDF.

The general EXAFS integral equation for a one-component system, i.e. for a system with only one type of backscattering atoms, has the form

where the kernel

results from the matrix element of the electron scattering according to Fermi's Golden Rule (see Stern, 1974; Lee & Pendry, 1975; Teo, 1986). Here χ(k) is the normalized oscillating part of the measured EXAFS spectrum, k is the electron wavevector, and r is the distance from the central to the backscattering atoms. The function's backscattering amplitude F(k,r), the central atom phase shift δ(k), the backscattering atom phase shift φ(k,r) and the mean free path of the photoelectron λ(k) are obtained from the curved wave approximation using the FEFF code (Ankudinov et al., 1998). The reduction factor is set to S₀² = 1 in this article.

Equation (2) may be simplified in spherical coordinates due to the fact that the kernel depends only on the radius r. In this case the volume element is dV = 4πr²dr. For EXAFS data analysis it is more demonstrative and more adapted to introduce a radial particle distribution n(r) instead of g(r),

with

so that N is the number of backscattering atoms in the spherical shell between r₁ and r₂, i.e. the analogue to the coordination number.

The EXAFS equation (2) now takes the form

The function n(r), depending on the radius, is unknown and contains all the information about the EXAFS spectrum.

The `shell fitting' procedure is the most frequently applied method to evaluate EXAFS spectra in practice. In this procedure it is usually assumed that the particle density is represented by separated Gaussian peaks where the half width of the Gauss curve is the Debye–Waller factor σ²,

Thus, the integral in the EXAFS equation (6) can be resolved analytically and, assuming some well justified approximations, will be reduced to the EXAFS equation for a number of coordination shells,

Here, N_m is the coordination number and is the Debye–Waller factor of the backscattering atoms in the mth shell.

While the shell fitting procedure approximates n(r) by symmetrical Gaussians, knowledge of the exact RPDF, or equivalent n(r), is of greatest importance, e.g. for the investigation of anharmonic corrections caused by thermal vibrations (Yang & Joo, 1998; Dalba & Fornasini, 1997; Crozier et al., 1988). This is the main physical reason to determine the complete RPDF independently.

Equation (6) may be rewritten in a condensed form as

This is a Fredholm integral equation of the first type connecting the known spectrum χ(k) with the unknown function n(r). The inversion of equation (9), i.e. the determination of n(r), is a so-called ill-posed problem. From the three conditions for a well posed problem, suggested by Hadamard (1932), i.e. (i) the existence, (ii) the uniqueness and (iii) the stability of the solution depending continuously on the data, condition (iii) is most often violated owing to the uncertainties in the data, i.e. experimental error. The two most important well established methods for the solution of this type of integral equations, like equation (9), are Tikhonov's regularization (Tikhonov & Arsenin, 1981) and Landweber's iteration (Landweber, 1951) procedures. The principal difference between the methods is the following. The integral equation, like equation (9) or (10), is solved using the Tikhonov regularization by multiplication with the inverse matrix (i.e. the inverse integral operator), whose irregularities are avoided by addition of a regularization matrix. The Landweber iteration (see details below) does not use the inverse matrix and therefore does not contain irregularities. Both methods are independent and different.

In the mathematical literature some comparisons of the two methods exist (Kirsch, 1996; Latham, 1998). They stated that for non-noisy model spectra both methods lead to the same result. A general comparison of both methods for the case of realistic noisy data is not possible because of the different character of data errors for each individual case.

It has been confirmed that for the inversion of equation (9), and for data with experimental error, more stable results are supplied by the Landweber iteration (Kirsch, 1996), which on the other hand needs more iteration steps. The latter statement agrees with our initial experiences of the analysis of EXAFS spectra using the Landweber iteration. However, no direct comparison of both methods, applied to real noisy EXAFS spectra, exists. In §3 a qualitative noise analysis for the Landweber iteration using model spectra is presented.

Babanov and co-workers (Babanov et al., 1981, 2007; Ershov et al., 1981; Kunicke et al., 2005; see also Yang & Bunker, 1996; Khelashvili & Bunker, 1999; Bunker et al., 2005; Ageev et al., 2007) have applied the Tikhonov regularization to the analysis of EXAFS spectra. Using special assumptions, adequate RPDF of various systems were reconstructed. In this context no systematic noise analysis exists.

Lee & Yang (2006) presented a new analytical method, the Matched EXAFS RPDF Projection (MEPP), to obtain the RPDF of an EXAFS spectrum via Fourier filtering.

An interesting article by Yamaguchi et al. (1999) deals with the application of the Landweber iteration to the inversion of the EXAFS equation. Their work focuses on the application of the wavelet Galerkin method for the reconstruction of the integral kernel A(k,r) in equation (9). The Landweber iteration method contains two parameters (see below): a convergence parameter α and a stopping number ν_opt, which indicates the optimal number of iterations. Both parameters were set by Yamaguchi et al. (1999) by `trial and error' to 5000.

In the present work we present a methodical introduction and some examples for the analysis of EXAFS spectra using the Landweber iteration method, including the convergence parameter α, which depends on the norm ||A^*A||², and the stopping parameter ν_opt, which depends on the error in the data. The calculation of the matrix A is performed by the use of the FEFF8.2 program (Ankudinov et al., 1998). We also give a physical interpretation of the iteration steps of the Landweber method for the inversion of the EXAFS equation. In favor of simplicity and clarity, the article does not contain analysis of multi-component spectra and of multiple scattering effects.

2. Landweber's iteration for a one-component backscattering system without multiple scattering effects

2.1. Landweber's iteration and constraints

2.1.1. Definition

The basic EXAFS equation (6), in operator form

when multiplied by αA^*, is identically rewritten as

Based on this equation, with initial condition n₀ = 0 and a parameter α > 0, Landweber (1951) suggested the iteration scheme

where ν is the iteration number. For mathematical details, see for example the textbooks by Kirsch (1996) and Louis (1989). This iteration procedure can be interpreted as the steepest descent algorithm to minimize the functional ||An − χ|| with α as convergence parameter. We remark that the acceleration term (χ − An_ν) is the difference spectrum between the original spectrum χ and the last iteration spectrum An_ν.

2.1.2. Convergence parameter

It has been proven, e.g. by Louis (1989), that for α in the range 0 < α < 2/||A^*A||, the Landweber algorithm is equivalent to a linear regularization procedure. This means that the iteration converges to an optimal solution n_νopt at an optimal number of iterations ν_opt. In the following examples, the convergence parameter α^EXAFS = 1/||A^*A|| will be chosen.

2.1.3. Stopping parameter

While α defines only the `speed' of the iteration, the optimal stopping parameter ν_opt must be defined for each problem separately. On one hand, iterating too long causes the experimental error in the data to increasingly disturb the solution. On the other hand, too few iterations will lead to a loss of resolution. This convergence behavior is called semi-convergence (see also Fig. 1a).

Figure 1 Definition of the stopping parameter for an idealized theoretical spectrum which contains artificial data error. (a) The difference between the iterated (nν) and the theoretical (ntheor) RPDF, demonstrating clearly the semi-convergent behavior of the Landweber iteration. (b) The difference between the iterated and the theoretical spectrum plotted over the iteration parameter ν. τδ is the total error.

Apparently, the stopping parameter ν_opt depends on the error in the data (δ), e.g. caused by noise, spline errors etc. All other sources of uncertainty, e.g. the choice of α and the model matrix A, are described in the literature (Elfving & Nikazad, 2007; Kirsch, 1996) by an additional parameter τ. There, it is proved that the iterative method is convergent for 0 < τ < 2. Our initial model calculations (see below) show best results for τ ≃ 1.

If the total data error τδ were known, the stopping rule defining ν_opt would simply be

This means that the solution is such that the norm of the residual ||An − χ|| is of the same order as the total data error τδ (Fig. 1b). We have used this condition several times with simulated noise, i.e. with known δ and τ ≃ 1, to validate the following procedure developed for realistic noisy data.

2.1.4. L-curve criterion

In the case of real EXAFS measurements, the total error τδ is not known. To determine the stopping parameter in the Landweber iteration we adapted the concept of the L-curve criterion, derived for the determination of an optimized regularization parameter for the Tikhonov regularization (Hansen, 1992, 2001; Hansen & Oleary, 1993; Kunicke et al., 2005). In general, the L-curve is a log–log plot of the norm of a regularized solution versus the norm of the corresponding residual. Both norms depend parametrically on ν_opt. The resulting curve shows a shape like the capital letter L (Fig. 2). The optimal regularization parameter is chosen as the maximum of the L-curve curvature, corresponding to the corner of the L. As ν increases up to the value ν_opt, the norm of the residual decreases much faster than the norm of the solution increases. A region of ν values follows where the norm of the solution starts to increase dramatically. This means that the value ν_opt, corresponding to the corner of the L-curve, gives an optimum in the effort to minimize simultaneously the norm of the residual and the norm of the solution.

Figure 2 A typical idealized L-curve. The function nν is the RPDF after the νth iteration, χ is the measured EXAFS spectrum, and A is the EXAFS integral operator. The iteration number ν increases from right to left.

In the ideal case the L-curve is defined by a smooth computable formula. Then the maximum of the curvature of the graph log||n_ν|| over log||An_ν − χ|| gives us the value of ν_opt.

The curvature c_ν of the L-curve is obtained using equation (14) where x = log||An_ν − χ|| and y = log||n_ν||. The prime denotes the derivative with respect to ν if both functions depend parametrically and continuously on the parameter ν,

Unfortunately, in the case of realistic data the L-curve does not have the clear form shown in Fig. 2, and is limited by the knowledge of only a finite set of points. As a first approximation for the numerical differentiations in equation (14) we used the integral values of the parameter ν as nodes. The corresponding stopping parameter is further denoted by ν_curv. The use of ν_curv resulted in satisfying curvature maxima for all of our examples.

A complete discussion of the numerical problems in locating the corner of the L-curve, e.g. by use of cubic spline curves, is given by Hansen & Oleary (1993).

2.1.5. Positivity constraint

Generally, different constraints are used for the construction of solutions of ill-posed problems in order to improve their stability and convergence. A common method is the utilization of positivity constraints (Kirsch, 1996). The positivity constraint [see equation (1)] is included in the Landweber iteration steps because it is an inherent part of the definition of the RPDF. It should be noted that the consideration of the positivity constraint leads to a slight degradation of the smoothness of the L-curve. Instead of the terminus `positivity constraint' the terminus `projected inversion' is equivalently used in the relevant literature (Bunker, 2009).

2.2. Physical interpretation of the transformation of the spectrum by the transposed matrix A^*

In the first step of the Landweber iteration (n₀ = 0) and with α = 1, equation (12) has the form n₁ = A^*χ, where χ is the EXAFS spectrum. In each further iteration step the operator A^* acts on a difference of two EXAFS spectra (χ − An_ν). We consider the effect of the transposed matrix A^* acting on the spectrum χ(k),³

and rewrite (15) according to (6) in its continuous form,

from which the following interpretation of transformation (16) is obvious. The oscillating part of the integration kernel, sin(2kr), is the imaginary part of the Fourier transform, which results usually in a very precise determination of the distances of the backscattering atoms. The argument of the sine function is phase corrected by the phases of the central and backscattering atoms, and an amplitude window function is formed by the atomic amplitudes and the mean free paths. Therefore the transform (16) may be interpreted as a windowed and phase-corrected Fourier transform, whereas the window and phase functions originate from the underlying model.

2.3. Construction of the matrix A

For the discretization of equation (6) by an equal-spaced step function, n(r_m) is approximated by a linear function in each Δr interval.

Consequently, the right-hand side of (8) represents the columns of the matrix A_im, which correspond to the k-dependence of the kernel for each distance r_m. If the coordination number in each interval is set to N_m = 1, we obtain the columns of A_im,

Therefore, for the distances r_m, the columns of A_im are calculated directly from the FEFF program (the chi.dat files without Debye–Waller factors) in the form of the theoretical spectra. We used step sizes of Δr = 0.02 Å or Δr = 0.002 Å in our calculations.

Most algorithms use for the description of an RPDF peak the theoretical scattering functions calculated for one radial distance supplied by the FEFF structural model. It is well known that the theoretical scattering functions depend on the radial distance. Therefore the theoretical scattering functions must be calculated for a narrow grid of radial distances in the r-interval of the RPDF. If the theoretical scattering functions are calculated only for one radial distance then, especially for strong disordered systems with a broad RPDF, the calculated RPDF might be incorrect.

2.4. Determination of the energy threshold shift ΔE₀

The k-axis, which is the basis for the computation of the matrix A, is defined by

where is Plank's constant, m is the electron mass and E is the photon energy. E₀ is the threshold energy, i.e. the minimum energy to free the electron. To define the k-axis of the real experiment, usually a correction term ΔE₀ is introduced. ΔE₀ is a priori unknown and usually it is determined by a shell fit of the spectrum before an inversion method is applied. If the shell fit model is incorrect then the determined ΔE₀ is erroneous and an inverse method would give a wrong RPDF. This leads in turn to the conclusion that the recent inversion procedures are only applicable if the structural model for the investigated system is already known.

In our approach, ΔE₀ is defined by the iteration procedure itself, without recourse to fit calculations or models. The only assumption is that for each ν an optimal ΔE₀ exists, and that over a broad range around ν_opt the best values of ΔE₀ are similar. By using this assumption for a given ν, the optimal ΔE₀ is determined simply by a one-parameter search, using the minimization of the standard deviation SD_ν, defined as

SD_ν is the norm of the residual of the EXAFS equation (10), and depends on the iteration number ν, normalized to the maximum number of data points in r-space (i_max).

We start with arbitrarily chosen values ν_start and ΔE₀. With this start energy threshold shift (e.g. ΔE₀ = 0) the columns in the matrix A_im are re-calculated on the new k-axis, based on equation (18). Then the Landweber iteration is performed up to ν_start and the new SD_ν is determined. In the next cycle, ΔE₀ is modified so that SD_ν decreases. If ΔE₀ converges to a constant value, the approximation procedure based on ν_start is stopped. This first optimal ΔE₀ is already near the true value, and we calculate with this ΔE₀ the new L-curve and the new ν_curv. Then the algorithm to find ΔE₀ is repeated by using ν_curv as the new ν_start. If the actual ν_start is equal to the next ν_curv, then ΔE₀ and SD_ν have reached their final values.

This procedure will be illustrated below, where ΔE₀ is determined for the noise contaminated spectrum of example #6.

3. Application of Landweber's iteration to model spectra

In this section, we examine theoretically how the Landweber iteration method responds to noisy `close-to-reality' spectra. Noise is added in ten steps with a noise generator to a noise-free theoretical spectrum (Table 1, example #1). Subsequently, the ten noise-contaminated spectra are analyzed. The change of the structural parameter (N, r, σ²) and the stopping parameter are inspected as a function of the amplitude of the added noise.

Table 1
Summary of the structural parameters N, r, σ², the theoretical (ν_theor) and L-curve defined stopping parameter (ν_curv) for the noise-contaminated model spectra

The sign ∞ for example #1 remarks that, for a spectrum without error, an a priori limitation rule for the iteration does not exist, as explained in §2.1.3. Standard deviations of the parameter are given in parentheses.

Example	Cnoise × 10^4	N	r (Å)	σ^2 × 10^3 (Å^2)	νtheor	νcurv
#1	0.00	1	2.5	5	∞	∞
#2	0.50	0.999 (6)	2.5000 (5)	4.99 (7)	84	93
#3	0.75	0.995 (7)	2.4999 (6)	4.93 (8)	75	79
#4	1.00	0.991 (8)	2.4999 (6)	4.88 (9)	70	72
#5	1.50	0.98 (1)	2.4998 (8)	4.8 (1)	62	64
#6	2.00	0.98 (1)	2.500 (1)	4.7 (1)	57	59
#7	2.50	0.97 (1)	2.500 (1)	4.6 (2)	53	56
#8	3.00	0.97 (2)	2.499 (1)	4.5 (2)	50	55
#9	3.50	0.97 (2)	2.499 (2)	4.4 (2)	48	53
#10	4.00	0.96 (2)	2.499 (2)	4.3 (2)	46	51

Additionally, the difference between the theoretical stopping parameter ν_theor (equivalent to ν_opt in Fig. 1a) and the stopping parameter ν_curv (equivalent to ν_opt in Fig. 2), determined from the curvature of the L-curve, will be discussed.

By using equation (7) we calculated for one oxygen shell a Gaussian-shaped model RPDF (n_theor) with N = 1, r = 2.50 Å and σ² = 0.005 Ų. The terms in equation (8) [F(k,r), λ(k), 2δ(k) and φ(k,r)] are calculated by FEFF8.20 (see §2.3), where the amplitude reduction factor is set to S₀ = 1. The FEFF structural model consists of a curium hydrate cluster based on the EXAFS results of Skanthakumar et al. (2007). The matrix A_im is then calculated using (17) with an energy shift of ΔE₀ = 5 eV [see (18)] and the model spectrum χ is calculated using (10). For the calculation of the k-axis we used (18). We added to the unweighted spectrum χ a noise equal to C_noise·rnd, where the random number generator (rnd) gives values in the range −0.5 < rnd < 0.5 (Table 1). Fig. 3 shows for examples #2, #6 and #10 (Table 1) the theoretical RPDF and the noise-contaminated spectra.

Figure 3 Theoretical RPDFs (ntheor) and the corresponding spectra with added artificial noise are shown in black. The RPDFs (nνcurv) reconstructed using the Landweber iteration and the corresponding spectra are shown in red. Residuals are blue.

For examples #2–#10 the L-curves are shown in Fig. 4 with an inset showing example #10 on a more convenient scale. For the presented theoretical examples, n_theor is known, hence for each example the theoretical optimum number of iterations ν_theor can be determined as the minimum in ||n_ν − n_theor|| normalized to the maximum number of data points in r-space (m_max). Figs. 5(a) and 5(b) illustrate the determination of ν_theor and ν_curv for the examples #2–#10.

Figure 4 L-curves for the noise-contaminated spectra of examples #2–#10 (Table 1).

Figure 5 Definition of the stopping parameters νtheor and νcurv for the noise-contaminated spectra of examples #2–#10 (Table 1). (a) The difference between the theoretical and the iterated RPDF, normalized to the maximum number of data points in r-space. (b) The curvature [see equation (14)] of the L-curves (Fig. 4). For graphs (a) and (b) the symbols correspond to the numbers of the examples.

The inspection of the shape of the L-curves (Fig. 4) leads to the conclusion that for examples with a high noise level the norm of the residual ||An_ν − χ|| decreases and the norm of the solution ||n_ν|| increases faster than for the examples with a lower noise level; hence the higher the experimental error in the data the lower is the number of iterations needed for the solution. This trend is in accordance with the shift of the minimum in ||n_ν − n_theor|| and c(ν) to smaller ν_theor and ν_curv, in case the noise level increases [Figs. 5(a) and 5(b)]. These observations are in line with the semi-convergent behavior of the Landweber iteration. The observed good accordance between ν_curv and ν_theor (Table 1) shows the reliability of the L-curve concept to estimate the number of required iterations for noise-contaminated spectra.

For all L-curves, the position of the maximum curvatures determined by equation (14) shows a very similar norm ||n_ν|| (open circles connected by line in Fig. 4). This is a first indication that for all calculated RPDFs the coordination number will be very similar as discussed later.

For all examples, the determined n_υcurv and the theoretical n_theor are in good accordance, as can be seen for examples #2, #6 and #10 in Fig. 3. Also, the spectra calculated using n_υcurv are in good agreement with the noise-contaminated spectra within the noise level (Fig. 3).

The sought structural data (N, r, σ²) in Table 1 are derived from a comparison of the RPDFs with a Gaussian-shaped RPDF [see equation (7)]. For all noise levels the structural parameters (N, r) are in good agreement with those of the initial model (Table 1, #1). The deviation of σ² from its true value increases systematically and significantly with the noise level. For example #2, no significant deviation is observed, whereas, for example #10, σ² is 14% below the true value (Table 1). It is interesting to note that the precision in the determination of the structural parameters follows the order r > N > σ². This trend is also visible in the deviations between n_υcurv and n_theor (Fig. 3). The position of the maxima in n_υcurv does not change with the noise level, while the left- and right-hand sides of n_υcurv show deviations from n_theor with increasing noise (Fig. 3). It is important to realise that these deviations influence the shape of the RPDF and that therefore σ² is the parameter which is most affected by noise, i.e. experimental error, of the EXAFS spectrum.

Note that for the preceding comparison of models the energy shift for the noise-free spectrum was fixed at ΔE₀ = 5 eV. In the following we demonstrate, using example #6, the procedure for calculating ΔE₀ as described in §2.4. For this spectrum, ν_theor = 57 and ν_curv = 59 (Table 1, #6).

The starting values are ν_start = 10 and ΔE₀ = 0 eV. At ν_start = 10 the resulting ΔE₀ is 5.08 eV (Table 2), hence already near the value of a noise-free initial spectrum (Table 1, #6). The curvature shows for this first value of ΔE₀ a maximum at ν_curv = 64 (Table 2). By taking this value as the next ν_start, a new improved ΔE₀ is determined. After four such steps, ν_start is equal to ν_curv and ΔE₀ has reached a final value of 5.14 eV. When starting at ν_start = 100, the method leads to the same result (Table 2). During the approximation cycles the value of SD_ν decreases and the difference between the theoretical (ΔE₀ = 5 eV) and the calculated ΔE₀ increases slightly. This demonstrates that a part of the noise in the spectrum influences ΔE₀ and consequently influences the RPDF, which reflects the normal behavior of error propagation.

Table 2
Determination of ΔE₀ for the noisy spectrum of example #6, arbitrarily starting at ν_start = 10 and ν_start = 100

Step	νstart	νcurv	ΔE0 (eV)	SDν × 10^2 (Å^−3)
1	10	64	5.08	6.286
2	64	67	5.13	6.284
3	67	68	5.14	6.283
4	68	68	5.14	6.283
1	100	66	5.11	6.284
2	66	68	5.14	6.283
3	68	68	5.14	6.283

Table 3 shows the structural parameters for example #6, determined either by the Landweber iteration with calculated ΔE₀ (`iteration') or by shell fit. Both results are in good agreement. Furthermore, they are in good agreement with the results obtained by keeping ΔE₀ constant during the Landweber iteration (Table 1, #6).

Table 3
Structural parameters for example #6 determined by the Landweber iteration including the calculation procedure of ΔE₀ and the structural parameters determined by shell fitting

Standard deviations of the parameter are given in parentheses. In the case of the shell fit, the standard deviations were calculated by using the software EXAFSPAK (George & Pickering, 1995).

Example #6	N	r (Å)	σ^2 (Å^2)	ΔE0 (eV)
Iteration	0.98 (1)	2.500 (1)	0.0047 (1)	5.14
Shell fit	0.96 (3)	2.494 (2)	0.0046 (2)	4.9 (3)

4. Application of Landweber's iteration to curium(III) hydrate

The coordination environment of the hydrated Cm³⁺ ion was recently investigated by Cm L₃ EXAFS and high-energy X-ray scattering (Skanthakumar et al., 2007). It was found that the Cm³⁺ ion is surrounded by nine coordinating water molecules at two different Cm—O distances. The EXAFS spectrum was fitted by six oxygen atoms at 2.47 Å and three oxygen atoms at 2.63 Å, while maintaining the coordination numbers fixed to stabilize the shell fit (Skanthakumar et al., 2007).

In order to test our full Landweber approach (including the procedure for determining ΔE₀) for deriving the RPDF of a real world spectrum, the authors of Skanthakumar et al. (2007) kindly provided their experimental Cm³⁺ hydrate spectrum [Fig. 6(a), #1].

Figure 6 (#1) Experimental data of Cm(III) hydrate. (#2) Theoretical example. Data (black), result of the Landweber iteration (red), residual between data and result of Landweber iteration (blue). (a) EXAFS spectra, (b) corresponding Fourier transforms, (c) RPDFs. Decomposition of the RPDFs into two Gaussians and the sum of the two Gaussians (green).

For the calculation of the matrix A an r-interval of 1.3–4.0 Å with a step size of Δr = 0.002 Å is used. To determine ΔE₀, the method described in §2.4 is applied and yields ΔE₀ = 6.42 eV at ν_curv ≃ 900 (Table 4). The determined RPDF shows an asymmetry at higher r-distances [Fig. 6(c), #1]. The small features at 3.1–3.8 Å in the RPDF [Fig. 6(c), #1] may arise from reproduction of small multiple-scattering contributions by the Cm—O functions or other sources. The determined asymmetric RPDF is only properly reproduced when two Gaussians are considered in equation (7) [Fig. 6(c), #1]. The resulting structural parameters for the two Cm—O shells are within the standard deviations in full agreement with those gained by a shell fit of the Cm(III) hydrate spectrum (Table 4). The coordination numbers of the two Cm—O shells were assumed to be 6 and 3 and were fixed to stabilize the shell fit. Note that in the case of the Landweber iteration the RPDF is calculated without such constraints and assumptions contrary to the shell fit.

Table 4
Summary of the structural parameters for Cm(III) hydrate and a theoretical example gained by the Landweber iteration and the shell fit, the theoretical (ν_theor) and L-curve (ν_curv) defined stopping parameters in the case of the Landweber iteration

Example	Cnoise × 10^4	N	r (Å)	σ^2 × 10^3 (Å^2)	νtheor	νcurv	ΔE0 (eV)
Cm(III) hydrate	–	6.0 (6)	2.469 (4)	6.0 (3)	–	900	6.42
(Landweber iteration)		3.1 (7)	2.61 (2)	10 (2)
Cm(III) hydrate	–	6†	2.464 (1)	5.8 (2)	–	–	6.51 (4)
(Shell fit)		3†	2.612 (4)	9.5 (8)
Theoretical example	5	6.0 (2)	2.464 (1)	5.6 (1)	1709	1922	6.51†
(Landweber iteration)		2.9 (2)	2.614 (6)	8.7 (5)

†Fixed structural parameter. Standard deviations of the parameter are given in parentheses. In the case of the shell fit, the standard deviations were calculated by using the software EXAFSPAK (George & Pickering, 1995).

In the following we tested the ability of the Landweber iteration in resolving these two overlapping Cm—O shells by constructing a theoretical model. Both the model RPDF and χ were constructed using the structural parameters gained from the shell fit of the Cm(III) hydrate spectrum (Table 4) and the procedure discussed in §3. To the unweighted spectrum χ, we added noise C_noise (Table 4) similar to the noise amplitude of the experimental spectrum of Cm(III) hydrate [compare #1 and #2 in Fig. 6(a)]. The optimum number of iterations, ν_curv = 1922, is near the theoretical optimum number of iterations of ν_theor = 1709 (Table 4). For Cm(III) hydrate, ν_curv is much higher than for the theoretical examples discussed in §3, which were analyzed in the smaller r-interval of 2–3 Å. The higher ν_curv can be explained by the decrease of the convergence parameter α owing to the increasing norm ||A^*A|| (see §2.1) in the case of the larger investigated r-interval of 1.3–4.0 Å.

The resulting RPDF [Fig. 6(c), #2] reproduces the noise-contaminated EXAFS spectrum [Fig. 6(a), #2] within the limits of the artificially added noise. Owing to the spatially distinct two oxygen shells, the RPDF shows an asymmetry at higher r-distances which is not visible in the corresponding Fourier transform [compare Figs. 6(b) and 6(c), #2], similar to the experimental data [compare Figs. 6(b) and 6(c), #1]. For the calculation of the structural parameters, the RPDF is reconstructed with two Gaussians according to equation (7). The resulting structural parameters (Table 4, theoretical example) are in good accordance with those taken for the construction of the theoretical RPDF (Table 4, shell fit), hence there is strong confidence in the correct calculation of the RPDF by the Landweber iteration in the case of the experimental Cm(III) hydrate spectrum. The structural parameter for Cm(III) hydrate, gained by the Landweber iteration and the shell fit (Table 4), agrees favorably with the results published by Skanthakumar et al. (2007) within the bounds of the errors in the determination of EXAFS structural parameters.

5. Conclusions

The analysis of model spectra with artificial noise, as well as of an experimental EXAFS spectrum of Cm(III) hydrate, demonstrates that the Landweber iteration approach is well suited to solving the EXAFS integral equation, even if the spectra contain experimental error. The robustness in relation to the error in the data represents the main advantage of this method.

All variables used in the iteration approach [equation (12)] have a well defined physical meaning. The determined RPDF n(r) is the density of coordination numbers, the key transformation with the matrix A^* is a windowed and phase-corrected Fourier transform of the spectrum, and the acceleration term (χ − An_ν) is the difference spectrum between the original spectrum χ and the last iteration spectrum An_ν. Convergence and stopping parameters are uniquely defined by the matrix A and the L-curve concept.

In the presented approach, the threshold energy shift ΔE₀ is defined by the iteration procedure itself, without recourse to fit-calculations or models.

The forthcoming investigations focus on the expansion of the stable Landweber iteration with the aim of considering more than one type of backscattering atoms and multiple-scattering effects for the inversion of the EXAFS equation.

Footnotes

¹This boundary condition is irrelevant for RPDFs, which are derived from EXAFS measurements. The distance of neighboring atoms measured by EXAFS spectra is smaller than the asymptotic value of r, for which g(r) → 1 is valid.

²A^* denotes the conjugate transpose. Here and in the following the norm of a matrix is the Frobenius norm.

³For the discretization of the EXAFS equation (6) we use the following. The indices i,j,… run from 1 to the maximum number i_max, j_max etc. of points in k-space. The indices m,n,… run from 1 to the maximum number m_max, n_max etc. of points in r-space.