1. INTRODUCTION

The so-called multi-slice (MS) approximation has been widely used to perform wave optical simulations of beams propagating through diffracting structures too large for the projection approximation to be valid (e.g., [1–6]) and is thus well established. The MS approximation instead divides the diffracting structure into a partition of slices to which the projection approximation may be applied individually. Algorithms employing the MS approximation generally use the theory of Fourier optics [7] to propagate beams across a slice between planes. Numerical implementations of such algorithms use the discrete Fourier transform (DFT), which makes such algorithms computationally tractable. Careful attention must be paid to sampling in DFT-based beam propagation techniques, which can become prohibitive when modeling divergent beams, since correct sampling must be maintained in both the spatial and reciprocal spaces throughout the entire simulated volume. This challenge can be exacerbated when modeling diffracting objects that may not be properly represented on a sampled grid. Despite the large number of published studies that employ the MS approximation, I am unaware of any that have analyzed, or proposed a solution to, the challenge posed by sampling divergent beams in a MS simulation containing diffracting objects. A method of transforming a divergent-beam geometry into a plane-wave geometry for single-slice simulation has been demonstrated [6,8,9]. In this paper, I provide a full description of how to integrate this transformation into a MS simulation in order to calculate how a coherent divergent beam propagates through a diffracting sample extended in the axial direction. Furthermore, I provide a thorough analysis of the sampling requirements of this new approach and show how it can be implemented in such a way that guarantees that correct sampling is maintained throughout the simulation.

I begin this paper with an overview of preliminary theory required to develop the simulation technique, including angular spectrum approaches to propagating beams within a system described by the paraxial wave equation. The sampling requirements of this technique, in both the spatial and reciprocal spaces, are reviewed before describing how diffracting objects are treated. I then introduce the divergent-wave-to-plane-wave transformation, which allows for significant relaxation of the sampling requirements. It is then shown how this transformation can be applied in the MS approximation, and I provide details of its numerical implementation and sampling requirements. I conclude the paper with a series of examples that provide verification of the new simulation method and illustrate its usefulness.

A. Preliminary Theory

I consider a system as depicted in Fig. 1, where a monochromatic point source is located at the origin of the global coordinate system. Extended sources that are spatially incoherent can be represented by collections of point sources, and polychromatic sources may be modeled by incoherently superimposing simulations performed at wavelengths throughout the source spectrum. Although beyond the scope of this paper, sources with more general states of coherence may be modeled by performing one simulation for each mode of the source’s coherent mode decomposition. I assume that free space exists for $0\le z\le \mathrm{\Delta }{z}_{\mathrm{so}}$, where I use the subscript so to denote source-to-object distance. I also assume that some degree of refractive index inhomogeneity exists in the region $z\ge \mathrm{\Delta }{z}_{\mathrm{so}}$ and that our ultimate goal is to calculate the field in the plane $z=\mathrm{\Delta }{z}_{\mathrm{so}}+\mathrm{\Delta }{z}_{\mathrm{od}}$, where $\mathrm{\Delta }{z}_{\mathrm{od}}$ is the object-to-detector distance. I thus refer to the space defined by $\mathrm{\Delta }{z}_{\mathrm{so}}\le z\le \mathrm{\Delta }{z}_{\mathrm{so}}+\mathrm{\Delta }{z}_{\mathrm{od}}$ as the computational region. I describe fields propagating in this system according to the form $u(x,y,z)\exp (ikz-i\omega t)$, where $k$ is the wavenumber, $\omega$ is the angular frequency of the radiation emitted by the source, and it is assumed that $u(x,y,z)$ varies only weakly with $z$. If $u(x,y,z)\exp (ikz-i\omega t)$ is substituted into the time-harmonic scalar wave equation, it is not difficult to show that the paraxial wave equation is obtained as

 

Fig. 1. Schematic diagram of the system studied in this paper. A point source is located at the origin of the Cartesian coordinate system. Refractive index inhomogeneities may exist for $z>\mathrm{\Delta }{z}_{\mathrm{so}}$, one example of which is illustrated in a layer bounded by the planes $z=z_i$ and $z=z_{i+1}=z_i+\mathrm{\Delta }{z}_i$.

