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.

## I. INTRODUCTION

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 ($\Delta \lambda /\lambda $. 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}

## II. METHODS

### A. Analytical simulations for neutron interferometry

^{14}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*(

*P*

_{s})

*.*The neutron field amplitude at a detector for a single-phase-grating system, such as the MPG (Fig. 1), can be expressed as

*U*(

*P*

_{s}) 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\pi \lambda $ is the wave number,

*λ*is the wavelength,

*r*

_{1}is the distance between the grating point

*y*

_{1}to the detector point

*y*given by $r1=Dgd2+y\u2212y12$, and

*θ*is the angle between $r1\u20d7$ and the normal of the MPG. The transmission function is given by

*A*

_{1}(

*y*

_{1}) is the amplitude transmission of the neutron beam through the grating due to attenuation and

*ϕ*(

*y*

_{1}) is the phase shift determined by the spatial heights of the MPG.

*y*

_{1}along the phase grating.

^{16}In all simulations, the neutron source wave function

*U*(

*P*

_{s}) is assumed to be a point source at point

*y*

_{0}; however, the fringe patterns from a line source are later acquired with a convolution method. Therefore, Eq. (1) was simplified to

*r*

_{0}is the distance between the point source at point

*y*

_{0}and the point

*y*

_{1}on the MPG given by $r0=Dsg2+y1\u2212y02$ and the cos

*θ*(between

*r*

_{1}and the

*z-*axis) has been replaced by

*D*

_{gd}/

*r*

_{1}, explaining the $r12$ term inside the SRDI shown in Eq. (4).

*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

In all simulations (MPG or standard-dual-grating systems), we obtained results with a point source [that is ** U**(

*P*_{s}) was a point source at point

*y*_{o}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.

### B. Fringe pattern blurring due to slit width and pixel size via convolution method

#### 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.

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 *Y*_{s}-axis of Fig. 3 with center separation multiples of *P*_{0}. The *P*_{0} (or source period *P*_{s} 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 *P*_{0} and slit opening are available in the literature.^{14,19} In the Discussion section, we calculated the pitch *P*_{0} 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 $\lambda =0.44nm$) exiting a 200 *µ*m slit source with the system parameters of *L*_{1} = 1.2 m, *D* = 12 mm, and *L*_{2} = 1.78 m (Fig. 2). The identical phase gratings had pitch $PG1=PG2=2.4\mu 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(\pi s/Ps)$, where *V*_{o} is the optimal fringe visibility from one point source, s is the width of the slit, and *P*_{s} 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.

### C. Evaluation of N-SRDI simulations for a standard-dual-grating system, comparisons with simulations and experiments

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 *D*_{sd}) and *L*_{1} (or *D*_{sg}) fixed (Fig. 2). The *L*_{2} (or *D*_{gd}) 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 G_{0} has a direct one-to-one effect on the fringe periods at the detector without a pinhole geometric magnification when P_{G1} = P_{G2}. Thus, the source period (or the pitch of G_{0}) is directly equal to the fringe period as shown in Refs. 8 or 9.

### D. Theory of modulated phase grating (MPG) and comparison with N-SRDI MPG simulations

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 L_{1} = 100 cm and MPG to detector distance is L_{2} = 200 cm, (L_{2}/L_{1} = 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, L_{2}/L_{1} = 2. In general, if the actual slit size is L_{s} *μ*m, the convolution kernel to be used is K_{s} = L_{s} (L_{2}/L_{1}).

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.

*y*-axis as the axis along the detector, the

*y*

_{1}-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

*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*(

*iφ*).

For the RECT MPG discussed in this paper, the envelope function is given by $gy1=expi\varphi 1recty1W/2+expi\varphi 2recty1\u2212W/2W/2\u22971W\u2211n=\u2212\u221e\u221e\delta (y1\u2212nW)$, the *ϕ*_{1} and *ϕ*_{2} being the phase heights in regions *h*_{1} and *h*_{2}, 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(*y*_{1}) from first principles outlined in Goodman.^{16}

*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,

*U*

_{1}, the frequencies are related to $mW+nP$ and in the second term

*U*

_{2}to $nP$.

*C*

_{ij}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.

#### 3. Conditions for maximum/minimum visibility

*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.

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 | 2 | 4.5 | 4.5 (for k = 3) | ⋯ | ⋯ |

120 | 1 | 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,* *g*_{m}*) 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 *h*_{2} 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 *h*_{1} such that it yields π-shift and varied the *h*_{2} to provide 0, π/4, π/2, and 3π/4 phase shifts (for $\lambda =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 *h*_{2} 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 *L*_{s} = 75 *µ*m with the magnified convolution kernel *K*_{s} = *L*_{s} × L_{2}/L_{1} were applied for each point.

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 *h*_{2} = 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 *h*_{2} = 0 for all MPG–detector distances, henceforth in this work we used (π, 0), as shown in the schematic diagram Fig. 8.

##### 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 *S*_{w}, source-to-MPG distances *D*_{sg}, and MPG-to-detector distances *D*_{gd}, as shown in Sec. III.

### E. N-SRDI simulations for a modulated phase grating (MPG) neutron interferometer illuminated by a polychromatic beam

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.

### F. Single-shot phase contrast recovery with MPG

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.

*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}

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 *D*_{od} 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/(2

*W*′) 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

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\lambda \u22121$. We show this by combining Eqs. (10) and (11) to obtain

