Modeling of beam hardening effects in a dual-phase X-ray grating interferometer for quantitative dark-field imaging

X-ray grating interferometry (XGI) can provide access to unresolved sub-pixel information by utilizing the so-called dark-field or visibility reduction contrast. A recently developed variant of conventional XGI named dual-phase grating interferometer, based only on phase-shifting structures, has allowed for straightforward micro-structural investigations over multiple length scales with conventional X-ray sources. Nonetheless, the theoretical framework of the image formation for the dark-field signal has not been fully developed yet, thus hindering the quantification of unresolved micro-structures. In this work, we expand the current theoretical formulation of dual-phase grating interferometers taking into account polychromatic sources and beam hardening effects. We propose a model that considers the contribution of beam hardening to the visibility reduction and accounts for it. Finally, the method is applied to previously acquired and new experimental data showing that discrimination between actual micro-structures and beam hardening effects can be achieved. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
X-ray grating interferometry (XGI) [1,2] is an imaging technique that simultaneously provides three complementary contrasts; absorption, differential phase, and visibility reduction also known as dark-field contrast. There has been a great interest, particularly in the dark-field contrast [3], as it provides information about underlying scattering structures, which are much smaller than the pixel size of the utilized detector [4][5][6]. The utility of the dark-field contrast has been demonstrated by applications both in the industrial and medical fields. For example, it can be used to obtain the orientation of unresolved carbon fibers in composite materials [7][8][9] which determine their mechanical properties. Furthermore, the observations of the time-evolving hydration process of mortars were enabled by the high sensitivity of the dark-field contrast to unresolved pores [10,11]. Recently, dark-field was also used to detect unresolved structural defects in composite materials [12]. In the medical research community, it was demonstrated that the dark-field signal could be used to differentiate between various types of microcalcifications which are associated with different levels of breast cancer malignancy [13,14]. Similarly, distinct types of kidney stones can also be discriminated [15]. Another highly promising medical application is lung imaging, by exploiting the strong scattering taking place in the microstructure of the lung, high contrast images can be obtained, opening the potential for improved diagnosis [16,17]. All of these applications are aiming at utilizing the high contrast of the dark-field signal compared to conventional transmission, and do not exploit the fact that quantitative information about the scattering structure can be obtained. The reason for this is that a GI setup has a specific length scale sensitivity which is defined by the autocorrelation length of the interferometer [6]. This effectively means that a single dark-field measurement corresponds to a single point of the autocorrelation function of the scattering structure which is insufficient for microstructural characterization. Therefore, to obtain quantitative information about the unresolved microstructures in a sample, multiple measurements over a range of autocorrelation lengths needs to be performed [18][19][20].
In general, the autocorrelation length ξ depends on the period of the interference pattern at the detector plane, the distance between the sample and the detector, and lastly, the wavelength. In a conventional XGI (Talbot-Lau interferometer, with three gratings) changing autocorrelation length can be achieved by either imaging the sample with phase gratings of different periods, an impractical and cost-inefficient solution, or by moving the sample along the optical axis, which in cone-beam geometries has the undesirable effect of varying magnification hindering further analysis. Changing the wavelength also has its disadvantages since, ideally, for each energy, a new set of gratings would be required. A so-called dual-phase interferometer was recently proposed to overcome these limitations of conventional XGI [21,22]. By using two phase gratings of small periods, enlarged fringes are generated at the detector plane, which can be recorded directly without the need for an analyzer grating [22]. Moreover, the relative distance between the two gratings can be tuned to change the fringe period, correspondingly changing the autocorrelation length. When XGI is performed with a polychromatic X-ray source, beam hardening effects from the sample can induce an apparent visibility reduction which does not correspond to small angle scattering [23][24][25][26], hence hindering any structural quantification. Nonetheless, in typical XGI setups, namely Talbot and Talbot-Lau interferometers, changing the autocorrelation length by moving the sample in the beam path neither changes the amount of beam hardening that takes place nor the spectral sensitivity (defined by the gratings and geometry) of the interferometer, hence leading to a constant bias across all measurements. On the contrary, in a dual-phase XGI, changing the relative distance between the gratings changes the geometry of the setup and thus the spectral acceptance of the system as well. Therefore, a varying effect of the sample beam hardening across different autocorrelation lengths is incorporated in the observed dark-field signal, which means that a usual imaging protocol where autocorrelation length is changed without accounting for beam hardening would not provide quantifiable structural information.
Although some recent works have presented the theoretical calculations for determining the fringe period seen at the detector and visibility [27] in the dual-phase XGI grating system, a complete theory incorporating the scattering and beam hardening effect of the sample is not available. This letter addresses the effect of the sample on the visibility reduction in the case of a polychromatic dual-phase XGI. We introduce a generalized formulation to calculate the interference fringe visibility at the detector when a sample is placed in the beam path using wave optics formulations. Furthermore, we demonstrate the varying effect of beam hardening on the visibility reduction across different autocorrelation lengths and eventually propose a method to account for it. Finally, we demonstrate the utility of this correction to quantify scattering structures and discriminate between visibility induced by beam hardening and actual scattering in unknown samples.

