Neutron grating interferometry provides information on phase and small-angle scatter in addition to attenuation. Previously, phase grating moiré interferometers (PGMI) with two or three phase gratings have been developed. These phase-grating systems use the moiré far-field technique to avoid the need for high-aspect absorption gratings used in Talbot–Lau interferometers (TLI) that reduce the neutron flux reaching the detector. We first demonstrate, through theory and simulations, a novel phase grating interferometer system for cold neutrons that requires a single modulated phase grating (MPG) for phase-contrast imaging, as opposed to the two or three phase gratings in previously employed PGMI systems. The theory shows the dual modulation of MPG with a large period and a smaller carrier pitch P, resulting in large fringes at the detector. The theory was compared to the full Sommerfeld–Rayleigh diffraction integral simulator. Then, we proceeded to compare the MPG system to experiments in the literature that use a two-phase-grating-based PGMI with best-case visibility of around 39%. The simulations of the MPG system show improved visibility in comparison to that of the two-phase-grating-based PGMI. An MPG with a modulation period of 300 µm, the pitch of 2 µm, and grating heights with a phase modulation of (π,0, illuminated by a monochromatic beam produces visibility of 94.2% with a comparable source-to-detector distance (SDD) as the two-phase-grating-based PGMI. Phase sensitivity, another important performance metric of the grating interferometer, was compared to values available in the literature, viz. the conventional TLI with the phase sensitivity of 4.5 × 103 for an SDD of 3.5 m and a beam wavelength of 0.44 nm. For a range of modulation periods, the MPG system provides comparable or greater theoretical maximum phase sensitivity of 4.1 × 103 to 10.0 × 103 for SDDs of up to 3.5 m. This proposed MPG system appears capable of providing high-performance PGMI that obviates the need for the alignment of two phase gratings.

Neutrons are a useful probing tool in measuring material properties and imaging bulk materials due to their dual particle and wave nature, the latter described quantum mechanically by de Broglie wave packets showing interference phenomena.1 As neutron waves pass through matter, they undergo phase-shifts due to interactions with local, spatially dependent potentials, the most common being neutron–nucleus interactions.2,3 Small-angle neutron scattering (SANS) of neutrons, due to nuclear or magnetic interaction potential variations in the sample, also locally degrades the coherence of a well-defined neutron wave front.4,5 Neutron beams are also attenuated by nuclear reactions and incoherent scattering. Thus, neutron grating interferometry can image variations in phase change (differential phase contrast image), small-angle scattering (dark-field image), and attenuation (transmission image).1,5

Currently, there are two grating interferometry methods at the forefront of neutron phase imaging: the Talbot–Lau interferometer (TLI)6,7 and the phase-grating moiré interferometer (PGMI).8–11 The TLI can operate in the full field of a cold neutron beam (contrary to typical Mach–Zehnder interferometers) and has flexible chromatic coherence requirements (Δλ/λ. However, the absorption gratings reduce the neutron flux reaching the detector by about a factor of 4 for similar source conditions.9,10 The PGMI also operates in the full-field of a neutron beam, but since it produces directly resolvable interference fringes in the far-field, only a source grating is required, thus reducing the intensity by about a factor of 2. In contrast to the TLI, the PGMI has relaxed grating fabrication requirements, had a broader wavelength acceptance, and permitted control of the fringe period by varying the separation of the phase modulating gratings.10 

The purpose of this study is to investigate by simulation a novel PGMI system for cold neutrons that requires only a single modulated phase grating (MPG). The phase grating has a rectangular modulation in spatial width with a period of W and a finer “carrier” pitch P. An advantage of using a single grating is reducing grating misalignment issues commonly found in multiple-grating systems, which can lead to fringe visibility or contrast loss.9–11 

A non-interferometric grating system that uses only a single attenuation grating was reported by Strobl et al.,12 which has the advantages of its fringe visibilities being independent of the wavelength (achromatic) and its ability to extend the range of autocorrelation lengths probed in materials (particularly in nanometer scale) for dark-field imaging (DFI). However, the system as demonstrated has an SDD of 7.26 m with the grating-to-detector distances of only 50–300 mm and fringe visibilities of 10%–70%, limiting the phase sensitivity with its small grating-to-detector distances. Fringe visibility also falls off with higher grating-to-detector distance. This system’s visibilities are highly dependent on the geometric blurring due to its pinhole source with the collimation ratio L/D (source-to-grating distance/pinhole size). Thus, the fringe visibilities are susceptible to dropping sharply if the source-to-grating distance is too small or the pinhole size is too large.

The MPG concept was originally investigated in simulations for x ray by our group.13,14 In Ref. 14, we showed two different components of the grating interfere to create a larger fringe pattern on the detector, directly observable without the analyzer.

In another development, the modulated phase grating was demonstrated in experiments by Ref. 15. The mathematical treatise in Ref. 15 was unclear and, to the best of our understanding, appears to consider the Fresnel kernel to operate directly on the intensity instead of amplitude. The Fresnel kernel approximation applies to the amplitude,8,16–18 and then the intensity may be found similar to Ref. 8.

In this paper, we derived the theory behind a rectangular modulated phase grating from first principles using Fresnel zone approximation (which is also valid for far-field systems) to show that the intensity pattern on the detector is indeed of a directly observable period of magnification times large modulation period W. We also estimate distance conditions for maximum visibility for a given W and P.

We also developed an analytical Sommerfield–Rayleigh simulator with which we investigated the advantages of performance, visibility, and sensitivity of the MPG system over the two-phase-grating-based PGMI, which we will hereinafter refer to as the “standard-dual-grating” system. Before thoroughly investigating the performance of the MPG, we compared some results of our simulation method with the experiments previously conducted with the standard-dual-grating system and reported by Pushin et al.9 

We have constructed an efficient wave-propagation simulator for neutron PGMI applications that can accept various designs (e.g., modulated phase gratings14 or standard gratings). The code, named “N-SRDI,” was written in the C++ programming language and can simulate systems, such as those shown in Figs. 1 and 2, using the Sommerfeld–Rayleigh diffraction integrals (SRDI).16 These predict the observed complex-valued field amplitude A(y) from a wave that has diffracted from a grating or aperture and which originated from a source wave function U(Ps). The neutron field amplitude at a detector for a single-phase-grating system, such as the MPG (Fig. 1), can be expressed as
A(y)=1jλUPsejkr1r1T1y1cosθ1dy1,
(1)
and the intensity is
I(y)=Ay2,
(2)
where U(Ps) is the neutron source wave function, T1y1 is the neutron transmission function though the MPG in a plane perpendicular to the neutron beam propagation direction along the z-axis, k=2πλ is the wave number, λ is the wavelength, r1 is the distance between the grating point y1 to the detector point y given by r1=Dgd2+yy12, and θ is the angle between r1 and the normal of the MPG. The transmission function is given by
T1y1=A1(y1)ejϕ(y1),
(3)
where A1(y1) is the amplitude transmission of the neutron beam through the grating due to attenuation and ϕ(y1) is the phase shift determined by the spatial heights of the MPG.
FIG. 1.

Schematic diagram of the simulated MPG system (not to scale). In our simulations, the neutron source is a point source emitting a spherical wavefront along the z-axis (fringe patterns from a line source are later acquired with a convolution method). The source-to-MPG distance Dsg was chosen to be 0.5 or 1.0 m (depending on the pitch), and the MPG-to-detector distance Dgd can range from 1.0 to 5.0 m. The grating has a rectangular functional form that modulates in spatial height with a slow-varying period of W and a smaller “carrier” pitch P. The mathematical form of T(y1) is given in the theory section. The duty cycle of the grating is 50%. The different heights (h1,h2 of the grating along the z-axis correspond to phase shifts, for example, (π,π/4 or (π,0, for neutrons of wavelength 0.44 nm. The grating simulated is an ideal phase grating so no attenuation of the beam is considered, i.e., A1 = 1. The fringe period at the detector W′ is determined by the geometric magnification Dsd/Dsg from the point source.14 

FIG. 1.

Schematic diagram of the simulated MPG system (not to scale). In our simulations, the neutron source is a point source emitting a spherical wavefront along the z-axis (fringe patterns from a line source are later acquired with a convolution method). The source-to-MPG distance Dsg was chosen to be 0.5 or 1.0 m (depending on the pitch), and the MPG-to-detector distance Dgd can range from 1.0 to 5.0 m. The grating has a rectangular functional form that modulates in spatial height with a slow-varying period of W and a smaller “carrier” pitch P. The mathematical form of T(y1) is given in the theory section. The duty cycle of the grating is 50%. The different heights (h1,h2 of the grating along the z-axis correspond to phase shifts, for example, (π,π/4 or (π,0, for neutrons of wavelength 0.44 nm. The grating simulated is an ideal phase grating so no attenuation of the beam is considered, i.e., A1 = 1. The fringe period at the detector W′ is determined by the geometric magnification Dsd/Dsg from the point source.14 

Close modal
FIG. 2.

Schematic diagram of the simulated standard-dual-grating system (not to scale). Simulations were done to compare to fringe visibility measurements done by Pushin et al.9 In our simulations, the neutron source is a point source emitting a spherical wavefront along the z-axis (fringe patterns from a line source are later acquired with a convolution method). L1 = 1.2 m, and L2 = 1.79 m. The inter-grating spacing D ranges from 7 to 16 mm. The gratings have a standard binary form with no modulation in spatial height. The periods of the gratings are PG1 = PG2 = 2.4 µm. The duty cycles of the gratings are 50%. The heights of the G1 and G2 gratings are h1 and h2, respectively. They both correspond to a phase shift of 0.27π for neutrons of wavelength 0.44 nm. The gratings simulated are ideal phase grating so there is no attenuation of the beam for either grating, i.e., A1 = A2 = 1. The fringe period at the detector Pd is given by the ratio LPG2/D as given in Ref. 8.

FIG. 2.

Schematic diagram of the simulated standard-dual-grating system (not to scale). Simulations were done to compare to fringe visibility measurements done by Pushin et al.9 In our simulations, the neutron source is a point source emitting a spherical wavefront along the z-axis (fringe patterns from a line source are later acquired with a convolution method). L1 = 1.2 m, and L2 = 1.79 m. The inter-grating spacing D ranges from 7 to 16 mm. The gratings have a standard binary form with no modulation in spatial height. The periods of the gratings are PG1 = PG2 = 2.4 µm. The duty cycles of the gratings are 50%. The heights of the G1 and G2 gratings are h1 and h2, respectively. They both correspond to a phase shift of 0.27π for neutrons of wavelength 0.44 nm. The gratings simulated are ideal phase grating so there is no attenuation of the beam for either grating, i.e., A1 = A2 = 1. The fringe period at the detector Pd is given by the ratio LPG2/D as given in Ref. 8.

Close modal
The SRDI shown in Eq. (1) is a representation of the Huygens–Fresnel principle since the observed complex-valued field amplitude at the detector is a superposition of diverging spherical waves ejkr1/r1 originating from secondary sources located at each point y1 along the phase grating.16 In all simulations, the neutron source wave function U(Ps) is assumed to be a point source at point y0; however, the fringe patterns from a line source are later acquired with a convolution method. Therefore, Eq. (1) was simplified to
Ay=Dgdjλejk(r0+r1)r0r12T1y1dy1,
(4)
where r0 is the distance between the point source at point y0 and the point y1 on the MPG given by r0=Dsg2+y1y02 and the cos θ (between r1 and the z-axis) has been replaced by Dgd/r1, explaining the r12 term inside the SRDI shown in Eq. (4).
To verify the accuracy of the N-SRDI code, we also simulated standard-dual-grating systems for neutron PGMI applications and compared our results to experimental data from Pushin et al.9 Since two separate gratings are involved for this interferometer system, we use two SRDIs to calculate the observed complex-valued field amplitude A(y) at the detector. Using the Sommerfeld–Rayleigh formulation for diffraction, the neutron field amplitude at a detector for a standard-dual-grating system (Fig. 2) was simplified, using the same logic as before, to
Ay=DL2λ2ejk(r0+r1+r2)r0r12r22T1y1T2y2dy1dy2.
(5)

In all simulations (MPG or standard-dual-grating systems), we obtained results with a point source [that is U(Ps) was a point source at point yo emitting spherical wave eikr0/r0] and then used a convolution method on the fringe pattern of a point source to simulate the line-source, as will be described next.

1. Slit width

In the first step of the simulation, the neutron source was assumed to be a point source emitting a spherical wave in the direction of the z-axis toward the imaging detector. A line source, such as that in Ref. 9, can be thought of as an array of point sources that are mutually incoherent. As shown in Fig. 3, each point along a line source creates a fringe pattern at the detector that is spatially displaced with respect to the fringe pattern from the central point source. Their superposition causes the observed composite fringe pattern to be washed out, reducing the overall visibility of the fringe pattern. Consequently, a line source will have lower fringe visibility compared to a single, monochromatic point source.

FIG. 3.

Schematic illustration demonstrating visibility loss of an approximated line source (not to scale). The red fringe pattern on the detector is due to the central point source, and the green fringe patterns are due to the endpoints of the approximated line source. The superposition of these fringe patterns will result in visibility loss due to the green fringe patterns being shifted away from the red fringe pattern. With more intermediate point sources added to approximate the line source, the visibility is improved since there are more fringe patterns that are closer in phase to the red fringe pattern. However, the visibility of this approximated line source will always be lower than that of the fringe pattern of a single point source because of the wash-out effect many slightly shifted fringe patterns create.

FIG. 3.

Schematic illustration demonstrating visibility loss of an approximated line source (not to scale). The red fringe pattern on the detector is due to the central point source, and the green fringe patterns are due to the endpoints of the approximated line source. The superposition of these fringe patterns will result in visibility loss due to the green fringe patterns being shifted away from the red fringe pattern. With more intermediate point sources added to approximate the line source, the visibility is improved since there are more fringe patterns that are closer in phase to the red fringe pattern. However, the visibility of this approximated line source will always be lower than that of the fringe pattern of a single point source because of the wash-out effect many slightly shifted fringe patterns create.

Close modal

Instead of repeating costly SDRI computations by taking a series of point sources along with the extent of a line source (Fig. 3), we modeled the fringe pattern of a line source by first simulating the fringe-pattern due to a center-point source and then convolving it with a rectangular window function of the same size as the line source. This has the equivalent “wash-out” effect on the fringe intensity pattern at the detector from simulating an array of point sources that make up a line source.

2. Note on simulation of source grating G0

The G0 simulations are similar to simulations of a slit. When a source grating is used, instead of a single slit, we have several slits adjacent to each other in the Ys-axis of Fig. 3 with center separation multiples of P0. The P0 (or source period Ps in Ref. 8) is chosen such that each slit sends a net fringe pattern co-registered with the fringes from the other slits, resulting in a simple scale factor of intensity. Therefore, simulating a single slit vs source grating simulation is identical in the procedure, except for a large-scale factor increasing the intensity at the detector. The calculations of P0 and slit opening are available in the literature.14,19 In the Discussion section, we calculated the pitch P0 of G0 and open-ratio requirement for the MPG system for an example geometry we investigated.

3. Pixel-size convolution and subsampling

Furthermore, to account for the fringe blurring due to the pixel resolution and to have the correct number of sampled data points given a certain pixel size, two steps are taken: (1) another convolution of the fringe pattern at the detector by a rectangular window function that is the same size as the pixel (2) subsampling of the fringe intensity pattern at pixel-size increments since the intensity pattern at the detector is originally sampled at 0.1 µm. This is followed by interpolation to reconstruct the signal at higher sampling rates.

We simulated a case from the experimental setups from Ref. 9 to assess our convolution method that used a standard-dual-grating system with a monochromatic beam λ=0.44nm) exiting a 200 µm slit source with the system parameters of L1 = 1.2 m, D = 12 mm, and L2 = 1.78 m (Fig. 2). The identical phase gratings had pitch PG1=PG2=2.4μm and a 0.27π phase shift for 0.44 nm wavelength neutrons. Both gratings had duty cycles of 50%. The detector pixel size resolution was reported to be ∼100 µm. In simulations, the grating function is sampled at every 2 nm and the intensity pattern at the detector is sampled at every 0.1 µm. We used the two independent convolutions previously mentioned on the fringe pattern to account for blurring due to the slit width and pixel size. The fringe pattern is subsampled at 100 µm increments and intensity values are interpolated between the sampled values by a factor of 1000 using MATLAB’s interp function to return to a 0.1 µm sampling rate.

We compared our convolution method against two other methods: (1) the expected closed-form visibility function of slit size given in Ref. 9: V=Vosinc(πs/Ps), where Vo is the optimal fringe visibility from one point source, s is the width of the slit, and Ps is the source period. The fringe visibility for a 200 µm source using our convolution method was 32.9%, which closely matched the expected close form 32.7% visibility (blue and orange bars in Fig. 4). (2) The brute force method of simulating many individual point sources on the slit, illustrated in Fig. 3, was also compared. The yellow bar in Fig. 4 represents the visibility for 21 points evenly spaced about the 200 µm width of the slit, yielding a 32.7% visibility. As more intermediate point sources were uniformly added as represented by the red bar (41 total point sources), the visibility remained at the expected visibility value of 32.7%. Thus, we verified that a central point source simulation followed by a convolution with a rectangular window function can adequately model the fringe pattern of a slit source, reducing computation cost considerably by avoiding brute force computations of several point sources.

FIG. 4.

Verification of convolution method visibility for a 200 µm slit source in a monochromatic beam. The fringe visibility using the convolution method is 32.9% closely matched the expected 32.7% visibility, theoretically calculated with Eq. (10) given in Ref. 9: V=Vosinc(πs/Ps), where Vo is the optimal fringe visibility given by one point source (39.6%), s is the slit size (200 μm, and Ps is the source period (598 μm. Another method tested was modeling the slit as a series of evenly spaced point sources. The visibility converged after 21 points, with 21 and 41 points each yielding the expected visibility of 32.7%.

FIG. 4.

Verification of convolution method visibility for a 200 µm slit source in a monochromatic beam. The fringe visibility using the convolution method is 32.9% closely matched the expected 32.7% visibility, theoretically calculated with Eq. (10) given in Ref. 9: V=Vosinc(πs/Ps), where Vo is the optimal fringe visibility given by one point source (39.6%), s is the slit size (200 μm, and Ps is the source period (598 μm. Another method tested was modeling the slit as a series of evenly spaced point sources. The visibility converged after 21 points, with 21 and 41 points each yielding the expected visibility of 32.7%.

Close modal

We evaluated the N-SRDI obtained fringe visibilities V against experimental results with standard-dual-grating systems for the neutron PGMI in Ref. 9 using the same methods detailed in Sec. II B. When the gratings were close to D = 7 or 8 mm, cos(Φ1) in Fig. 2 is assumed to be close to 1 to avoid numerical problems. This is akin to a Fresnel-zone approximation in this region. We also evaluated whether the simulated fringe periods agreed with the theoretically predicted values.8 The evaluations are shown in Results under the corresponding heading.

1. Comparison of visibility

The simulation parameters to match experiments are detailed in Sec. II B. The grating-to-grating separation distance D was varied from 7 to 16 mm while keeping the L (or Dsd) and L1 (or Dsg) fixed (Fig. 2). The L2 (or Dgd) was adjusted accordingly. Our convolution method previously described was applied to the fringe patterns, accounting for the 200 µm slit and 100 µm pixel used in the experiments. Again, evaluations are compiled in Sec. III.

It is important to note that the slit size or pitch of a source grating G0 has a direct one-to-one effect on the fringe periods at the detector without a pinhole geometric magnification when PG1 = PG2. Thus, the source period (or the pitch of G0) is directly equal to the fringe period as shown in Refs. 8 or 9.

This is not true for the single MPG. The pinhole factor Refs. 14 and 19 is present for the slit-opening as tested in Sec. II D with simulations and theory.

The idea behind modulated phase grating was originally demonstrated by our group in simulations using the Sommerfeld–Rayleigh simulations in Ref. 14. The modulated phase grating has an envelope function with two grating phase heights, with a large period W (order of 100 µm) and a small modulation of pitch P (order of few μm) (Fig. 1). The grating has two clear components with different phase heights. As shown in Ref. 14, these two components interact in amplitude to produce intensity patterns on the detector, which have periods with a magnified version of W.

In what follows, we show the (a) effect of slit size on MPG system fringe patterns, (b) the mathematical theory behind the MPG (with details in  Appendix A), and finally the (c) Z-condition for maximum and minimum visibility.

1. Effect of slit size on MPG system interference pattern

Consider a geometry where the center source to MPG distance is L1 = 100 cm and MPG to detector distance is L2 = 200 cm, (L2/L1 = 2). Figure 5 shows the N-SRDI simulated intensity pattern [Eqs. (1) and (2)] for this case with a point source at the center, and two others shifted at +100 and −100 µm (shift perpendicular to propagation direction). As shown, these create three beam-patterns shifted from each other by ±200 µm. This shows that for a 200 µm slit (with extremities at ±100 µm), one expects the patterns to shift over a net 400 µm range (±200 µm). This can be explained by the pinhole magnification factor of the slit, L2/L1 = 2. In general, if the actual slit size is Ls μm, the convolution kernel to be used is Ks = Ls (L2/L1).

FIG. 5.

N-SRDI simulations for point sources at the center (black, 0 µm) and ±100 µm (blue and red, respectively). The dashed vertical lines show that red and blue patterns are separated by a total of 400 µm, which is explained by point source separation of 200 µm multiplied by a pinhole-factor of 2 for the geometry setting.

FIG. 5.

N-SRDI simulations for point sources at the center (black, 0 µm) and ±100 µm (blue and red, respectively). The dashed vertical lines show that red and blue patterns are separated by a total of 400 µm, which is explained by point source separation of 200 µm multiplied by a pinhole-factor of 2 for the geometry setting.

Close modal

This result is mathematically proved in  Appendix B (after the main theory in  Appendix A, highlighted next).

2. Mathematical theory behind the MPG

In this section and  Appendix A, the field equation for MPG is derived using the Fresnel approximation of the Sommerfeld–Rayleigh equation.

We consider a 1D detector, the y-axis as the axis along the detector, the y1-axis as the axis of the grating, and the z-axis as the propagation axis. First, we consider the grating function such as in Fig. 1. This can be written as
Ty1=gy1n=δy1nP+n=δy1nP2P/2×recty1P/2.
(6)
Note that second term inside the square brackets, 1×n=δ(y1nP2P/2), is to account for the fact that in between the spikes of the comb-train +δ(y1np) sampling the basic envelope function (which is given below), there is a transmission factor of 1 (and not zero). The convolution by the small half-pitch width rect function is to take into account the small post-size (P/2) of the grating (assuming a duty cycle of 50%).20 The equation is similar to a standard phase grating as shown in Wilde and Hesselink,17 except with the envelope function is allowed to be more complex than a single constant phase transmission factor of exp().

For the RECT MPG discussed in this paper, the envelope function is given by gy1=expiϕ1recty1W/2+expiϕ2recty1W/2W/21Wn=δ(y1nW), the ϕ1 and ϕ2 being the phase heights in regions h1 and h2, each of width W/2 for the MPG of interest. The first curly bracket shows one period of the envelope, which is repeated via the comp-convolution in the second curly bracket.

Since P/2 is effectively the opening of the grating, the coherence requirement should be of the order of P (not W). One of the periodicities of T(y1) is W from the envelope function. But this alone would necessitate high coherence requirement. But since the function is also multiplied by a carrier “sampling” period P, we expect the T(y1) also creates multiple periodic frequencies in amplitude, a combination of harmonics of 1W and 1P and not just 1W. We explicitly derive this in  Appendix A.

Note, the Fresnel approximation holds where the incident and generated Huygens’s spherical waves in Eq. (4) can be approximated as paraboloid.16 In  Appendix A, the field amplitude is derived explicitly in terms of the Fourier coefficients of the grating transmission factor T(y1) from first principles outlined in Goodman.16 

The field amplitude, U(y, z), for the case of the RECT MPG in Eq. (6) is derived in  Appendix A. The key findings are represented here. The amplitude field may be divided into (W, P)-dependent and P-dependent parts,
Uy,z=AexpiB1L1+zn=+m=+gmsinc12mPW+n×expjπλL1zL1+zmW+nP2expj2πL1yL1+zmW+nP+n=+expjπnsincn2expjπλL1zL1+znP2×expj2πL1yL1+znP=AexpiB1L1+zU1y,z,W,P+U2y,z,P.
(7)
These terms are reminiscent of Refs. 18 and 21 except here the frequency is related to a combination of harmonics of 1W and 1P. More specifically, in the first term U1, the frequencies are related to mW+nP and in the second term U2 to nP.
The intensity as both mW+nP modulation terms ( Appendix A) but since n/P will be hardly visible directly in the detector, it is approximated by the n = 0 harmonic,
Iy,zAL1+z2C22z,0+m=+C11z,m,0+C12z,m,0+C12*z,m,0expj2πL1yL1+zmW=AL1+z2I0+m=1+Im(z)cos2πL1yL1+zmW+θm,
(8)
where Cij is autocorrelation terms shown in  Appendix A. Therefore, the intensity produces fringe patterns of period mMW, where the magnification M = L1+zL1. In our case of an even grating, θm = 0.

Figure 6 shows an example where the Cij [in Eq. (8) above] were evaluated for the MPG with phase heights (π, 0) and W = 300 µm in order to calculate the intensity. The N-SRDI simulations for the same MPG parameters and geometry settings are also plotted in Fig. 6, showing an excellent match.

FIG. 6.

Intensity from Eq. (8) compared to the N-SRDI simulation, for MPG W = 300 µm, phase heights (π, 0) with geometry Dgd (or L2) = 200 cm and Dsg (or L1) = 100 cm. The pixel size is 100 µm, and the actual slit size is Ls = 200 µm, with the slit convolution kernel Ks = Ls × L2/L1 = 400 µm.

FIG. 6.

Intensity from Eq. (8) compared to the N-SRDI simulation, for MPG W = 300 µm, phase heights (π, 0) with geometry Dgd (or L2) = 200 cm and Dsg (or L1) = 100 cm. The pixel size is 100 µm, and the actual slit size is Ls = 200 µm, with the slit convolution kernel Ks = Ls × L2/L1 = 400 µm.

Close modal

3. Conditions for maximum/minimum visibility

As shown in  Appendix A, the visibility of the first harmonic will oscillate in Z with extrema occurring at
Z=L1WPk2λL1WPk=L12λL1WPk1,
(9)
where k′ is an integer. For the Z = 1–5 m of interest, we found some peaks and troughs as reported in Table I, which was corroborated by the N-SRDI simulations. For example, for W = 120 µm and P = 1 µm, the visibility has a local maximum at Z = 2.14 m and a local minimum at Z = 4.5 m.
TABLE I.

Distance for maximum or minimum visibility with simulation vs Eq. (9).

W (μm)P (μm)Z@Vmax (N-SRDI) (m)Z@Vmax [Eq. (9)] (m)Z@Vmin (N-SRDI) (m)Z@Vmin [Eq. (9)] (m)
120 4.5 4.5 (for k = 3) ⋯ ⋯ 
120 2.14 2.14 (for k = 5) 4.5 4.5 (for k = 6) 
W (μm)P (μm)Z@Vmax (N-SRDI) (m)Z@Vmax [Eq. (9)] (m)Z@Vmin (N-SRDI) (m)Z@Vmin [Eq. (9)] (m)
120 4.5 4.5 (for k = 3) ⋯ ⋯ 
120 2.14 2.14 (for k = 5) 4.5 4.5 (for k = 6) 

Note that the maximum Z condition (which is independent of the envelope function Fourier parameters, gm) should be similar for all MPG height pairs. But the envelope function or the height difference in the RECT MPG will affect the visibility.

To confirm this, we applied the N-SRDI to simulate different MPGs with h2 for the distance ranges of interest here (Dgd = 1–5 m or Dsd = 2–6 m from the source) and the wavelength 0.44 nm. We kept the maximum height of grating h1 such that it yields π-shift and varied the h2 to provide 0, π/4, π/2, and 3π/4 phase shifts (for λ=0.44nm. The results for W = 120 µm, P = 2 µm and P = 1 µm are shown in Fig. 7. We, therefore, confirmed that the Z at maximum visibility is independent of the MPG-type, but the visibility will drop as the h2 is increased. Note for the simulations, the location of peaks was very close whether we applied the slit smoothing or not. In Fig. 7, a pixel size of 100 µm and slit of Ls = 75 µm with the magnified convolution kernel Ks = Ls × L2/L1 were applied for each point.

FIG. 7.

The visibility plots for N-SRDI simulations for different phase height pairs, keeping the larger height fixed at π and changing the smaller height from 0 to 3π/4. (Left panel) W = 120 µm, P = 2 µm and (Right panel) W = 120 µm, P = 1 µm. Note maximum or minimum visibility position is independent of the phase heights. The (π, 0) yielded the highest visibility. The maxima and minima correspond to what is predicted by Eq. (9).

FIG. 7.

The visibility plots for N-SRDI simulations for different phase height pairs, keeping the larger height fixed at π and changing the smaller height from 0 to 3π/4. (Left panel) W = 120 µm, P = 2 µm and (Right panel) W = 120 µm, P = 1 µm. Note maximum or minimum visibility position is independent of the phase heights. The (π, 0) yielded the highest visibility. The maxima and minima correspond to what is predicted by Eq. (9).

Close modal

We observed that the peak visibility and minimum visibility matched for the N-SRDI and the values from Eq. (9). In Table I, we tabulate the distances and the corresponding k-values in the distance range of interest Z from 1 to 5 m.

The fringe visibility was highest for h2 = 0 for all MPG–detector distances. This is expected because it would provide the maximum difference between the two components of the grating. This was also the grating that was investigated for x rays in Ref. 15.

Since the fringe visibility was highest for h2 = 0 for all MPG–detector distances, henceforth in this work we used (π, 0), as shown in the schematic diagram Fig. 8.

FIG. 8.

Schematic of the MPG system (not to scale) with a slit source with width Sw, MPG with grating heights with phase modulation of (0, π), and detector with a phase object at an object-to-detector distance Dod.

FIG. 8.

Schematic of the MPG system (not to scale) with a slit source with width Sw, MPG with grating heights with phase modulation of (0, π), and detector with a phase object at an object-to-detector distance Dod.

Close modal
a. Evaluations of N-SRDI for MPG with (π, 0) case.

Fringe visibility, fringe period, and maximum phase sensitivity were evaluated for different modulation periods W, pitches P, slit widths Sw, source-to-MPG distances Dsg, and MPG-to-detector distances Dgd, as shown in Sec. III.

We investigated how a polychromatic beam may degrade the visibility of our MPG system. The polychromatic beam was modeled after the one at the NG6 Cold Neutron Imaging Facility (NCI) at the NCNR, which is approximately given by a Maxwell–Boltzmann distribution with a peak wavelength of λc = 0.5 nm.9 The same MPG is used in both the monochromatic and polychromatic simulations. Since the polychromatic beam can be thought of as the aggregate of incoherent sources, the fringe patterns of each wavelength in the spectrum can be added in intensity according to their weight in the spectrum. The evaluation results are in Sec. III.

Phase objects, that is, objects that introduce a pure phase shift on the wavefront on the path of the beam, were simulated mathematically, and the effect on the interference fringe at the detector was analyzed.

The shift in the interference pattern Δy at the detector is related to the refractive angle α imposed by an object on the wave field, and therefore the object’s differential phase-shift dΦ/dy is as follows:14,22
Δy=DodtanαλDod2πdΦdy.
(10)

Estimating Δy (which itself is a function of y in general) can therefore yield the object’s differential phase shift dΦ/dy. Note that Δy also depends on system parameters: object-to-detector distance Dod and the wavelength λ, which have to be corrected for a true estimate of the object’s differential phase.

Since we have a large interference pattern period W′ (usually greater than the pixel size), and the detector sampling rate is above the Nyquist sampling rate, we use a single-shot recovery using Fourier transforms to demonstrate our system properties, following the method in Ref. 23. The steps followed are shown below. Before showing the steps, we like to remind you that the object’s differential phase at the detector space is denoted by dΦ/dy (this is to be calculated) and the measured Δϕ is the actual phase difference calculated via delay Δy observed at the detector.

Step 1. (a) Take fast Fourier transform (FFT) of each simulated fringe pattern with and without an object. Isolate the first harmonic by windowing the Fourier transforms on one side. A window-width of total width 1/W′ is chosen symmetrically around the first harmonic peak located at 1/W′ [i.e., the window extends ±1/(2W′) around 1/W′]. (b) Perform Inverse Fast Fourier transform (IFFT) of the one-sided windowed FFTs and take the angle difference of the two complex spatial domain signals. The difference Δϕ is related to the spatial shift between the fringe pattern signals Δy as
Δϕ=ϕϕb=2πWΔy.
(11)

Note: the IFFT of the windowed FFT yields complex exponentials in the spatial domain, ideally the first harmonic signals. These have the frequency 2π/W′. The difference of phase angles Δϕ between the blank fringe pattern phase angle ϕb and the with-object fringe pattern phase ϕ is related to the spatial shift Δy between these two signals as given by Eq. (11).

Implementation details: Note: The N-SRDI code outputs detector data with a detector pixel size of 0.1 µm. We convolve the fringe patterns by independent rectangular window functions to compensate for the fringe blurring due to the slit width and pixel size. Then we subsample the data at a rate equal to the pixel size which is 100 µm.

We interpolate the subsampled pixel-size detector data up by a factor of 1000 by using MATLAB’s interp function to resample the function at 0.1 µm before performing the FFT and IFFT (these are the same processing steps used for the standard-dual-grating simulations). The FFT, IFFT, and angle differences of the fringe patterns are all performed in MATLAB. Each phase angle ϕb and ϕ is obtained modulo ±π. The resultant phase angle difference has to be unwrapped before integration. Note for some small objects, the Δy shift can be sub-pixel. However, since we sample the fringe pattern above the Nyquist rate (i.e., pixel size < W′/2), we can always reconstruct the signals with interpolation and obtain the phase difference. The interpolation used is the default one used by MATLAB’s interp function, which uses a symmetric special FIR filter that passes the original samples unchanged while minimizing the sum-squared error between the original and the ideal values.

The interpolation step is optional but adds to the robustness of the recovery for cases where the Nyquist edge (frequency) is too close to the first harmonic window. The subpixel interpolation yields a higher sampling rate and moves the Nyquist edge away from the first harmonic in Fourier domain, which makes it easier to place an automated symmetric window in the Fourier domain to retrieve the first harmonic.

Step 2(a). To obtain the object’s differential phase dΦ/dy at the detector space (y) from the measured phase difference Δϕ, we multiply by Sλ1. We show this by combining Eqs. (10) and (11) to obtain
Δϕ=λDodWdΦdy=SλdΦdy,
(12)
where S = Dod/W′ is the phase sensitivity for the interferometer. Expressed in terms of the object’s differential phase profile dΦ/dy,
dΦdy=Sλ1Δϕ.
(13)
Thus, the measured phase difference Δϕ must be “corrected” by 1/(λS) to obtain Φ/dy.

Step 2(b). To obtain the object’s differential phase dΦ/dyob at the object space (yob), we have to scale (de-magnify) the spatial variation of the detector (y) by 1/Mobj, where Mobj = object magnification.

Step 3. We integrate dΦ/dyob to obtain the object phase profile. We use MATLAB’s cumtrapz function for the integration. The scaled grid (object grid) is provided in MATLAB’s cumtrapz function as the grid argument. Note that we could have also magnified the true object’s grid and compared both phase profiles in the detector space (y), but we wanted to avoid any processing with the true object’s phase profile.

Figure 9 shows an example simulated blank and with-object fringe pattern. A phase object is a triangular object with a maximum phase shift of 8π rad, ramping up/down over ±800 µm spatial extent. It is placed at an object-to-detector distance of 2.5 m. In reality, this phase profile could be a silicon wedge sample with a maximum height of 275.6 µm at the center and falling off to zero over ±800 µm on either side of its peak. Figure 9(a) shows both fringe patterns just after being convolved with two independent rectangular filters to account for the blurring of the slit width of 50 µm and the pixel size of 100 µm. Figure 9(b) shows the fringe patterns after being subsampled at every 100 µm and then interpolating between sampled values to return to a 0.1 µm sampling rate. In the N-SRDI code, the objects are represented simply as mathematical phasors eiΦyob, which are applied to the wave field at the object-to-detector distance, where Φ(yob) is the object phase profile.

FIG. 9.

(a) Interference pattern with and without a triangle phase object with a peak phase of 8π rad using the system parameters W = 50 µm, pitch = 2.0 µm, Dsd = 5 m, Dsg = 1 m, and Dod = 2.5 m. Monochromatic beam (wavelength 0.44 nm) used Intensity values are originally sampled along the detector plane at 0.1 µm increments and are filtered by two independent rectangular window functions for blurring effects due to slit width and pixel size. (b) shows both fringe patterns subsampled at 100 µm increments and interpolated back to a sampling rate of 0.1 µm. The dark lines on the left and right are locations where the Δy shift of the fringe pattern was checked manually.

FIG. 9.

(a) Interference pattern with and without a triangle phase object with a peak phase of 8π rad using the system parameters W = 50 µm, pitch = 2.0 µm, Dsd = 5 m, Dsg = 1 m, and Dod = 2.5 m. Monochromatic beam (wavelength 0.44 nm) used Intensity values are originally sampled along the detector plane at 0.1 µm increments and are filtered by two independent rectangular window functions for blurring effects due to slit width and pixel size. (b) shows both fringe patterns subsampled at 100 µm increments and interpolated back to a sampling rate of 0.1 µm. The dark lines on the left and right are locations where the Δy shift of the fringe pattern was checked manually.

Close modal

The expected spatial shifts in the fringe pattern due to this triangular phase object are given by Eq. (10). Since the object phase consists of an up and down ramp, the spatial shift Δy is a positive constant on the left side of the fringe pattern center (detector position 0 μm and a negative constant on the right side of the fringe pattern center. As a manual check, we calculated the average Δy in the two regions over three cycles each at the 90% normalized intensity line in Fig. 9(b). In Table II, we compare the expected Δy and our measurements showing that the values are very close to each other.

TABLE II.

Measured fringe patterns shift from left and right windows in Fig. 9.

Fringe pattern ΔyLeft window Δy(μm)Right window Δy(μm)
Expected 5.5 −5.5 
Measured from plot 5.0 ± 1.4 −5.1 ± 1.3 
Fringe pattern ΔyLeft window Δy(μm)Right window Δy(μm)
Expected 5.5 −5.5 
Measured from plot 5.0 ± 1.4 −5.1 ± 1.3 

Figure 10 shows the Steps 1–3 recovery with the blank and with-object fringe patterns shown in Fig. 10(b). First, we show the measured Δϕ after Step 1 in Fig. 10(a). This approximately shows the rectangular pattern expected from the differential phase of the triangular object consisting of two ramps. The differential phase is shown in Fig. 10(b) after correction by Sλ1 as in Step 2. In Fig. 10(c), the detector grid is scaled by the object magnification to switch to the object’s grid space (yob) instead of the detector grid space (y). Then, Fig. 10(d) shows the integrated phase of the object (Step 3).

FIG. 10.

(a) The measured phase difference Δϕ at the detector from Fig. 9(b). Obtained via the single-shot method explained in Step 1. (b) The measured phase difference Δϕ is corrected by a Sλ1 factor to obtain the object differential phase dΦ/dy at the detector space (y). (c) The grid of dΦ/dy is scaled by the inverse object magnification factor Dso/Dsd to obtain the object differential phase dΦ/dyob in the object space (yob). (d) The object phase Φ retrieved from dΦ/dyob by integration (Step 3). For this instance, no detector noise was added.

FIG. 10.

(a) The measured phase difference Δϕ at the detector from Fig. 9(b). Obtained via the single-shot method explained in Step 1. (b) The measured phase difference Δϕ is corrected by a Sλ1 factor to obtain the object differential phase dΦ/dy at the detector space (y). (c) The grid of dΦ/dy is scaled by the inverse object magnification factor Dso/Dsd to obtain the object differential phase dΦ/dyob in the object space (yob). (d) The object phase Φ retrieved from dΦ/dyob by integration (Step 3). For this instance, no detector noise was added.

Close modal

Finally, Fig. 11 elaborates the Step 1 and Step 2(b) of the phase recovery process for the triangular and parabolic phase objects for object-to-detector distances varied from 0.5 to 3.5 m. Figures 11(a) and 11(c) show the measured phase Δϕ (Step 1) at the detector grid. Figures 11(b) and 11(d) show the object differential phase dΦ/dyob [Step 2(b)] on the object grid space (yob). They both approximately show the rectangular and linear dΦ/dyob expected for triangular and parabolic phase profiles, respectively.

FIG. 11.

(a) Measured phase difference ΔΦ at the detector for object-to-detector distances of 0.5–3.5 m for the 8π rad triangle phase object. (b) The corresponding object differential phase dΦ/dyob at the object plane. (c) Measured phase shifts ΔΦ at the detector for object-to-detector distances of 0.5–3.5 m for the 8π rad parabolic phase object. (d) The corresponding object differential phase dΦ/dyob at the object plane. Δϕ is corrected by a Sλ1 factor to obtain dΦ/dy, and the detector grids are also corrected by the inverse object magnification factor Dod/Dsd.

FIG. 11.

(a) Measured phase difference ΔΦ at the detector for object-to-detector distances of 0.5–3.5 m for the 8π rad triangle phase object. (b) The corresponding object differential phase dΦ/dyob at the object plane. (c) Measured phase shifts ΔΦ at the detector for object-to-detector distances of 0.5–3.5 m for the 8π rad parabolic phase object. (d) The corresponding object differential phase dΦ/dyob at the object plane. Δϕ is corrected by a Sλ1 factor to obtain dΦ/dy, and the detector grids are also corrected by the inverse object magnification factor Dod/Dsd.

Close modal

While these concepts were demonstrated in Sec. II without adding detector noise, we added realistic Poisson noise to the detector in Sec. III before phase recovery.

In Sec. III, we show quantitative phase recovery error analysis of triangular, parabolic, and trapezoidal objects with object-to-detector distances of 0.5–3.5 m and maximum object phases of 0.8π rad and 8π rad. This was done after adding realistic Poisson noise to the detector pixels, which we explain in detail in Sec. III.

To test our simulator N-SRDI, we simulated the standard-dual-phase grating system and compared our results to theory and experiments. We varied the grating-to-grating distance from 7 to 16 mm in increments of 1 mm and obtained fringe periods and fringe visibilities for the monochromatic experiments performed in Ref. 9 with λ = 0.44 nm. The fringe periods obtained via the simulations were compared to the theoretical closed-form prediction given in Eq. (11) in Ref. 8, as shown in Fig. 12. We observed an excellent agreement of the simulated fringe periods to the corresponding theoretical values.

FIG. 12.

Fringe period from standard-dual-grating simulations compared to theory. N-SRDI simulations were done with the monochromatic configuration parameters used in Ref. 9. The theoretical fringe period is given by Eq. (11) in Ref. 8.

FIG. 12.

Fringe period from standard-dual-grating simulations compared to theory. N-SRDI simulations were done with the monochromatic configuration parameters used in Ref. 9. The theoretical fringe period is given by Eq. (11) in Ref. 8.

Close modal

Figure 13 plots the N-SRDI with convolution fringe visibility results and the experimental fringe visibility results for the monochromatic experiments in Ref. 9. Not all but every other grating-to-grating distance was simulated and shown below since we varied the grating-to-grating distance in spacings of 1 mm. The general trend was captured by the N-SRDI with convolution simulations with the visibility peaking around 11–12 mm. The curve also fell within the statistical agreement of experiments (indicated by the error bars) except for the grating-to-grating distances D of 7, 8, and 11 mm.

FIG. 13.

Fringe visibility comparison of N-SRDI with convolution simulations to experiments performed in Ref. 9 The dashed lines show the closed-form expected visibility given in Eq. (12) in Ref. 8, which was plotted in Ref. 9.

FIG. 13.

Fringe visibility comparison of N-SRDI with convolution simulations to experiments performed in Ref. 9 The dashed lines show the closed-form expected visibility given in Eq. (12) in Ref. 8, which was plotted in Ref. 9.

Close modal

We simulated the MPG system with grating heights modulation of (0, π), the same slit-source in a monochromatic beam (λ = 0.44 nm) used for the standard-dual-grating system in Sec. III A, different modulation periods W, and the two different pitches of 2.0 and 1 µm. The pitch of P = 2.0 µm is approximately the period of some of the gratings currently available at the National Institute of Standards and Technology Center for Neutron Research (NCNR).

We show an excellent performance for an equivalent geometry as the set up in Ref. 9. This is shown in Fig. 14(a) of interference fringe “carpet” for W = 300 µm and P = 2 µm, where the intensity carpet is obtained by placing the detector at different distances from the MPG. The carpet shows a diverging self-image of the modulation pattern. At the source-to-MPG distance of 1 m and MPG-to-detector distance of 2.0 m (SDD of 3 m), we are at a nearly equivalent SDD of 2.99 m used in the setup in Ref. 9 for their monochromatic configuration. A fringe visibility analysis of the whole carpet is shown in Fig. 14(b). The normalized fringe pattern with the maximum visibility of 94.2% at the MPG-to-detector distance of 2 m is shown in Figs. 14(b) and 14(c).

FIG. 14.

(a) Interference fringe “carpet” generated with N-SRDI and pixel and source convolution, for the source-to-MPG distance of 1 m, W = 300 µm, and pitch P = 2.0 µm. The pixel size is 100 µm, and the slit width Ls = 200 µm (kernel Ks = Ls × L2/L1). (b) Fringe visibility analysis of the carper shows a maximum visibility of about 94.2% at the MPG-to-detector distance of about 2 m. (c) The normalized fringe pattern at the MPG-to-detector distance of 2 m, with maximum visibility of about 94.2%.

FIG. 14.

(a) Interference fringe “carpet” generated with N-SRDI and pixel and source convolution, for the source-to-MPG distance of 1 m, W = 300 µm, and pitch P = 2.0 µm. The pixel size is 100 µm, and the slit width Ls = 200 µm (kernel Ks = Ls × L2/L1). (b) Fringe visibility analysis of the carper shows a maximum visibility of about 94.2% at the MPG-to-detector distance of about 2 m. (c) The normalized fringe pattern at the MPG-to-detector distance of 2 m, with maximum visibility of about 94.2%.

Close modal

The example shows that the MPG system may be expected to yield more than twice the visibility compared to the standard-dual-grating system’s maximum visibility of ∼39% in a nearly equivalent set-up geometry to the one used in Ref. 9. Note that at this magnification the fringe period is about 900 µm; therefore, the slit or pixel size effects are relatively small.

We summarize other cases in Figs. 15 and 16 with high visibility. In Fig. 16, we also point out in which cases the geometry is more compact than the system in Ref. 9. In Figs. 15(a) and 15(b), we show the visibility and fringe period for different modulation periods, W = 50–600 µm (all with a pitch of 2.0 µm), at varying MPG-to-detector distances. We make sure to not include cases where the pixel size > W′/2, which would lead to insufficient sampling of the fringe pattern (Nyquist theorem) and cause aliasing. We observed several operating points with fringe visibility V > 40% in Figs. 15(a) and 15(b) for different fringe periods. The ability to control the fringe period is important as the fringe period determines the autocorrelation probing lengths in dark-field imaging (DFI),24 and it is also inversely proportional to the phase sensitivity.14 

FIG. 15.

Measured (a) fringe visibility and (b) fringe periods from simulations with varying MPG-to-detector distances of 1.0–5.0 m. The fringe period W′ is given by the geometric magnification Dsd/Dsg of W. The source-to-MPG distance is 1.0 m, and the pitch is 2.0 µm. The pixel size is 100 µm, and the slit width Ls = 200 µm (kernel Ks = Ls × L2/L1). The fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2).

FIG. 15.

Measured (a) fringe visibility and (b) fringe periods from simulations with varying MPG-to-detector distances of 1.0–5.0 m. The fringe period W′ is given by the geometric magnification Dsd/Dsg of W. The source-to-MPG distance is 1.0 m, and the pitch is 2.0 µm. The pixel size is 100 µm, and the slit width Ls = 200 µm (kernel Ks = Ls × L2/L1). The fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2).

Close modal
FIG. 16.

Fringe visibility vs period for a source-to-MPG distance of 1.0 m, varying MPG-to-detector distances from 1.0 to 5.0 m, and pitch of 2.0 µm. The pixel size is 100 µm and the slit Ls = 200 µm (kernel Ks = Ls × L2/L1). The cases with an SDD of 2.99 m (≈3.0 m) like the monochromatic experimental configuration in Ref. 9 are shown with the red arrows.

FIG. 16.

Fringe visibility vs period for a source-to-MPG distance of 1.0 m, varying MPG-to-detector distances from 1.0 to 5.0 m, and pitch of 2.0 µm. The pixel size is 100 µm and the slit Ls = 200 µm (kernel Ks = Ls × L2/L1). The cases with an SDD of 2.99 m (≈3.0 m) like the monochromatic experimental configuration in Ref. 9 are shown with the red arrows.

Close modal

To visualize the information in a meaningful way, we plot the visibility vs fringe-period in Fig. 16, with the MPG-to-detector distance implicitly varied. The nine circles in each curve correspond to the visibilities at nine different MPG-to-detector distances (100–500 cm, with intervals of 50 cm), i.e., the SDDs range from 200 to 600 cm. We partitioned into a green and light blue zone where the MPG provided high visibility, V > 39%, and acceptable 20%–39%, respectively. In order to compare to the standard-dual-grating system (in Fig. 13), the peak visibility from those experiments is shown as the asterisk with an error bar (falling on the “good visibility” zone). The dashed (purple) vertical line (through the asterisk case) shows that for the same fringe period for the asterisk case (from the monochromatic setup in Ref. 9 with the best-case visibility of 39%), there is the W = 300 µm MPG design that can yield better visibility of ∼75% visibility for MPG-to-detector distance of 100 cm (SDD of 200 cm). Other designs, such as W = 200 µm, yield 40%–65% visibility for a fringe period 0.5 mm or lower.

If one is looking for high visibility V > 39% with a compact setup geometry, we point out the operating points marked with yellow stars in Fig. 16 where the SDD ≤ 300 cm.

We, therefore, show several operating points of the MPG system with higher visibility than the highest reported visibility of 39% from the standard-dual-grating system in Ref. 9, with similar, smaller or larger fringe periods and the same or more compact geometry.

We investigated a lower pitch for the lower fringe period region in the dashed square in Fig. 16. Lowering pitch also enabled the Dsg to be smaller (0.5 m) and the coherence requirement to hold. The smaller fringe period zone is shown in Fig. 17 for the P = 1 µm cases. We make sure not to include cases where the pixel size >W′/2, which would lead to insufficient sampling of the fringe pattern (Nyquist theorem). Unfortunately, the slit kernel convolution effect dominated at these lower fringe periods and the visibility was poor for Ls = 200 µm with convolution kernel varying as Ks = Ls × L2/L1.

FIG. 17.

Fringe visibility vs period for Dsg = 0.5 m, Dgd = (1.0–5.0 m), and pitch = 1.0 µm. The fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2). The pixel size is 100 µm, and the slit Ls = 200 µm (kernel Ks = Ls × L2/L1). The slit kernel convolution dominated at the lower fringe periods and visibility suffered.

FIG. 17.

Fringe visibility vs period for Dsg = 0.5 m, Dgd = (1.0–5.0 m), and pitch = 1.0 µm. The fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2). The pixel size is 100 µm, and the slit Ls = 200 µm (kernel Ks = Ls × L2/L1). The slit kernel convolution dominated at the lower fringe periods and visibility suffered.

Close modal

We also investigated if the visibility for the smaller W cases (50, 70, and 120 µm) where the fringe period is smaller, could be improved by using narrower slit widths Sw (albeit knowing that the flux reaching the detector would be reduced proportionately and typically a G0 grating would be required instead of a slit). We kept the Sw (or Ls) = 200 µm as a reference. As expected in general, we see in Figs. 18(a)18(c), using a narrower slit width Sw improved the visibility for the W = 50, 70, and 120 µm for these fringe patterns with smaller periods.

FIG. 18.

(a)–(c) Fringe visibility for slit sizes Sw (or Ls) of 50, 110, and 200 µm, respectively. The source-to-MPG distance is 1.0 m, and the pitch is 2.0 µm. The MPG-to-detector distances range from 1.0 to 5.0 m. The pixel size is 100 µm. The slit convolution kernel is calculated as Ks = Sw × L2/L1).

FIG. 18.

(a)–(c) Fringe visibility for slit sizes Sw (or Ls) of 50, 110, and 200 µm, respectively. The source-to-MPG distance is 1.0 m, and the pitch is 2.0 µm. The MPG-to-detector distances range from 1.0 to 5.0 m. The pixel size is 100 µm. The slit convolution kernel is calculated as Ks = Sw × L2/L1).

Close modal

Finally, we calculated the maximum phase sensitivities of our MPG system. The maximum phase sensitivity is given by Smax = Dgd/W′ (where the object is closest to the MPG, i.e., Dod = Dgd). This is analyzed for different MPG-to-detector distances as shown in Fig. 19(a) for different grating modulation periods with source-to-MPG distance kept at 100 cm. There is also an increase in phase sensitivity with lower W′ (lower fringe period) as shown in Fig. 19(b) for a grating-to-detector distance of 5.0 m. Phase sensitivity information is sparse for the standard-dual-grating system but available in the literature for the conventional neutron Talbot–Lau interferometer (TLI). For MPG modulation periods W between 50 and 120 µm, we note that the theoretical maximum phase sensitivities are between 4.1 × 103 to 10.0 × 103 (0.41–1 × 104) when using SDDs up to 3.5 m (MPG-to-detector distance up to 2.5 m). These sensitivities are comparable or greater than that of the conventional neutron TLI of 4.5 × 103 with a SDD of 3.5 m and a beam wavelength of 0.44 nm.7,25

FIG. 19.

(a) Maximum phase sensitivity for different MPG-to-detector distances. (b) Phase sensitivity vs fringe period for the MPG-to-detector distance of 5 m. The pixel size is 100 µm and the slit width is 200 µm. Fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2).

FIG. 19.

(a) Maximum phase sensitivity for different MPG-to-detector distances. (b) Phase sensitivity vs fringe period for the MPG-to-detector distance of 5 m. The pixel size is 100 µm and the slit width is 200 µm. Fringe sampling frequency is ensured to be greater than the Nyquist rate (pixel size < W′/2).

Close modal

We used a polychromatic beam for the setup which produced the best visibility of 94.2% in the monochromatic simulations (λ=0.44nm to see how the fringe visibility would degrade. The polychromatic beam is approximately described by Maxwell–Boltzmann spectrum with a peak wavelength λc = 0.5 nm. The MPG used in both monochromatic and polychromatic simulations were the same and had fixed phases of (0, π) at λ = 0.44 nm. The system parameters include a grating with modulation period W of 300 µm, pitch of 2 µm, source-to-MPG distance of 1 m, and MPG-to-detector distance of 2.0 m. The pixel size is 100 µm, and the slit width is 200 µm. While we note that the visibility dropped from 94.2% to 68.9% (as calculated from intensities shown in Fig. 20), this polychromatic visibility is still significantly higher than the standard-dual-grating system with a best-case monochromatic visibility of around 39% in Ref. 9.

FIG. 20.

Fringe visibility degradation due to polychromatic beam for our MPG system. We consider the configuration that produced the best visibility of 94.2% using a monochromatic beam of λ = 0.44 nm (blue trace). The polychromatic beam described by a Maxwell–Boltzmann spectrum with peak wavelength λc = 0.5 nm produces a fringe pattern with visibility of 68.9% (orange trace). The system parameters include a grating with modulation period W of 300 µm, pitch of 2 µm, source-to-MPG distance of 1 m, and MPG-to-detector distance of 2.0 m.

FIG. 20.

Fringe visibility degradation due to polychromatic beam for our MPG system. We consider the configuration that produced the best visibility of 94.2% using a monochromatic beam of λ = 0.44 nm (blue trace). The polychromatic beam described by a Maxwell–Boltzmann spectrum with peak wavelength λc = 0.5 nm produces a fringe pattern with visibility of 68.9% (orange trace). The system parameters include a grating with modulation period W of 300 µm, pitch of 2 µm, source-to-MPG distance of 1 m, and MPG-to-detector distance of 2.0 m.

Close modal

Figure 21 shows examples of linear, trapezoid, and quadratic phase objects recovered using the single-shot method described in Sec. II F. The recovery is shown for the noiseless case and with realistic Poisson noise. The noise addition step is explained as follows: First, Poisson statistics with a realistic flux level to our detector intensity was added. To determine a realistic flux, we have taken the average intensity for experiments shown in Ref. 9, shown to be about 10 000. Taking into account the scintillation gain camera settings, the actual neutron counts are about a factor of ten smaller for the case shown in Ref. 9. This makes the neutron counts realistically 1000. We add the noise to each pixel after pixel and slit convolutions and then perform the subsampling and interpolation. Instead of phase stepping acquisitions, we take several single shot acquisitions and add them. In Refs. 8 and 9, phase steps and three frames were taken for each step. For our case, about 16 lines of noisy acquisitions were added together—this would be similar to eight steps at two frames. We compared a typical noisy recovery with the noiseless case. The parameters of the MPG system used are W = 50 or 120 µm, pitch = 2.0 µm, Dsg = 1 m, and Dsd = 3 m. The pixel size is 50 or 100 µm, and the slit width is 37.5 or 45 µm (convolution kernel of Ls × L2/L1). A 50% duty cycle G0 grating may be used instead of a single slit (see Sec. IV). The main detriments to the recovery were at the transitions due to the phase-sensitivity effects. For each phase object, the noisy recovery process described above was repeated 100 individual times (100 different seeds) and the root-mean-square-error (RMSE) and the standard deviation are reported. For system parameters of W = 120 µm, pixel size 50 or 100 µm, and a slit width of 45 µm, the average RMSE error was higher as expected from the smaller phase sensitivity (larger fringe period 360 µm at W = 120 µm as opposed to 150 µm for W = 50 µm).

FIG. 21.

Example phase recovery for a ramp, trapezoid, and a quadratic shaped object for noiseless and a typical with-noise case. Monochromatic beam(wavelength 0.44 nm) used. The acquisition parameters are MPG W = 50 µm, MPG pitch = 2.0 µm, Dsd = 3 m, and Dsg = 1 m. The object distance to detector Dod = 1.7 m. The pixel size is 50 µm, and the slit size 37.5 µm (slit kernel Ks = 37.5 µm × L2/L1 = 75 µm). A G0 with a 37.5 µm opening and 50% duty cycle may be used instead of a single slit. Realistic Poisson noise was added to the detector simulation. 16 noisy lines were added before phase recovery.

FIG. 21.

Example phase recovery for a ramp, trapezoid, and a quadratic shaped object for noiseless and a typical with-noise case. Monochromatic beam(wavelength 0.44 nm) used. The acquisition parameters are MPG W = 50 µm, MPG pitch = 2.0 µm, Dsd = 3 m, and Dsg = 1 m. The object distance to detector Dod = 1.7 m. The pixel size is 50 µm, and the slit size 37.5 µm (slit kernel Ks = 37.5 µm × L2/L1 = 75 µm). A G0 with a 37.5 µm opening and 50% duty cycle may be used instead of a single slit. Realistic Poisson noise was added to the detector simulation. 16 noisy lines were added before phase recovery.

Close modal

Table III shows the phase recovery RMSE over 100 noise realizations.

TABLE III.

MPG phase recovery errors for triangular, trapezoid, and paraboloid phase objects.

TriangularTrapezoidParaboloid
RMSE (±std)RMSE (±std)RMSE (±std)
MPGSlitPixelOver 100 seedsOver 100 seedsOver 100 seeds
WPSizeaSizeNoWithNoWithNoWith
(μm)(μm)Ls (μm)(μm)NoiseNoiseNoiseNoiseNoiseNoise
50 37.5 50 1.44 2.22 ± 0.86 1.88 2.48 ± 0.91 2.55 2.91 ± 0.93 
120 45 50 1.64 3.04 ± 1.20 2.82 3.84 ± 1.70 2.80 3.74 ± 1.86 
120 45 100 1.68 4.08 ± 1.74 2.89 4.89 ± 2.30 2.66 4.69 ± 2.37 
TriangularTrapezoidParaboloid
RMSE (±std)RMSE (±std)RMSE (±std)
MPGSlitPixelOver 100 seedsOver 100 seedsOver 100 seeds
WPSizeaSizeNoWithNoWithNoWith
(μm)(μm)Ls (μm)(μm)NoiseNoiseNoiseNoiseNoiseNoise
50 37.5 50 1.44 2.22 ± 0.86 1.88 2.48 ± 0.91 2.55 2.91 ± 0.93 
120 45 50 1.64 3.04 ± 1.20 2.82 3.84 ± 1.70 2.80 3.74 ± 1.86 
120 45 100 1.68 4.08 ± 1.74 2.89 4.89 ± 2.30 2.66 4.69 ± 2.37 
a

The convolution kernel used is Ks = 2Ls due to the pinhole magnification, L2/L1.

Our wave-propagation simulator for neutron PGMI applications can accept various designs (e.g., modulated phase gratings and standard-dual-grating). These analytical simulations use Sommerfeld–Rayleigh diffraction integrals (SRDI) to calculate the intensity of the diffracted neutron wave at the detector plane. We derived the theory behind the MPG and compared intensity and minima/maxima conditions with the SRDI. Fringe pattern blurring at the detector due to the spatial coherence loss of a slit source and pixel size resolution was accounted for by convolving the fringe pattern with independent appropriate rectangular window functions. Fringe patterns are also subsampled according to the pixel size of the detector and interpolated to increase the sample rate of their signal. Simulations of the standard-dual-grating system with a monochromatic configuration produced visibility results in good agreement with experimental data from Ref. 9 and fringe periods, which also agreed well with the theory given in Ref. 8. Note also that experiments include gratings with imperfect square profiles of gratings, a slight misalignment between gratings, coherence loss from air scattering over long distances, and mechanical vibrations.10,11 While our analytical SRDI simulations and fringe processing techniques do not account for these sources of noise, the fringe visibilities of our simulations still matched well to experimental values. The SRDI and fringe pattern processing techniques were used to investigate our MPG system, which demonstrated it to be a promising type of PGMI system by producing higher visibilities (50%–94.2%) than the standard-dual-grating system for a range of possible modulation period W parameters (Fig. 16). At the lower end of the modulation periods, a narrower slit can aid in yielding higher visibilities (Fig. 18). The margin of improvement the MPG system displays is high for most cases and appreciable from what is currently reported in the literature for the standard-dual-grating system. For example, W = 300 µm and pitch = 2 µm produce visibility of near 94.2% with a comparable SDD as the standard-dual-grating system operating in a monochromatic configuration. When that same MPG configuration was illuminated by a polychromatic beam, the visibility dropped to 68.9%, which is still higher than the best visibility case in Pushin et al.,9 and it demonstrates that our MPG system is potentially robust to polychromatic neutron sources. Operating in the smaller W range also improves the phase sensitivity, yielding as much as 10.0 × 103 for SDDs up to 3.5 m.

Our mathematical treatise of intensity and the equation for maximum visibility was borne out with N-SRDI simulations, as shown in Fig. 6 and Table I.

For a small fringe period, a smaller slit size is important due to the pinhole magnification amplifying the effects of slit blurring. For a smaller slit, a G0 grating will be useful. Therefore, one aspect to consider is the design of the G0 for the MPG system as shown in Ref. 14. If P0 = actual G0 pitch, the projected pitch at the detector will have the pinhole magnification P0 = Dgd/Dsg × P0 (or in other notation L2/L1 × P0). This should equate the fringe period W= M × W = Dsd/Dsg × W for coherence.14  So, Dgd/Dsg × P0 = Dsd/Dsg × W. Thus, P0 can be calculated from geometry for a given W as P0 = Dsd/Dgd × W. For clarity, it is emphasized that the magnification M = Dsd/Dsg is a different factor than the pinhole magnification effect, which is L2/L1 (that is Dgd/Dsg).

For a sample calculation, we consider the phase recovery example in Fig. 21. For this case, W = 50 µm, P = 2 µm, slit size = 37.5 µm for Dsg = 1 m, Dsd = 3 m (or Dgd = 2 m). The required P0 = Dsd/Dgd × W = 3/2 × 50 µm = 75 µm. The derivation can also be shown step by step. For W = 50 µm, the fringe period at the detector is W′ = M × W = Dsd/Dsg × W = 3 × 50 µm = 150 µm. The projected P0 has to be therefore 150 µm. Therefore, P0 = 75 µm. This makes the duty cycle Ls/P0 = 37.5/75 or 50%.

The performance improvement with MPG from the dual phase grating system comes with the added benefit of the MPG system being a single-phase-grating system that does not require precise alignment of two gratings.

A potential advantage of the MPG system compared to TLI (for x rays or neutron interferometry) is that the fluence absorbing analyzer grating between the object and detector that is necessary for the TLI is not necessary for MPG. The advantage is potentially more than just one less grating element and less alignment. For the same fluence at the detector for both systems (and for similar object positions relative to the source), the dose to the object will be less for the MPG system than for the TLI, where the analyzer absorbs half the x rays or neutrons. This is particularly important when imaging bio-samples in x ray or neutrons. The caveat is that source brilliance, exposure time, and source grating G0 have to be considered to get the same fluence at the detector. If the duty-cycle or open-ratios are the same, say 50% for both systems, then exposure time may be kept similar for similar fluence-level at the detector. However, for example, if the G0 has a 50% open ratio for the TLI and 40% for the MPG system, then the exposure time has to be increased by about 25% to obtain the same fluence at the detector.

The gratings can be manufactured using existing manufacturing techniques according to the lead manufacturer Microworks, GmbH.

For x-ray imaging (energies around 25 keV) we have obtained from Microworks GmbH rectangular MPGs with W = 120 µm, pitch P = 1 µm, as well as triangular MPGs with W = 120 µm and P = 1.8 µm. We have conducted successful experiments with both and are preparing a manuscript for submission of those results.26 

For neutrons, Si is the material of choice with deep UV gray tone lithography being one option to modulate etch depth in the silicon with locally different heights of the resist serving as etch stops.

Other modulation shapes such as triangular or sinusoidal MPG is possible.14 We observed in the simulation that the fringe pattern may have ringing and harmonics for the RECT modulated grating at some geometries. The quality of fringes may be better for an ideal sinusoidal modulation in the sense that there are fewer harmonics for sinusoidal MPG due to the smoother fall off. However, the RECT modulated grating is easier to manufacture and maintain quality control and is of lower cost. For these reasons, in this work, we limited our investigation to RECT MPG.

The 1D N-SRDI simulations were performed using the LSU LONI QB3 cluster. A single compute node was utilized, which has two Intel Xeon Platinum 8260 Processors each with 24 cores (2.4 GHz). OpenMP, an application program interface (API) for thread-based parallelism, is utilized to generate multiple threads that run in parallel in order to speed up calculation times when running simulations. On average, a single simulation can take about 30 core-hours to complete. We expect more complicated/2D geometries to be amenable to an MPI (Message Passing Interface) implementation for parallel computing since it would allow the use of multiple compute nodes.

We have proposed a novel phase-grating moiré interferometer system for cold neutrons that only requires a single modulated phase grating (MPG) for phase-contrast imaging, as opposed to the two or three phase gratings in previously employed PGMI systems. The MPG system promises to deliver significantly better fringe visibility relative to previous experiments in the literature that use a two-phase-grating-based PGMI. When operating in the smaller W range and using SDDs of up to 3.5 m, the MPG system could provide comparable or greater theoretical maximum phase sensitivity when compared to the conventional Talbot–Lau interferometer. The single MPG reduces the precise alignment requirements needed for multi-grating systems. Like other PGMI systems, the MPG system does not require the high-aspect-ratio absorption gratings used in Talbot–Lau interferometers, which are challenging to manufacture and reduce the neutron flux reaching the detector.

The work presented here is, in part, from the thesis work of the first author Ivan Hidrovo27 and from ongoing thesis work of Hunter Meyer28 (both advised by co-author J.D.), toward their medical physics M.Sc. degrees, Department of Physics and Astronomy, LSU, Baton Rouge, LA. This work was funded, in part, by NSF EPSCOR RII Track 4, Award No. 1929150. The simulations were conducted with high-performance computational resources provided by Louisiana State University and the Louisiana Optical Network Infrastructure.

Certain commercial equipment, instruments, or materials (or suppliers, software, …) are identified in this paper to foster understanding. Such identification does not imply recommendation or endorsement by the National Institute of Standards and Technology, nor does it imply that the materials or equipment identified are necessarily the best available for the purpose.

The authors have no conflicts to disclose.

I. Hidrovo: Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). J. Dey: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal). H. Meyer: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (supporting). D. S. Hussey: Methodology (supporting); Writing – review & editing (equal). N. N. Klimov: Writing – review & editing (supporting). L. G. Butler: Methodology (supporting); Writing – review & editing (supporting). K. Ham: Writing – review & editing (supporting). W. Newhauser: Writing – review & editing (supporting).

Data used in this work is available upon emailing the corresponding author, Dr. J. Dey.

1. Field amplitude derivation

See Fig. 1 for a general set up of the system. The Fresnel condition is derived by Patorski et al.21 for a standard periodic grating with period d in terms of the discrete Fourier coefficients. Since our grating has W and P components, we expect to have complicated periodic frequencies, with harmonic combinations of mW+nP. Hence, we derived the Fresnel conditions from first principles.28 First, we derive the amplitude of the field under Fresnel condition for MPG grating T(y) [see Eq. (6), reproduced here for convenience],
Ty1=gy1n=δy1nP+n=δy1nP2P/2×recty1P/2,
where gy1=expiϕ1recty1W/2+expiϕ2recty1W/2W/21Wδ(y1nW), the ϕ1 and ϕ2 being the phase heights in regions h1 and h2, each of width W/2 for the MPG of interest. The first curly bracket shows one period of the envelope that is repeated via the comp-convolution in the second curly bracket.
Here, the field in the Fresnel region is derived in Goodman [Eqs. (4)–(10) in Ref. 16],
Ux,y,z=exp(ikz)iλzexpik(x2+y2)2z×FTUx,y,zexpik(x2+y2)2zfx=xλz,fy=yλz,
(A1)
where Ux,y,z=expikr0r0TyexpikL1L1expik(x2+y2)2L1T(y), making a parabolic approximation of incident spherical wave from point source, as done in Patorski.21 
Substituting Ux,y,z, the FT-term in Eq. (A1) becomes
FTUx,y,zexpik(x2+y2)2z=FTT(y)FTexpikx2+y221L1+1z.
(A2)
The FT(Ty) of the particular MPG Ty [Eq. (6)] is given by
FTTy=P2sincPfy2F(gy)1Pn=δfynP+1Pexpj2πPfy2n=δfynP.
The P and 1P scaling terms cancel. Note the grating is assumed to be constant lines along the x-direction, which just introduces an overall δfx and leaves the function to convolute intact in the x-direction.
Then, carrying out the convolution in Eq. (A2) and plugging in fy=yλz and fx=xλz and simplifying, we obtain the amplitude in Eq. (A1) as
Uy,z=AexpiB1L1+zn=+m=+gmsinc12mPW+n×expjπλL1zL1+zmW+nP2expj2πL1yL1+zmW+nP+n=+expjπnsincn2expjπλL1zL1+znP2×expj2πL1yL1+znP,
(A3)
=AexpiB1L1+zU1y,z,W,P+U2y,z,P.
(A4)
Here, U1y,z,W,P is the first dual sum (depends on W and P) and U2y,z,P is the second single sum. Also, we lumped some real (and z-independent) and exponential phasor terms in Aexp(iB) as they will result in z-independent scaling term A2 in intensity Iy,z=Uy,zU*y,z=|Uy,z|2. Note that the gm are the Fourier coefficients of the envelope function, with period W,
gy=m=gmexpj2πmyWor,Fgy=m=gmδfymW,
where gm is found by Fourier transform integral of gy. Equation (A3) is reminiscent of Patorski21 or Guigay,18 except with a double summation due to the dual-harmonics. We are interested in Iy,z=|U(y,z)|2, which is what we observe in detector.

2. Field intensity and visibility

From Eq. (A4), we can write
Iy,z=AL1+z2U1U1*+U2U1*+U1U2*+U2U2*.
Writing each of U1 and U2 in terms of z- and y-dependent terms,
U1y,z=n=+m=+b1m,n,zexpj2πL1yL1+zmW+nP,
U2y,z=n=+b2m,n,zexpj2πL1yL1+znP,
where
b1m,n,z=gmsinc12mPW+nexpjπλL1zL1+zmW+nP2
and
b2m,n,z=expjπnsincn2expjπλL1zL1+znP2.
We can show
U1U1*=n=+m=+C11m,n,zexpj2πL1yL1+zmW+nP,
(A5)
where C11m,n,z=b1m,n,zb1*(m,n,z), which is the auto-correlation of the b1m,n,z. Note immediately that this component of the intensity (and, in fact, for the others) has expj2πL1yL1+znP, which is a modulation of the harmonics of P and cannot really be detected by the practical detector of 50100 µm without an analyzer grating as is typically used in Talbot-Lau Interferometers.
Similarly, U1U2*, U2U1*, and U2U2* can be written in equivalent forms to get
I(y,z)=AL1+z2n=+m=+C11z,m,n+C12z,m,n+C12*z,m,nexpj2πL1yL1+zmW+nP+C22(z,n)expj2πL1yL1+znP,
where C12z,m,n is the cross-correlation of b1m,n,z and b2n,z and C22(z, n) is the auto-correlation of b2n,z.

It is noted again that since this is the final intensity at the detector, the P-harmonics, i.e., nP will be blurred by pixel size. This is the reason why for a Talbot-Lau system with standard phase grating with pitch P (of the order of 1–4 µm), a second grating is needed (the G2, absorption grating) to observe the fringes with a practical detector of pixel size 50100 µm.

Therefore, we can consider only n = 0 and the intensity can be shown to be
Iy,zAL1+z2C22z,0+m=+C11z,m,0+C12z,m,0+C12*z,m,0expj2πL1yL1+zmW=AL1+z2I0+m=1+Im(z)cos2πL1yL1+zmW+θm.
(A6)
Therefore, the intensity produces fringe patterns of period mMW where magnification M = L1+zL1. In our case of an even grating, θm = 0. The Im=κm=2C11z,m,0+C12z,m,0+C12*z,m,0. The average is contributed by I0=C220+κ0.
The auto-correlations and cross-correlations are given by (for n = 0, i.e., not considering standalone n/P harmonic modulations, as explained earlier)
C11z,m,0=mnb1mm,nb1*m,n=mngmmgm*sincπ2mmPWn×sincπ2mP2Wn×expjπL1zL1+z2mnWP+2mmm2W2,
and
C12z,m,0=gmnsincπ2mPWnsincπn2expjπn×expjπλL1zL1+z2mnWPm2W2.
The visibility is guided by the first harmonic component that is 1MW that is for m = 1.
In other words, visibility is given by
Visibility=A1A0=2|C11z,1,0+C12z,1,0+C12*z,1,0|C22z,0+|C11z,0,0+C12z,0,0+C12*z,0,0|.
The visibility will maximize and minimize according to the exponentials in C11z,1,0andC12z,1,0. This is difficult to pin down, in general; however, setting m = 1 in the equations for C11z,m,0andC12z,m,0 ignoring 1W2 in comparison to 1WP terms, we can see both terms have expjπλL1zL1+z2nWP.
This may lead to a condition, 2πλL1zL1+z1WP=πk,
z=L1WPk2λL1WPk,
(A7)
where k′ is an integer.

We also observe that the location of Z is independent of the grating envelope function (though the visibility is).

Here, the effect of applying a shift to a parabolic wavefront source will be shown to follow pinhole-type magnification under the Fresnel approximation. First, the diffracted field amplitude as derived by Goodman (ignoring constant multipliers),16 
Ux,y,zUx,y,0expjk2z(x2+y2)×expjkz(xx+yy)dxdy.
(B1)
For a source located on the optical axis, the field at the grating under parabolic wavefront illumination is
Ux,y,0=T(y)ejkL1L1expjk2L1x2+y2,
where T(y’) is the grating transmission function. Shifting the source by y0 yields
Ux,y,0=TyejkL1L1expjk2L1x2+yy02.
Plugging this into (B1) and extracting the constants yields
Ux,y,zTyexpjk2L1x2+yy02×expjk2z(x2+y2)expjkz(xx+yy)dxdy.
Expanding and combining the exponentials yields
Ux,y,zexpjky022L1Tyexpjk2L1x2+y2+jk2z×x2+y2jkyy0L1jkz(xx+yy)dxdy,
which can then be arranged to
Ux,y,zexpjky022L1Tyexpjk2L1x2+y2+jk2zx2+y2jkzxx+yy+y0zL1dxdy.
It is immediately evident that (besides multiplicative factors), the field for a grating illuminated by a source shifted by y0 will result in a field that is equivalent but shifted (in the opposite direction) by y0zL1. This is corroborated by Fig. 5. It should be noted that Fig. 5 is generated using a true spherical wave (point source) for illumination, and the corroboration further validates the parabolic wave approximation used in  Appendix A.
1.
F.
Pfeiffer
,
C.
Grünzweig
,
O.
Bunk
,
G.
Frei
,
E.
Lehmann
, and
C.
David
,
Phys. Rev. Lett.
96
,
215505
(
2006
).
2.
K. W.
Herwig
, “
Introduction to the neutron
,” in
Neutron Imaging and Applications. A Reference for the Imaging Community
, edited by
H. Z.
Bilheux
,
R.
McGreevy
, and
I.
Anderson
(
Springer
,
Boston, MA
,
2009
), pp.
3
12
.
3.
H.
Rauch
and
S. A.
Werner
,
Neutron Interferometry: Lessons in Experimental Quantum Mechanics
(
Oxford University Press
,
Oxford
,
2000
).
4.
C.
Grünzweig
,
C.
David
,
O.
Bunk
,
M.
Dierolf
,
G.
Frei
,
G.
Kühne
,
J.
Kohlbrecher
,
R.
Schäfer
,
P.
Lejcek
,
H. M. R.
Rønnow
et al,
Phys. Rev. Lett.
101
,
025504
(
2008
).
5.
B.
Betz
,
R. P.
Harti
,
M.
Strobl
,
J.
Hovind
,
A.
Kaestner
,
E.
Lehmann
,
H.
Van Swygenhoven
, and
C.
Grünzweig
,
Rev. Sci. Instrum.
86
,
123704
(
2015
).
6.
T.
Neuwirth
,
A.
Backs
,
A.
Gustschin
,
S.
Vogt
,
F.
Pfeiffer
,
P.
Böni
, and
M.
Schulz
,
Sci. Rep.
10
,
1764
(
2020
).
7.
Y.
Kim
,
J.
Kim
,
D.
Kim
,
D. S.
Hussey
, and
S. W.
Lee
,
Rev. Sci. Instrum.
90
,
073704
(
2019
).
8.
H.
Miao
,
A.
Panna
,
A. A.
Gomella
,
E. E.
Bennett
,
S.
Znati
,
L.
Chen
, and
H.
Wen
,
Nat. Phys.
12
,
830
(
2016
).
9.
D. A.
Pushin
,
D.
Sarenac
,
D. S.
Hussey
,
H.
Miao
,
M.
Arif
,
D. G.
Cory
,
M. G.
Huber
,
D. L.
Jacobson
,
J. M.
LaManna
,
J. D.
Parker
et al,
Phys. Rev. A
95
,
043637
(
2017
).
10.
D.
Sarenac
,
D. A.
Pushin
,
M. G.
Huber
,
D. S.
Hussey
,
H.
Miao
,
M.
Arif
,
D. G.
Cory
,
A. D.
Cronin
,
B.
Heacock
,
D. L.
Jacobson
et al,
Phys. Rev. Lett.
120
,
113201
(
2018
).
11.
B.
Heacock
,
D.
Sarenac
,
D. G.
Cory
,
M. G.
Huber
,
D. S.
Hussey
,
C.
Kapahi
,
H.
Miao
,
H.
Wen
, and
D. A.
Pushin
,
AIP Adv.
9
,
085115
(
2019
).
12.
M.
Strobl
,
J.
Valsecchi
,
R. P.
Harti
,
P.
Trtik
,
A.
Kaestner
,
C.
Gruenzweig
,
E.
Polatidis
, and
J.
Capek
,
Sci. Rep.
12
,
3461
(
2019
).
13.
J.
Dey
,
N.
Bhusal
,
L.
Butler
,
J. P.
Dowling
,
K.
Ham
, and
V.
Singh
, “
Phase contrast x-ray interferometry
,” U.S. patent 10872708 (
22 December 2020
).
14.
J.
Xu
,
K.
Ham
, and
J.
Dey
, “
X-ray interferometry without analyzer for breast CT application: A simulation study
,”
J. Med. Imaging
7
(
2
),
023503
(
2020
).
15.
A.
Pandeshwar
,
M.
Kagias
,
Z.
Shi
, and
M.
Stampanoni
, “
Envelope modulated x-ray grating interferometry
,”
Appl. Phys. Lett.
120
,
193701
(
2022
).
16.
J. W.
Goodman
,
Introduction to Fourier Optics
,
2nd ed.
(
McGraw-Hill
,
New York
,
1996
).
17.
J. P.
Wilde
and
L.
Hesselink
, “
Statistical optics modeling of dark-field scattering in x-ray grating interferometers: Part 1. Theory
,”
Opt. Express
29
,
40891
40916
(
2021
).
18.
J. P.
Guigay
, “
On Fresnel diffraction by one-dimensional periodic objects, with application to structure determination of phase objects
,”
Opt. Acta: Int. J. Opt.
18
(
9
),
677
682
(
1971
).
19.
T.
Weitkamp
,
C.
David
,
C.
Kottler
,
O.
Bunk
, and
F.
Pfeiffer
, “
Tomography with grating interferometers at low-brilliance sources
,”
Proc. SPIE
6318
,
63180S
(
2006
).
20.
A.
Macovski
,
Medical Imaging Systems
(
Prentice Hall, Inc.
,
NJ
,
1983
).
21.
K.
Patorski
, “
I the self-imaging phenomenon and its applications
,”
Prog. Opt.
27
,
1
108
(
1989
).
22.
T.
Donath
,
M.
Chabior
,
F.
Pfeiffer
,
O.
Bunk
,
E.
Reznikova
,
J.
Mohr
,
E.
Hempel
,
S.
Popescu
,
M.
Hoheisel
,
M.
Schuster
et al,
J. Appl. Phys.
106
,
054703
(
2009
).
23.
N.
Bevins
,
J.
Zambelli
,
K.
Li
,
Z.
Qi
, and
G. H.
Chen
,
Med. Phys.
39
,
424
428
(
2012
).
24.
Y.
Kim
,
J.
Valsecchi
,
J.
Kim
,
S. W.
Lee
, and
M.
Strobl
,
Sci. Rep.
9
,
18973
(
2019
).
25.
J.
Kim
,
S. W.
Lee
, and
G.
Cho
,
Nucl. Instrum. Methods Phys. Res., Sect. A
746
,
26
32
(
2014
).
26.
H.
Meyer
,
J.
Dey
S.
Carr
,
K.
Ham
,
L.
Butler
,
I.
Hidrovo
,
T.
Varga
,
J.
Schulz
,
T.
Beckenbach
, and
K.
Kaiser
, “
X-ray interferometry with a modulated phase grating: Theory and experiments
”, (unpublished).
27.
I. J.
Hidrovo Giler
, “
Neutron interferometry using a single modulated phase grating
,” M.S. thesis, (
Louisiana State University
,
2022
).
28.
H.
Meyer
, “
The theory of x-ray and neutron interferometry using a modulated phase grating
,” M.S. thesis, (
Louisiana State University
2023
).