*S*=

*D*

_{od}/

*W*′ is the phase sensitivity for the interferometer. Expressed in terms of the object’s differential phase profile

*d*Φ/

*dy*,

*ϕ*must be “corrected” by 1/(

*λS*) to obtain Φ/

*dy*.

**Step 2(b)**. To obtain the object’s differential phase *d*Φ/*dy*_{ob} at the object space (*y*_{ob}), we have to scale (de-magnify) the spatial variation of the detector (y) by 1/*M*_{obj}, where *M*_{obj} = object magnification.

**Step 3.** We integrate *d*Φ/*dy*_{ob} 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\Phi yob$, which are applied to the wave field at the object-to-detector distance, where Φ(*y*_{ob}) is the object phase profile.

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 $\mu 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.

Fringe pattern Δy
. | Left 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 Δy
. | Left 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\lambda \u22121$ 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 (*y*_{ob}) instead of the detector grid space (y). Then, Fig. 10(d) shows the integrated phase of the object (**Step 3**).

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*Φ/*dy*_{ob} [**Step 2(b)**] on the object grid space (*y*_{ob}). They both approximately show the rectangular and linear *d*Φ/*dy*_{ob} expected for triangular and parabolic phase profiles, respectively.

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.

## III. RESULTS

### A. Evaluation of N-SRDI simulations for a standard-dual-grating system, comparisons with theory and experiments

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.

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.

### B. N-SRDI simulations for a modulated phase-grating (MPG) neutron interferometer illuminated by a monochromatic beam

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).

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}

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 D_{sg} 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 L_{s} = 200 *µ*m with convolution kernel varying as K_{s} = L_{s} × L_{2}/L_{1.}

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 *S*_{w} (albeit knowing that the flux reaching the detector would be reduced proportionately and typically a G_{0} grating would be required instead of a slit). We kept the *S*_{w} (or *L*_{s}) = 200 *µ*m as a reference. As expected in general, we see in Figs. 18(a)–18(c), using a narrower slit width *S*_{w} improved the visibility for the W = 50, 70, and 120 *µ*m for these fringe patterns with smaller periods.

Finally, we calculated the maximum phase sensitivities of our MPG system. The maximum phase sensitivity is given by *S*_{max} = *D*_{gd}/*W*′ (where the object is closest to the MPG, i.e., *D*_{od} = *D*_{gd}). 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 × 10^{3} to 10.0 × 10^{3} (0.41–1 × 10^{4}) 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 × 10^{3} with a SDD of 3.5 m and a beam wavelength of 0.44 nm.^{7,25}

### C. N-SRDI simulations for a modulated phase grating (MPG) neutron interferometer illuminated by a polychromatic beam

We used a polychromatic beam for the setup which produced the best visibility of 94.2% in the monochromatic simulations ($\lambda =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.

### D. Single-shot phase contrast recovery with MPG

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, D_{sg} = 1 m, and *D*_{sd} = 3 m. The pixel size is 50 or 100 *µ*m, and the slit width is 37.5 or 45 *µ*m (convolution kernel of L_{s} × L_{2}/L_{1}). 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).

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

. | . | . | Triangular . | Trapezoid . | Paraboloid . | ||||
---|---|---|---|---|---|---|---|---|---|

. | . | . | RMSE (±std) . | RMSE (±std) . | RMSE (±std) . | ||||

MPG . | Slit . | Pixel . | Over 100 seeds . | Over 100 seeds . | Over 100 seeds . | ||||

W . | P . | Sizea . | Size . | No . | With . | No . | With . | No . | With . |

(μm)
. | (μm)
. | Ls (μm)
. | (μm)
. | Noise . | Noise . | Noise . | Noise . | Noise . | Noise . |

50 | 2 | 37.5 | 50 | 1.44 | 2.22 ± 0.86 | 1.88 | 2.48 ± 0.91 | 2.55 | 2.91 ± 0.93 |

120 | 2 | 45 | 50 | 1.64 | 3.04 ± 1.20 | 2.82 | 3.84 ± 1.70 | 2.80 | 3.74 ± 1.86 |

120 | 2 | 45 | 100 | 1.68 | 4.08 ± 1.74 | 2.89 | 4.89 ± 2.30 | 2.66 | 4.69 ± 2.37 |

. | . | . | Triangular . | Trapezoid . | Paraboloid . | ||||
---|---|---|---|---|---|---|---|---|---|

. | . | . | RMSE (±std) . | RMSE (±std) . | RMSE (±std) . | ||||

MPG . | Slit . | Pixel . | Over 100 seeds . | Over 100 seeds . | Over 100 seeds . | ||||

W . | P . | Sizea . | Size . | No . | With . | No . | With . | No . | With . |

(μm)
. | (μm)
. | Ls (μm)
. | (μm)
. | Noise . | Noise . | Noise . | Noise . | Noise . | Noise . |

50 | 2 | 37.5 | 50 | 1.44 | 2.22 ± 0.86 | 1.88 | 2.48 ± 0.91 | 2.55 | 2.91 ± 0.93 |

120 | 2 | 45 | 50 | 1.64 | 3.04 ± 1.20 | 2.82 | 3.84 ± 1.70 | 2.80 | 3.74 ± 1.86 |

120 | 2 | 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 = 2*L*_{s} due to the pinhole magnification, L_{2}/L_{1}.

## IV. DISCUSSION

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 × 10^{3} 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 *G*_{0} grating will be useful. Therefore, one aspect to consider is the design of the *G*_{0} for the MPG system as shown in Ref. 14. If P_{0} = actual *G*_{0} pitch, the projected pitch at the detector will have the pinhole magnification $P0\u2032$ = D_{gd}/D_{sg} × P_{0} (or in other notation L_{2}/L_{1} × P_{0}). *This should equate the fringe period* *W*′ *= M* × *W = D*_{sd}*/D*_{sg} × *W for coherence*.^{14} *So,* D_{gd}/D_{sg} × P_{0} = D_{sd}/D_{sg} × W*. Thus, P0 can be calculated from geometry for a given W as* *P*_{0} **= D**

