1. Introduction

Grating-based X-ray interferometry is an important imaging technology finding use in multiple application areas, predominantly medical, security, and materials characterization [1–9]. The Talbot-Lau geometry allows a laboratory X-ray tube source to be used in conjunction with three microfabricated gratings to acquire images of the sample transmission (attenuation), differential phase contrast (DPC), and small-angle scattering (dark field). All three imaging modalities can provide important information about the sample. While attenuation and DPC imaging probe the atomic properties of a material, the dark-field effect is distinctly different in that it arises from ultra-small-angle coherent scattering due to micron-scale density fluctuations, like those associated with a material’s granularity. This makes dark-field imaging particularly well suited for detecting powders, exploring grain formation, identifying small voids, or looking at fibrous structures [10–17]. The dark-field signal is found from changes in the X-ray fringe visibility, and while other physical mechanisms such as beam hardening and sharp-edge refraction can also cause localized visibility variations [9,18], we assume here that “dark-field imaging” refers solely to a measurement of visibility contrast arising from small-angle scattering. There have been several studies centered on using dark-field imaging for biomedical applications (see [19–21] and references therein), in particular for examination of lung tissue where healthy alveoli create small-angle scattering that reduces fringe visibility, so that any observed increase in fringe visibility may indicate the presence of disease that degrades the alveoli structure [21–24].

Existing models of dark-field imaging include convolutional blurring [25–27] and wave-optics formulations [28–30]. In this two-part paper we develop a statistical optics model of the dark-field signal that rigorously captures the key features from a new, more physically intuitive perspective, and demonstrate through various examples that it is also amenable to efficient numerical implementation. Moreover, this model, unlike the classical model often used to describe X-ray scattering (e.g., see [31], Sec. 1.1), is not subject to the Born approximation (weak-scattering) in the far-field limit. In contrast, the statistical optics model can accommodate strong scattering with subsequent propagation subject only to the much less restrictive Fresnel approximation.

In the literature on X-ray dark-field grating interferometry, reference is often made to the so-called “correlation length” of a Talbot-Lau interferometer, which depends on its construction parameters as well as the location of the sample [21,27,29,32]. In physical terms this correlation length is simply an effective lateral separation that is imparted to the dominant diffraction orders (at the design energy), so the resulting fringe visibility can be estimated by evaluating the field autocorrelation function in the detection plane at this particular separation distance. The width of the field autocorrelation function also depends on the scattering object location. However, this approach neglects other diffraction orders that may in fact play an important role in determining the Talbot fringe visibility. In addition, when considering the dark-field signal associated with particle suspensions, the autocorrelation function of the sample’s refractive index distribution (which is included in the expression for the field autocorrelation function) is often assumed to be given by that for a single particle (e.g., a microsphere), but doing so is only valid in the regime of low particle concentration. In this paper we overcome both of these limitations.

We begin by taking a more general approach that includes evaluation of the ensemble-average mutual coherence of the beam upon passing through a Gaussian random-phase scattering object, combined with the change to the mutual coherence as the beam propagates to the detector plane. The role of the beam splitter grating (typically denoted G$_1$) is then incorporated, with all of the relevant diffraction orders included, leading to a more accurate model for the Talbot pattern and its visibility. The result demonstrates that there are in fact multiple “interferometer correlation lengths” involved, one for each fringe harmonic of the Talbot pattern. More correctly, we should say there are multiple separations at which the field correlation function must be evaluated in the detection plane. We then apply this model to a sample consisting of random dielectric microspheres having an arbitrary packing fraction (by means of a generalized autocorrelation analysis, although subject to a near-field or geometric projection approximation). The flexibility of our model allows us to readily determine the visibility as a function of various object and system parameters in a very computationally efficient manner. Our overall goal with this new framework is to provide a more comprehensive and unified theoretical picture of the dark-field effect, one that may prove helpful in quantifying experimental data or in optimizing system construction to enhance the dark-field signal for a particular class of objects. The primary focus here centers on a monochromatic point-source configuration that can serve as the basis for modeling a more general polychromatic extended-source setup.

Part 1 of the paper is organized as follows. Section 2 provides a brief review of the conventional X-ray grating-based system layout upon which our model is based. Section 3 then outlines a simple method for calculating the Talbot fringe intensity pattern that takes into account the spatial coherence of the beam. Section 4 provides a detailed development of the statistical optics model. In Appendix A we provide details associated with formulating the 2D projected autocorrelation or autocovariance functions of the scattered beam’s phase profile (as needed for use in the statistical optics model) from known 3D autocorrelation functions of the scattering medium’s random index distribution. In Supplement 1 we provide a derivation of the index-distribution correlation function for an ensemble of monodisperse hard microspheres based on an approach first put forth by Torquato and Stell [33]. Lastly, in Supplement 1, we also review the convolutional blurring dark-field model that has been utilized in the literature, address some of its deficiencies, and propose a modified version that exhibits close agreement with the statistical optics model.

2. System layout

Figure 1 shows a schematic representation of a typical three-grating Talbot-Lau shearing interferometer using a standard X-ray tube source. The system is nominally designed at a specific beam energy $E_0$ (or wavelength $\lambda _0$), but it will operate with a polychromatic source over some usable bandwidth centered about $E_0$. The gratings are generally fabricated from etched silicon, although other materials can be utilized. G$_1$ is a diffraction grating, typically a phase grating, that creates multiple diffraction orders, thereby providing a beam-splitter function central to operation of the shearing interferometer. Upon propagation over a distance $d$, the overlapping G$_1$ diffraction orders create a Talbot fringe pattern, or interferogram. The relationships between the design energy and the locations and periods of the gratings are well-established in the literature and are not reviewed here (e.g., see [34]).

 

Fig. 1. Diagram of a Talbot-Lau X-ray interferometer. The object is shown in front of G$_1$, but it can also be placed after G$_1$.

Download Full Size | PDF

Placing objects in the beam path alters the interferogram, and a measurement of the perturbed pattern relative to the unperturbed reference pattern yields the object information of interest. G$_0$ and G$_2$ are absorption gratings, made for example by electroplating gold on top of the corresponding silicon structures, that function as attenuating masks. The G$_0$ grating spatially modulates the source intensity so that the fringes formed from radiation passing through each of the transmitting slots reinforce one another; therefore this grating effectively modifies the spatial coherence of the extended source to allow fringe formation to occur. The G$_2$ analyzer grating, which has a period equal to that of the X-ray Talbot fringe, is used for detection of the fringe properties (average intensity, phase, and visibility) by translating it across the fringe pattern and measuring the average transmittance per pixel as a function of the grating position. This allows use of an imaging detector that has pixels with a linear dimension much larger than the fringe period. Because we assume a point source, the G$_0$ grating properties do not come into play, although a collection of uncorrelated point sources could be used to mimic an extended incoherent source modulated by G$_0$. For detection, we assume that G$_2$ is an ideal square-wave amplitude grating having complete modulation of its intensity transmittance from zero to one.

3. Talbot fringe formation

3.1 G$_1$ phase grating diffraction analysis

We assume that G$_1$ is a binary phase grating with period $p_1$ as depicted in Fig. 2. It has a single-period unit cell comprised of two adjacent rect functions, one having a width of $\alpha p_1$ and an amplitude of $e^{i\theta }$, and the other a width of $(1-\alpha )p_1$ and an amplitude of one. The complex transmittance of the grating can be represented by a convolution of the unit cell with an array of delta functions separated by the grating period (in the form of an appropriately scaled comb function, e.g., see [35], p. 83). Thus,

Here $\theta = 2\pi (n-1)\Delta _1/\lambda$ is the amplitude of the G$_1$ phase grating, $n$ being the refractive index of the grating substrate and $\Delta _1$ its etch depth. The diffraction amplitudes are found from a Fourier transform of the complex transmittance ([36], Sec. 2.4.1),

 