Download Full Size | PDF

The spherical wave emitted by the point source has complex amplitude in the free-space region of $\exp (ikr)/r$, where $r=\sqrt{x^2+y^2+z^2}$. If I make the Fresnel approximation to this field and drop the $\exp (ikz)$ dependence, I obtain the expression

Our objective in this paper is to employ DFT techniques to determine the complex amplitude that emerges from the plane $z=\mathrm{\Delta }{z}_{\mathrm{so}}+\mathrm{\Delta }{z}_{\mathrm{od}}$ when refractive index inhomogeneities may exist in the region $z\ge \mathrm{\Delta }{z}_{\mathrm{so}}$. The principal limitation when using DFT-based techniques is that the field must at all times be sampled in a manner that satisfies the Nyquist criterion. When considering spherical waves, impractically dense sampling requirements may result. Since the main aim of this paper is to avoid such dense sampling, I first show how these dense sampling requirements arise from Eq. (5). In particular, if one writes $u_{\mathrm{inc}}$ as

The final aspect of preliminary theory deals with inhomogeneous refractive index distributions, which I treat using the so-called projection approximation [6,10]. I explain this with the aid of Fig. 1, which contains a position-dependent refractive index distribution, within the region $z_i\le z\le z_{i+1}$, described by

 

Fig. 2. Diagram illustrating the principal coordinate system notation used in this paper.

Download Full Size | PDF

B. Divergent-Beam-to-Plane-Wave Transformation

The sampling requirements expressed in Eqs. (18) and (19) for modeling divergent spherical waves cannot be achieved simultaneously. However, they can be circumvented by applying a coordinate system transformation that transforms the divergent spherical wave into a plane wave [6,8,9]. Although this subject is treated in detail by Paganin [6], I follow here the notation introduced by Sziklas and Siegman [8,9]. I introduce this transformation by considering the problem depicted in Fig. 1, whereby the complex amplitude incident upon the plane $z=z_i$ is known, which I denote by $u(x,y,z_i^{-})$. Consider for now that the space $z_i\le z\le z_{i+1}$ is composed only of homogeneous space, i.e., there is no refractive index inhomogeneity present within the slice. I note that two coordinate systems are depicted in this diagram: the global coordinate system $(x,y,z)$ and a coordinate system $(x_i^{\prime },y_i^{\prime },z_i^{\prime })$, valid only for $z_i\le z\le z_{i+1}$. In order to evaluate $u(x,y,z_{i+1}^{-})$, I first perform a transformation into a primed coordinate system defined by [8,9]

Having obtained $v(x_i^{\prime },y_i^{\prime },z_i^{\prime }(z_i)=0)$ from Eq. (22), it remains to calculate $v(x_i^{\prime },y_i^{\prime },z_i^{\prime }(z_{i+1}))$ using angular spectrum propagation. I begin by calculating the angular spectrum of $v(x_i^{\prime },y_i^{\prime },0)$ using Eq. (2), thus giving $V(a_i,b_i,0)$, which can be propagated a distance $z_i^{\prime }(z_{i+1})$ according to Eq. (4) as

C. Simulation of an Inhomogeneous Refractive Index Distribution

I use the projection approximation discussed in Section 1.A to model refractive index inhomogeneities present within the slice $z_i\le z\le z_{i+1}$. The projection approximation can be applied in either the global or primed coordinate systems due to the linearity of Eq. (22). I opt to apply the projection approximation directly in the primed coordinate system, i.e., I apply it directly to $v(x_i^{\prime },y_i^{\prime },z_i^{\prime }(z_{i+1}))$, where I have dropped the dependence of $x_i^{\prime }$ and $y_i^{\prime }$ on $(x,z)$ and $(y,z)$, respectively, for brevity. It is important to note, however, that I must transform the argument of the projection function [Eq. (21)] from the global to the primed coordinate system. In particular, the projection approximation must be applied as $v(x_i^{\prime },y_i^{\prime },z_i^{\prime }(z_{i+1}))\exp (i\varphi (M_i{x}_i^{\prime },M_i{y}_i^{\prime }))$ when operating in the primed coordinate system. While this is evident from the transformation expressed in Eqs. (22)–(25), this may be understood intuitively by noting that the lateral cross section of an object must effectively be de-magnified in the primed coordinate system to compensate for the lack of divergence in the incident wavefront. However, despite this de-magnification, its optical thickness must remain the same in both coordinate systems.