Theoretical description for visibility reduction in a dual-phase grating interferometer
Let us consider a standard dual-phase grating setup, as shown in Fig. 1. The first phase grating G 1 with a period of p 1 is placed at a distance of L 1 from the X-ray source. A second phase grating G 2 with a period of p 2 is placed at a distance d from the first grating. Both gratings have a duty cycle of 0.5. The detector is placed at a distance L from the source. The transmission functions of these phase gratings can be described by Fourier series as following where a n and b m are the Fourier components, and x 1 and x 2 the spatial coordinates at each grating plane respectively. By assuming that the gratings are pure phase gratings and that they introduce a phase shift of π at given wavelengths λ 1 and λ 2 respectively, the Fourier components a n for a given wavelength λ can be written as , when n = 0. 0, when n is even.
, when n is odd.
(2) Similar expressions can be written for b m by changing the design energy λ 1 in the above expression to λ 2 . The design energy can be calculated from the depth of the etching in the phase grating. For gratings which do not have a 0.5 duty cycle, generalized parameters can be obtained by calculating the Fourier transform of the grating transmission function. To calculate the fringe at the detector plane, we perform a series of wave propagations (Fresnel propagation) of the complex wave amplitude with and without a sample in the beam path. For the sample and the gratings, we utilize the projection approximation, thus taking into account their complex transmission functions. For each energy, the complex transmission of the sample is multiplied by G 1 and propagated for a distance d, and the resulting complex wavefield is then multiplied by G 2 and propagated for a distance L 2 to the detector where the intensity is calculated as the squared modulus. Eventually, the final fringe is calculated by incoherently summing the monochromatic fringes over the whole energy spectrum. The complete derivation is given in the Appendix. The visibility reduction, which corresponds to the dark-field image, is the ratio of the sample visibility V over the flat field visibility V 0 . It is shown in the Appendix that this visibility reduction can be written as the product of two terms. A spectral term V V 0 spectral which is the result of the beam hardening due to the sample, and a scattering term V V 0 SAXS originating from small angle X-ray scattering (SAXS) due to locally unresolved variations of the refractive index decrement δ of the sample The spectral term depends on the geometry of the setup (L 1 , L 2 , L, and d), the grating parameters a n and b m , the spectrum of the X-ray tube s (λ), the pixel size of the detector D, and the absorption coefficient of the sample α s (x, y, λ). Eventually, the spectral term is written as where the functions H and T are derived in the Appendix (see Eq. (19)). The spectral term changes when the distance between the gratings d is changed since the geometrical parameter T (p 1 , p 2 , L 1 , L 2 , d, λ) changes. Hence, the visibility reduction arising from the spectral term depends on the distance between the gratings when a sample is placed in the beam. The scattering based term is linked to the real space correlation function of the sample as shown in multiple publications [4][5][6]. Specifically, it corresponds to the value of the real space correlation function B at a distance ξ which is the autocorrelation length of the imaging setup.
Assuming weak dependence of phase signal on wavelength and given that phase variations are much smaller compared to pixel size of detector , it is shown that the autocorrelation length for a dual-phase grating interferometer is given by (see Appendix for complete derivation) where Therefore, the claim that the autocorrelation length can be changed by changing the distance between the gratings is established. To perform the microstructural analysis only the scattering term is of interest, hence an appropriate correction method is required to account for the spectral term. To estimate this term, we need the following parameters: the geometric parameters of the system, the grating parameters such as period, material, and depth, the spectrum of the X-ray source, and finally, estimation of the absorption coefficient of the sample or knowledge of the material composition. Most of these parameters can be either measured, estimated or simulated. The most complex is the energy dependent absorption coefficient of an unknown sample which can be estimated by the method presented in the Appendix (subsection (5.3)).