Fig. 2. Binary G$_1$ phase grating having period $p_1$, duty cycle $\alpha$, and amplitude $\theta$.

Download Full Size | PDF

 

Fig. 3. Diffraction-order amplitudes for a binary G$_1$ $\pi$-phase grating with 50% duty cycle. All even orders are zero.

Download Full Size | PDF

The diffraction efficiencies, on an intensity basis, are given by $|\mathbf {F}_m|^2$. Going forward, we will focus on a $\pi$-phase grating with 50% duty cycle ($\alpha = 0.5$), for which

3.2 Talbot fringe analysis

As the diffraction orders propagate away from G$_1$, they interfere with each other and produce grating images at integral multiples of the Talbot distance $Z_T$ (assuming we are in the region of spatial overlap behind the grating). For an amplitude grating illuminated by a plane wave $Z_T = p_1^2/\lambda$ ([37], Sec. 4.5.3). A similar phenomenon occurs for phase gratings, and in particular for a binary $\pi$-phase grating illuminated by a plane wave, frequency-doubled subimages arise at odd multiples of $Z_T/8$ [38,39]. If the grating is illuminated by a spherical wave emanating from a point source at distance $L$, the fundamental Talbot distance changes to $Z_{T}^{\prime } = L Z_T/(L - Z_T)$, thereby shifting the locations of the various self-images and scaling their periods by a factor of $Z_{T}^{\prime }/Z_T$ [39].

For the geometry of Fig. 1, we can calculate the Talbot intensity pattern for a monochromatic point source located at the origin (no G$_0$ present). The diffraction orders of G$_1$ can then be viewed, as shown in Fig. 4, as originating from a set of virtual point sources distributed along a circular arc with a radius of curvature equal to $L$. These virtual sources are weighted by their corresponding diffraction amplitudes, so the fields are given by

 

Fig. 4. Talbot fringe formation arising from the interference of overlapping G$_1$ diffraction orders. For a mononchromatic point source, the diffraction orders can be viewed as virtual point sources distributed along a circular arc of radius $L$. The amplitudes of the point sources are given by the Fourier coefficients of the grating as shown in Fig. 3.

Download Full Size | PDF

Here $r_m$ is the distance from the $m^{th}$ point source to an arbitrary observation point having coordinates $(x,y,z)$. From straightforward geometrical considerations, the virtual point source coordinates are

In practice, we also need to account for the spatial coherence of the interfering orders when a scattering object is placed in the beam, which can modify the observed Talbot fringe contrast when measured by pixels that are large compared to the fringe period. Accordingly, the apparent degree of coherence between any two diffraction orders depends on some measure of their relative spatial separation in the detection plane. More detail is provided in Section 4.4, but at this juncture it suffices to note that for our purposes here, the coherence can be cast in terms of a real-valued function, $-1 \leq \gamma (q) \leq 1$, where $q$ is related to the spatial frequency of an elemental fringe. While mathematically we can treat $\gamma (q)$ as a continuous function of $q$, our primary interest lies in a discrete version $\gamma _q$ indexed by the difference between order integers, $q = |m-n| = 0, 1, 2, \ldots 2M$. In vector form we have $\overline {\gamma } = [\gamma _0, \: \gamma _1, \: \ldots \gamma _{2M}]^{T}$. All pairs of diffraction orders having the same $q$ will interfere to produce elemental fringes of the same period. We can now use $\overline {\gamma }$ to form a symmetric Toeplitz coherence matrix,

It is instructive to manipulate Eq. (9) in order to isolate the various spatial-frequency terms comprising the Talbot pattern. We begin by noting that the elements of each diagonal of $\overline {I}_T(x,y)$ correspond to elemental fringes of the same period, all multiplied by the same coherence scaling value $\gamma _q$. Therefore, the objective is to recast Eq. (9) in terms of diagonal sums indexed by $q = |m - n|$ as follows

As a simple example, consider a discrete Gaussian coherence function of the form

 

Fig. 5. (a) Sample Gaussian coherence functions ($q_c$ = 20, 3.0, and 1.5). (b) Corresponding Talbot fringe patterns (normalized to unit average intensity), calculated using Eq. (9) with $M=25$, yielding visibility values of 1.0, 0.82, and 0.22. The Talbot fringes are based on a 28-keV system configuration as described in Sec. 2.

Download Full Size | PDF

4. Statistical optics model for the dark-field signal

To assess the impact of a scattering object on the observed fringe visibility, it is helpful to first visualize the system without any gratings as shown in Fig. 6. Here a phase-only volume scattering object is placed a distance $z_s$ from the point source. We employ the projection approximation and assume that the phase perturbation introduced by the object can be represented as a thin random phase screen having a 2D phase profile found from the geometric projection of accumulated phase along undeviated ray paths. Even though the object is illuminated by a spherical wave, we utilize a plane-wave projection by making two reasonable assumptions. First, the thickness of the medium, $T$, is taken to be much smaller than $z_s$. Second, for simulation purposes, we can restrict attention to a small paraxial field-of-view patch on the detector. It is instructive to examine these approximations in more detail, which we do next.

 

Fig. 6. Scattering of a spherical wave by a random phase object located at $z = z_s$. The thickness of the medium, $T$, is assumed to be much smaller than the distance from the source to the medium $z_s$ (i.e., $T << z_s$).

Download Full Size | PDF

Our analysis is applicable in the Fresnel regime, which is equivalent to a paraxial approximation [40]. In the limit of a single on-axis detection pixel it is reasonable to treat that portion of the beam propagating through the object, and illuminating only this one pixel, as a small quasi-collimated bundle of rays parallel to the z-axis (for $T << z_s$). This “beamlet” will then acquire a random phase profile in accordance with the plane-wave projection approximation; therefore, upon exiting the medium the local wavefront will be characterized by its phase statistics, namely a mean value, $\overline {\phi }$, a variance, $\sigma ^2_{\phi }$, and a normalized autocovariance function, $c_{\phi }(\Delta \alpha, \Delta \beta )$. Now consider a uniform scattering medium displaying homogeneous statistics (i.e., the medium is spatially stationary). For a paraxial off-axis pixel, the corresponding collimated ray bundle propagates at a small angle. The local wavefront for this off-axis beamlet will have a different mean phase by virtue of the slightly longer propagation path, particularly in the X-ray regime owing to the very short wavelength. However, the phase variance and autocovariance function are much less sensitive to slight off-axis propagation through the medium, so they will be similar to those of the wavefront received by the on-axis pixel. As we will see, it is the phase variance that dictates the scattering magnitude, and ultimately when combined with the autocovariance function will determine the Talbot fringe visibility. Therefore, with a uniform scattering object, all of the paraxial pixels surrounding the optical axis should have comparable visibility signals stemming from a plane-wave projection analysis.

We should mention that it may be possible to circumvent this approximation of a spherical wave as an angular collection of quasi-collimated beamlets. Doing so could remove the constraint that $T << z_s$. In the study of optical wave propagation through the turbulent atmosphere, the case of a point source embedded in the scattering medium has been analyzed ([41], and [42], Sec. 8.5.5). This earlier work includes a general formulation for the wavefront statistical properties that is not tied to the specific structure of atmospheric inhomogeneity, but instead only assumes the power spectral density of the refractive index distribution is known. Therefore, in principle, this analysis could be applied to scattering objects of interest in dark-field X-ray imaging. We only mention this possibility here as further investigation along these lines is outside the scope of this paper.

More detail related to the plane-wave projection approximation approach is provided in Appendix A of both Part 1 and Part 2 [49], with a focus on application to a random distribution of microspheres. Again, we leverage existing mathematical analysis for optical propagation through random turbulent media, including conditions under which the geometric projection approach, also referred to as a near-field analysis, is valid ([42], Ch. 8). Of particular note is the fact that the expression for the field autocorrelation function found using the projection approximation remains valid for much thicker media in which diffraction cannot be neglected ([42], Sec. 8.5.6).

