Abstract
The simulation of the propagation of divergent beams using Fourier-based angular spectrum techniques can pose challenges for ensuring correct sampling in the spatial and reciprocal domains. This challenge can be compounded by the presence of diffracting objects, as is often the case. Here, I give details of a method for robustly simulating the propagation of beams with divergent wavefronts in a coordinate system where the wavefronts become planar. I also show how diffracting objects can be simulated, while guaranteeing that correct sampling is maintained. These two advances allow for numerically efficient and accurate simulations of divergent beams propagating through diffracting structures using the multi-slice approximation. The sampling requirements and numerical implementation are discussed in detail, and I have made the computer code freely available.
© 2019 Optical Society of America
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 , 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 and that our ultimate goal is to calculate the field in the plane , where is the object-to-detector distance. I thus refer to the space defined by as the computational region. I describe fields propagating in this system according to the form , where is the wavenumber, is the angular frequency of the radiation emitted by the source, and it is assumed that varies only weakly with . If is substituted into the time-harmonic scalar wave equation, it is not difficult to show that the paraxial wave equation is obtained as
The angular spectrum, , associated with a wavefield, , can be used to propagate such a field over some axial distance in a homogeneous space. It has been shown that and are related according to [7] It is then relatively straightforward to show that a component of an angular spectrum may be propagated a distance in homogeneous space according to [7] , where is the well-known free-space propagation operator.The spherical wave emitted by the point source has complex amplitude in the free-space region of , where . If I make the Fresnel approximation to this field and drop the dependence, I obtain the expression
which satisfies the paraxial wave equation, Eq. (1).Our objective in this paper is to employ DFT techniques to determine the complex amplitude that emerges from the plane when refractive index inhomogeneities may exist in the region . 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 as
the local spatial frequencies may be defined as [7] However, a further constraint on sampling arises after calculating the angular spectrum of using Eq. (2), which gives which must also be sampled according to the Nyquist criterion. By following the same procedure as was used to derive Eqs. (7) and (8), I can obtain the following expressions for the local frequency content of the angular spectrum: In any DFT-based field propagation technique, discrete spatial coordinates and propagation vectors must be defined. I denote these by and , respectively. Furthermore, for simplicity, in the remainder of this paper, I will assume that and . Then, if I have (assumed even without loss of generality) sample points for each such quantity, these discrete coordinates may be defined as For this analysis only, I shall neglect the additional complication arising from the implicit periodicity in underlying the DFT. This discretization must satisfy the Nyquist criterion in both the spatial and reciprocal spaces. In particular, and must satisfy where is taken to mean the maximum value of . Equations (12) and (13) reveal that and , which, with the aid of Eqs. (7) and (10), allows us to write In order to proceed, it is assumed that is fixed and that may be varied to achieve correct sampling. Then, by substituting Eqs. (16) and (17) into Eqs. (14) and (15), respectively, and using the relationship , the following inequalities are obtained, which must be satisfied if the Nyquist criterion is to be satisfied in both the spatial and reciprocal spaces: which are unable to be satisfied simultaneously. Equation (18) states that since is fixed, the transverse width of the simulation must not exceed a particular value if the field incident upon the computational volume is to be sampled correctly in the spatial domain. Equation (19), however, states that in order to correctly sample the angular spectrum in the reciprocal space for fixed , at some plane beyond the computational volume entrance plane, the transverse width of the simulation must exceed a particular value. The problem arises because these two constraints cannot be satisfied simultaneously.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 , described by
Under the projection approximation, it is assumed that if the field exiting the region , in the absence of a refractive inhomogeneity, is given by , the perturbed field is given as , where and is the free-space wavenumber. In obtaining this result, it is assumed that outside of the support of the refractive index inhomogeneity, and . Typical implementations of the MS method work by partitioning the volume containing a diffracting structure into slices, as depicted in Fig. 2, such that the th slice is defined by . The field exiting the plane is first calculated assuming the slice is composed of free space. This is performed by calculating the angular spectrum of the field according to Eq. (2), propagating the angular spectrum to the end of the slice using the free-space propagation operator [Eq. (4)], before evaluating the field at , , using Eq. (3). Refractive index inhomogeneity in the slice is then accounted for by multiplying by , where is evaluated according to Eq. (21). This procedure is then repeated until the field at is obtained.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 is known, which I denote by . Consider for now that the space 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 and a coordinate system , valid only for . In order to evaluate , I first perform a transformation into a primed coordinate system defined by [8,9]
where . It is clear from Eqs. (22)–(25) that at , I have , , and ; thus, the transverse coordinates are identical. Furthermore, it is straightforward to show that satisfies the paraxial wave equation [Eq. (1)] in the primed coordinate system. This means that techniques such as angular spectrum propagation, developed for calculating the propagation of fields that satisfy the paraxial wave equation, can be applied to instead of . This overcomes the sampling problem expressed by Eqs. (18) and (19), since a divergent spherical wave in the global coordinate system becomes a plane wave in the primed coordinate system.Having obtained from Eq. (22), it remains to calculate using angular spectrum propagation. I begin by calculating the angular spectrum of using Eq. (2), thus giving , which can be propagated a distance according to Eq. (4) as
where, from Eq. (25), and . The complex amplitude in the primed coordinate system can thus be obtained by applying Eq. (3) to , thus giving . The complex amplitude in the global coordinate system is found by inverting the coordinate system transformation expressed in Eqs. (22)–(25), which will result in which illustrates how returning from the primed (i.e., local) coordinate system to the global coordinate system entails a magnification of the transverse coordinates. Having obtained , the procedure outlined above may be repeated to propagate the complex amplitude from the plane to the plane and so on.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 . 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 , where I have dropped the dependence of and on and , 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 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 will not, in general, be band limited. Thus, if a complex amplitude in the absence of a diffracting object, , is band limited, the complex amplitude after the application of the projection approximation 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 , I define the band-limited version of as
where is the band-limited version of , is a windowing function, , is the DFT operator, and is the continuous Fourier transform operator. I note here that it is convenient to introduce the spatial frequency parameters . Variables under the tilde sign are assumed to be discretized. A variety of windowing functions could be used; however, in this work, a Tukey window was employed.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 is divided into slices. Propagation from the plane to the plane can be achieved by performing 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 is denoted . Propagation through the first slice (i.e., the slice defined by ) is achieved by evaluating according to Eqs. (22)–(25), giving the field in the primed coordinate system at . Under this transformation, the primed and global coordinate systems coincide at , i.e., and . Propagation to would be achieved by multiplying the angular spectrum associated with (i.e., ) according to
where . Returning to the global coordinate system at entails the change of coordinates An important detail emerges when propagating the field through the second slice beginning at . In particular, following the same strategy as for propagation through the first slice, I apply the transformation into the primed coordinate system , where and , meaning that After propagation through this slice, the transformation becomes I can thus write the following general expressions for the coordinate system transformation at the entrance plane of a slice: and for the exit plane .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 , which is immediately transformed into the primed coordinate system according to Eqs. (22)–(25), yielding . 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
It is also necessary to make the following assignment when transferring from the end of one slice into the beginning of an adjacent slice: At the end of the final slice located at , the field in the primed coordinate system, , must be transformed back to the global coordinate system according to Eqs. (22)–(25). It is worth noting that at each iteration of Eq. (43), is simply a matrix that is continually updated by the DFT and multiplication operations contained within Eq. (43). In particular, no additional operations such as resampling are necessary. Instead, the underlying real-space coordinate system vectors are continually being rescaled according to where and are established at the beginning of the calculation, and it is understood that and correspond to the beginning of the slice . Furthermore, as is outlined in Section 2.B below, the reciprocal space coordinate system vectors are also constantly being rescaled according to where and are established at the beginning of the simulation.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 . 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 I have , the sampled grid in the primed coordinate system is defined as
and the reciprocal space propagation vectors as The only non-trivial aspect of this algorithm from a numerical point of view is that each time the field is propagated through a slice, the sampled grid [Eq. (49)] is magnified after transformation back to the global coordinate system. For example, after performing angular spectrum propagation and transforming back to the global coordinate system after the first slice, the sampled complex amplitude will now be implicitly defined on a global coordinate system sampled grid defined by and , respectively. In general, using Eqs. (41) and (42), the sampled grid in the global coordinate system at exit plane will be given by and , respectively. This has two principal consequences for the algorithm, the first of these being that at the exit plane , the sampled reciprocal space propagation vectors are given by and , respectively. The second of these implications is that any band-limited projection functions [see Eq. (31)] must be evaluated on a different spatially sampled grid for every slice. In some cases, this may not entail significant computational cost. However, in the case of diffracting spheres, e.g., where the Fourier transform of the projection function must be evaluated via numerical integration (see Appendix A), direct evaluation of this Fourier transform on a different grid for each slice may be too computationally costly. In this case, however, use can be made of the Fourier transform relationship meaning that must be resampled onto the reciprocal space frequency grid , which may be done efficiently by one-dimensional interpolation of a lookup table. As an additional step, prior to performing inverse Fourier transformation to obtain , it is convenient to multiply by to allow for modeling spheres located at an arbitrary position , where it is assumed that lies within the slice being considered.C. Sampling Requirements
I shall define the sampling requirements with reference to a cutoff frequency , such that the windowing function , assumed to be separable in and , takes the value of 0 for and . From this point on, I shall assume that the spatial, and therefore reciprocal, computational grids are entirely isotropic. The first requirement that must satisfy is
The next requirement is that propagation of the angular spectrum from to using Eq. (4) must be correctly sampled in the reciprocal space. Following Matsushima and Shimobaba [11], I find that the sampling period in the reciprocal space, must satisfy but, since , this can also be written as The criteria expressed in Eqs. (52) and (54) do not consider the sampling requirements imposed by diffracting objects. In the case of modeling a single diffracting object, I shall assume that is chosen such that the band-limited projection function, , has satisfactory agreement with the actual projection function, . I do not propose criteria for making this assessment; however, one can assess this by considering the difference between scattered fields predicted by both treatments. In order to avoid aliasing when using the projection approximation, a guard band of width should be employed as illustrated in Fig. 3. Whenever a complex amplitude, , encounters a diffracting object, , imposed via the projection approximation, aliasing will not occur if both and do not contain spectral components within the guard band. The field after diffraction, , will in general contain spectral components within the guard band. In order to prevent aliasing at a subsequent diffracting object, the spectral components of the field within the guard band can be filtered out. The energy lost due to such filtering can be monitored to ensure it is kept below a tolerable level, dependent upon the application. In practice however, it has been found in the examples considered in this paper that the amount of energy entering the guard band is negligible, meaning that artifacts arising from aliasing are also negligible.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 to have a background refractive index . If the slice contains a refractive index homogeneity as in Eq. (20), this should be redefined as
recalling that this applies to inside the inhomogeneity only. The free-space propagation operator [Eq. (4)] must be redefined as An additional final step must then be applied where the total field exiting the slice must be multiplied by which is necessary because the field is considered within the paraxial approximation, meaning that a term is factored out of the field. It is useful to note that the sampled reciprocal space parameters do not require rescaling due to the redefinition of the free-space propagation operator in Eq. (56).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 . 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.
The objective of this test is to demonstrate how the choice of , the number of sample points along each Cartesian direction, impacts upon the accuracy of the simulation method. In particular, for a given choice of , it remains only to choose either (the isotropic spatial sampling period), or equivalently, (half the total spatial width of the simulation; see Section 1.A). I choose the value of to ensure that, within the primed coordinate system, a spherical wave that has its source in the plane is correctly sampled after propagating distance [see Eq. (25)], resulting in the requirement given by Eq. (14):
where was chosen to have a value of 0.95 of its maximum allowable value. Note that Eq. (58) represents the worst-case scenario in which a diffracting object located at leads to the creation of a spherical wave within the primed coordinate system. The cutoff frequency is thus dictated by Eq. (52). Increasing the value of allows for to be increased, which reduces the difference between the band-limited and true projection functions of the aperture, and later spheres, as explained in Section 1.D.The results of this simulation are shown in Figs. 5 and 6. Figure 5 shows the magnitude of the diffracted field at 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 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 , the diffracted field at is spatially truncated as a result of increased filtering of the aperture’s projection function. These plots show clearly that increasing increases the spatial extent over which the calculated diffracted field remains accurate. To further probe the effect that has upon accuracy, I introduce an error metric, defined as
to quantify the error between complex amplitude and a reference complex amplitude . This error metric is plotted in Fig. 7, where the field calculated using Fresnel–Kirchhoff diffraction theory is used as the reference field. This plot shows that the error in the field calculated using the primed coordinate system angular spectrum model reduces approximately linearly with the log of .B. Diffracting Aperture and Single Sphere
In this example, I introduce a single sphere of radius 5 μm with , 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.
C. Diffracting Aperture and Two Spheres
In this example, I include an additional sphere with the same radius and 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 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 is obtained, thus illustrating the accuracy of the primed coordinate system angular spectrum model.
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 . In any one simulation, each sphere had the same 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 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 , the complex amplitude can be considered to be a fully developed speckle pattern.
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).
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 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 the spheres were assumed to have. Each computation such as are displayed in Fig. 11 took an average of 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.
APPENDIX A: BAND-LIMITED PROJECTION FUNCTION OF A SPHERE
The projection function of a sphere, located at the origin, with refractive index and radius is given as
and its Fourier transform by where is Dirac’s delta function. I can simplify Eq. (A1) by substituting in polar coordinates for both the spatial and frequency coordinates as allowing us to write Eq. (A1) as where is the zero-order Bessel function of the first kind. As discussed in Section 1.D, the band-limited version of is then obtained as recalling that and are the discretized spatial and reciprocal space sample grids, respectively.APPENDIX B: FRESNEL–KIRCHHOFF DIFFRACTION THEORY
I consider a complex amplitude incident upon a diffracting object, which I define as . I assume that the diffracting object is contained within the planes and . Then, I apply the projection approximation, as outlined in Section 1.A, to obtain the field in the plane as , where is evaluated according to Eq. (21) with and . The field at some arbitrary observation point, downstream of the diffracting object, may then be found according to
where is the Fresnel–Kirchhoff integral kernel given by [13]The integral in Eq. (B1) must be evaluated over an infinite plane, which is unsuitable for numerical integration. Instead, an often-used technique is to rewrite the integral as where is the field incident directly on in the absence of the diffracting aperture, and is the transverse extent of the diffracting object. Using this approach, the concept of the scattered field, , can be defined by decomposing as If now a second diffracting object is introduced within the planes and , with , the field at may be found by evaluating Eq. (B3), by noting that the field incident upon the second sphere is given as , which allows the field on the observation plane to be evaluated as Up until this point, I have used the notation to refer to the field scattered by the first diffracting object only. In order to avoid confusion, I now introduce and , which refer to the fields scattered by the two diffracting objects, respectively, each in isolation. This allows the field due to two diffracting objects to be decomposed as where is the field directly incident in the absence of any diffracting objects, and is the multiply scattered field.Funding
Royal Society (UF130304); Engineering and Physical Sciences Research Council (EPSRC) (EP/P005209/1).
Acknowledgment
P. M. is supported by a Royal Society University Research Fellowship.
REFERENCES
1. J. M. Cowley and A. F. Moodie, “The scattering of electrons by atoms and crystals. I. A new theoretical approach,” Acta Crystallogr. 10, 609–619 (1957). [CrossRef]
2. A. Hare and G. Morrison, “Near-field soft x-ray diffraction modelled by the multislice method,” J. Mod. Opt. 41, 31–48 (1994). [CrossRef]
3. Y. Wang, “A numerical study of resolution and contrast in soft x-ray contact microscopy,” J. Microsc. 191, 159–169 (1998). [CrossRef]
4. A. Malecki, G. Potdevin, and F. Pfeiffer, “Quantitative wave-optical numerical analysis of the dark-field signal in grating-based x-ray interferometry,” Europhys. Lett. 99, 48001 (2012). [CrossRef]
5. K. Li, M. Wojcik, and C. Jacobsen, “Multislice does it all—calculating the performance of nanofocusing x-ray optics,” Opt. Express 25, 1831–1846 (2017). [CrossRef]
6. D. Paganin, Coherent X-Ray Optics (Oxford University, 2006).
7. J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 2005).
8. E. A. Sziklas and A. E. Siegman, “Mode calculations in unstable resonators with flowing saturable gain. 2: fast Fourier transform method,” Appl. Opt. 14, 1874–1889 (1975). [CrossRef]
9. E. A. Sziklas and A. E. Siegman, “Diffraction calculations using fast Fourier transform methods,” Proc. IEEE 62, 410–412 (1974). [CrossRef]
10. K. S. Morgan, K. K. W. Siu, and D. M. Paganin, “The projection approximation and edge contrast for x-ray propagation-based phase contrast imaging of a cylindrical edge,” Opt. Express 18, 9865–9878 (2010). [CrossRef]
11. K. Matsushima and T. Shimobaba, “Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields,” Opt. Express 17, 19662–19673 (2009). [CrossRef]
12. P. R. T. Munro, “Multi-slice x-ray beam propagation code,” http://prtmunro.net.
13. P. R. T. Munro, K. Ignatyev, R. Speller, and A. Olivo, “The relationship between wave and geometrical optics models of coded aperture type x-ray phase contrast imaging systems,” Opt. Express 18, 4103–4117 (2010). [CrossRef]
14. J. W. Goodman, Laser Speckle and Related Phenomena (Springer-Verlag, 1975).