Sample phase gradient and fringe phase shift in dual phase grating X-ray interferometry

One of the key tasks in grating based x-ray phase contrast imaging is to accurately retrieve local phase gradients of a sample from measured intensity fringe shifts. To fulfill this task in dual phase grating interferometry, one needs to know the exact mathematical relationship between the two. In this work, using intuitive analysis of the sample-generated fringe shifts based on the beat pattern formation mechanism, the authors derived the formulas relating sample's phase gradients to fringe phase shifts. These formulas provide also a design optimization tool for dual phase grating interferometry.


Introduction
In x-ray interferometry based phase contrast imaging, one of the key tasks is to retrieve local phase gradients of a sample. More specifically, the sample's local refraction angle, α(x, y), is proportional to the local phase gradient of the sample [1][2][3][4]: where λ is x-ray wavelength and Φ s (x, y) denotes local sample phase shift, which is equal to Φ s (x, y) = −λr e ∫ ρ e (x, y, s)ds. In this integral ρ e denotes the sample's electron density and r e is the classical electron density which equals to 2.82 × 10 −15 m. In Eq. (1), ∂Φ s (x, y)/∂ x denotes sample's phase gradient along the direction perpendicular to grating slits. In x-ray interferometry one measures the fringe shift ∆φ generated by sample refraction. In order to retrieve local phase gradients of the sample, one must find out functional relationship between intensity fringe shifts and sample phase gradients. It turns out that the mathematical relation between them depends not only on the interferometer's geometrical configuration and phase grating periods, but also on the position of the sample. So one important task in x-ray interferometry is to rigorously derive this functional relationship which will enable retrieval of sample phase gradients from measured interference fringe shifts. Combined with tomography, the retrieved sample gradients from angular projections can be used to reconstruct 3D maps of sample electron densities. In Talbot-Lau interferometry, where only a single phase grating is employed as the beam splitter, thereby high-contrast interference fringes form at certain distances downstream ( Fig. 1(a)). Sample refraction distorts the intensity fringe pattern, and the mathematical relations between fringe shifts and sample phase gradient is well described in literature, and the specific formulas depend on whether the sample is placed upstream or downstream of the phase grating [1][2][3][4].
∆φ(x, y) =      λR 1 (R 1 +R 2 )p eff L D · ∂Φ s (x,y) ∂x , if Sample at downstream of G 1 , where R 1 is the source-to-phase grating (G 1 ) distance, and R 2 the distance from G 1 to detector entrance as is shown in Fig. 1(a). Note that the absorbing grating G 2 in Fig. 1(a) serves as an analyzer for resolving the intensity fringes. In Eq. (2), p eff denotes the effective grating period of the phase grating G 1 , p eff = p 1 /2 if it is a π-grating, and p eff = p 1 otherwise. In Eq. (2) L D denotes sample-to-detector distance, and L S the sample-to-source distance. However, to measure fringe phase shifts with common image detectors, one usually has to use a fine-pitch absorbing grating placed at the detector's entrance. One indirectly detects fringe pattern through grating scanning, which is also called phase stepping procedure [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. The absorbing grating blocks more than half of transmitting x-ray, and will significantly increase radiation dose in imaging exams. This is a serious disadvantage of Talbot-Lau interferometry for radiation-dose sensitive imaging applications such as medical imaging.
Recently, dual phase grating x-ray interferometry emerges as an attractive alternative [16][17][18][19]. A typical setup of dual phase grating interferometers employs two phase gratings G 1 and G 2 as the beam splitters, as is shown in Fig. 1(b). The split waves transmitting through the phase gratings create different diffraction orders interfering with each other. The intensity fringe pattern includes a beat pattern [18]. The imaging detector D has a pixel size much larger than periods of both the phase gratings ( Fig. 1(b)). The detector just resolves the beat patterns of large periodicities, and renders other fine patterns to a constant background. Hence, different from Talbot-Lau interferometry [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15], the dual phase grating interferometers directly detect interference fringes without the need of absorbing analyzer grating. This advantage brings significant radiation dose reduction as compared to Talbot-Lau interferometry.
However, to retrieve local phase gradients of a sample, those formulas in Eq. (2) are not applicable to dual phase grating interferometers, in which refracted x-rays may pass through two phase-gratings rather than a single phase grating. In this work we set out to derive the corresponding formulas for dual phase grating x-ray interferometry. In section 2, using intuitive analysis of the sample-generated fringe shifts and the intensity beat pattern formation, we derived formulas that relate the sample's phase gradients to fringe phase shifts. In section 3 we perform wave propagation simulations to validate the formulas derived in section 2. We also show how the angular sensitivity of a dual phase grating interferometer is determined in dual phase grating interferometry. We conclude the paper in section 4.