Fig. 1.
The standard dual-phase interferometer setup consisting of two phase gratings G 1 and G 2 . The first grating is placed at a distance L 1 from the X-ray source, the second grating at a varying distance d from the first grating and the detector with a pixel size of D at a distance L from the source. Finally, the sample is placed right before the first grating.

Results and discussion
In order to evaluate our model, we conducted a series of experiments demonstrating that the spectral visibility reduction can be predicted and accounted for. We started with imaging a water phantom, where visibility reduction would solely be due to beam hardening effects which our model can predict. For our second experiment, we showed that the already known real space correlation function of microspheres could be measured after the proposed correction is applied. Finally, we imaged a sample containing both areas with and without microstructure and demonstrated that with the proposed model, we could successfully distinguish those areas. The setup we used for all the reported experiments is the one proposed in [21]. Specifically, we utilized two identical Si gratings with a period of 1.2 µm and etched to a depth of 23 µm, introducing a π phase shift at 17 keV. The gratings were fabricated by interference Talbot lithography and deep reactive ion etching [28]. The source we utilized was a HAMAMATSU L0101 microfocal X-ray tube, operated at 70 kVp and 100 µA which resulted in a focal spot size of 9.5 µm. The detection of the generated fringe was performed with a Princeton Instruments PI-SCX:4300 X-ray detector. The pixel size of the detector was 24 µm allowing to record fringes with a period down to approximately 100 µm. The PI-SCX:4300 X-ray detector is equipped with a Gd2O2S:Tb (Gadox) scintillator, which is very efficient at energies up to 20 keV (50% efficiency). The first phase grating G 1 was placed at a distance of 0.5 m from the X-ray source. The second grating G 2 was placed at a distance of d downstream G 1 , and this distance was varied in order to change the autocorrelation length. Specifically, for all imaging examples, we changed d from 2 mm to 18 mm in steps of 1 mm. This resulted in a range of autocorrelation lengths between 104 nm and 998 nm. The X-ray detector was placed at a distance of 0.5 m from G 1 . Finally, the samples were placed right before the first phase grating G 1 resulting in an effective magnification of approximately two, which determined the effective imaging field of view (FOV) to 2.5 × 2.5 cm 2 .
For all the corrections, we utilized the measured X-ray spectrum at the operating current and acceleration voltage as reported in the Appendix. The grating parameters were estimated from cross sectional scanning electron images (SEM) of the gratings (SEM images are shown in Appendix).
The system can be operated both in phase stepping and in single shot mode. Since in the single shot mode the image resolution is defined by the fringe period, which in turn varies across measurements at different autocorrelation lengths, we decided to work with the phase stepping mode. In this mode, one of the two gratings is laterally shifted by one or more periods in multiple steps which are smaller than the period. At each step, an image is recorded which eventually allows recording a period of the interference fringe at each pixel. For all the imaging experiments we utilized five phase steps with an exposure time of 60 sec each. The retrieval of the three complementary signals is identical to the methods used in standard grating interferometry when phase stepping is used [29].

