Improving the spatial and statistical accuracy in X-ray Raman scattering based direct tomography

An algorithm to simultaneously increase the spatial and statistical accuracy of X-ray Raman scattering (XRS) based tomographic images is presented. Tomography that utilizes XRS spectroscopy signals as a contrast for the images is a new and promising tool for investigating local atomic structure and chemistry in heterogeneous samples. The algorithm enables the spatial resolution to be increased based on a deconvolution of the optical response function of the spectrometer and, most importantly, it allows for the combination of data collected from multiple analyzers and thus enhances the statistical accuracy of the measured images.

1. Introduction

X-ray Raman scattering (XRS) spectroscopy or non-resonant inelastic X-ray scattering (NRIXS) from core levels is a rapidly emerging technique that utilizes hard X-rays to probe soft X-ray absorption edges (Schülke, 2007; Sahle et al., 2015). Recently, Huotari et al. (2011) presented a fundamentally new tomographic imaging technique that can use the inelastic X-ray scattering signal as contrast for the three-dimensional images. This imaging technique was coined direct tomography (DT) by the authors and is based on the fact that the energy analysis of the scattered hard X-rays in NRIXS is performed with spherically bent analyzer crystals in the Johann geometry, which yield a point-to-point focus. Using these analyzer crystals in combination with a pixelated area detector, one can construct three-dimensional images from simple translation scans of the sample. Complicated algorithms to reconstruct images from different projections are not necessary in DT.

The DT technique bears tremendous potential for the study of heterogeneous sample systems such as catalysts, batteries and fuel cells (Inkinen et al., 2015; Sahle et al., 2016). Most notably, the hard X-rays used in NRIXS (usually in the 5–15 keV range) render DT compatible with complicated in situ cells and reaction chambers (Sahle et al., 2013, 2016).

Current drawbacks for DT include the limited spatial resolution when compared with standard absorption- or phase-contrast computed tomography (CT), the limited penetration depth for the normally used incident photon energies, and the long acquisition times that are due to the small overall cross section of NRIXS. The dominating factors for the limit in spatial resolution are the quality of the bent analyzer crystals and their focusing properties, the size of the incident X-ray beam, and the detector pixel size. The former two are defined by the instrument used, typical beam sizes are of the order of 10 µm, whereas the detector pixel size is of the order of 50–100 µm.

In this contribution, we introduce an algorithm to simultaneously enhance the spatial resolution and the voxel-wise signal-to-noise-ratio for DT imaging. With this approach, we extend the capabilities and range of applications for DT into regimes that notably improve its appeal for the study of the named applications by approaching micrometer-resolution in millimeter-sized systems. Some theoretical background of XRS is given in §2, the NRIXS-based imaging (or DT) technique and the algorithm are explained in detail in §3 using an illustrative example. We conclude briefly and give an outlook in §4.

2. Non-resonant inelastic X-ray scattering

In non-resonant inelastic X-ray scattering (NRIXS), hard X-ray photons are scattered inelastically from the sample system. In the scattering process, the energy = − and momentum q = k₁ − k₂ are transferred to the system. Here, () and k₁ (k₂) are the energy and momentum of the incident (scattered) X-ray photons, respectively (see inset of Fig. 1). If the energy transfer ω is of the order of an electronic binding energy of the system, one can measure a core edge spectrum similar to an X-ray absorption spectrum. At low transferred momentum q, these core edge spectra are directly proportional to soft X-ray absorption spectra (Mizuno & Ohmura, 1967). At higher exchanged momenta, channels allowing excitations of higher than dipolar order open up and NRIXS becomes sensitive to multipolar transitions (Gordon et al., 2007).

Figure 1 Schematic NRIXS spectrum that exemplifies the different electronic excitations that can be measured using this inelastic X-ray scattering technique. The inset shows the scattering geometry used in a typical NRIXS experiment.

The measured quantity in an NRIXS experiment is the double differential scattering cross section,

with the Thomson scattering cross section

where r_e is the classical electron radius and () is the polarization vector of the incident (scattered) X-rays. is the so-called dynamic structure factor,