Methods
We start our derivation from explaining the fringe formation mechanism in dual phase grating interferometry. Figure 1(b) shows the geometrical configuration of a dual phase grating interferometer. It consists of a source S, two phase gratings G 1 and G 2 , and an imaging detector D. The periods of the first and second phase gratings are p 1 and p 2 respectively, R s is the source-to-G 1 distance. Note that R g is the spacing between the two phase gratings, and R d denotes the G 2 -to-detector distance ( Fig. 1(b)). For sake of convenience in discussion, we define several magnification factors as follows: where M g 1 represents the geometric magnification factor from G 1 to detector plane, and M g 2 is the geometric magnification factor from G 2 to detector plane. Obviously, in absence of G 2 -grating, a diffraction order of the intensity pattern generated by G 1 phase grating alone would be represented by exp i2π(l · x)/(M g 1 · p 1 ) , where l is an integer and indexes the diffraction order. Similarly, in absence of G 1 -grating, a diffraction order of the intensity pattern generated by G 2 grating alone would be represented by exp i2π(m · x)/(M g 2 p 2 ) , where m indexes the diffraction order of the fringe in absence of G 1 -grating. As is shown from our theory of dual phase grating interferometry [18], x-ray irradiance at detector entrance is a result of cross-modulation between the fringe patterns generated by phase gratings G 1 and G 2 respectively. Hence the irradiance pattern at the detector entrance is a weighted sum of different diffracted orders. Each of the orders in dual phase grating interferometry is indexed by two integers (l, m) and represented by a product as: However, among these diffracted orders, there are beat patterns formed by those diffraction orders characterized by l = −m. With proper setup such that 1/(M g 1 p 1 ) close to 1/(M g 2 p 2 ) (for example R g R s , R d and p 1 ≈ p 2 ), the period of these beat patterns can be made much larger than p 1 and p 2 . As long as the detector pitch p D p 1 , p 2 , the detector-resolved intensity pattern is just the beat pattern, since the detector renders all other fine fringes of those l −m orders to a constant background. The beat pattern consists of different harmonics resolved by the imaging detector. The period of the fundamental mode of the beat patterns is [16,18]: Obviously, the period of the m-th harmonics is p fr /m. According to Eq. (5), the period of the beat pattern can be much larger than grating periods p 1 and p 2 , provided that the inter-grating spacing R g is set much smaller than the distances R s and R g . With large fringe periods the beat pattern may be resolved by a common imaging detector of pixels in few tens of micrometers. Once the fringe formation mechanism is understood, we are now ready to derive the mathematical relationship between fringe phase shift and sample phase gradient. First, we consider a relatively simple case of sample being downstream of the second phase grating G 2 . Let a sample's local refraction angle be α(x, y), and the sample-detector distance be L D . Referring to Fig. 2(a), the fringe lateral shift caused by the refraction is equal to αL D . Since the period of the m-th diffracted order in the beat pattern is p fr /m, using Eq. (1) and Eq. (5) we found the corresponding fringe phase shift ∆φ down (x, y, m) as In above equation the subscript in ∆φ down (x, y; m) indicates that the sample is positioned downstream of the second grating. Equation (6) shows that the fringe phase shift scales with sample phase gradient, and the proportional constant is given by mλL D /p fr . In practice, the dominant diffracted order for π/2-gratings is m = 1 and that for π-gratings is m = 2.
On the other hand, there is a different way to place the sample in the beam. In dual phase grating interferometry, one does not place sample between two phase gratings, since the spacing between them can be as narrow as few millimeters. Therefore, another way of sample placement is to place the sample upstream of the first phase grating G 1 . Assume that the sample is placed upstream of G 1 and the source-sample distance is L S . Let us firstly consider how the sample phase gradient would affect G 1 -associated fringe in absence of G 2 grating. Consider a ray refracted by the sample with an angle α. Referring to Fig. 2(b), the refracted ray propagates to G 1 grating, and the refracted ray makes an angle α eff,1 to the locally undisturbed ray. We call α eff,1 the effective refraction angle. In practice, sample refraction angles are much smaller than 1 radian. Based on the geometry depicted in Fig. 2(b), under assumption of α 1, it is easy to find that α eff,1 = (L S /R s )α. That is, compared to the sample refraction angle, the effective refraction angle is reduced by a factor, which is the ratio of source-to sample distance over the sample to grating distance. Propagating over a distance R g + R d to the detector, in absence of G 2 grating, this refracted ray would cause a lateral shift of the G 1 -associated fringe by α eff, It is easy to see from Eq. (4) that in absence of G 2 grating the period of the l-th order G 1 -fringe is M g 1 p 1 /l, hence the G 1 -fringe phase shift caused by the sample refraction is equal to: A similar derivation can be applied to G 2 -fringe phase shift, assuming absence of G 1 grating. Referring to Fig. 2(c), the effective refraction angle is α eff,2 . As compared to the refraction angle α, it should be reduced by the ratio of source-to-sample distance over the sample to grating distance, thereby α eff,2 becomes α eff,2 = L S /(R s + R g )α. Propagating over a distance R d to the detector, this refracted ray would cause a lateral shift of the G 2 -associated fringe by α eff,2 · R d = (L S /(R s + R g )) · R d α. Equation (4) shows that in absence of G 1 grating the period of the m-th order G 2 -fringe is M g 2 p 2 /m. Hence, the sample refraction generates a phase shift for the m-th order G 2 -fringe: In dual phase grating interferometry the interference fringe pattern interferometry is indexed by two integers (l, m) and represented by a product as exp i2π(l · x)/(M g 1 · p 1 ) × exp i2π(m · x)/(M g 2 · p 2 ) , thus the refraction-generated fringe phase shift of the order (l, m) is In dual phase grating x-ray interferometry, one usually adopts a common imaging detector of pixels much larger than grating periods. As is discussed in section 2 on fringe formation mechanism, due to the pixel averaging effects, the detector resolves only the beat patterns of diffracted orders (l = −m, m). This being so, Eq. (9) gives the sample-refraction generated fringe phase shift for the m-th order beat patterns as follows: where Comparing Eqs. (5) and (11), one can see p 0 is always equal to p fr if the two gratings have the same period (p 1 = p 2 ). But for p 1 p 2 , p 0 can be greater or less than p fr , depending on the geometric setup. More explicitly, we rewrite the sample-refraction generated fringe phase shift for the m-th order beat patterns as: In above equation the subscript in ∆φ up (x, y; m) indicates that the sample is positioned upstream of the first grating G 1 . Equation (12) shows that, when a sample is positioned upstream of the 1 st phase grating, the fringe phase shift scales again with sample phase gradient, and the proportional constant is −mλL S /p 0 . Especially, the closer to the 1 st phase grating the sample is, the larger magnitude the fringe phase shift is. Obviously, this proportional constant in Eq. (12) is different to that in Eq. (6), which corresponds to the setups with sample positioned downstream of the 2 nd phase grating.