The scattering object and its corresponding thin phase screen proxy can be located either before or after the G$_1$ grating. In either case, we are interested in understanding the impact of the object on the spatial coherence of the beam as this will determine the fringe contrast. Because the diffusing screen is a static structure placed in a spatially coherent beam, on a time-average basis the radiation passing through the screen will remain spatially coherent. However, as will become apparent, the Talbot fringe contrast averaged over a pixel dimension that is large relative to the fringe period depends on the ensemble-average coherence of the beam ([42], Sec. 5.9). Here the ensemble consists of macroscopically similar but microscopically different scattering objects. Therefore, our initial objective in this section is to calculate the ensemble-average mutual intensity for the beam leaving the phase screen and to determine how the mutual intensity changes upon propagation to the detection plane. This portion of the analysis largely follows that of Goodman (see Ref. [43], Ch. 5). We then subsequently take the G$_1$ splitter grating into account to find an expression for the measured fringe visibility as a function of the scattering object’s statistical phase properties as well as its placement within the interferometer.

Use of a monochromatic point source model may seem overly simplified because in an actual Talbot-Lau X-ray interferometer a polychromatic extended tube source is typically used. However, this more general configuration can be modeled as an incoherent superposition of the monochromatic point-source intensity results, albeit subject to the paraxial approximation. So for this reason the work presented here centers on the basic monochromatic point-source arrangement.

4.1 Scattered-field mutual intensity

Referring again to the propagation path illustrated in Fig. 6, the scattering object is placed with its output surface at the $z = z_s$ plane where the transverse coordinates are denoted $(\alpha, \beta )$. In accordance with the projection approximation, the scattering object is represented by a thin transmissive screen with an average amplitude transmission $\mathbf {t}_0$ and a spatially varying random phase $\phi (\alpha,\beta )$,

To assess the impact of the random screen on the wavefront, we construct the ensemble-average mutual intensity (i.e., the statistical autocorrelation function of the field) immediately behind the screen,

4.2 Propagation of the mutual intensity

Our next task is to propagate the mutual intensity of the radiation field from the scattering medium, $\mathbf {J}_{s}$, to the detection plane, $\mathbf {J}_{d}$. Afterwards we will turn our attention to determining a specific mathematical form for $\boldsymbol {\mu }_{s}(\Delta \alpha, \Delta \beta )$, but for now we treat it as a spatially localized 2D function centered on the origin. In the Fresnel regime, the mutual intensity propagation law is (see, e.g., [44] Eq. A1, or [45] Eq. 4.30)

As mentioned previously, spatial coherence effects come into play because diffraction from the G$_1$ grating causes portions of the wavefront to split and shift with respect to one another in the detection plane. However, this shift only occurs in one direction, say along the $x$-axis, so in the detection plane we have $\Delta y = 0$. As a consequence, we can restrict the analysis from 4D to 2D and use the following simplified version of the propagation equation,

For example, if we use the following numerical values: $\lambda = 4.4$ x $10^{-8}$ mm ($E = 28$ keV), $\ell = 10$ mm, $z_s = 650$ mm, and $z = L + d - z_s = 650$ mm, then $\Delta \alpha _0 = 1.4$ x $10^{-6}$ mm = 1.4 nm. This very small value of $\Delta \alpha _0$ implies that for all intents and purposes, the sinc function behaves as a delta function (at least relative to the $\boldsymbol {\mu }_{s}$ functions of interest to us here). By using the scaling law for delta functions, $\delta (ax) = \delta (x)/|a|$, Eq. (30) becomes

To more closely examine this scaling, we can explicitly write the autocorrelation function of the field in the detection plane as

4.3 Gaussian random-phase scattering object

We now further investigate the autocorrelation function of the scattered field as given by Eq. (20),

Given that $\Delta \phi$ is a zero-mean Gaussian random variable, it is well established that the scattered field autocorrelation function is given by,

To close this section, we note that the interferometer as described here is constructed with a one-dimensional G$_1$ grating, so the wavefront shearing between the various diffraction orders occurs only along one direction. Therefore, even if the scattering medium is anisotropic with a correlation length that varies with angle in the $\alpha \beta$-plane, we need only consider the direction in which the system is sensitive, which we take to be the $\alpha$-axis. Accordingly, going forward, we will drop the $\beta$ dependence and use the following 1D scattered field correlation function,

4.4 Talbot fringe visibility (dark-field signal)

We are now in a position to calculate the Talbot fringe and its detected visibility, or dark-field signal, as a function of the phase properties of the scattering object as well as the location of the object along the optical path. We continue to assume the source is a monochromatic point source, again recognizing that a more general polychromatic extended source can be modeled as a spectral and spatial superposition of uncorrelated point sources. Regarding the scattering object location, there are two distinct arrangements to consider as shown in Fig. 7, namely: (a) the object placed before the G$_1$ grating, and (b) the object located behind G$_1$. In Fig. 7 the diffraction angles are grossly exaggerated for illustration purposes. In reality, the angles are exceedingly small and the orders have significant spatial overlap in the detection plane.

 

Fig. 7. Illustration of the lateral shearing between diffraction orders of G$_1$ for the purposes of calculating the relative spatial coherence of the orders in the detection (Talbot) plane. There are two distinct cases: (a) sample before the G$_1$ grating, and (b) sample after G$_1$. The various relative shearing values ($\Delta x_{m,n}$) are given by Eq. (45). For simplicity, only orders 0-3 are shown here; however, for analysis all relative pairwise shearing separations for orders from $-M$ to $+M$ are utilized.

Download Full Size | PDF

In the first case of Fig. 7(a), with the scattering sample placed before G$_1$, the spherical source wavefront first illuminates the object, and the scattered wavefront (which in general contains a specular component) propagates to the grating. Then the diffraction orders, which are tilted replicas of the wavefront incident on the grating, propagate to the detection plane and arrive laterally sheared with respect to one another. The shearing occurs over a propagation distance of $d$, the spacing between the grating and the detection plane. So, for any sample position in front of G$_1$, the lateral shearing between orders remains constant. However, as just discussed in the previous section, the width of the field autocorrelation function (which is a measure of ensemble-average spatial coherence) depends on the total propagation distance from the sample to the detection plane, $z = L + d - z_s$. The greater the total propagation distance, the more the coherence of the scattered wavefront improves. Therefore, we see there are two separate effects coming into play that together determine the fringe visibility; one is the change in the field correlation function width with propagation distance from the scattering object to the detector, and the second is the lateral shearing in the detector plane between the various diffracted copies of the scattered wavefront (i.e., the separation values at which the correlation function is evaluated).

In the second case of Fig. 7(b), in which the sample is placed after the grating, the diffracted spherical wavefronts begin shearing as they leave the grating, but they are still spatially coherent with one another. Once the orders arrive at the front surface of the object they begin to scatter. If the object is sufficiently thin, then all of the orders will exit the medium with essentially the same phase perturbation, so as a group they lie within an isoplanatic angular patch. At this point shearing for ensemble-average coherence purposes begins as these scattered diffraction orders now start to laterally separate upon free-space propagation to the detection plane (depicted in Fig. 7(b)).

However, as the medium (still located after G$_1$) becomes thicker, the diffraction orders begin probing slightly different regions of the scattering volume, and as a consequence they start acquiring different random phase profiles and decorrelating. The sensitivity to thickness increases with propagation angle, so therefore the higher orders are the first to suffer this change and fall outside the isoplanatic patch. In this circumstance, the medium can be modeled as a stack of statistically independent phase plates, spaced far enough apart to ensure their associated phase profiles are uncorrelated (i.e., a spacing greater than the medium’s average particle size), yet close enough together to adequately approximate the continuous nature of the medium ([48], Sec. 17.2). For $N$ Gaussian random-phase plates, the correlation coefficient between diffraction orders $m$ and $n$, with propagation angles $\varphi _m$ and $\varphi _n$, is