Water phantom
To demonstrate that beam hardening due to the sample can result in an observed visibility reduction even in lack of small angle X-ray scattering events, we imaged an absorbing sample (non-scattering) at different inter-grating distances. As a sample, we chose a Poly(methyl methacrylate) (PMMA) cuvette filled with water; the width of the cuvette was 0.5 cm. In Figs. 2(a), (b), and (c) the retrieved visibility reduction images for autocorrelation lengths of 104 nm, 525 nm, and 998 nm are shown respectively. A signal both for the PMMA container and the water is present. Nevertheless, no microstructure is present in the sample, and therefore this varying signal is due to the beam hardening as explained by our theory. Given the spectrum of the X-ray source (shown in the Appendix), the detector efficiency, the grating parameters, and the attenuation coefficient of water and PMMA, we can theoretically calculate the expected visibility reduction of the sample. In Fig. 2(d) the average visibility reduction for the marked area is compared to the theoretical calculation for the range of probed autocorrelation lengths. In general, a strong agreement is observed validating our theoretical model. Nonetheless, small discrepancies are observed due to uncertainties regarding the exact parameters utilized such as the X-ray spectrum and grating parameters.

Glass microspheres
As a next step, we wanted to demonstrate that quantitative visibility reduction due to small angle X-ray scattering can be retrieved by correcting for beam hardening effects. As samples, we utilized glass spheres of known diameters and volume fractions whose small angle scattering signals can be estimated. Two different microsphere powders [30] with diameters of 0.26 µm and 1.7 µm were placed on top of each other in a cylindrical container as shown in Fig. 3(a). These samples were also utilized for the demonstration of the tunability of the dual-phase grating interferometer in [21]. In light of our theoretical analysis, we will reinterpret the previous experimental results in order to demonstrate quantitativeness. In Figs. 3(b) and (c) representative visibility reduction images at two different autocorrelation lengths are shown. The two powders have the same absorption coefficient which would lead to similar beam hardening; therefore, the difference in the observed visibility reduction between the two powders is due to their difference in the size of microstructures. In Fig. 3(d) the average visibility reduction for the measured autocorrelation length range in the two marked areas is shown. Given that the material of the samples is known, we can estimate the spectral term V V 0 spectral which we use to correct the experimental data. This results in the curves shown in Fig. 3(e) which are compared to the expected visibility reduction curves from the numerical calculation of the small angle scattering from SASview [31]. In the Appendix, the scattering intensity profiles for both powders are shown, these scattering intensities can be translated to visibility reduction by applying the Hankel transform as explained here [18]. The calculated and experimental curves show a high accuracy which corroborates our model. As explained in the previous experiment a few deviations are present. In particular, the high frequency oscillations in the theoretical model of the 0.26 µm powder are blurred out in the experimental curve. This is due to the polychromaticity of the utilized source which is not modeled in the SAXS calculations.

Generalised sample
In this section, we show the applicability of the model for a sample whose sample properties (absorption coefficient, microstructured or not) are unknown. As samples we used, a glass block and polydispersed glass powder of unknown granularity (see Fig. 4(a) for a sketch). Assuming, the sample is made of single material, we can calculate the absorption coefficient. The absorption coefficient of the samples was estimated as explained in the Appendix in subsection 5.3. Both samples had a thickness of 0.5 cm. The visibility reduction images of the samples at two different autocorrelation lengths are shown in Figs. 4(c) and (e). The corresponding images obtained after correction are shown in Figs. 4(b) and (d). The average visibility reduction within the marked areas of Fig. 4(b) of the uncorrected and corrected data is shown in Figs. 5(a) and (b) for the glass powder and the glass block respectively. Interpreting the uncorrected visibility   Fig. 4(b), the errorbars correspond to the standard deviation. reduction signals would lead us to conclude that both samples contain subpixel structures since they both show a visibility reduction. Furthermore, one would possibly claim that the glass powder contains structures finer than 104 nm since the visibility reduction signal seems to have saturated. On the other hand, one would conclude that the glass block has structures in the range of the autocorrelation length since a signal fluctuation is observed. Alas, both statements are erroneous. After correction, the visibility reduction curves dramatically change as shown in Figs. 5(a) and (b). The corrected visibility reduction for the glass block seems to be constant and equal to one, meaning that no unresolved microstructures are present which is the ground truth. However, we see a strong visibility reduction for the glass spheres which is changing as the autocorrelation length is increased. Eventually, this curve could be fitted with an appropriate model in order to retrieve structural information.