D. Mitigating Aliasing Due to Use of Discrete Fourier Transforms

Even after applying the divergent-to-plane-wave transformation, aliasing, due to the use of DFTs, may still result due to the presence of diffracting refractive index inhomogeneities. In particular, the projection function $\exp (i\varphi (x,y))$ will not, in general, be band limited. Thus, if a complex amplitude in the absence of a diffracting object, $v(x^{\prime },y^{\prime },z^{\prime })$, is band limited, the complex amplitude after the application of the projection approximation $v(x^{\prime },y^{\prime },z^{\prime })\exp (i\varphi (M_i{x}_i^{\prime },M_i{y}_i^{\prime }))$ will not be, in general. I suggest an approach to overcoming this problem by calculating a band-limited projection function. This approach is possible only for diffracting objects for which the Fourier transform of its projection function can be calculated either analytically or numerically to high precision. Supposing I define a projection function $t(x,y)=\exp (i\varphi (x,y))$, I define the band-limited version of $t(x,y)$ as

It is worth noting that for many applications, this step may be omitted without significantly perturbing the calculation of the transmitted complex amplitude, as many other implementations of the MS method do. This approach is presented here as a means of eliminating aliasing for applications where this is important.

2. CALCULATION OF A FIELD PROPAGATING THROUGH AN ARBITRARY NUMBER OF SLICES

Consider the geometry depicted in Fig. 2, where the space $z_0\le z\le z_{N_s}$ is divided into $N_s$ slices. Propagation from the plane $z=z_0$ to the plane $z=z_{N_s}$ can be achieved by performing $N_s$ individual propagation calculations between the planes indicated in Fig. 2. When using the divergent-beam-to-plane-wave transformation, it is important to note that a different transformation is required for each slice. The details of these transformations are indicated in Fig. 2. I denote the primed coordinate system in a particular slice by the subscript associated with the plane closest to the source. For example, the primed coordinate system corresponding to the region $z_i\le z\le z_{i+1}$ is denoted $(x_i^{\prime },y_i^{\prime },z_i^{\prime })$. Propagation through the first slice (i.e., the slice defined by $z_0\le z\le z_1$) is achieved by evaluating $v(x_0^{\prime },y_0^{\prime },0)$ according to Eqs. (22)–(25), giving the field in the primed coordinate system at $z=z_0$. Under this transformation, the primed and global coordinate systems coincide at $z=z_0$, i.e., $x_0^{\prime }(x,z_0)=x$ and $y_0(y,z_0)=y$. Propagation to $z=z_1$ would be achieved by multiplying the angular spectrum associated with $v(x_0^{\prime },y_0^{\prime },0)$ (i.e., $V(a_0,b_0,0)$) according to

A. Outline of Algorithm

All of the elements of the MS algorithm introduced in this paper have now been introduced, and the algorithm can be explained in its entirety. Assuming a partition of slices as illustrated in Fig. 2, the complex amplitude due to a monochromatic point source is evaluated analytically on the sampled grid at $z=z_0$, which is immediately transformed into the primed coordinate system according to Eqs. (22)–(25), yielding $\widehat{v}(x_0^{\prime },y_0^{\prime },z_0^{\prime }(z_0))$. This field is the base case of a recursive definition of the algorithm, which propagates the field in the primed coordinate system to the end of the final slice given as

B. Numerical Implementation

The simulation methodology outlined in this paper is designed to be implemented using DFTs. I assume that a divergent spherical wave is incident upon the plane $z=z_0$. After applying the coordinate system transformation into the primed coordinate system, the spatially sampled grid must be defined as indicated in Eq. (12). By noting that at $z=z_0$ I have $(x_0^{\prime },y_0^{\prime })=(x,y)$, the sampled grid in the primed coordinate system is defined as