As an example, assume the medium has phase fluctuations corresponding to $\sigma ^2_{\phi } = 2.0$ (total scattering fraction = $1 - \exp (-\sigma ^2_{\phi }) = 0.86$) that are distributed uniformly across $N = 10$ equally spaced layers, and each layer has the same Gaussian autocovariance function $c_{\phi }(\xi ) = \exp [-(\xi /R)^2]$ with a characteristic particle size of $R = 3\mu$m. Let the wavelength equal 4.4 x $10^{-5}\mu$m ($E = 28$ keV) and the G$_1$ grating period $p_1$ = 9.2 $\mu$m (these values correspond to our example simulation system in Part 2 [49]), from which we can determine the diffraction-order angles ($\varphi _m = m \lambda /p_1$ for the $m^\textrm {th}$ order). We can now inquire as to what object thickness causes the value of $\rho _q$ to drop to say 0.8 for the largest $q$ of interest and denote this as a critical thickness, $T_c$. This is the maximum thickness for which all diffraction orders of interest remain in an isoplanatic patch. For simulation purposes, we typically include orders up to $\pm$15, so $q_\textrm {max} = 30$. A simple numerical implementation of Eq. (44) reveals, for this particular example, $T_c = 12$ mm. For a thickness less than $T_c$, decorrelation of the diffraction orders upon propagation through the object can be neglected; otherwise, the $\rho _q$ terms should be included in the dark-field analysis. Even for thicker objects that require inclusion of the $\rho _q$ terms, we can still employ the projection approximation for the zeroth-order beam and treat the object as planar by ascribing the appropriate statistical properties that mimic the actual thicker medium.

After the diffraction orders leave the scattering object and propagate toward the detector, they will continue to decorrelate because of lateral shearing. Therefore, in the configuration with the object placed after G$_1$, as the effective scattering plane moves closer to the detection plane the relative ensemble-average coherence between the orders improves by virtue of the smaller degree of lateral shearing. On the flip side, a shorter gap means less coherence restoration via propagation, but regardless we find that the smaller lateral shearing dominates and thus the fringe visibility improves as the scattering plane moves closer to the detection plane.

With these two configurations in mind, we designate $\Delta x_{m,n}$ to be the ensemble-average-coherence lateral shearing between the scattered wavefronts of diffraction orders $m$ and $n$ in the Talbot plane,

For a $\pi$-phase G$_1$ grating, the $\pm 1$ orders are dominant at the design energy. If we neglect the other orders, then an approximate expression for the Talbot fringe is found by keeping only the $q = 2$ term in Eq. (48), so that

The signal $S(x)$ recorded by a detector pixel is found from a convolution of the local Talbot X-ray fringe with the intensity transmission function of the G$_2$ analyzer grating. The corresponding detected visibility equals $(S_\textrm {max}-S_\textrm {min})/(S_\textrm {max}+S_\textrm {min})$. As is customary, we calculate the normalized detected visibility, $V_d$, given by a ratio of the detected visibilities with and without the scattering object in the beam path,

5. Conclusion

Dark-field X-ray imaging with a Talbot-Lau interferometer has the ability to detect and image small-angle scattering from a sample with high sensitivity. As its importance for a variety of applications continues to grow, further theoretical understanding of the dark-field effect will help advance this relatively new technology. In Part 1 of this paper we utilize the basic properties of grating diffraction, combined with the concept of ensemble-average coherence as it relates to a scattering object, in order to develop a statistical optics model for the dark-field signal. A convenient and analytically tractable model for the scattering object is based on a random collection of dielectric microspheres (although more general refractive index distributions, including anisotropic distributions, can be handled numerically). This approach allows us to calculate the detected fringe visibility for an object having arbitrary particle concentration placed at any location within the interferometer. The model is based on fundamental optical principles and therefore helps provide insight into the underlying physical phenomena. In Part 2 [49] of the paper we present the results of several example monochromatic simulations based on a 28-keV system configuration.

In formulating our model, we make use of the projection approximation in which the 2D phase profile for a beam passing through the object is found from the accumulated phase along straight ray paths. The autocovariance function of this random phase profile, along with the related autocorrelation function of the radiation field, figure prominently in the analysis. Therefore, important details related to the mathematical derivations of the phase covariance functions for both an individual sphere and for a random collection of microspheres are provided here in Appendix A of Part 1. While it might seem that the projection approximation constrains the scattering object to be very thin, we provide detailed arguments in Part 2 [49] (Appendix A) as to why the results of this approach can actually extend to much thicker lossless media. In this two-part paper we focus solely on transparent media with phase-only variations. However, the projection formulation in Appendix A can likely be modified to include particles that have both refraction and attenuation, thereby allowing for a more general polychromatic simulation that includes beam hardening, albeit potentially restricted to a thin object now that loss is involved.

Lastly, in Supplement 1, we: (a) explore the convolutional blurring model of the dark-field effect as found in the literature and propose a revised version that is consistent with the statistical optics model, and (b) provide a derivation of the 3D autocorrelation function for a random collection of hard spheres. In summary, this paper provides a comprehensive theoretical formulation of the dark-field effect from a statistical optics perspective, allowing for a quantitative analysis of the dark-field signal that is versatile, amenable to efficient numerical implementation, and helps provide additional insight into the relevant basic physical phenomena.

Appendix A. Computation of projected phase autocovariance functions

Here we review the projection approximation and consider its use in calculating the 2D autocorrelation and autocovariance functions of the phase distribution imparted to the beam by a volume scattering object.

A.1. Projection approximation

The basic idea of the projection approximation is illustrated in Fig. 8. Parallel rays propagate through an inhomogeneous scattering medium having a random spatial distribution of refractive index, for example a suspension of small randomly located microspheres. The scale sizes associated with the inhomogeneities are assumed to be much larger than the wavelength (similar to the assumption typically made for turbulent media, in contrast to turbid media in which the scale sizes are comparable to or smaller than the wavelength). In this approximation the rays incur no bending, but each ray does acquire an accumulated phase delay along its path. Therefore, this approach comprises a purely geometrical model with diffraction being neglected during propagation within the scattering medium itself.

 

Fig. 8. Illustration of the plane-wave projection approximation as applied to a random distribution of microspheres. Rays propagating parallel to the $z$-axis yield a 2D projected phase distribution at the output plane. We generally neglect the uniform background phase ($\phi _0 = kn_0T$, where $k$ is the free-space wavenumber) and retain only the spatial phase modulation $\phi (x,y)$ at $z = T$.

Download Full Size | PDF

We assume that the scattering medium is statistically isotropic, described by a 3D autocorrelation function of the refractive index variations, $\Gamma _{\Delta n}(r)$, that has spherical symmetry. For simplicity, as well as to provide a specific example, we will restrict our subsequent analysis to a volume scatterer comprising a random ensemble of microspheres having a fixed radius. However, before addressing the problem of an ensemble of microspheres, we first consider the projection approximation applied to an individual spherical particle, with the result being representative of an ensemble in the limit of low concentration (where the spheres are non-interacting and randomly positioned with free arrangement [52]).

A.2. Phase autocorrelation function for a single spherical particle

Assume beam propagation along the $z$ axis, which becomes the direction along which we seek to evaluate the projection. Therefore, our objective is to find the 2D autocorrelation function of the wavefront phase projected onto the output (x,y) plane as shown in cross section in Fig. 9.

 

Fig. 9. Projection of a sphere having radius $R$.

Download Full Size | PDF

For this exercise we take the sphere of radius $R$ to be represented by its so-called indicator function, $m(r)$, expressed in terms a radial Heaviside step function,

