Calculation of the Kramers-Kronig transform of X-ray spectra by a piecewise Laurent polynomial method

Abstract: An algorithm is presented for the calculation of the KramersKronig transform of a spectrum via a piecewise Laurent polynomial method. This algorithm is demonstrated to be highly accurate, while also being computationally efficient. The algorithm places no requirements on data point spacing and is capable of integrating across the full spectrum (i.e. from zero to infinity). Further, we present a computer application designed to aid in calculating the Kramers-Kronig transform on near-edge experimental X-ray absorption spectra (extended with atomic scattering factor data) in order to produce the dispersive part of the X-ray refractive index, including near-edge features.


Introduction
Knowledge of the complex index of refraction of materials is required for the analysis of a wide range of X-ray experiments, especially those making use of near-edge regions of the spectrum where interactions between atoms become significant, causing strong variations in both the real and imaginary parts of the scattering factors. Resonant soft X-ray scattering (RSoXS) [1,2] and reflectivity (RSoXR) [3,4] as well as multiple-wavelength anomalous dispersion (MAD) [5][6][7][8] and diffraction anomalous fine structure (DAFS) [8][9][10] routinely require calculations based on both the real, dispersive and the imaginary, absorptive parts of the refractive indices of the samples component materials. Indeed, a recently published spectroscopic database [11] of conjugated polymer materials presented the data in terms of scattering factors in order to provide easy accessibility to both the dispersive and absorptive parts of the index of refraction.
While the imaginary, absorptive part of the refractive index can be easily determined via simple experiments, the real, dispersive part is far more difficult to measure. Fortunately, the real and imaginary parts of the refractive index are related to each other by the Kramers-Kronig relations [12] and so the full refractive index can be calculated from measurements of only the imaginary part. Let us describe the complex index of refraction following the convention used by Henke et al. [13] in defining their atomic scattering factors, f 1 and f 2 , which are published for the first 92 elements: where n represents the refractive index (with dispersive, δ , and absorptive, β , components) of a material consisting of a set of q elements, each with a number density of N q , and r 0 is the classical electron radius, λ is the wavelength and E is the photon energy of the incident light. Then the Kramers-Kronig transform can be written in terms of the Henke scattering factors as: where E is the photon energy at which we wish to calculate the real part of the atomic scattering factor, Z * is the relativistic correction, and P denotes the Cauchy principal value. Equation (1) is problematic in two ways. Firstly, the integral spans zero to infinity, meaning that in order to calculate f 1 at a single photon energy, the entire f 2 spectrum is required. The second issue is that the the integral contains an undefined point at x = E and so must be treated as a Cauchy principal value.
The Cauchy principal value is central to many calculations that occur in a wide range of physics and engineering and a number of quadrature strategies have been applied to the problem [14][15][16], such as the FFT [17,18], Euler-MacLaurin [14,[19][20][21], Gaussian [22], Chebyshev [23] and variable transformation [24] (e.g. TANH [25] and Iri-Moriguti-Takasawa [26]) methods. However, these all have significant disadvantages: the FFT and Euler-Maclaurin methods, for example, require evenly separated data-points, which is especially cumbersome over wide ranges, while the others are very CPU intensive and in practice require compromising on accuracy to offset the computational expense. Cross et al. [19,27] demonstrated a creative way of using the linearity of the Kramers-Kronig transform to only require calculating the transform of the difference between two spectra, which serves to limit the integral to a manageable range. They use the theory of Cromer and Liberman [28] to generate the scattering factor spectra ( f 1 (E) and f 2 (E)) for a collection of isolated atoms corresponding to the stoichiometry of the material of interest such that the difference between the measured near-edge imaginary scattering factor spectrum and the calculated one is only significant in a narrow near-edge region. The Kramers-Kronig transform of this difference spectrum is then added to the calculated real scattering factor spectrum to obtain the desired result.
Henke et al. [13] demonstrated a simple method involving manual piecewise polynomial representation and symbolic evaluation of the integrals that is accurate, but labour-intensive. The method described in this work is an extension of that used by Henke et al. to allow computer automation and so can provide computationally cheap calculations whose accuracy is limited only by the representation of the absorption spectrum by a piecewise set of Laurent polynomials.
If the input spectrum f 2 for Eq. (1) is expressed (or approximated) as a continuous piecewise Laurent polynomial: Now we can integrate each interval separately and add the areas together to get the total integral. So taking the first integral in Eq. (2), for some polynomial function p j (x) with terms of degree M ∈ Z to N ∈ Z that defines F (x) in the jth interval: x j+1 Similarly, the integral having an x − E denominator can be evaluated via: x j+1 While Eq. (3) can be easily evaluated for all j, Eq. (4) contains a problematic ln term that cannot be evaluated when either x j or x j+1 is equal to E. (Note that cases in which x j < E < x j+1 pose no issue.) In such cases, we can avoid evaluating Eq. (4) at these points by integrating across the two intervals adjacent to x = E at once. Since F (x) is defined piecewise, the adjacent intervals are described by separate polynomials -let's call them p j−1 (x) = ∑ N n=M a j−1,n x n , and p j (x) = ∑ N n=M a j,n x n and avoid the case where x j = E by integrating across the double interval x j−1 → x j+1 . F (x) must be continuous in order to satisfy the Kramers-Kronig relations and so the polynomials of the adjacent intervals must intersect at the shared point (E, F (E)) between them. Therefore F (E) = p j−1 (E) = p j (E) and so ∑ N n=M a j−1,n E n = ∑ N n=M a j,n E n and the problematic terms cancel to give: x j+1 The double interval result above in Eq. (5) is valid for all j except the endpoints (x j = x 1 , x i ) and so the x − E integral can be evaluated at any energy, E, between the integral limits (x 1 < E < x i ) by using a combination of Eqs. (4) and (5). Further, instead of fully switching to Eq. (5) to integrate across the singularity, we can simply take the result for the problematic logarithm term given in Eq. (5) and incorporate it into Eq. (4) using the Kronecker delta function: where δ is the Kronecker delta function, defined as: and is easily implemented in computer code using logic functions. Equation (6) contains three instances of the Kronecker delta function which provide three different versions of the problematic logarithm term. In the case where E is equal to neither x j nor x j+1 , the deltas all evaluate to zero and the standard form seen in Eq. (4) is provided. In the case where E = x j+1 , the term is set to zero and 1 is added to the numerator of the logarithm in order to avoid evaluating the logarithm of zero. Finally, in the case where E = x j , the x j is replaced with x j−1 in order to provide the two-interval form seen in Eq. (5).
Finally, assembling Eqs. (2), (3) and (6) gives: It is a common strategy to restrict the integration limits of Eq. (1) to non-zero and finite values via the assumption that f 2 (x) is not significant outside of the chosen integration range . However, integration from zero to infinity is feasible with a piecewise polynomial approach so long as the input data is appropriately described. Examining Eq. (4) in the limit where x 1 → 0, it is clear that the polynomial describing the imaginary scattering factor in the photon energy interval adjacent to 0 eV must have its coefficients for all even terms of degree -2 or lower equal to zero (i.e. a n = 0 for all even n ≤ −2). However, F (0) must be finite in order to satisfy the Kramers-Kronig relations, which provides the stricter requirement that the coefficients for all terms of degree less than zero be equal to zero (i.e. a n = 0 for all n < 0). Similarly, in the limit where x i → ∞, the polynomial describing the imaginary scattering factor in the photon energy interval extending to infinity must have its coefficients for all terms of degree -1 or greater equal to zero (i.e. a n = 0 for all n ≥ −1) in order for the integrals to converge.

Methods
Following the method of Henke et al. [13], values for the complex atomic scattering factors between 10 eV and 30 keV were downloaded from the website of the center for X-ray optics [29] and supplemented with higher photon energy data (up to 500 keV) from Biggs and Lighthill [30]. The Henke data is provided as point values and was converted to a piecewise polynomial function by applying simple linear functions between successive data points. The Biggs and Lighthill data is provided in the form of a four-term piecewise polynomial approximation of photoelectric cross-sections in reciprocal powers of the photon energy. In translating the Biggs and Lighthill data into the atomic scattering factor format, the polynomials increase in degree by one such that the terms are of degree 0 to degree -3. While Biggs and Lighthill provide data extending to infinity, the polynomials given in the range extending from 500 keV to infinity have non-zero coefficients for terms of degree 0 and -1 and so do not allow the integrals to converge, as discussed in the previous section. Combining the two data sets produces a set of piecewise polynomial functions representing the imaginary part of the atomic scattering factors for the first 92 elements of the form: where f 2 is the complex part of the atomic scattering factor, x is the photon energy and a n is the coefficient of the term of degree n. Again following Henke et. al, the relativistic correction, Z * , was calculated via: 2.37 (9) where Z q and n q represent the atomic number and stoichiometric quantity of the qth element in the material. Applying Eq. (7) to a piecewise polynomial function restricted to terms of degree 1 to -3, as shown in Eq. (8), produces: The combination of coefficients and terms in Eq. (10) are summarised in table 1.
Calculating the full scattering factors for a material involves the following steps: 1. Look up the imaginary scattering factors in the database for each element included in the material and then add them together according to the material stoichiometry.
2. If a section of spectrum is provided, convert it to scattering factor format and then to piecewise polynomial form, then merge with the database data. 6. Compile the set of real an imaginary scattering factors along with the set of photon energies at which they were calculated and present them to the user. These methods have been implemented as a Python-based library and application named "KKcalc". It is cross-platform, includes a graphical user interface (as shown in Fig. 1) and can be downloaded from https://bitbucket.org/benajamin/kkcalc.

Results and discussion
In order to prove the accuracy of the algorithm, the model extinction function of Ohta and Ishida [14] was applied and the output compared to the analytical result. Figure 2 plots one peak of the Ohta and Ishida antisymmetric Lorentzian pair absorption function k (ν) and its analytically derived Kramers-Kronig transform, ∆n (ν). k (ν) was approximated with a piecewise polynomial function, k PP (ν) with terms ranging from 4 th to 0 th order, via values of k (ν), k ′ (ν) and k ′′ (ν) evaluated at a set of 101 points distributed logarithmically (in terms of distance from the peak frequency) between 5 cm −1 and 10 10 cm −1 . The error in this approximation (magnified by 10 4 ) is shown in Fig. 2. k PP (ν) was then transformed by the Kramers-Kronig algorithm presented in this work to produce ∆n PP (ν) (evaluated at 1001 evenly spaced points from 1000 cm −1 to 3000 cm −1 ) and its difference from the analytical result, ∆n (ν), is also plotted (again magnified by 10 4 ) in Fig. 2. It is clear from Fig. 2 that the piecewise polynomial Kramers-Kronig algorithm is very precise and calculating the root mean square (rms) values of the difference from the analytical result across the interval 1000 cm −1 to 3000 cm −1 provides a value of 2.4 × 10 −7 . This rms error is nearly two orders of magnitude better than the best method tested by Ohta and Ishida [14], despite using far fewer integration points. Part of the reason for such a strong improvement in the accuracy of the algorithm is due to its ability to easily integrate over a very wide spectral range, eliminating artefacts that commonly occur near the edge of the integration range when the effective input function drops suddenly, and significantly, to zero. ∆n(ν) (k(ν)−k PP (ν)) ×10 4 (∆n(ν)−∆n PP (ν)) ×10 4 Fig. 2. The model extinction coefficient spectrum, k (ν), and the analytic solution to its Kramers-Kronig transform, ∆n (ν)), plotted as described by Ohta and Ishida [14]. Their differences from a piecewise polynomial approximation of the input function, k PP (ν), and its transform using the algorithm presented here, ∆n PP (ν), are also plotted to demonstrate a close fit to the analytical functions. . While difference between the extended (green) and limited (yellow) integrations is negligible in the case of light elements such as carbon (Fig. 3), it becomes just noticeable with iron ( Fig. 4) and significant for heavy elements like lead (Fig. 5), where there is significant absorption of X-rays above 30 keV.
The methods described in this work allow easy merging of highly detailed, material-specific data into the coarse, long-range elemental data from a database without sacrificing accuracy or committing significant computing power. Recently, Yan et al. [31] demonstrated an experimental method for measuring the complex index of refraction of polymer films at the carbon K-edge, along with a doubly-subtracted Kramers-Kronig calculation (utilising a fast Fourier transform algorithm) for calculating the real part from absorption measurements (i.e. the imaginary part). Figure 6 reproduces Fig. 3 from Yan et al. [31] overlayed with a plot of a recalculation of δ (E) using Eq. (10) (orange line in left panel) and integrating over the energy range 10 eV to 500 keV, as discussed in the previous section. Conversion between β , δ and the corresponding scattering factors was performed using a density of 1.20 gcm −3 , as determined experimentally by Yan et al. While the two δ (E) calculations match well at the "subtraction energies" of 270 eV and 320 eV, the doubly-subtracted δ values show some distortion that exaggerate the significance of the local structure in the β spectrum due to the limited integration range used by Yan et al. Note that carbon K-edge spectroscopic data, as shown in Fig. 6, is prone to distortions due largely to the significant carbon contamination typically found on X-ray optics, whose effects are diffi-   cult to fully remove [32]. Such distortions may explain the systematic differences between the NEXAFS (blue line in Fig. 6 (right)) and reflectivity (green points in Fig. 6) β measurements and hence also the differences between the reflectivity δ measurements and the two δ functions calculated from the NEXAFS β measurement. The dashed black lines in Fig. 6 show spectra derived from the data published by Henke et al. [13], which describe isolated atoms and thus display none of the near-edge resonance features that are observed in the experimental data. The success of the presented Kramers-Kronig algorithm stems largely from the aptness of piecewise Laurent polynomial representation for X-ray spectra, involving occasional highresolution features separated by wide stretches of gradual decay. Because the spectra can be represented compactly yet accurately, the Kramers-Kronig transform algorithm presented here can process broad sections of spectrum at once, while maintaining accuracy and precision. Also, the presented method allows for complete freedom in the choice of data points in both the input f 2 (E) spectrum (i.e. near-edge data can be measured in any way) and the calculated f 1 (E) spectrum. Further, the lack of any approximation around the integration singularity and the ability to represent the entire spectrum (0 ≤ E ≤ ∞) makes it possible to calculate a guaranteed Kramers-Kronig conjugate pair.

Conclusion
An automated method for calculating the real, dispersive part of the scattering factors (or refractive index) from the imaginary, absorptive part is presented along with some sample calculations using published X-ray data [13,29,30]. The presented algorithm utilises a piecewise Laurent polynomial representation of the input data and is accurate, computationally inexpensive (plus is not memory-bound) and straightforward to implement. It requires no especial subdivision of data points and can operate over the full energy range from zero to infinity. A Python-based library and application named "KKcalc" that implements this algorithm is available from https://bitbucket.org/benajamin/kkcalc.
The piecewise Laurent polynomial algorithm would also be easily adaptable to other problems related to the Cauchy principle value and the Kramers-Kronig relations. An interesting extension of this work would be to express the algorithm in spline representation, in order to leverage the spline approximation and evaluation functions commonly available in modern programming libraries. This work further suggests that a piecewise continuous representation of spectral data in published databases, such as Henke et al. [13], would be more accurate and can be easier to manipulate than the simple tabulated format currently used.