C. Sampling Requirements

I shall define the sampling requirements with reference to a cutoff frequency $f_{\mathrm{co}}$, such that the windowing function $W(f_x,f_y)$, assumed to be separable in $f_x$ and $f_y$, takes the value of 0 for $|f_x|>f_{\mathrm{co}}$ and $|f_y|>f_{\mathrm{co}}$. From this point on, I shall assume that the spatial, and therefore reciprocal, computational grids are entirely isotropic. The first requirement that $\mathrm{\Delta }x=\mathrm{\Delta }y$ must satisfy is

 

Fig. 3. Illustration of the region in reciprocal space (the intersection of $|f_x|<f_{\mathrm{co}}$ and $|f_y|<f_{\mathrm{co}}$) where the angular spectrum of the field is permitted to reside, along with a guard band where angular spectrum content is filtered after propagation through each layer.

Download Full Size | PDF

D. Non-Uniform Background Refractive Indices

Up until this point, for clarity, it has been assumed that all scatterers are contained within free space. It is, however, relatively simple to allow each slice defined by $z_i\le z\le z_{i+1}$ to have a background refractive index $n_i=1-{\delta }_i+i{\beta }_i$. If the slice contains a refractive index homogeneity as in Eq. (20), this should be redefined as

3. EXAMPLES AND ANALYSIS

Code for performing each of the following examples may be freely downloaded [12].

A. Diffracting Aperture

I consider first the simple example of a square diffracting aperture and a monochromatic point source as depicted in Fig. 4, which shows a point source placed 1.6 m from a square aperture of width $W=20\mathrm{\mu m}$. The field is observed on the observation plane, which is located a further 0.3 m from the aperture. The attenuating part of the aperture is assumed to be perfectly attenuating, while the transmitting region of the aperture is assumed to be perfectly transmitting. The point source is assumed to be monochromatic with photon energy of 20 keV.

 

Fig. 4. Schematic diagram of the diffracting aperture and point source upon which all examples in Section 3 are based. Values of $W=20\mathrm{\mu m}$, $z_a=1.6\mathrm{m}$, and $z_o=1.9\mathrm{m}$ were chosen. Note that the diagram is not drawn to scale.

Download Full Size | PDF

The objective of this test is to demonstrate how the choice of $N$, the number of sample points along each Cartesian direction, impacts upon the accuracy of the simulation method. In particular, for a given choice of $N$, it remains only to choose either $\mathrm{\Delta }x$ (the isotropic spatial sampling period), or equivalently, $X$ (half the total spatial width of the simulation; see Section 1.A). I choose the value of $\mathrm{\Delta }x$ to ensure that, within the primed coordinate system, a spherical wave that has its source in the plane $z=z_a$ is correctly sampled after propagating distance $z^{\prime }=(z_o-z_a)z_a/z_o$ [see Eq. (25)], resulting in the requirement given by Eq. (14):

The results of this simulation are shown in Figs. 5 and 6. Figure 5 shows the magnitude of the diffracted field at $z=z_o$ directly evaluated using Fresnel–Kirchhoff diffraction theory as per Eq. (1) of Ref. [13], without using the primed coordinate system. The field magnitudes as evaluated using the primed coordinate system angular spectrum model for $N=1024,3072,5120$ and 7168 appear very similar to that plotted in Fig. 5. As a result, the magnitude of the differences between the complex amplitudes evaluated using the angular spectrum approach and those evaluated directly are plotted in Fig. 6. These plots show that at lower values of $N$, the diffracted field at $z=z_o$ is spatially truncated as a result of increased filtering of the aperture’s projection function. These plots show clearly that increasing $N$ increases the spatial extent over which the calculated diffracted field remains accurate. To further probe the effect that $N$ has upon accuracy, I introduce an error metric, defined as

 

Fig. 5. Magnitude of the field diffracted by the aperture found by direct evaluation of the Fresnel–Kirchhoff diffraction integral.