and holds all information about the sample system obtainable by NRIXS; is the initial state (typically a shallow or low-energy core level) and the summation is over all reachable final states . E_i and E_f are the according energies of the initial and final state. Expansion of the transition operator in a Taylor series (or, complementary, in terms of spherical harmonics) leads to the mentioned q-dependence, i.e. sensitivity to dipolar transitions for small and higher-order transitions (2^k, k = 2, 3, 4…) for larger . Fig. 1 shows a schematic spectrum that exemplifies the different excitations that can be measured using NRIXS and, consequently, can serve as a contrast mechanism in DT.

3. NIXS-based imaging (direct tomography)

The direct tomography technique was introduced by Huotari et al. (2011). It is based on the point-to-point focusing properties of bent analyzer crystals that are frequently used in NRIXS spectrometers for the energy analysis of the inelastically scattered radiation from the sample in combination with pixelated area detectors. Fig. 2 shows a schematic sketch of the DT geometry. The incident X-ray beam impinges on a heterogeneous sample, and the radiation scattered off the sample at different positions along the X-ray beam through the sample is focused onto different positions on a pixelated area detector.

Figure 2 Schematic of the direct tomography imaging principle.

In principle, any one of the excitations shown in Fig. 1 can be used as a contrast for the direct imaging technique simply by tuning the energy loss to the appropriate value. In the most elaborate approach, an entire NRIXS spectrum can be recorded for each voxel.

For a given sample position denoted by y,z and for a given energy loss ω, each analyzer denoted by a provides, if we neglect multiple scattering, a one-dimensional image of the beam path. These images provide a view, along the x-axis, of the superposed constituent densities of scatterers , where c is the index of the respective constituent. We show in Fig. 3(a) an example of a typical detector image from the multiple-analyzer spectrometer at ID20 where all images for the 72 different analyzers a and for a given translational position y,z of the sample are recognizable as 72 elongated spots.

Figure 3 Images of the beryllium gasket at a given y,z sample position (a) and point spread function for analyzer a = 17, measured with a 25 µm Kapton foil (b). In this paper we use the ensemble of the point spread functions to deconvolve the signal in the x direction.

The measured images can be represented by

where is the image produced by analyzer a, for a single emitting point located along the beam at coordinate x (the point spread function, in other words), the factor takes into account the emission physics which depends on the values of the energy loss and the exchanged momentum q. Finally the transmission factors t take into account the absorption from the source to the sample and from the sample to the detector. These transmission factors are given by the exponential of the absorption integral along the photon path. In the above formula we neglect the angular aperture of the analyzer crystals as seen from the sample, so that in this approximation the value of q_a and the absorption integrals are well defined.

We are able to retrieve the density of scatterers if, for a given sample position y,z, we are able to extract the profile along x of from the measure of the images . These measures are the images of the beam portion which crosses the sample, as imaged by the curved crystal analyzers. One of these images can be considered as the result of a summation of all images [integral in equation (4)] generated by a set of point sources located along the beam path.

An important ingredient enabling the inversion of equation (4) is thus the measurement of such point-generated images. We realise these measurements by placing a thin foil, thick, of a constituent in the beam, and translating it along the x axis (imaging direction), while measuring the reference images , given by

These images must be acquired at the same energy loss ω as the sample images to avoid slight translations of the incident X-ray beam due to mechanical errors when its energy is changed. It is important that the reference images are well aligned with the sample images. The algorithm for inverting equations (4) and (5) must be adapted to the specific case under study, as a sample may contain different constituents, each one with a particular spatial distribution, spectral properties and absorption coefficients.

We illustrate in this paper the solution of this inversion problem for a simple case with low absorbance made up of a single constituent. As an outlook, we provide a further development of the method, which aims at the faithful imaging of low-concentration impurity profiles exploiting two (or several) measurements around the onset of an electronic excitation of one of the constituents of a heterogeneous sample to subtract out, from the formula, the absorption (t_x,y,z ^a, t_source ^x,y,z) of the main constituent of the sample.

3.1. Imaging of a single constituent in a low-absorption sample