Conclusion
Like conventional grating interferometers, the dual-phase grating interferometer suffers from beam hardening effects when operated with polychromatic sources. These effects cause visibility reduction to be detected even in the absence of small angle scattering. Through theoretical modeling and experimental investigations, we have demonstrated that this effect strongly depends on the distance between the two phase gratings, consequently affecting measurements across different autocorrelation lengths which are performed by changing the inter-grating distance. Therefore, having a deleterious effect on any statements made on alleged micro-structural properties of the sample. With a reasonable pre knowledge of the system parameters (such as spectrum and grating geometry), a rigorous theoretical model allows for mitigation of this effect. Specifically, we have demonstrated that the measured visibility reduction can be written as the product of two terms. The first one purely depends on the beam hardening introduced by the sample, and the second corresponds to the small angle X-ray scattering of the microstructure. Our model allows us to predict the spectral term and eventually retrieve the true scattering signal. Our description was validated with a series of experiments ranging from basic proof of principle of the beam hardening induced visibility reduction to the quantitative retrieval of the real space correlation function of microparticles. To conclude, we have demonstrated that quantitative retrieval of the underlying real space autocorrelation function can be achieved with polychromatic dual-phase X-ray grating interferometers. This comprises a significant step towards the application of the method for the quantitative investigation of the (sub)microstructure of materials and devices. Taking into account the advances in grating fabrication technologies [32,33], the method could be translated to higher energy ranges, hence allowing the investigation of devices such as Polymer Electrolyte Fuel Cells (PEFCs) under operation. Such fuel cells contain multiple layers whose pore sizes distributions affect water transport and eventually their performance. With our method, we could get a close look at the water propagation through the submicrometer porous structures which could lead to the design optimization of the fuel cell.