Download Full Size | PDF

 

Fig. 6. Plots of the magnitude of the error in each calculated complex amplitude relative to that calculated using Fresnel–Kirchhoff diffraction theory for a diffracting aperture.

Download Full Size | PDF

 

Fig. 7. Plot of the error $ϵ$ as a function of $N$, where the reference field is that calculated using Fresnel–Kirchhoff diffraction theory.

Download Full Size | PDF

B. Diffracting Aperture and Single Sphere

In this example, I introduce a single sphere of radius 5 μm with $\delta =2\times {10}^{-6}$, located a distance 0.1 m downstream of the aperture and situated at transverse position (5,0) μm relative to the center of the aperture. Following the same method of presentation as in the previous example, the diffracted field found by direct evaluation of the Fresnel–Kirchhoff diffraction integral is plotted in Fig. 8. For the case of a single sphere, a two-dimensional extension of Eq. (8) in Ref. [13] was used to find the directly evaluated field. The difference between the primed coordinate system angular spectrum model and the directly evaluated field is also plotted in Fig. 9, which is seen to be very similar to the case where a sphere was not present. I do not plot the error metric in this case, as it was seen to be nearly identical to that which arose in the absence of a sphere, as is plotted in Fig. 7. This example serves the purpose of demonstrating that the primed coordinate system angular spectrum model converges to the directly evaluated field for a reasonably general example.

 

Fig. 8. Magnitude of the field diffracted by the aperture and a single sphere found by direct evaluation of the Fresnel–Kirchhoff diffraction integral.

Download Full Size | PDF

 

Fig. 9. Plots of the magnitude of the error in each calculated complex amplitude relative to that calculated using Fresnel–Kirchhoff diffraction theory for the case of a single sphere and diffracting aperture.

Download Full Size | PDF

C. Diffracting Aperture and Two Spheres

In this example, I include an additional sphere with the same radius and $\delta$ value as in the previous example, however, located 0.15 m downstream of the aperture and at transverse location (5,5) μm with respect to the center of the aperture. This example was calculated for the case $N=7168$ only. The field at the observation plane, calculated using the primed coordinate system angular spectrum model is shown in the left of Fig. 10. For validation purposes, a field was evaluated directly using Fresnel–Kirchhoff diffraction theory. The Fresnel–Kirchhoff field was evaluated by first considering each of the two spheres in isolation. The field resulting from this calculation, as if Born’s first-order approximation holds, is shown in the right hand of Fig. 10, which demonstrates that the first-order approximation is not appropriate in this case. The field can be evaluated correctly by taking into account the field that is scattered by both spheres, as is detailed in Appendix B. Once the multiply scattered field is taken into account, an error of $ϵ=1.2\times {10}^{-3}$ is obtained, thus illustrating the accuracy of the primed coordinate system angular spectrum model.

 

Fig. 10. Magnitude of the field obtained when two spheres are located between the aperture and observation plane (left) and the magnitude of the field obtained by considering the two spheres in isolation (right).

Download Full Size | PDF

D. Diffracting Aperture and an Ensemble of Spheres

This final example illustrates the type of application where this model may be particularly useful. Using the same combination of source and aperture as in the previous example, this example considers an ensemble of non-overlapping spheres. Spheres were arranged within a cuboid with a square transverse cross section of width 100 μm and bounded by the planes located 5 cm and 10 cm, respectively, downstream of the diffracting aperture. A total of 47,746 spheres were arranged within this volume resulting in a sphere density of 5% by volume. Although the computation time could have been reduced by considering slices containing multiple spheres, the simulation was performed with one sphere per slice. This was done for two main reasons: (1) to demonstrate how computationally efficient the simulation is, and (2) to ensure the validity of the projection approximation. Simulations were performed for a single ensemble of spheres, but each simulation considered a different value of $\delta$. In any one simulation, each sphere had the same $\delta$ value.