Results
Combining Eqs. (6) and (12) together, we found the measured fringe phase shift ∆φ is proportional to sample's phase gradient with a proportional constant ξ given by where m = 2 if the phase gratings are π gratings and m = 1 otherwise. In Eq. (13), p 0 and p fr are given in Eqs. (11) and (5) respectively. Equation (13) shows that this formula lays foundation of quantitative imaging for sample phase gradients with dual phase grating x-ray interferometry. Based on our understanding of fringe formation mechanism [18], we developed above intuitive and heuristic method to derive Eqs. (13).
To validate the fringe shift formula of Eq. (13), we compare the calculated fringe phase shift of formula Eq. (13) to that determined through numerical simulations. We employed four different interferometer setups in the simulation study. In the first simulation, a point x-ray source of 20keV design energy and dual-π phase gratings of period p 1 = p 2 = 1µm are employed. The geometric setup is R s = R d = 450mm and R g = 5mm. With this geometric setup, the fringe can attain a high visibility pattern and large fringe period p fr = p 0 = 181µm as well [18]. So with a detector of p D = 36.2µm pitch, one period of the resolved fringe will take up five pixels. To validate the fringe shift formula of Eq. (13), a sample is placed half way between the source and G 1 -plane, i.e. L S = R s /2 = 225mm. The sample is assumed a sphere of diameter 4 mm filled with 100% adipose tissue. We employed ray-tracing to get the projected sample attenuation A 2 (x, y), and the phase map of the sample. From sample phase map we computed sample phase gradients.  (13). In the simulation, the setup consists of two π-gratings and a 20keV point source. A sphere shaped sample of adipose tissue is placed upstream of the 1 st grating and downstream of the 2 nd grating respectively. We employed ray-tracing to compute the theoretical map of sample's attenuation A 2 (x, y) and fringe phase shift ∆φ(x, y). We then simulated Fresnel wave propagation in the interferometer to generate the intensity fringe patterns. The fringe phase shift ∆φ(x, y) and attenuation map A 2 (x, y) are retrieved through the Fourier method. The ray-tracing generated theoretical maps of attenuation A 2 (x, y) and fringe phase shift ∆φ(x, y), when sample is place half way between the source and the 1 st grating, are shown in Fig. 3(a) and (d); while the retrieved A 2 (x, y) and ∆φ(x, y) are shown in Fig. 3(b) and (e) respectively. For the purpose of comparison, the profiles of the attenuation and fringe phase shift along the central line across grating are shown in Fig. 3(c) and (f). The blue curves correspond to the theoretical values, while the red curves correspond to the values retrieved from intensity fringes. Figure 3(g) is the profiles of the fringe phase shift ∆φ(x) along the central line across grating, when the sample is placed to the plane half way between G 2 grating and detector. The good agreement between the theoretical values and the fringe-pattern retrieved values provide a validation of Eq. (13). For details, see text.
Using Eq. (13) and the sample phase gradients we theoretically predict a map of the fringe shift ∆φ(x, y). The theoretical map of sample attenuation A 2 (x, y) and map of fringe shift ∆φ(x, y) are shown in Fig. 3(a) and (d) respectively. We then employed Fresnel propagation down stream from sample plane through the gratings G 1 , G 2 and finally to the detector plane to get the projected image. We then employed the Fourier method for retrieval of the fringe phase shift ∆φ(x, y) and attenuation map A 2 (x, y). The maps of retrieved attenuation and phase shift are shown in Fig. 3(b) and (e). As comparison, profiles of the attenuation A 2 (x) and phase shift ∆φ(x) along the central line across grating are shown in Fig. 3(c) and (f) respectively. The dashed blue curves represent the theoretical value, while the solid red curves are the retrieved values from the intensity fringes. The good match between the theoretical and retrieved values of ∆φ(x) validates the fringe phase shift formula of Eq. (13) for sample upstream G 1 .
In the second simulation, we keep the same setup as in the previous simulation but placed sample to the half way between the G 2 -grating and detector plane, i.e, L D = R d /2 = 225mm. The diameter of the sample sphere is also increased to 12mm to increase the resolution. The profiles of the phase fringe shifts ∆φ(x) along the central line across grating are shown in Fig. 3(g), in which the dashed blue line represents the theoretical value and the solid red line is the retrieved one. The closeness of the two profiles validates the fringe phase shift formula of Eq. (13) for sample downstream G 2 . Comparing Fig. 3(f) and (g) in previous page one finds the profiles of ∆φ(x) are reflected, in other words, the values of ∆φ(x) change signs when the sample moves from upstream to downstream, as is predicted by the sign difference in Eq. (13) for the two scenarios. Noting this interesting feature is important in practice for accurately retrieving sample phase gradients from measured fringe phase shifts. Fig. 4. In this simulation, the second grating is replaced with a p 2 = 1.1µm period π grating but the first grating is kept the same as that in Fig. 3. The geometric setup is also changed to R s = 450mm, R g = 40mm, and R d = 380mm and detector pixel pitch p D = 39.1µm. By placing the sample to plane L S = R s /2 = 225mm and to plane L D = R d /2 = 190mm, one gets the retrieved plots of ∆φ(x) for sample upstream of G 1 (Fig. 4(a)) and for sample downstream of G 2 (Fig. 4(b)).
In the third and fourth simulations, we replace the second phase grating G 2 with a p 2 = 1.1µm period π-phase grating and keep the first grating G 1 the same. Corresponding to this change, the geometric setup is changed accordingly for attaining good fringe visibility [18]. So in the third and fourth simulations we set the interferometer with R s = 450mm, R d = 380mm, and R g = 40mm.
In this configuration, the fringe period, as is determined by Eq. (5), is p fr = −191.4µm and p 0 = 11.67µm [18]. In addition, the detector pixel size is changed to p D = 31.9µm, thereby one fringe period can take up 6 pixels. We then simulate the fringe formation processes by placing the sample upstream to plane L S = R s /2 = 225mm in the third simulation, and downstream to plane L D = R d /2 = 190mm in the fourth simulation. The results are shown in Fig. 4. The profiles of ∆φ(x) along the central line for sample upstream of G 1 are shown in Fig. 4(a) and those for sample downstream of G 2 are shown in Fig. 4(b), in which the dashed blue lines are the theoretical values and the solid red lines correspond to the retrieved ones. The close match between the red and blue lines again validated Eq. (13).
In literature the ratio ∆φ/(2πα) is called the angular sensitivity. It is a measure of fringe phase shift generated per unit sample-refraction angle. Apparently, the larger the ∆φ/(2πα) is, the larger the magnitude of the fringe phase shift is. Equation (13) can be used to compute the angular sensitivity of a given setup for design optimization of an interferometer. For example, Tab. 1 lists the angular sensitivity values computed by using Eq. (13) for the two setups investigated in the 3 rd and 4 th simulations. The notations used for describing the setups are explained in the text.
This table shows clearly that the angular sensitivity of the setup with the sample placed upstream of the 1 st phase grating is about nineteen times that of the setup with sample positioned downstream of the 2 nd phase grating. With a given geometric setup, the closer the sample is placed to the first grating for sample upstream or the closer the sample is to the second grating for sample downstream, the larger the angular sensitivity is. For given sample, the setup with a higher angular sensitivity will generate larger magnitude of fringe phase shifts for measurement. Carefully examining Fig. 4(a) and (b), one may notice that the values of ∆φ(x) do not change signs when the sample moves from upstream to downstream, as their features differ from that exhibited in Fig. 3(e)-(f). In fact, this difference is exactly predicted by Eq. (13). In fact, because of the different interferometer configuration, p fr and p 0 in these simulations are changed to p fr = −191.4µm < 0 and p 0 = 11.67µm > 0. Equation (13) dictates that ∆φ(x) will not change sign when sample is moved from upstream of G 1 to downstream of G 2 . This example demonstrates again that Eq. (13) gives accurate mathematical relation between sample phase gradient and fringe phase shift, accurate not only for predicting their magnitudes but also for their signs. Hence Eq. (13) provides a useful tool in quantitative phase contrast imaging based on dual phase grating interferometry.