Wave propagation
We consider that the sample is placed right in front of grating G 1 as shown in Fig. 1. The complex refractive index of the sample can is written as 1 − δ + iβ. The electric field for a given wavelength λ after the sample is given by E 0 exp[α s (x 1 , y 1 ) /2] exp[iφ (x 1 , y 1 )] according to the projection approximation, where E 0 is the amplitude of the incident X-rays, and α s and φ are given by 2πδ(x,y) λ and 4πβ(x 1 ,y 1 ) λ respectively, where (x 1 , y 1 ) are the coordinates at the sample and G 1 plane. Additionally, we assume that the sample has only fine phase variations that cannot be resolved by the detector. Initially, we propagate the multiplied transmission functions of the G 1 and the sample. By utilising the Fresnel propagator. This results in the field at the second grating G 2 being equal to where M is the magnification of G 1 in relation to the source given by Since the beam has a cone geometry we need to rescale the coordinates (x 1 , y 1 ) at the plane of the second grating G 2 with coordinates (x 2 , y 2 ). Eventually resulting in the following expression (10) After the second grating we propagate for a distance L 2 and rescale the coordinates with N = L 1 +d L which results in the field at the detector plane expressed in coordinates (x, y) For the sake of simplicity, we will consider only one dimension x. We can now calculate the intensity of a single wavelength λ at the detector plane as follows.
Which can be further written as The above equation can be regarded as a series of trigonometric functions with periods that depend on the specific combination of the n, m, q, and l. In particular, the trigonometric terms can be written as At this point, we assume that the utilized detector has a pixel size of D/2 which is much larger than the fine periods p 1 and p 2 . Given this, it is reasonable to assume that a resolvable interference pattern will only occur when q = −l. Therefore, we can simplify the intensity as where the interference periods at the detector P det is given by and The terms H and T are given by H (a n , a m , l) = exp (−α s (MNx, MNy, λ)) a n (λ T (n, m, l, p 1 , To calculate the recorded intensity from the detector we need to take into account the point spread function P(x), the finite source size S, the spectral acceptance of the detector sa (λ), and the spectrum of the source ss (λ). Without loss of generality, sa (λ) and ss (λ) can be combined as a single spectral distribution, given by s (λ) = ss (λ) sa (λ). All of these, result in the following intensity within one pixel at the detector plane where P(x) is The term H (a n , a m , l) is inversely proportion to l 2 , therefore for larger l the diffraction efficiency of the gratings is dropping fast, which in turn allows us to consider only the following terms l = 0, l = 1 and l = −1. Meaning that the recorded intensity can be written as the following sum When l = 0, there is no dependence on x which leads to a DC term that corresponds to the transmission. The terms from l = 1 and l = −1 contribute to the interference pattern that is seen on the detector which carries the scattering and phase information.

Absorption contrast
To calculate the transmission of a sample, we first acquire a flat image without a sample in the beam, the DC term of the fringe is then given by When a sample is introduced in the beam, the DC term becomes Which leads to the following sample transmission

Visibility reduction contrast
The dark-field signal is retrieved from the visibility of the interference pattern which corresponds to the modulation of the oscillating term which is given by Since, we use phase gratings with a duty cycle of 0.5, a n , b m = 0 when m, n are even numbers. In order to calculate the visibility of the above equation, a number of assumptions are made.
First, we assume a weak dependence of the phase function φ(x) with regard to the wavelength. This is achieved by introducing an effective wavelength which is discussed in detail in the next section. Furthermore, we assume that the exp (−α s (MNx , λ)) is slow varying within a pixel, and therefore this term can be taken out of the spatial integral. These assumptions allow us to simplify the integral to the following expression.
The attenuation, spectrum and geometry related parameters can be grouped together in a single term Which allows us to write the oscillating part of the fringe as In the case of no sample in the beam, the above expression simplifies to Which eventually leads to the following visibility reduction that can be written in the following form where and Hence, we are able to show that the visibility reduction in a dual-phase interferometer arises from two terms: one spectral induced and one scattering based term (SAXS).

Deriving effective wavelength and autocorrelation function
To calculate the autocorrelation length, we first need to calculate the effective wavelength. By starting from Eq. (40) and assuming that the sample is weakly absorbing we get (36) From the above integral we group together the terms which vary with x as I partial (x) As the detector pixel size D is small or comparable to λL 2 p 2 we get We further simplify the above expression by assuming that the spectrum can be written as a discreet function where Using the sifting property of the Dirac delta function, the integral can be simplified as follows Under the assumption that φ(x) can be modelled as Gaussian random process [34] we get where and τ is the lateral length of correlation and σ = φ 2 (0). In the case of q<τ, we use the first order Taylor expansion on γ(q) Therefore the autocorrelation length and the effective wavelength can be written as where, When q ≥ τ the autocorrelation function saturates and therefore Which means that the approximation of λ eff holds in this case as well.

Estimating absorption contrast
In order to estimate the absorption coefficient µ of the sample, we assume that it varies with the third power of the wavelength λ, such that µ = rλ 4 + qλ 3 is a valid approximation [35]. Eventually, we would like to find a way to estimate the parameters: r and q. For a given spectrum s(λ) and a sample thickness T the transmission of the sample can be calculated as By conduction measurements with different spectra, the spatially modulated parameters q(x, y) and r(x, y) can be estimated by least-square fitting of the measured transmission A(x, y). For our experiments, we used spectra at 35 kVp, 50 kVp, and 70 kVp, which we measured using the Ametek X-123 Complete X-Ray Spectrometer with Si-PIN Detector. The measured spectrum (70 kVp and 100 µA) which was used for the imaging experiments is shown in Fig. 6. When there are one or more absorption edges the estimation of attenuation coefficient may not be very accurate in the absorption edge regions (e.g. K edge). Other attenuation coefficient models could be used for accurate estimation around the edge region. However, such models require additional parameters like atomic number. Given that additional information about the sample like atomic number and density are known, such models could be used to estimate wavelength dependent absorption coefficient [35].

Theoretical SAXS curves
We were able to obtain the SAXS curves shown in Fig. 7 for Silicon dioxide microsphere samples using SASview. These curves were transformed into visibility reduction curves, as described in [18].

SEM images of the gratings used in the experiments
We used two phase gratings of period 1.2 µm, with duty cycle of 0.5 and etched to depth of 23µm. The SEM image of the gratings is shown in Fig. 8.