(57)$$m(r) = H(R-r) = \left\{ \begin{matrix} 1, \quad & r \leq R \\ \, 0, & r \;>\; R \end{matrix} \right.$$

with $r = \sqrt {x^2 + y^2 + z^2}$. In other words, the indicator function assumes a value of one inside the sphere and zero outside. With this representation, the 3D autocorrelation function of the sphere (with volume $V_s$) is given by ([52], Eq. 10),

(58)$$A(r) = \left\{\begin{matrix} V_s \left[ 1 - \frac{3}{4}\frac{r}{R} + \frac{1}{16}\left( \frac{r}{R} \right)^3 \right] & \textrm{for } \: r \leq 2R, \\ 0 & \textrm{otherwise}. \end{matrix} \right.$$

Multiplying this result by $\Delta n^2$ yields the index-contrast correlation function, $\Gamma _{\Delta n}(r) = \Delta n^2 A(r)$. Following the approach outlined by Krouglov et al. in Ref. [53] (Eq. 21), we can write the corresponding 2D autocorrelation function of the projected phase as

(59)$$\Gamma_{\phi}(r_\perp) = \left\{\begin{matrix} 2 (k \Delta n)^2 \int_{r_{{\perp}}}^{2R} \frac{ A(r) \, r}{(r^2 - r_{{\perp}}^2)^{1/2}} \, dr & \textrm{for } \: |r_{{\perp}}| \leq 2R, \\ \, 0 & \textrm{otherwise}. \end{matrix} \right.$$

Substituting $A(r)$ from Eq. (58) into Eq. (59) and carrying out the integration yields the following normalized version of the phase correlation function ([53], Eq. 32),

(60)$$\begin{aligned} \gamma_{\phi}(r_\perp) = \frac{\Gamma_{\phi}(r_\perp)}{\Gamma_{\phi}(0)} & = \left[ 1 - \frac{1}{4} \left( \frac{r_{{\perp}}}{R} \right)^2 \right]^{1/2} \left[ 1 + \frac{1}{8} \left( \frac{r_{{\perp}}}{R} \right)^2 \right]\\ & \hspace{-3 mm} + \frac{1}{2} \left( \frac{r_{{\perp}}}{R} \right)^2 \left[ 1 - \frac{1}{16} \left( \frac{r_{{\perp}}}{R} \right)^2 \right] \ln \left[ \frac{|r_{{\perp}}|/R}{2 + \sqrt{4 - (r_{{\perp}}/R)^2}} \right] \hspace{2 mm}\textrm{for } \: |r_{{\perp}}| \leq 2R, \end{aligned}$$

where $\Gamma _{\phi }(0) = 3R(k \Delta n)^2 V_s/2$. This result has been used by various researchers for dark-field fringe visibility analysis of dilute microsphere suspensions [12,27–30]. While Eq. (60) is the exact analytical solution for the spherical-particle phase autocorrelation function, a simple Gaussian function provides a reasonable approximation:

(61)$$\gamma_{\phi}(r_\perp) \simeq \exp \left[ -\left( \frac{r_{{\perp}}}{R} \right)^2 \right].$$

The exact and approximate versions for $\gamma _{\phi }$ are plotted in Fig. 10. We see that for increasing radial displacement the correlation function falls to zero (which implies the mean phase associated with the random index distribution is negligible, separate from the constant background phase shift), so for purposes of calculating the field correlation function of Eq. (42) we can equate the phase covariance function $c_{\phi }(\Delta \alpha )$ to the phase correlation function $\gamma _{\phi }(\Delta \alpha )$. While these results are useful when analyzing dilute monodisperse spheres for which the volume fraction $f \ll 1$, we next turn our attention to the more general case of arbitrary $f$ (up to a maximum value of $\sim 0.64$ based on the limit for random hard-sphere packing). Before doing so, we close this section by mentioning that Eq. (60), and its approximation of Eq. (61), are also applicable to random distributions of penetrable spheres (i.e., spheres that can physically overlap with one another) over a much larger range of volume fractions (see [54], Fig. 6).

 

Fig. 10. Exact (Eq. (60)) and approximate (Eq. (61)) forms of the normalized phase autocorrelation function for a homogeneous spherical particle after projection onto a plane.

Download Full Size | PDF

A.3. Phase autocovariance function for an ensemble of spherical particles

For a dilute random distribution of microspheres, the 3D and projected 2D correlation functions for the ensemble are well approximated by their geometrical counterparts for a single microsphere. However, as the concentration of hard microspheres goes up, the 3D ensemble autocorrelation function, which we denote as $\gamma _e(r)$, becomes increasingly dependent on the volume fraction of microspheres. This correlation function for monodisperse hard spheres (i.e., spheres of the same radius) having a volume fraction $f$ can be found using the approach outlined by Torquato and Stell [33]. The details are provided in Supplement 1; here we simply state the result:

(62)$$\Gamma_e(r) = f \frac{A(r)}{V_s} + f^2 \left[ 1 + \frac{h(r) \ast m(r) \ast m(r)}{V_s^2} \right],$$

where $A(r)$ is the autocorrelation function of a single sphere (Eq. (58)), $h(r) = g_2(r)-1$ where $g_2(r)$ is the radial distribution function (see [55], Fig. 3.2), $m(r)$ is the indicator function (Eq. (57)), and the asterisks represent 3D convolutions. For computational purposes it is convenient to evaluate the convolution term as a product in the Fourier domain followed by an inverse transformation. Also, it can be shown that $\Gamma _e(0) = f$, so that we arrive at the following final expression for the normalized ensemble correlation function,

(63)$$\gamma_e(r) = \frac{\Gamma_e(r)}{f} = \frac{A(r)}{V_s} + f \left[ 1 + \frac{1}{V_s^2} \mathcal{F}^{{-}1}_{3D} \left\{ \tilde{h}(s) \widetilde{m}^2(s) \right\} \right].$$

Here $\tilde {h}(s)$ and $\widetilde {m}(s)$ are 3D Fourier transforms of $h(r)$ and $m(r)$, respectively; both of these transforms have reasonably straightforward analytical expressions. Therefore, evaluation of the inverse transform term above is readily amenable to numerical computation (via a simple sum) when using a judicious choice of sampling in both the spatial and spatial-frequency domains as described by Lado [56]. From Eq. (63) we see that $\gamma _e(r) \simeq A(r)/V_s$ when $f$ is small, on the order of few percent. Example plots of $\gamma _e(r)$ are provided in Supplement 1 (Fig. S4).

Our primary goal now is to use $\gamma _e(r)$ to determine an expression for the corresponding projected phase correlation function. In this case, we take a statistical approach, similar to that used to describe propagation through the turbulent atmosphere (e.g., see [42], Sec. 8.5.2). Figure 11 illustrates the geometry for the analysis. As mentioned previously, we restrict our analysis to a volume scatterer comprising a random ensemble of microspheres with fixed radius $R$. The refractive index distribution is therefore an isotropic random process described by the following 3D autocorrelation,

(64)$$\Gamma_{\Delta n}(r) = \Delta n^2 \Gamma_e(r) = \Delta n^2 f \gamma_e(r),$$

where $\gamma _e(r)$ is given by Eq. (63) and $r = |\vec {r}|$ represents a radial displacement in 3D space. We further assume that the scattering medium is spatially ergodic, so the statistical and spatial autocorrelation functions are equivalent,

(65)$$\Gamma_{\Delta n}(r) = \overline{\Delta n(\vec{u}) \, \Delta n(\vec{u} + \vec{r})} = \frac{1}{V} \iiint_{V} \Delta n(\vec{u}) \Delta n(\vec{u} + \vec{r}) \, d\vec{u}.$$

 

Fig. 11. Geometry for computing the 2D projected phase autocorrelation function of a random scattering medium characterized by a spherically-symmetric 3D autocorrelation function of the refractive index perturbation, $\Gamma _{\Delta n}(r)$.

Download Full Size | PDF

Referring again to Fig. 11, consider the phase delays incurred by the two dashed rays (neglecting the overall constant phase delay associated with the background index $n_0$),

(66)$$\begin{aligned} \phi_1 (\vec{r}_{{\perp} 1}) & = k \int_0^T \Delta n (\vec{r}_1) dz \\ \phi_2 (\vec{r}_{{\perp} 2}) & = k \int_0^T \Delta n (\vec{r}_2) dz, \end{aligned}$$

where $k=2\pi /\lambda$ is the free-space wavenumber, with $\lambda$ being the wavelength in vacuum. Let the two position vectors be expressed in terms of their transverse and longitudinal coordinates, $\vec {r}_1 = (\vec {r}_{\perp 1}, z_1)$ and $\vec {r}_2 = (\vec {r}_{\perp 2}, z_2)$. Doing so lets us construct the projected phase autocorrelation function as,

(67)$$\begin{aligned}\Gamma_{\phi}(\vec{r}_{{\perp} 1}, \vec{r}_{{\perp} 2}) & = \overline{\phi_1 (\vec{r}_{{\perp} 1}) \, \phi_2 (\vec{r}_{{\perp} 2})}\\ & = k^2 \int_0^T \int_0^T \overline{ \Delta n(\vec{r}_{{\perp} 1}, z_1) \, \Delta n(\vec{r}_{{\perp} 2}, z_2) } \, dz_1 dz_2\\ & = k^2 \int_0^T \int_0^T \Gamma_{\Delta n} \! \left( \sqrt{ \left| \vec{r}_{{\perp} 1} - \vec{r}_{{\perp} 2} \right|^2 + \left( z_1 - z_2 \right)^2} \, \right) dz_1 dz_2. \end{aligned}$$

The last step assumes the index fluctuations are isotropic, so the index contrast autocorrelation function $\Gamma _{\Delta n}(r)$ has spherical symmetry and depends only on the magnitude $r = |\vec {r}_1 - \vec {r}_2|$ of its displacement argument. We see $\Gamma _{\Delta n}$ is now an even function of $z_1 - z_2$, allowing the double integral to be reduced to a single integral ([42], Eqs. 3.7-26 and 8.5-33),

(68)$$\Gamma_{\phi}(r_{{\perp}}) = 2k^2 \int_0^T (T - z) \, \Gamma_{\Delta n} \! \left( \sqrt{r_{{\perp}}^2 + z^2} \, \right) dz,$$

where $r_{\perp } = |\vec {r}_{\perp 1} - \vec {r}_{\perp 2}|$ and $z = z_1 -z_2$. This is an important basic result; it expresses the 2D projected autocorrelation function of the wavefront phase, $\Gamma _{\phi }(r_{\perp })$, as an integral (along the propagation direction) of the 3D autocorrelation function of the underlying refractive index variation, $\Gamma _{\Delta n}(r)$ with $r = (r_{\perp }^2 + z^2)^{1/2}$.

By substituting Eq. (64) into Eq. (68), the expression for the projected phase autocorrelation function becomes

(69)$$\Gamma_{\phi}(r_{{\perp}}) = 2f (k \Delta n)^2 \int_0^T (T - z) \, \gamma_{e} \! \left( \sqrt{r_{{\perp}}^2 + z^2} \, \right) dz.$$

We can express the radial displacement arguments of the correlation functions in units of particle radius, so $\Gamma _{\phi }(r_{\perp }) \rightarrow \widehat {\Gamma }_{\phi }(r_{\perp }/R)$ and $\gamma _e(r) \rightarrow \widehat {\gamma }_e(r/R)$. In this way, the correlation functions accept normalized displacements as input. Making this change allows us to write

(70)$$\widehat{\Gamma}_{\phi}(r_{{\perp}}/R) = 2f (k \Delta n)^2 \int_0^T (T - z) \, \widehat{\gamma}_{e} \! \left( \sqrt{r_{{\perp}}^2 + z^2} \, /R \, \right) dz.$$

For numerical evaluation purposes, it is convenient to recast this result into a slightly different form by using a coordinate change $z = uR$, leading to

(71)$$\widehat{\Gamma}_{\phi}(r_{{\perp}}/R) = 2f (kR \Delta n)^2 \int_0^{T/R} \left( \frac{T}{R} - u \right) \, \widehat{\gamma}_{e} \! \left( \sqrt{(r_{{\perp}}/R)^2 + u^2} \, \right) du.$$

So for any desired value of normalized displacement, $r_{\perp }/R$, the integral now evaluates to a dimensionless value for a given ratio of sample thickness to particle radius, $T/R$.

In Fig. 12(a) we use Eq. (71) in combination with Eq. (63) to calculate and plot the normalized version of the phase autocorrelation function,

(72)$$\widehat{\gamma}_{\phi}(r_{{\perp}}/R) = \frac{\widehat{\Gamma}_{\phi}(r_{{\perp}}/R)}{ \widehat{\Gamma}_{\phi}(0)},$$

for a scattering region thickness of $T = 10R$. Even for such a relatively thin scatterer, a noticeable asymptote exits for large $r_{\perp }$, with its value increasing as $f$ becomes larger. This asymptote arises because the projected phase is not a zero-mean random variable. As the quantity of scattering material that is being integrated in the projection process increases, so does the average projected phase. Recall the particle index differs from the background index by a fixed amount $\Delta n$, which can be either positive or negative, meaning the phase projection is a random variable biased by the same sign as $\Delta n$. Note this mean value for the phase projection is over and above the contribution from the uniform background refractive index (which is already neglected in the analysis). We can easily quantify this effect by noting from Eq. (63) that the large-$r$ limit of $\gamma _e(r)$ equals $f$ (the first and third terms in Eq. (63) go to zero); therefore, in this limit Eq. (69) yields

(73)$$\left( \overline{\phi} \right)^2 = \lim_{r_{{\perp}} \rightarrow \infty} \Gamma_{\phi}(r_{{\perp}}) = 2f^2 (k \Delta n)^2 \int_0^T (T - z) \, dz = (k \Delta n f T)^2.$$

Just as we neglected the uniform phase accumulation due to the background index, we can likewise neglect the average projected phase because it only imparts an additional constant phase offset to the beam and does not impact scattering. So it is more convenient to instead use the phase autocovariance function,

(74)$$C_{\phi}(r_{{\perp}}) = \Gamma_{\phi}(r_{{\perp}}) - \left( \overline{\phi} \right)^2,$$

which in normalized form is $c_{\phi }(r_{\perp }) = C_{\phi }(r_{\perp }) / C_{\phi }(0)$. For numerical evaluation of this function based on Eq. (71), we let $c_{\phi }(r_{\perp }) \rightarrow \widehat {c}_{\phi }(r_{\perp }/R)$. Figure 12(b) shows plots of $\widehat {c}_{\phi }(r_{\perp }/R)$, giving a better view of the way in which the phase correlations change with particle concentration.

 

Fig. 12. (a) Normalized phase autocorrelation functions for a collection of monodisperse microspheres with a sphere radius $R$, and (b) corresponding normalized autocovariance functions. The sample thickness is $T = 10 R$.

Download Full Size | PDF

It is also of interest to determine how the variance of the projected phase, $\sigma _{\phi }^2$, depends on both scattering thickness and particle concentration. This is an important parameter that governs the degree of scattering. In general we know

(75)$$\Gamma_{\phi}(0) = \overline{\phi^2} = \left( \overline{\phi} \right)^2 + \sigma_{\phi}^2 \quad \Rightarrow \quad \sigma_{\phi}^2 = \Gamma_{\phi}(0) - \left( \overline{\phi} \right)^2.$$

So, by using Eq. (69) for $\Gamma _{\phi }(0)$, combined with Eq. (73) for $(\overline {\phi })^2$, the phase variance becomes

(76)$$\sigma_{\phi}^2 = 2(k \Delta n f)^2 \int_0^T (T - z) \left(\frac{\gamma_{e}(z)}{f} - 1 \right) dz.$$

Once again, letting $\gamma _e(z) \rightarrow \widehat {\gamma }_{e}(z/R)$ and then invoking the coordinate change $z = uR$, we have

(77)$$\sigma_{\phi}^2 = 2(kR \Delta n f)^2 \int_0^{T/R} \left( \frac{T}{R} -u \right) \left(\frac{\widehat{\gamma}_{e}(u)}{f} - 1 \right) du.$$

Upon inspection of the integrand we note that the term $(\widehat {\gamma }_{e}(u)/f - 1)$ quickly approaches zero with increasing $u$, so as a practical matter the upper limit of the integral can be terminated at $T/R \simeq 10$, which is helpful when evaluating this integral for the typical case of an object having large $T/R$. Moreover, for $T/R \gg 1$ in Eq. (77), we have $(T/R - u) \simeq T/R$, and thus

(78)$$\sigma_{\phi}^2 \simeq 2(kR \Delta n f)^2 \, \frac{T}{fR} \int_0^{T/R} \left(\widehat{\gamma}_{e}(u) - f \right) du.$$

If we neglect interaction between the microspheres (valid for small $f$), we can use a modified version of the autocorrelation function for a single sphere (Eq. (58)), with $A(r) \rightarrow \widehat {A}(r/R)$, to write $\widehat {\gamma }_{e}(u) - f \simeq (1-f) \widehat {A}(u)/V_s$, so

(79)$$\begin{aligned}\sigma_{\phi}^2 & \simeq 2(kR \Delta n f)^2 \, \frac{T}{R} \frac{(1-f)}{f} \int_0^2 \frac{\widehat{A}(u)}{V_s} du\\ & = 2f (1-f)TR(k \Delta n)^2 \int_0^2 \left( 1 - \frac{3}{4}u + \frac{1}{16}u^3 \right) du\\ & = 1.5 f(1-f) TR (k \Delta n)^2, \end{aligned}$$

which provides a simple expression for an estimate of the phase variance.

Funding

Transportation Security Administration (HSTS04-17-C-CT7224).

Acknowledgments

We would like to gratefully acknowledge review of a preliminary manuscript by Prof. Joseph W. Goodman, as well as many helpful discussions with members of the research team at Stanford including Max Yuen, Bill Aitkenhead, Yao-Te Cheng, Paul Hansen, Ludwig Galambos, Niharika Gupta, Ching-Wei Chang, and George Herring. In addition, we thank the anonymous reviewers for their thoughtful insight and suggestions that helped improve the paper.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). [CrossRef]  

2. A. Momose, “Recent advances in x-ray phase imaging,” Jpn. J. Appl. Phys. 44(9A), 6355–6367 (2005). [CrossRef]  

3. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]  

4. T. Donath, F. Pfeiffer, O. Bunk, C. Grünzweig, E. Hempel, S. Popescu, P. Vock, and C. David, “Toward clinical x-ray phase-contrast CT: Demonstration of enhanced soft-tissue contrast in human specimen,” Invest. Radiol. 45(7), 445–452 (2010). [CrossRef]  

5. M. Willner, J. Herzen, S. Grandl, S. Auweter, D. Mayr, A. Hipp, M. Chabior, A. Sarapata, K. Achterhold, I. Zanette, T. Weitkamp, A. Sztrókay, K. Hellerhoff, M. Reiser, and F. Pfeiffer, “Quantitative breast tissue characterization using grating-based x-ray phase-contrast imaging,” Phys. Med. Biol. 59(7), 1557–1571 (2014). [CrossRef]  

6. E. A. Miller, T. A. White, B. S. McDonald, and A. Seifert, “Phase contrast x-ray imaging signatures for security applications,” IEEE Trans. Nucl. Sci. 60(1), 416–422 (2013). [CrossRef]  

7. J. Herzen, T. Donath, F. Pfeiffer, O. Bunk, C. Padeste, F. Beckmann, A. Schreyer, and C. David, “Quantitative phase-contrast tomography of a liquid phantom using a conventional x-ray tube source,” Opt. Express 17(12), 10010–10018 (2009). [CrossRef]  

8. S. Zabler, ed., Phase-Contrast and Dark-Field Imaging (J. Imaging, 2018). Special Issue.

9. J. P. Wilde and L. Hesselink, “Modeling of an x-ray grating-based imaging interferometer using ray tracing,” Opt. Express 28(17), 24657–24681 (2020). [CrossRef]  

10. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, C. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7(2), 134–137 (2008). [CrossRef]  

11. F. Pfeiffer, M. Bech, O. Bunk, T. Donath, B. Henrich, P. Kraft, and C. David, “X-ray dark-field and phase-contrast imaging using a grating interferometer,” J. Appl. Phys. 105(10), 102006 (2009). [CrossRef]  

12. V. Revol, I. Jerjen, C. Kottler, P. Schütz, R. Kaufmann, T. Lüthi, U. Sennhauser, U. Straumann, and C. Urban, “Sub-pixel porosity revealed by x-ray scatter dark field imaging,” J. Appl. Phys. 110(4), 044912 (2011). [CrossRef]  

13. A. Velroyen, M. Bech, A. Malecki, A. Tapfer, A. Yaroshenko, M. Ingrisch, C. C. Cyran, S. D. Auweter, K. Nikolaou, M. Reiser, and F. Pfeiffer, “Microbubbles as a scattering contrast agent for grating-based x-ray dark-field imaging,” Phys. Med. Biol. 58(4), N37–N46 (2013). [CrossRef]  

14. C. Hannesschläger, V. Revol, B. Plank, D. Salaberger, and J. Kastner, “Fibre structure characterisation of injection moulded short fibre-reinforced polymers by x-ray scatter dark field tomography,” Case Stud. Nondestruct. Test. Eval. 3, 34–41 (2015). [CrossRef]  

15. S. Senck, M. Scheerer, V. Revol, B. Plank, C. Hannesschläger, C. Gusenbauer, and J. Kastner, “Microcrack characterization in loaded cfrp laminates using quantitative two- and three-dimensional x-ray dark-field imaging,” Compos. Part A: Appl. Sci. Manuf. 115, 206–214 (2018). [CrossRef]  

16. M. S. Nielsen, K. B. Damkjær, and R. Feidenhans’l, “Quantitative in-situ monitoring of germinating barley seeds using x-ray dark-field radiography,” J. Food Eng. 198, 98–104 (2017). [CrossRef]  

17. F. Yang, F. Prade, M. Griffa, R. Kaufmann, J. Herzen, F. Pfeiffer, and P. Lura, “X-ray dark-field contrast imaging of water transport during hydration and drying of early-age cement-based materials,” Mater. Charact. 142, 560–576 (2018). [CrossRef]  

18. W. Yashiro, P. Vagovič, and A. Momose, “Effect of beam hardening on a visibility-contrast image obtained by X-ray grating interferometry,” Opt. Express 23(18), 23462–23471 (2015). [CrossRef]  

19. A. Yaroshenko, K. Hellbach, M. Bech, S. Grandl, M. F. Reiser, F. Pfeiffer, and F. G. Meinel, “Grating-based x-ray dark-field imaging: a new paradigm in radiography,” Curr. Radiol. Reports 2(7), 57 (2014). [CrossRef]  

20. A. Velroyen, A. Yaroshenko, D. Hahn, A. Fehringer, A. Tapfer, M. Müller, P. Noël, B. Pauwels, A. Sasov, A. Yildirim, O. Eickelberg, K. Hellbach, S. Auweter, F. Meinel, M. Reiser, M. Bech, and F. Pfeiffer, “Grating-based X-ray dark-field computed tomography of living mice,” EBioMedicine 2(10), 1500–1506 (2015). [CrossRef]  

21. V. Ludwig, M. Seifert, C. Hauke, K. Hellbach, F. Horn, G. Pelzer, M. Radicke, J. Rieger, S.-M. Sutter, T. Michel, and G. Anton, “Exploration of different x-ray Talbot–Lau setups for dark-field lung imaging examined in a porcine lung,” Phys. Med. Biol. 64(6), 065013 (2019). [CrossRef]  

22. J. Vignero, N. Marshall, G. Vande Velde, K. Bliznakova, and H. Bosmans, “Translation from murine to human lung imaging using x-ray dark field radiography: A simulation study,” PLoS One 13(10), e0206302 (2018). [CrossRef]  

23. K. Willer, A. Fingerle, W. Noichl, F. De Marco, M. Frank, T. Urban, R. Schick, A. Gustschin, B. Gleich, J. Herzen, T. Koehler, A. Yaroshenko, T. Pralow, G. Zimmermann, B. Renger, A. Sauter, D. Pfeiffer, M. Makowski, E. Rummeny, P. Grenier, and F. Pfeiffer, “X-ray dark-field chest imaging can detect and quantify emphysema in COPD patients,” pre-publication in medRxiv (2021).

24. K. Taphorn, K. Mechlem, T. Sellerer, F. De Marco, M. Viermetz, F. Pfeiffer, D. Pfeiffer, and J. Herzen, “Direct differentiation of pathological changes in the human lung parenchyma with grating-based spectral X-ray dark-field radiography,” IEEE Transactions on Medical Imaging (2021). Accepted for publication.

25. Z.-T. Wang, K.-J. Kang, Z.-F. Huang, and Z.-Q. Chen, “Quantitative grating-based x-ray dark-field computed tomography,” Appl. Phys. Lett. 95(9), 094105 (2009). [CrossRef]  

26. M. Bech, O. Bunk, T. Donath, R. Feidenhans, C. David, and F. Pfeiffer, “Quantitative x-ray dark-field computed tomography,” Phys. Med. Biol. 55(18), 5529–5539 (2010). [CrossRef]  

27. M. Strobl, “General solution for quantitative dark-field contrast imaging with grating interferometers,” Sci. Rep. 4(1), 7243 (2015). [CrossRef]  

28. W. Yashiro, Y. Terui, K. Kawabata, and A. Momose, “On the origin of visibility contrast in x-ray Talbot interferometry,” Opt. Express 18(16), 16890–16901 (2010). [CrossRef]  

29. S. K. Lynch, V. Pai, J. Auxier, A. F. Stein, E. E. Bennett, C. K. Kemble, X. Xiao, W.-K. Lee, N. Y. Morgan, and H. H. Wen, “Interpretation of dark-field contrast and particle-size selectivity in grating interferometers,” Appl. Opt. 50(22), 4310–4319 (2011). [CrossRef]  

30. A. Malecki, G. Potdevin, and F. Pfeiffer, “Quantitative wave-optical numerical analysis of the dark-field signal in grating-based x-ray interferometry,” Europhys. Lett. 99(4), 48001 (2012). [CrossRef]  

31. L. Feigin and D. Svergun, Structure Analysis by Small-Angle X-Ray and Neutron Scattering (Springer, 1987).

32. S. Gkoumas, P. Villanueva-Perez, Z. Wang, L. Romano, M. Abis, and M. Stampanoni, “A generalized quantitative interpretation of dark-field contrast for highly concentrated microsphere suspensions,” Sci. Rep. 6(1), 35259 (2016). [CrossRef]  

33. S. Torquato and G. Stell, “Microstructure of two-phase random media. V. The n-point matrix probability functions for impenetrable spheres,” J. Chem. Phys. 82(2), 980–987 (1985). [CrossRef]  

34. T. Donath, M. Chabior, F. Pfeiffer, O. Bunk, E. Reznikova, J. Mohr, E. Hempel, S. Popescu, M. Hoheisel, M. Schuster, J. Baumann, and C. David, “Inverse geometry for grating-based x-ray phase-contrast imaging,” J. Appl. Phys. 106(5), 054703 (2009). [CrossRef]  

35. R. N. Bracewell, The Fourier Transform and Its Applications (McGraw-Hill, 2000), 3rd ed.

36. D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather, Diffractive Optics: Design, Fabrication, and Test (SPIE, 2003), vol. TT62, chap. 5, pp. 84–87.

37. J. W. Goodman, Introduction to Fourier Optics (W. H. Freeman & Company, 2017), 4th ed.

38. T. J. Suleski, “Generation of Lohmann images from binary-phase Talbot array illuminators,” Appl. Opt. 36(20), 4686–4691 (1997). [CrossRef]  

39. T. Weitkamp, C. David, C. Kottler, O. Bunk, and F. Pfeiffer, “Tomography with grating interferometers at low-brilliance sources,” in Developments in X-Ray Tomography V, vol. 6318U. Bonse, ed., International Society for Optics and Photonics (SPIE, 2006), pp. 249–258.

40. The Fresnel approximation is valid for small propagation angles θ < (4λ/z)¹/4 (see, e.g., [37] Eq. 4-23), which is equivalent to a paraxial approximation. In our case, with λ = 4.4e-05 µm and z ∼ 1 m, the angular spread of a cone beam over which our analysis is applicable should be less than about 0.4 degrees. If needed, a larger field of view can be constructed from a collection of conical pencil beams, each pointed in a different nominal direction but still satisfying the local paraxial condition with respect to this direction. In this case, the system benefits from use of curved or segmented gratings with a center of curvature located at the source.

41. D. L. Fried, “Limiting resolution looking down through the atmosphere,” J. Opt. Soc. Am. 56(10), 1380–1384 (1966). [CrossRef]  

42. J. W. Goodman, Statistical Optics (John Wiley & Sons, 2015), 2nd ed.

43. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (SPIE, 2020), 2nd ed.

44. B. McMorran and A. D. Cronin, “Model for partial coherence and wavefront curvature in grating interferometers,” Phys. Rev. A 78(1), 013601 (2008). [CrossRef]  

45. J. Mertz, Introduction to Optical Microscopy (Roberts & Company, 2010).

46. E. Wolf, Introduction to the Theory of Coherence and Polarization of Light (Cambridge, 2007).

47. W. Yashiro, S. Harasse, K. Kawabata, H. Kuwabara, T. Yamazaki, and A. Momose, “Distribution of unresolvable anisotropic microstructures revealed in visibility-contrast images using x-ray Talbot interferometry,” Phys. Rev. B 84(9), 094106 (2011). [CrossRef]  

48. T. S. McKechnie, General Theory of Light Propagation and Imaging Through the Atmosphere (Springer, 2016).

49. J. P. Wilde and L. Hesselink, “Statistical optics modeling of dark-field scattering in X-ray grating interferometers: Part 2. Simulation,” Opt. Express 29(25), 40917–40993 (2021). [CrossRef]  

50. F. Prade, A. Yaroshenko, J. Herzen, and F. Pfeiffer, “Short-range order in mesoscale systems probed by x-ray grating interferometry,” Europhys. Lett. 112(6), 68002 (2015). [CrossRef]  

51. K. Taphorn, F. De Marco, J. Andrejewski, T. Sellerer, F. Pfeiffer, and J. Herzen, “Grating-based spectral x-ray dark-field imaging for correlation with structural size measures,” Sci. Rep. 10(1), 13195 (2020). [CrossRef]  

52. C. Mätzler, “Autocorrelation functions of granular media with free arrangement of spheres, spherical shells or ellipsoids,” J. Appl. Phys. 81(3), 1509–1517 (1997). [CrossRef]  

53. T. Krouglov, I. M. de Schepper, W. G. Bouwman, and M. T. Rekveldt, “Real-space interpretation of spin-echo small-angle neutron scattering,” J. Appl. Crystallogr. 36(1), 117–124 (2003). [CrossRef]  

54. H. Lim, M. Veysoglu, S. Yueh, R. Shin, and J. Kong, “Random medium model approach to scattering from a random collection of discrete scatterers,” J. Electromagn. Waves Appl. 8(7), 801–817 (1994). [CrossRef]  

55. S. Torquato, Random Heterogeneous Materials, Microstructure and Macroscopic Properties (Springer, 2002).

56. F. Lado, “Numerical Fourier transforms in one, two, and three dimensions for liquid state calculations,” J. Comput. Phys. 8(3), 417–433 (1971). [CrossRef]