Discussion and conclusions
One important task in x-ray interferometry is to rigorously derive the functional relationship between measured interference fringe shifts and sample phase gradients. Combined with tomography, the retrieved sample gradients from different angular projections can be used to reconstruct quantitative volumetric images of sample electron densities. The mathematical relationship between fringe shifts and sample phase gradients depends not only on an interferometer's grating periods and geometric configuration, but also on the sample position within the interferometer. The formulas for retrieving sample phase gradient in Talbot-Lau interferometry are well known as summarized in Eq. (2). However, there is a need to provide corresponding formulas applicable to dual phase grating interferometry with various sample position. In this work we fulfill this need using a novel intuitive analysis based on the fringe pattern formation mechanism. The derived formulas of Eq. (13) provides a useful tool for retrieving sample phase gradients in dual phase grating interferometry with any sample position. These formulas have been validated by extensive simulations, not only in the magnitudes but also in the signs of the retrieved sample phase gradients, as is shown in section 3. In addition Tab. 1 shows that sample position has a significant effect on angular sensitivity of a setup.
Moreover, the formulas of Eq. (13) can also be used to compute the angular sensitivity ∆φ/(2πα), which represents the fringe phase shift generated per unit sample-refraction angle in a given setup. As is discussed in section 3, although it is generally desirable to have a large ∆φ/(2πα) for a setup, but it is only a nominal sensitivity measure. The true sensitivity of a interferometer setup is the minimum detectable refraction angle, which depends on the fringe visibility and photon quantum noise with the setup. Detail discussion of this topic is out the scope of this paper.
In Eq. (13), the proportion constant ξ between fringe phase shift and sample phase gradient is also called the auto-correlation length of an interferometer setup. This quantity ξ characterize the setup's sensitivity for detecting small angle x-ray scattering from the sample's fine structures. For detail readers are referred to dark field imaging literature [ref]. We found that the auto-correlation length ξ is inversely proportional to fringe period p fr of Eq. (5), if the sample is positioned downstream of G 2 grating. On the other hand, when the sample is positioned upstream of G 1 , the auto-correlation length ξ is inversely proportional to a length-quantity p 0 of Eq. (11). The size of p 0 has a physical meaning in a different context. Namely, if a source grating were incorporated, then the periodicity of the source grating should be set to p 0 for achieving good fringe visibility [20].
In conclusion, in this work, using intuitive analysis of the sample-generated fringe shifts based on the beat pattern formation mechanism, the authors derived the formulas relating sample's phase gradients to fringe phase shifts. These formulas provide also a design optimization tool for dual phase grating interferometry.