The first method that we present assumes that either the sample is made up of only one constituent c or that one can select an energy loss ω, which enhances the contribution of the component of interest c to equation (4) to such a degree that one can neglect all other constituents. We include in our method the possibility that the reference images are measured for a foil constituent different from c. It is important to introduce this flexibility in the method since it may be difficult to always find appropriate foils for which c holds.

In the current demonstration, we apply the method to the imaging of a beryllium gasket, frequently used in high-pressure X-ray experiments to contain a micrometer-sized sample between two diamond anvils, using a thin Kapton foil as a reference. As the q_a dependence of the factors may vary from one component to another, we accommodate these variations by introducing the adjustment factors f_a in the following equation, which expresses the recorded images as a function of and f_a,

where the set x_i are the foil positions at which we have probed the point response functions, while the purpose of the factors f_a is to absorb our disregard of the complexity of the physics underlying the factor e_c.

A practical method to extract consists of iterating the two following equations,

and

where β is an optional regularization factor, which multiplies the total 3D variation of the solution ρ. It is used to stabilize the inverse problem when the sampling points x_i are too dense. In such a case, the measured reference images may not be linearly independent and the inverse problem may become singular, if no regularization is applied. The solution is obtained by first solving equation (7), with all of the f_a set to unity, and then alternating the solutions of equations (8) and (7).

To solve this inverse problem, we used the Fast Iterative Shrinkage-Thresholding Algorithm (Beck & Teboulle, 2009); total variation regularization was included as described by Mirone et al. (2014) and Van der Walt et al. (2014). For the construction of the DT image shown in Fig. 4, the algorithm converged within less than 100 iterations in approximately 15 min on a single CPU. The regularization parameter β was chosen manually such that the resulting image provides the best compromise between high variance and high bias.

Figure 4 Left: phase-contrast computed tomography reconstruction of the specimen (red). Middle: isosurface plot of the DT volume, where the colormap applied to the isosurface represents spatial distances between the CT and DT isosurfaces. Right: zoomed-in detail, near the central hole of the specimen, of the DT reconstruction with different numbers of analyzers crystals. In all panels, CT images are red and DT − CT difference maps are in blue.

All DT images were taken at the inelastic scattering beamline ID20 of the ESRF, France (Huotari et al., 2017). Photons from four U26 undulators were monochromated by a Si(111) high-heat-load double-crystal monochromator and a consecutive Si(311) channel-cut post-monochromator to an incident bandwidth of approximately 0.4 eV at a primary photon energy of 9.7 keV. The X-ray beam was focused to a spot size of 10 µm × 20 µm (V × H) at the sample position using a Kirkpatrick–Baez mirror system. We utilized the multi-analyzer spectrometer with 72 spherically bent Si(660) analyzer crystals arranged in six separate analyzer modules each housing 12 analyzer crystals and a Maxipix detector (pixel size 55 µm) (Ponchut et al., 2011). We show images of all six detectors tiled in a 3×2 manner for a given sample position y,z and = 0 and = 9.7 keV in Fig. 3(a). Each detector records the images produced by 12 focusing analyzers. Therefore, there are 72 line-shaped spots in the figure, each line being the image of the scattering from the same sample region (y,z), intersected by the X-ray beam, but from different scattering angles, i.e. momentum transfers q_a. The overall time spent recording the DT image presented in Fig. 4 was 160 min, where most of the time was spent for movements of the motor, and the integrated counting time was only 9 min, leaving space for vast improvements. For comparison we measured the same sample at beamline ID17 of the ESRF using phase-contrast CT where a collimated X-ray beam produced by a 21-pole wiggler and monochromated by a double Si(111) crystal system was used. The X-ray energy was set to 40 keV. Images were acquired by using a PCO Edge 5.5 camera with a pixel size of 3.05 µm × 3.05 µm. The propagation distance was 2.5 m. This setup is designed for an optimal phase propagation enhancement for the medical studies which are routinely performed on this beamline and the overall accumulation time was approximately 10 min. The comparison volume was reconstructed using the PyHST2 code (Mirone et al., 2014) using phase retrieval (Paganin et al., 2002).