_{sd}

*/D*_{gd}×

**. For clarity, it is emphasized that the magnification M = D**

*W*_{sd}/D

_{sg}is a different factor than the pinhole magnification effect, which is L

_{2}/L

_{1}(that is D

_{gd}/D

_{sg}).

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 D_{sg} = 1 m, D_{sd} = 3 m (or D_{gd} = 2 m). The required **P**_{0} = **D**_{sd}**/D**_{gd} × **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 = D_{sd}/D_{sg} × W = 3 × 50 *µ*m = 150 *µ*m. The projected $P0\u2032$ has to be therefore 150 *µ*m. Therefore, P_{0} = 75 *µ*m. This makes the duty cycle *L*_{s}/*P*_{0} = 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 G_{0} 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 G_{0} 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.

## V. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

The work presented here is, in part, from the thesis work of the first author Ivan Hidrovo^{27} and from ongoing thesis work of Hunter Meyer^{28} (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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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 AVAILABILITY

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

### APPENDIX A: DERIVATION OF INTENSITY AND VISIBILITY FOR A MODULATED PHASE GRATING

#### 1. Field amplitude derivation

*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],

^{21}

*x*-direction, which just introduces an overall $\delta fx$ and leaves the function to convolute intact in the

*x*-direction.

*Aexp*(

*iB*) as they will result in z-independent scaling term

*A*

^{2}in intensity $Iy,z=Uy,zU*y,z=|Uy,z|2$. Note that the

*g*

_{m}are the Fourier coefficients of the envelope function, with period W,

*g*

_{m}is found by Fourier transform integral of $gy$. Equation (A3) is reminiscent of Patorski

^{21}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

*U*

_{1}and

*U*

_{2}in terms of z- and y-dependent terms,

**of the $b1m,n,z$.**

*auto-correlation**Note immediately that this component of the intensity (and, in fact, for the others) has*$expj2\pi L1yL1+znP$,

*which is a modulation of the harmonics of P and cannot really be detected by the practical detector of*50

*–*100

*µ*

*m without an analyzer grating as is typically used in Talbot-Lau*Interferometers.

*C*

_{22}(

*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 50*–*100 *µ*m*.*

*n =*0 and the intensity can be shown to be

*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=\kappa m=2\u2223C11z,m,0+C12z,m,0+C12*z,\u2212m,0\u2223$. The average is contributed by $I0=C220+\kappa 0$.

*k*′ is an integer.

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

### APPENDIX B: SOURCE SHIFT DERIVATION

^{16}

*y*

_{0}yields

*y*

_{0}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.

## REFERENCES

*Neutron Imaging and Applications. A Reference for the Imaging Community*

*Neutron Interferometry: Lessons in Experimental Quantum Mechanics*

*.*