Abstract
When calculating the focusing properties of cylindrically symmetric focusing optics, numerical wave propagation calculations can be carried out using the quasi-discrete Hankel transform (QDHT). We describe here an implementation of the QDHT where a partial transform matrix can be stored to speed up repeated wave propagations over specified distances, with reduced computational memory requirements. The accuracy of the approach is then verified by comparison with analytical results, over propagation distances with both small and large Fresnel numbers. We then demonstrate the utility of this approach for calculating the focusing properties of Fresnel zone plate optics that are commonly used for x-ray imaging applications and point to future applications of this approach.
© 2015 Optical Society of America
1. INTRODUCTION
X-ray nanofocusing optics can be used in x-ray imaging and spectroscopy techniques to provide new insights into the structure and functioning of cells and materials [1,2]. While there have been impressive advances in the development of x-ray mirrors [3,4], compound refractive lenses [5–7], and multilayer Laue lenses [8,9], Fresnel zone plates are used for the majority of applications requiring sub-100 nm beam spots [10,11] due to their high-quality focusing and simplicity of alignment. Therefore, efficient simulations of zone plate focusing are useful for developing new approaches and improvements in x-ray microscopy and spectromicroscopy.
Zone plate focusing represents one example of a wide variety of optics calculations involving the propagation of wavefields with wavelength through free space a distance from a plane to a plane . The general topic is well treated in textbooks (see, for example, [12]) and various papers describing numerical wave propagation in cylindrical coordinates using near-field expansions [13], Hankel transforms though without beam focusing examples [14], Helmholtz equation solutions for near-field propagation into waveguides [15], and mode propagation within optical fibers [16–18] as well as in free space. These approaches include mode expansions using Lanczos [15,16] and Arnoldi [17,18] methods with great utility for those applications. Our goal here is to describe numerical Hankel transform methods that can later be applied to optics which are less easy to characterize in terms of mode structures, such as Fresnel zone plates with various errors in zone placement [19,20] when cylindrical symmetry still applies, and departures from the standard zone plate formulation [21]. We also discuss criteria for choosing near-field versus far-field computational approaches.
In 2D Cartesian coordinates and within the Fresnel approximation [12], wave propagation can be carried out using forward and inverse 2D Fourier transform pairs of
to describe forward propagation either with or as a 2D convolution using In both cases a uniform phase shift term has been dropped. These expressions employ dimensionless propagator functions of with as spatial frequencies [12]. Note that for the convolution approach of Eq. (4), we have incorporated the prefactor into the propagator function of Eq. (6).In cases where cylindrical symmetry applies, the input plane and a plane located a distance away can be represented in cylindrical coordinates. In this case, the transforms change from the 2D Fourier transform pair of Eqs. (1) and (2) to a zeroth-order Hankel transform pair of
where is a Bessel function of the first kind. Using the Hankel transform, the expressions for wavefield propagation in cylindrical symmetry become or, in the equivalent convolutional approach, Here, the propagator functions in cylindrical coordinates become with representing a radial spatial frequency, similar to Eq. (7). As with Eqs. (4) and (6), the prefactor for the cylindrically symmetric convolution approach of Eq. (11) has been incorporated into Eq. (13).The two analytically equivalent methods of Eqs. (10) and (11) are useful in different regimes for optics with cylindrical symmetry. With the Hankel transform approach of Eq. (10), the input plane wavefield is multiplied by the real space propagator given in Eq. (12). In the convolution approach, the Hankel-transformed input plane wavefield is multiplied by the reciprocal space propagator of Eq. (13). It then becomes important to consider the nature of oscillations in the two propagator functions when deciding which method to use. This is shown in Fig. 1, which indicates that the reciprocal space propagator is slowly varying at shorter propagation distances (so that it minimizes aliasing artifacts when applying to discretely sampled functions), while the real space propagator is slowly varying at longer propagation distances.
To find the crossover point between the two approaches, consider the problem of propagating a monochromatic, coherent plane wave from a Fresnel zone plate to a plane a distance away. In real space, the argument of the real space propagator of Eq. (12) [appropriate for propagation using Eq. (10)] is , so the total number of Fresnel zones (number of phase shifts) within a radius is
where is the number of sampling points and is an averaged sampling pixel size of . Similarly, the argument of the Fourier space propagator of Eq. (13) [appropriate for propagation using Eq. (11)] is . A Nyquist sampling interval of in real space corresponds to a maximum spatial frequency of in Fourier space. Therefore, the maximum number of Fresnel zones in Fourier space is The matching distance at which we arrive at an identical number of Fresnel zones in real and Fourier space, or , is found from Eqs. (15) and (16) to be The Fresnel number for propagation from an aperture of radius over a distance is given by If we solve Eq. (17) for , we obtain where we see that the number of sampling points required at the matching distance is equal to twice the Fresnel number if the aperture spans over the whole space .The propagation distance does not set a hard boundary between the two propagation approaches of Eqs. (10) or (11); instead, both approaches are valid. However, it does suggest an approximate boundary for which approach will work with fewer Fresnel zones or and thus more sampling points per phase shift. For propagation over distances , we prefer to use the convolutional approach of Eq. (11) which can be written as
which has as . For , we prefer to use the Hankel transform approach of Eq. (10) written as which has tending toward the Hankel transform of as , which is the usual result for Fraunhofer diffraction. Consider the example of a Fresnel zone plate with radius and outermost zone width ; its focal length is found from in the paraxial approximation. One must have enough sampling points, or , to have good sampling of the wavefield within the outermost zone. When propagating a monochromatic, coherent plane wave to its focal plane, using Eq. (17) we have which indicates that the Hankel transform approach is preferred.2. QUASI-DISCRETE HANKEL TRANSFORM AND THE SAMPLING THEOREM
We now wish to consider the discrete form of the Hankel transforms of Eqs. (8) and (9). In this case the integration limits will be set to and for real and reciprocal space, respectively, and the wavefield will be sampled over discrete values. By using a Fourier–Bessel series to approximate the Hankel transform over a finite integral region, quasi-discrete Hankel transform (QDHT) methods have been developed by Yu [22] in zeroth order and Guizar–Sicairos [23] at higher orders. The QDHT uses discrete sampling points at
to yield with where and are the radial integration limits of and , respectively, and where is the th root of the zeroth-order Bessel function .A disadvantage of the QDHT is that one cannot calculate values at arbitrary radii or ; instead, one must use the sampling points of the QDHT shown in Fig. 2. Because no sampling points are close to , one cannot directly calculate wavefields at axial positions, where (for example) the intensity of a focused beam is strongest. To address this limitation, Norfolk used a generalized sampling theorem [24] which we restate as follows. The functions sampled on the grid , where [22], can be reconstructed at an arbitrary point as
with the sampling kernel We may therefore reconstruct the function of Eq. (26) as at any arbitrary point in , including with which is important for preserving overall wave intensity. The generalized sampling theorem of Eqs. (28) and (29) is also very useful for obtaining a smooth intensity distribution from coarser sampling.3. RAPID CALCULATION USING PARTIAL TRANSFORMS
To perform the Hankel transform (HT), the zeroth-order HT and inverse HT can be rewritten with and sampled over ranges up to in real space, and in reciprocal space, at the roots of the Bessel function:
By defining an transform matrix whose th element is we can simplify Eqs. (31) in terms of matrix multiplications on and of However, calculations of the full wavefield have large computational demands. Consider the case of simulating focusing from a Fresnel zone plate with 300 μm diameter and , where one might want in order to get a smooth intensity profile at the focus. Though the QDHT uses a nonuniform grid especially at low radial values, when is large the sampling point spacings approach a mean value of . If done with the same fine sample spacing and total field of view at both input and output planes, this would require determination of a matrix with which would occupy about 670 GBytes as single precision floating point complex numbers and require considerable computation time. As a result, previous work has involved simulations of small diameter zone plates with high resolution [25] or large diameter zone plates with low resolution [26].In order to overcome this limit, we have implemented the QDHT using partial transform matrices. Consider the case of calculating the focus spot profile of a Fresnel zone plate with a focal length , where the Hankel product approach of Eq. (21) is preferred. To propagate an entire wavefield within a radial distance from the zone plate to the focal plane, we would require a transform matrix which could be prohibitively large as noted above. However, in many cases what we are interested in is the detailed focal profile near the optical axis, with less need for a detailed calculation of the wavefield at larger radii. In this case we can use a matrix with , as shown in Fig. 3. For example, we might need only points to see the detailed focal profile (including several Airy fringes) of our example zone plate at calculation grid size. This saves a factor of 3000 in computation time and required storage over the example given above. The choice of depends on the radius range of the output plane that we want to see; it should be at least large enough to see the focus. We show such an example calculation in Fig. 4 which uses ; this involved a time of less than 0.1 s for a single propagation distance when using a server with dual Intel Xeon X5550 processors and 48 GB RAM. We note that for the convolution method, a full transform matrix is required as it involves both forward and inverse transforms.
4. COMPARISONS AGAINST ANALYTICAL CALCULATIONS
In order to check the accuracy of the QDHT propagation method described above, we have compared it against a situation with a well-known analytical result: the Airy pattern that results from far-field diffraction of light by a pinhole. This comparison was done with the single transform approach of Eq. (21), using a pinhole with a diameter of on a calculation grid of spacing extending to a distance of so that samples were involved. The far-field diffraction pattern for x rays at involved a Fresnel number of 0.1 and was calculated using calculation points at a spacing of . The results shown in Fig. 5 show that the maximum error in the far-field diffraction intensity was about 0.12%. At a farther distance of , the maximum error goes down to 0.08% as it is even farther away from having any non-Fraunhofer terms contributing.
We have shown in Fig. 4 an example calculation of the intensity distribution produced around the focal point of a Fresnel zone plate, which again should follow an Airy intensity profile [27]. In Fig. 6, we show both the intensity distribution, integral of intensity with radius, and percentage difference from the analytical result for a binary, fully absorptive Frensel zone plate with a diameter of 45 μm and outermost zone width of used to focus (10 keV) x rays. The position of the first minimum of the Airy pattern for such a zone plate is at a multiple of the first root of divided by times the outermost zone width, or , which is 30.5 nm in this example, while the integrated energy fraction should approach [28]. For the numerical QDHT calculation of Eq. (21), sampling points were used within a maximum radius , while in the output plane points were used at a spacing of . Again, the numerical results are in good agreement with the expected values.
The calculations of Figs. 5 and 6 show the fraction of incident energy present near the center of the Airy pattern (the diffraction pattern from a pinhole in Fig. 5 and the first-order Fresnel zone plate focus in Fig. 6). To check the accuracy of calculating the total energy leaving an input plane, one must integrate out to larger radii and compare it with the incident energy. Two such calculations are shown in Fig. 7. For the calculation of diffraction from a pinhole shown on the left, a wavelength was propagated a distance of 50 cm downstream, using a calculation grid with points at spacing on the input plane, and points up to a maximum radius of 500 μm at the output plane. As can be seen, this captures 99.7% of the energy, which is in excellent agreement with the analytical result. For the calculation for a binary absorption Fresnel zone plate, one can see that about 10% of the total energy is located on the optical axis in the form of the first focal order. About 25% of the energy is captured within the positive focal orders near the optical axis, and 50% of the beam energy is transmitted over all radii, as expected.
5. COMBINED PROPAGATION WITHIN AND BEYOND
While both the near-distance convolution approach of Eq. (20) and the far-distance single Hankel transform approach of Eq. (21) are valid at all distances, as described above they offer different sampling properties on either side of the distance of Eq. (17). Consider the case of first-order focusing from a Fresnel zone plate, where one can write , where is the focal length of the zone plate, is the sampling interval in real space, and is the width of the finest, outermost zone of the zone plate. Because one desires to have be much smaller than in order to have many samples within the width of the outermost zone, the focal length is always in the far-distance location (that is, is large compared to ). In Fig. 8, we show the radial and longitudinal focus intensity profile of a Fresnel zone plate calculated using both approaches. As can be seen, the far-distance method leads to a smooth intensity profile, whereas the near-distance method leads to irregularities in the intensity profile on the optical axis due to the fast oscillations in the reciprocal space propagator function, as shown in Fig. 1. Even so, the near-distance approach gives the correct result for the integrated intensity and indeed for intensities away from the optical axis.
Sometimes it is desirable to do calculations with multiple propagations over various distances. One might wish to propagate a wavefield to over a distance beyond from a lens to a cylindrically symmetric object near the focus, and then carry out a series of multislice propagations [29] over short distances through the object; or propagate a wavefield through several zone plates stacked a short distance from each other and then on to their common focus point located some farther distance away [25]. In these cases, one will want to use a mixture of both the near-distance propagation approach of Eq. (20) as well as the longer-distance propagation approach of Eq. (21), so the sampling grid and range must be considered. While in the continuous case we wrote the input and output radii as and , respectively, we will write their discrete counterparts as and .
For propagation over shorter distances , the convolution approach of Eq. (20) samples an input wavefield , transforms it to reciprocal space as where it is multiplied by a real space propagator [Eq. (13)], and inverse transforms it back to real space as on the same calculational grid. If we sample the input wavefield in points over a radial extent , the real space sampling interval is , and the maximum value of the Hankel function argument is [Eq. (27)]. When is large, can be chosen as , so
becomes the maximum spatial frequency, which when divided by the number of samples yields an interval in reciprocal space of Multiplication with the propagator of Eq. (13) modifies but does not change its sample positions. The inverse QDHT brings the wavefield back to real space, and the discretely sampled, propagated wavefield extends to a maximum radius with interval . Therefore, the new array is in real space with unchanged sampling and extent.For propagation over longer distances , the single Hankel transform method of Eq. (21) is preferred. We start with the discrete 1D array in real space and multiply it by the real space phase propagator of Eq. (12). We then perform the QDHT to bring the wavefield into reciprocal space , with the same extent as Eq. (34) and interval as Eq. (35). In this approach there is no second QDHT; instead, the reciprocal space array is multiplied by where the real space positions are found from . The propagated wavefield extends to a maximum radius of
with sampling interval Obviously, the output plane real space sampling interval might differ from the input plane interval . As an example, consider a zone plate with diameter and outermost zone width , where again the paraxial approximation gives . For propagation of a wavefield from the zone plate to the focus, or , Eq. (37) becomes If we want to keep the sampling interval constant, or , we need to choose the number of sampling points to be Equation (37) can be used to choose for maintaining at other propagation distances as well.6. CONCLUSION
We have described here an approach for the numerical propagation of cylindrically symmetric wavefields with increased speed and reduced array size requirements and demonstrated its accuracy by comparison with analytical results. For propagation over distances less than , as given by Eq. (17), the convolution approach of Eq. (20) which involves two QDHTs is preferred, while for longer distances the single QDHT approach of Eq. (21) is freer of aliasing artifacts. When simulating the focusing properties of Fresnel zone plates or other cylindrically symmetric optics, one can either choose a small number of output sampling points on a fine sampling interval to calculate the detailed profile of the beam focus near the optical axis, or one can use a coarser sampling yet still recover the efficiency of the optic as demonstrated in Fig. 7.
Calculations of this sort play an important role in the prediction of the performance of optics such as Fresnel zone plates, which are commonly used for high-resolution x-ray focusing [1,27]. In order to improve focusing efficiency within the limits of high aspect ratio nanofabrication approaches, several zone plates can be aligned onto successive axial positions, either in the near field [30] or at greater separation distances [25]. While this approach has recently been used to achieve 19% diffraction efficiency for focusing 25 keV x rays [31], there are several questions on the optimization of these approaches that we plan on addressing in future work. As zone plate thickness is increased further, one must begin to adjust the zones to meet the Bragg condition for volume diffraction [32,33], and QDHT multislice propagation calculations may provide a way of rapidly estimating focusing properties with subsequent validation using rigorous coupled-wave theory.
Funding
Basic Energy Sciences (BES) (DE-AC02-06CH11357).
Acknowledgment
We thank Michael Wojcik of Argonne National Laboratory for many helpful comments.
REFERENCES
1. A. Sakdinawat and D. T. Attwood, “Nanoscale x-ray imaging,” Nat. Photonics 4, 840–848 (2010). [CrossRef]
2. G. E. Ice, J. D. Budai, and J. W. Pang, “The race to x-ray microbeam and nanobeam science,” Science 334, 1234–1239 (2011). [CrossRef]
3. P. Kirkpatrick and A. V. Baez, “Formation of optical images by x-rays,” J. Opt. Soc. Am. 38, 766–774 (1948). [CrossRef]
4. H. Mimura, S. Handa, T. Kimura, H. Yumoto, D. Yamakawa, H. Yokoyama, S. Matsuyama, K. Inagaki, K. Yamamura, Y. Sano, K. Tamasaku, Y. Nishino, M. Yabashi, T. Ishikawa, and K. Yamauchi, “Breaking the 10 nm barrier in hard-x-ray focusing,” Nat. Phys. 6, 122–125 (2009). [CrossRef]
5. B. X. Yang, “Fresnel and refractive lenses for x-rays,” Nucl. Instrum. Methods Phys. Res. A 328, 578–587 (1993).
6. A. Snigirev, V. Kohn, I. Snigireva, and B. Lengeler, “A compound refractive lens for focusing high energy x-rays,” Nature 384, 49–51 (1996). [CrossRef]
7. C. Schroer, M. Kuhlmann, U. Hunger, T. Gunzler, O. Kurapova, S. Feste, F. Frehse, B. Lengeler, M. Drakopoulos, A. Somogyi, A. Simionovici, A. Snigirev, I. Snigireva, C. Schug, and W. Schroder, “Nanofocusing parabolic refractive x-ray lenses,” Appl. Phys. Lett. 82, 1485–1487 (2003). [CrossRef]
8. J. Maser, G. Stephenson, S. Vogt, W. Yũn, A. Macrander, H. Kang, C. Liu, and R. Conley, “Multilayer Laue lenses as high-resolution x-ray optics,” in Design and Microfabrication of Novel X-ray Optics II, A. Snigirev and D. Mancini, eds. (SPIE, 2004), Vol. 5539, pp. 185–194.
9. H. Yan, R. Conley, N. Bouet, and Y. S. Chu, “Hard x-ray nanofocusing by multilayer Laue lenses,” J. Phys. D 47, 263001 (2014).
10. A. G. Michette, “X-ray microscopy,” Rep. Prog. Phys. 51, 1525–1606 (1988). [CrossRef]
11. M. Howells, C. Jacobsen, and T. Warwick, “Principles and applications of zone plate x-ray microscopes,” in Science of Microscopy, P. W. Hawkes and J. C. H. Spence, eds., 1st ed. (Springer, 2006), Chap. 13.
12. J. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company, 2005).
13. M. D. Feit and J. A. Fleck, “Simple spectral method for solving propagation problems in cylindrical geometry with fast Fourier transforms,” Opt. Lett. 14, 662–664 (1989). [CrossRef]
14. D. Lemoine, “Highly accurate discrete Bessel representation of beam propagation in optical fibers,” J. Opt. Soc. Am. A 14, 411–416 (1997). [CrossRef]
15. R. P. Ratowsky, J. A. Fleck, and M. D. Feit, “Helmholtz beam propagation in rib waveguides and couplers by iterative Lanczos reduction,” J. Opt. Soc. Am. A 9, 265–273 (1992). [CrossRef]
16. Q. Luo and C. T. Law, “Nonparaxial propagation of a cylindrical beam with Lanczos reduction,” Opt. Lett. 25, 869–871 (2000). [CrossRef]
17. Q. Luo and C. T. Law, “Propagation of nonparaxial beams with a modified Arnoldi method,” Opt. Lett. 26, 1708–1710 (2001). [CrossRef]
18. Q. Luo and C. T. Law, “Discrete Bessel-based Arnoldi method for nonparaxial wave propagation,” IEEE Photon. Technol. Lett. 14, 50–52 (2002). [CrossRef]
19. M. J. Simpson and A. G. Michette, “The effects of manufacturing inaccuracies on the imaging properties of Fresnel zone plates,” Opt. Acta 30, 1455–1462 (2010). [CrossRef]
20. C. Pratsch, S. Rehbein, S. Werner, and G. Schneider, “Influence of random zone positioning errors on the resolving power of Fresnel zone plates,” Opt. Express 22, 30482–30491 (2014). [CrossRef]
21. K. Jefimovs, J. Vila-Comamala, T. Pilvi, J. Rabbe, M. Ritala, and C. David, “Zone-doubling technique to produce ultrahigh-resolution x-ray optics,” Phys. Rev. Lett. 99, 264801 (2007). [CrossRef]
22. L. Yu, M. C. Huang, M. Z. Chen, W. Z. Chen, W. D. Huang, and Z. Z. Zhu, “Quasi-discrete Hankel transform,” Opt. Lett. 23, 409–411 (1998). [CrossRef]
23. M. Guizar-Sicairos and J. C. Gutiérrez-Vega, “Computation of quasi-discrete Hankel transforms of integer order for propagating optical wave fields,” J. Opt. Soc. Am. A 21, 53–58 (2004). [CrossRef]
24. A. W. Norfolk and E. J. Grace, “Reconstruction of optical fields with the quasi-discrete Hankel transform,” Opt. Express 18, 10551–10556 (2010). [CrossRef]
25. J. Vila-Comamala, M. Wojcik, A. Diaz, M. Guizar-Sicairos, C. M. Kewish, S. Wang, and C. David, “Angular spectrum simulation of x-ray focusing by Fresnel zone plates,” J. Synchrotron Radiat. 20, 397–404 (2013). [CrossRef]
26. P. Srisungsitthisunti, O. K. Ersoy, and X. Xu, “Beam propagation modeling of modified volume Fresnel zone plates fabricated by femtosecond laser direct writing,” J. Opt. Soc. Am. A 26, 188–194 (2009). [CrossRef]
27. A. G. Michette, Optical Systems for Soft X Rays (Plenum, 1986).
28. J. Kirz, “Phase zone plates for x rays and the extreme UV,” J. Opt. Soc. Am. 64, 301–309 (1974). [CrossRef]
29. 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).
30. S. D. Shastri, J. M. Maser, B. Lai, and J. Tys, “Microfocusing of 50 keV undulator radiation with two stacked zoneplates,” Opt. Commun. 197, 9–14 (2001). [CrossRef]
31. S.-C. Gleber, M. Wojcik, J. Liu, C. Roehrig, M. Cummings, J. Vila-Comamala, K. Li, B. Lai, D. Shu, and S. Vogt, “Fresnel zone plate stacking in the intermediate field for high efficiency focusing in the hard x-ray regime,” Opt. Express 22, 28142–28153 (2014). [CrossRef]
32. J. Maser and G. Schmahl, “Coupled wave description of the diffraction by zone plates with high aspect ratios,” Opt. Commun. 89, 355–362 (1992). [CrossRef]
33. S. Werner, S. Rehbein, P. Guttmann, and G. Schneider, “Three-dimensional structured on-chip stacked zone plates for nanoscale x-ray imaging with high efficiency,” Nano Res. 7, 1–8 (2014).