The images are extracted by selecting a region of interest for each analyzer a, and acquiring images at a minimum of one energy loss ω and a set of sample positions (y,z). The same procedure is employed when measuring for the extraction of the point spread function . In the current case, we used a Kapton foil of 25 µm thickness and measured reference images at different x in steps of 40 µm. The extracted response function for analyzer a₁₇ at x = 0 and = 0 is shown exemplarily in Fig. 3(b).

The resulting DT image of the Be gasket is shown in the center of Fig. 4 (blue isosurface plot) together with a phase-contrast CT image recorded at beamline ID17 of the ESRF (red isosurface plot on the left of the figure). We obtained the DT image using signals from 69 analyzer crystals and a scanning step size of 20 µm in the horizontal (y) and 20 µm in the vertical (z) direction.

Instead of plotting the isosurface of the DT volume in a single color, we applied a colormap to the image that represents the spatial accuracy of the constructed DT image. It expresses, in micrometers, the distance between the isosurface of the DT image and the isosurface of the CT image, which we assume to be the ultimate benchmark in terms of spatial accuracy. These distances between the isosurface of the CT image and the isosurface of the DT image are calculated as follows: for each point of the DT isosurface we find the nearest point in the CT isosurface and calculate the distance between these two points. The isosurfaces are then aligned to minimize the sum of the squares of the distances. We plot these distances for the optimally aligned grids as color codes in the center part of Fig. 4.

In the right-hand part of Fig. 4, a comparison of DT images constructed using signals from different numbers of analyzer crystals with the phase-contrast CT image is shown for detailed views of the central hole in the gasket. The quality of the images clearly improves with an increasing number of used analyzer crystals as a proof of the increased information content when data from several analyzer crystals are combined. This improvement is quantified by the quadratic mean of the distances, which is 45, 19 and 16 µm for 1, 19 and 69 analyzer crystals, respectively.

3.2. Outlook for quantitative imaging of minor constituents

We consider two measurements taken at two slightly different energy losses and , below and above an absorption edge of a minor (approximately 1–10 at%) constituent . We neglect the energy dependence, between and , of the physical properties, except for the response of the minor constituent. We neglect further the effect of the constituent in question on the absorption coefficients t_x,y,z ^a, t_source ^x,y,z and their energy dependence.

Given these assumption, we can rewrite equation (4), at the two energies, in a form which allows for the quantitative extraction, up to a global factor, of the concentration of impurities in the probed volume,

and

From these two equations, when the minor constituent is embedded in a homogeneous matrix we can extract the following distribution functions: f(x,y,z), which is proportional to the concentration of the component in question, and , which is a gauge of , the latter being proportional to the off-resonant electron density. They are defined as

where u_a and are unknown factors, which depend only on the detector index a. The rationale, which justifies the factorization of u_a and , is that we assume a low concentration of constituent inside a homogeneous matrix. Given the above assumptions and definitions, we propose to extract the concentration profile f(x,y,z) as the solution of the following system of equations,

and

which can be solved by first solving equation (14) for the set of , given an arbitrary choice of u_a, then solving equation (15) by alternatively optimizing f(x,y,z) and .

4. Conclusions and outlook

We presented an algorithm to increase the spatial resolution of direct tomography images by means of deconvolution of the instrument's optical response function and, at the same time, increasing the statistical accuracy of the images by combining signals from different analyzer crystals of a multi-analyzer crystal spectrometer. We further outlined an approach to use information from direct tomography to extract quantitative information about the concentration of a minor constituent inside a homogeneous matrix. These developments of the direct tomography technique will facilitate greatly the application of this technique by reducing the overall measurement time and by making accessible spatial inhomogeneities of a few tens of micrometers. Direct tomography is a promising tool for studying the local structure and chemical bonding of low-Z-element-containing samples contained in complicated sample environments, such as in operando batteries and fuel cells. For example, we envision this technique to be capable of monitoring the spatial distribution of Li in a specific oxidation state across the anode–electrolyte–cathode assembly of an in operando Li-ion battery.

As an outlook, a more faithful modeling can be performed given a knowledge of the absorption coefficients at the incoming and outcoming energies. In this case the absorption paths t_x,y,z ^a and t_source ^x,y,z can be estimated. Tomography computed volumes at relevant absorption energies could be used.