The results of the simulation are illustrated in Fig. 11. Each of the four columns in Fig. 11 corresponds to a different value of $\delta$ of the spheres, including the free-space case. The first row of images shows a somewhat zoomed-out view of the field magnitude on the observation plane, while the middle row shows a zoomed-in view. The lower row of images shows a histogram of the proportion of pixels having field magnitude within a certain range. The bar plots are derived from the pixels within the magnified views of the middle row of images. The line plot shows the distribution of field magnitudes that would be expected from a Rayleigh distribution [14] having mean equal to that of the simulated field distributions shown in the middle row of images. This comparison with the Rayleigh distribution is performed purely to demonstrate that for $\delta ={10}^{-6}$, the complex amplitude can be considered to be a fully developed speckle pattern.

 

Fig. 11. Magnitude of the field that results when an ensemble of spheres is located between a square diffracting aperture and the observation plane. Each group of three images corresponds to a different value of $\delta$ for the spheres, including the top left group, which represents free space. The middle image in each group is a magnified view of the region denoted by the outline box in the top image. The lowest image shows a histogram of the magnitudes of the pixels within the magnified region along with a Rayleigh distribution with mean equal to that of the magnified region.

Download Full Size | PDF

E. Comparison of Computational Complexity

The primed coordinate system angular spectrum method is several orders of magnitude more computationally efficient than direct evaluations of the Fresnel–Kirchhoff integral for the case of a single sphere. For multiple spheres, direct evaluations of the Fresnel–Kirchhoff integral become unfeasible from a computational point of view. This is illustrated in Fig. 12, which plots the computation time per observation point against the number of observation points. This timing information was obtained for the case of a single sphere and a square diffracting aperture, as discussed in Section 3.B. Simulations were run on a computer containing 512 MB of RAM and two Intel Xeon Gold 6148 Processors each possessing 20 physical cores running at 2.40 GHz. I note, however, that only a small proportion of the RAM was used for either simulation. Both simulation techniques were implemented in MATLAB (R2017b), and all 40 physical cores were used. The Fresnel–Kirchhoff method was parallelized trivially by using a parfor loop to compute the field at different points in the observation plane in parallel. The primed coordinate system angular spectrum method was parallelized using the parallelization intrinsic to MATLAB’s implementation of the fast Fourier transform (i.e., fft2).

 

Fig. 12. Plots of computation time per sample point in the observation plane for the Fresnel–Kirchhoff and primed coordinate system angular spectrum methods.

Download Full Size | PDF

Both simulation techniques are seen to exhibit a reduction in the computation per sample point as the number of sample points increases as simulation overheads are amortized over an increasing number of points. The computation time per point for the primed coordinate system angular spectrum method is on the order of $5\times {10}^4$ times lower than that of the Fresnel–Kirchhoff method.

The simulations considered in Section 3.D each required approximately the same computation time, irrespective of the value of $\delta$ the spheres were assumed to have. Each computation such as are displayed in Fig. 11 took an average of $63\times {10}^4$ seconds to compute. While this may be considered a substantial computation time, such a calculation is intractable using the Fresnel–Kirchhoff formalism. Also, as discussed in Section 3.D, this simulation could have been evaluated significantly faster by considering multiple spheres per slice, thus substantially reducing the number of evaluations of fft2.

4. CONCLUSION

I have shown how the simulation of divergent beams propagating through axially extended diffracting objects can pose sampling challenges when using Fourier-based MS propagation techniques. Two principle solutions to these challenges have been demonstrated. The first of these is to use a divergent-wave-to-plane-wave transformation that significantly reduces the sampling requirement when simulating the propagation of such beams through homogeneous volumes. The second of these is the calculation of band-limited object projection functions, which guarantees that aliasing will not occur when using Fourier theory to propagate a beam that has been perturbed by a diffracting object under the projection approximation. Both of these solutions have been combined into a MS simulation technique. I have provided validation of the technique for an example employing a square aperture only and for the case of a square aperture with one or two spheres, respectively. I have also used the technique to simulate beam propagation through an ensemble of spheres dispersed throughout a macroscopic scale volume. I expect this technique to be useful in the study of modalities such as x-ray dark field imaging.