Rapid calculation of paraxial wave propagation for cylindrically symmetric optics

Kenan Li; Chris Jacobsen

doi:10.1364/JOSAA.32.002074

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 $z$ from a plane $(x_{0}, y_{0})$ to a plane $(x, y)$ . 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

F {g (x, y)} = \int_{- \infty}^{\infty} g (x, y) e^{i π (x f_{x} + y f_{y})} d x d y,

F^{- 1} {G (f_{x}, f_{y})} = \int_{- \infty}^{\infty} G (f_{x}, f_{y}) e^{- i π (x f_{x} + y f_{y})} d f_{x} d f_{y}

to describe forward propagation either with

ψ_{z} (x, y) = \frac{i}{λ z} h (x, y) F {ψ_{0} (x_{0}, y_{0}) h (x_{0}, y_{0})}

or as a 2D convolution using

ψ_{z} (x, y) = \frac{i}{λ z} [ψ (x_{0}, y_{0}) * h (x, y)] = F^{- 1} {F {ψ_{0} (x_{0}, y_{0})} H (f_{x}, f_{y})} .

In both cases a uniform phase shift term

e^{- i 2 π z / λ}

has been dropped. These expressions employ dimensionless propagator functions of

h (x, y) = e^{- i \frac{2 π}{λ} \sqrt{z^{2} + x^{2} + y^{2}}} \approx e^{- i π (x^{2} + y^{2}) / (λ z)},

H (f_{x}, f_{y}) = \frac{i}{λ z} F {h (x, y)} = e^{- i \frac{2 π}{λ} \sqrt{z^{2} - λ^{2} z^{2} (f_{x}^{2} + f_{y}^{2})}} \approx e^{i π λ z (f_{x}^{2} + f_{y}^{2})},

with

f_{x} = \frac{x}{λ z} and f_{y} = \frac{y}{λ z}

as spatial frequencies [12]. Note that for the convolution approach of Eq. (4), we have incorporated the

i / (λ z)

prefactor into the propagator function of Eq. (6).

In cases where cylindrical symmetry applies, the input plane $r_{0}$ and a plane $r$ located a distance $z$ 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

H {g (r)} = 2 π \int_{0}^{\infty} G (ρ) J_{0} (2 π ρ r) ρ d ρ,

H^{- 1} {G (ρ)} = 2 π \int_{0}^{\infty} g (r) J_{0} (2 π ρ r) r d r,

where

J_{0}

is a Bessel function of the first kind. Using the Hankel transform, the expressions for wavefield propagation in cylindrical symmetry become

ψ_{z} (r) = \frac{i}{λ z} h (r) H {ψ_{0} (r_{0}) h (r_{0})}

or, in the equivalent convolutional approach,

ψ_{z} (r) = H^{- 1} {H {ψ_{0} (r_{0})} H (ρ)} .

Here, the propagator functions in cylindrical coordinates become

h (r) = e^{- i \frac{2 π}{λ} \sqrt{z^{2} + r^{2}}} \approx e^{- i π r^{2} / (λ z)},

H (ρ) = \frac{i}{λ z} H {h (r)} = e^{- i \frac{2 π}{λ} \sqrt{z^{2} - λ^{2} z^{2} ρ^{2}}} \approx e^{i π λ z ρ^{2}},

with

ρ = \frac{r}{λ z}

representing a radial spatial frequency, similar to Eq. (7). As with Eqs. (4) and (6), the

i / (λ z)

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 $ψ_{0} (r_{0})$ 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.

Fig. 1. Real space [Eq. (12)] and reciprocal space [Eq. (13)] propagators shown for two different propagation distances $z$ with $λ = 1 nm$ . In each case, the real part is shown in blue and the imaginary part in green. The dots are the positions of $N = 1000$ sampling points over a radius of $R = 50 μm$ , for which Eq. (17) gives $z_{0} = 5 mm$ . The reciprocal space propagator is more slowly varying at short propagation distances, while the real space propagator is more slowly varying at longer distances.

Download Full Size | PDF

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 $z$ away. In real space, the argument of the real space propagator of Eq. (12) [appropriate for propagation using Eq. (10)] is $π r^{2} / (λ z)$ , so the total number of Fresnel zones $N_{real}$ (number of $π$ phase shifts) within a radius $R = r_{\max}$ is

N_{real} = \frac{R^{2}}{λ z} = \frac{N^{2} {〈 Δ r 〉}^{2}}{λ z},

where

N

is the number of sampling points and

〈 Δ r 〉

is an averaged sampling pixel size of

R / N

. Similarly, the argument of the Fourier space propagator of Eq. (13) [appropriate for propagation using Eq. (11)] is

π λ z ρ^{2}

. A Nyquist sampling interval of

〈 Δ r 〉

in real space corresponds to a maximum spatial frequency of

P = ρ_{\max} = 1 / (2 〈 Δ r 〉)

in Fourier space. Therefore, the maximum number of Fresnel zones in Fourier space

N_{Hankel}

is

N_{Hankel} = λ z P^{2} = \frac{λ z}{4 {〈 Δ r 〉}^{2}} .

The matching distance

z_{0}

at which we arrive at an identical number of Fresnel zones in real and Fourier space, or

N_{real} = N_{Hankel}

, is found from Eqs. (15) and (16) to be

z_{0} = \frac{2 R 〈 Δ r 〉}{λ} = \frac{2 R^{2}}{λ N} .

The Fresnel number

F

for propagation from an aperture of radius

a

over a distance

L

is given by

F = \frac{a^{2}}{λ L} .

If we solve Eq. (17) for

N

, we obtain

N = 2 \frac{R^{2}}{λ z_{0}} = 2 F_{0},

where we see that the number of sampling points

N

required at the matching distance

z_{0}

is equal to twice the Fresnel number

F_{0}

if the aperture

a

spans over the whole space

R

.

The propagation distance $z_{0}$ 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 $N_{real}$ or $N_{Hankel}$ and thus more sampling points per $π$ phase shift. For propagation over distances $z ≲ z_{0}$ , we prefer to use the convolutional approach of Eq. (11) which can be written as

ψ_{0} (r_{0}) \to H {} \to \times e^{i π λ z ρ^{2}} \to H^{- 1} {} \to ψ_{z} (r) (z ≲ z_{0}),

which has

ψ_{z} (r) \to ψ_{0} (r_{0})

as

z \to 0

. For

z ≳ z_{0}

, we prefer to use the Hankel transform approach of Eq. (10) written as

ψ_{0} (r_{0}) \to \times e^{- i π r_{0}^{2} / (λ z)} \to H {} \to \times \frac{i}{λ z} e^{- i π r^{2} / (λ z)} \to ψ_{z} (r) (z ≳ z_{0}),

which has

ψ_{z} (r)

tending toward the Hankel transform of

ψ_{0} (r_{0})

as

z \to \infty

, which is the usual result for Fraunhofer diffraction. Consider the example of a Fresnel zone plate with radius

R

and outermost zone width

d r_{zp}

; its focal length is found from

2 R d r_{zp} = λ f

in the paraxial approximation. One must have enough sampling points, or

N ≫ R / d r_{zp}

, 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

f = \frac{2 R d r_{zp}}{λ} > \frac{2 R^{2}}{λ N} = z_{0},

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 $R$ and $P$ for real and reciprocal space, respectively, and the wavefield will be sampled over $N$ 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

r \to α_{n} / (2 π P) in real space,

ρ \to α_{m} / (2 π R) in reciprocal space

to yield

G (α_{m} P / S) = \frac{1}{π P^{2}} \sum_{n = 1}^{N} \frac{g (α_{n} R / S)}{J_{1}^{2} (α_{n})} J_{0} (α_{n} α_{m} / S),

g (α_{n} R / S) = \frac{1}{π R^{2}} \sum_{m = 1}^{N} \frac{G (α_{m} P / S)}{J_{1}^{2} (α_{m})} J_{0} (α_{m} α_{n} / S),

with

S = 2 π R P,

where

R

and

P

are the radial integration limits of

g (r)

and

G (ρ)

, respectively, and where

α_{m}

is the

m

th root of the zeroth-order Bessel function

J_{0} (α)

.

A disadvantage of the QDHT is that one cannot calculate values at arbitrary radii $r$ or $ρ$ ; instead, one must use the sampling points of the QDHT shown in Fig. 2. Because no sampling points are close to $α = 0$ , 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 $r_{n}$ , where $r_{n} = α_{n} / 2 π P = (α_{n} / α_{N + 1}) R$ [22], can be reconstructed at an arbitrary point $r$ as

g (r) = \sum_{n = 1}^{\infty} g (r_{n}) K_{n} (r),

with the sampling kernel

K_{n} (r) = \frac{r_{n} J_{0} (2 π P r)}{π P (r_{n}^{2} - r^{2}) J_{1} (2 π P r_{n})} .

We may therefore reconstruct the function

g (r_{n})

of Eq. (26) as

g (r)

at any arbitrary point in

r

, including

r = 0

with

g (r = 0) = \sum_{n = 1}^{N} \frac{g (r_{n})}{π P r_{n} J_{1} (2 π P r_{n})},

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.

Fig. 2. Sampling grid of the QDHT (red cross marks) and FFT (blue square marks). The QDHT sampling points are nonequally spaced and do not include $α = 0$ .

Download Full Size | PDF

3. RAPID CALCULATION USING PARTIAL TRANSFORMS

To perform the Hankel transform (HT), the zeroth-order HT and inverse HT can be rewritten with $r_{n}$ and $ρ_{m}$ sampled over ranges up to $R$ in real space, and $P$ in reciprocal space, at the roots of the Bessel function:

{\begin{matrix} G (ρ_{m}) = \sum_{n = 1}^{N} \frac{g (r_{n})}{π P^{2} J_{1}^{2} (α_{n})} J_{0} (\frac{α_{n} α_{m}}{S}) \\ g (r_{n}) = \sum_{m = 1}^{N} \frac{G (ρ_{m})}{π R^{2} J_{1}^{2} (α_{m})} J_{0} (\frac{α_{m} α_{n}}{S}) \end{matrix} .

By defining an

N \times N

transform matrix

C_{N, N}

whose

(m, n)

th element is

C_{m, n} = \frac{J_{0} (\frac{α_{n} α_{m}}{α_{N + 1}})}{J_{1}^{2} (α_{m})},

we can simplify Eqs. (31) in terms of matrix multiplications on

g

and

G

of

{\begin{matrix} G {(ρ)}_{1, N} = \frac{1}{π P^{2}} [g {(r)}_{1, N} \cdot C_{N, N}] \\ g {(r)}_{1, N} = \frac{1}{π R^{2}} [G {(ρ)}_{1, N} \cdot C_{N, N}] \end{matrix} .

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

d r_{zp} = 20 nm

, where one might want

〈 Δ r 〉 = 1 nm

in order to get a smooth intensity profile at the focus. Though the QDHT uses a nonuniform grid especially at low radial values, when

N

is large the sampling point spacings approach a mean value of

〈 Δ r 〉 = R / N

. 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

C_{N, N}

with

N = 300, 000

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 $f > z_{0}$ , where the Hankel product approach of Eq. (21) is preferred. To propagate an entire wavefield within a radial distance $R$ from the zone plate to the focal plane, we would require a transform matrix $C_{N, N}$ 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 $C_{N, M}$ with $M ≪ N$ , as shown in Fig. 3. For example, we might need only $M = 50$ points to see the detailed focal profile (including several Airy fringes) of our example zone plate at $〈 Δ r 〉 = 2 nm$ calculation grid size. This saves a factor of 3000 in computation time and required storage over the example given above. The choice of $M$ 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 $M = 100$ ; 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.

Fig. 3. Schematic of the partial transform matrix $C_{N, M}$ . If one wishes to calculate the wavefield over only a subset of $M$ points on the output plane, a nonsymmetric matrix $C_{N, M}$ [Eq. (32)] can be used to reduce computational time and memory requirements in Eq. (33).

Download Full Size | PDF

Fig. 4. Rapid QDHT calculation of the focusing of $λ = 0.124 nm$ (10 keV) x rays by a fully absorptive Fresnel zone plate with 320 μm diameter and $d r_{zp} = 20 nm$ outermost zone width. The image shows the intensity as a function of radial distance $r$ and axial and defocus distance $Δ z$ , with a grid spacing of $〈 Δ r 〉 = 2 nm$ in radius and $Δ z = 1 μm$ in defocus positions (intensities at negative radii are simply copied from the positive radii calculation points). The wavefield exiting the zone plate is used as the input wavefield to generate output wavefields at various distances.

Download Full Size | PDF

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 $a = 5 μm$ on a calculation grid of spacing $〈 Δ r 〉 = 10 nm$ extending to a distance of $R = 270 μm$ so that $N = 27, 000$ samples were involved. The far-field diffraction pattern for $λ = 0.124 nm$ x rays at $z = 50 cm$ involved a Fresnel number of 0.1 and was calculated using $M = 500$ calculation points at a spacing of $〈 Δ r^{'} 〉 = 115 nm$ . 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 $Z = 250 cm$ , the maximum error goes down to 0.08% as it is even farther away from having any non-Fraunhofer terms contributing.

Fig. 5. Comparison of the numerical QDHT calculation of far-field diffraction intensity of a pinhole using Eq. (21), versus the analytical result of the Airy pattern. The image on the left shows the radial intensity distribution along with the integral of intensity with radius, while the image on the right shows the intensity distribution along with the percentage difference from the analytical result. As is shown, the QDHT calculation with parameters as described in the text is accurate to within 0.12% of the expected analytical result.

Download Full Size | PDF

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 $d r_{zp} = 25 nm$ used to focus $λ = 0.124 nm$ (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 $J_{0}$ divided by $π$ times the outermost zone width, or $(3.83 / π) d r_{zp} = 1.22 d r_{zp}$ , which is 30.5 nm in this example, while the integrated energy fraction should approach $1 / π^{2} = 10.1 %$ [28]. For the numerical QDHT calculation of Eq. (21), $N = 27,000$ sampling points were used within a maximum radius $R = 54 μm$ , while in the output plane $M = 100$ points were used at a spacing of $〈 Δ r 〉 = 2 nm$ . Again, the numerical results are in good agreement with the expected values.

Fig. 6. Comparison of the numerical QDHT calculation of the radial intensity distribution near the focus of a Fresnel zone plate ( $D = 45 μm$ diameter, $d r_{zp} = 25 nm$ outermost zone width using $λ = 0.124 nm$ or 10 keV x rays, yielding a focal length of $f = D d r_{zp} / λ = 9 mm$ ) using Eq. (21), versus the analytical result of the Airy pattern. The image on the left shows the radial intensity distribution along with the integral of intensity with radius, while the image on the right shows the intensity distribution along with the percentage difference from the analytical result.

Download Full Size | PDF

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 $λ = 0.124 nm$ wavelength was propagated a distance of 50 cm downstream, using a calculation grid with $N = 27,000$ points at $〈 Δ r 〉 = 10 nm$ spacing on the input plane, and $M = 440$ 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.

Fig. 7. Verification that the QDHT preserves energy over sufficiently large radii of integration. The image on the left shows the intensity calculated from propagating a $λ = 0.124 nm$ wavefield through a 5 μm diameter pinhole to a distance 50 cm downstream (Fresnel number $F = 0.1$ ). The image on the right shows the calculation for a binary absorption Fresnel zone plate with a diameter of 45 μm and outermost zone width of $d r_{zp} = 25 nm$ , at a distance of 9 mm which is one focal length away. As can be seen, about $1 / π^{2} = 10 %$ of the total energy is located on the optical axis in the form of the first focal order, and about half of the transmitted energy is located near the optical axis in the form of the positive (converging) focal orders. The integrated intensity increases up to the point of the radial extent of the zone plate at 22.5 μm, with the remaining energy arriving in the -1 focal order (extending to twice the radius or 45 μm) until nearly all of the 50% nonabsorbed energy is contained within a 60 μm radius, as expected.

Download Full Size | PDF

5. COMBINED PROPAGATION WITHIN AND BEYOND $z_{0}$

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 $z_{0}$ of Eq. (17). Consider the case of first-order focusing from a Fresnel zone plate, where one can write $z_{0} = f (〈 Δ r 〉 / d r_{zp})$ , where $f$ is the focal length of the zone plate, $〈 Δ r 〉$ is the sampling interval in real space, and $d r_{zp}$ is the width of the finest, outermost zone of the zone plate. Because one desires to have $〈 Δ r 〉$ be much smaller than $d r_{zp}$ in order to have many samples within the width of the outermost zone, the focal length $f$ is always in the far-distance location (that is, $f$ is large compared to $z_{0}$ ). 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.

Fig. 8. Comparison of the focused intensity profile of a Fresnel zone plate (which is a far-distance calculation) as calculated using the near-distance convolution approach of Eq. (20), and the far-distance single Hankel transform approach of Eq. (21). The radial intensity distribution and energy integral is shown on the left, while the longitudinal intensity distribution about the focal point is shown on the right. These calculations assumed a zone plate with diameter $D = 45 μm$ and outermost zone width $d r_{zp} = 25 nm$ , and an x-ray wavelength of $λ = 0.124 nm$ , yielding a focal length of $f = D d r_{zp} / λ = 9 mm$ . The input plane sampling was done with $N = 27,000$ points over a radius of $R = 54 μm$ , so that the distance $z_{0}$ for preferring near-distance or far-distance approaches [Eq. (17)] was $z_{0} = 1.7 mm$ . 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. The far-distance calculation approach works better for propagating by a distance of 9 mm when $z_{0} = 1.7 mm$ , but the near-distance approach still gives the correct overall intensity distribution at points away from the optical axis.

Download Full Size | PDF

Sometimes it is desirable to do calculations with multiple propagations over various distances. One might wish to propagate a wavefield $ψ_{0} (r_{0})$ to $ψ_{z} (r)$ over a distance beyond $z_{0}$ 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 $r_{0}$ and $r$ , respectively, we will write their discrete counterparts as $r_{n}$ and $r_{n}^{'}$ .

For propagation over shorter distances $z < z_{0}$ , the convolution approach of Eq. (20) samples an input wavefield $ψ_{0} (r_{n})$ , transforms it to reciprocal space as $Ψ (ρ_{n})$ where it is multiplied by a real space propagator [Eq. (13)], and inverse transforms it back to real space as $ψ_{z} (r_{n}^{'})$ on the same calculational grid. If we sample the input wavefield in $N$ points over a radial extent $R$ , the real space sampling interval is $〈 Δ r 〉 = R / N$ , and the maximum value of the Hankel function argument is $P = α_{N + 1} / 2 π R = S / 2 π R$ [Eq. (27)]. When $N$ is large, $S$ can be chosen as $α_{N} ≃ π N$ , so

P = \frac{π N}{2 π R} = \frac{N}{2 R}

becomes the maximum spatial frequency, which when divided by the number of samples

N

yields an interval

Δ ρ

in reciprocal space of

Δ ρ = \frac{1}{2 R} .

Multiplication with the propagator

\exp [i π λ z ρ_{n}^{2}]

of Eq. (13) modifies

Ψ (ρ_{n})

but does not change its sample positions. The inverse QDHT brings the wavefield back to real space, and the discretely sampled, propagated wavefield

ψ_{z} (r_{n}^{'})

extends to a maximum radius

R^{'} = \frac{π N}{2 π P} = R

with interval

〈 Δ r^{'} 〉 = 〈 Δ r 〉

. Therefore, the new array

ψ_{z} (r_{n}^{'})

is in real space with unchanged sampling and extent.

For propagation over longer distances $z > z_{0}$ , the single Hankel transform method of Eq. (21) is preferred. We start with the discrete 1D array $ψ_{0} (r_{n})$ in real space and multiply it by the real space phase propagator $\exp [- i π r_{n}^{2} / (λ z)]$ of Eq. (12). We then perform the QDHT to bring the wavefield into reciprocal space $Ψ (ρ_{n})$ , with the same extent $P$ as Eq. (34) and interval $Δ ρ$ as Eq. (35). In this approach there is no second QDHT; instead, the reciprocal space array $Ψ (ρ_{n})$ is multiplied by $(i / (λ z)) \exp [- i π r_{n}^{' 2} / (λ z)]$ where the real space positions are found from $r_{n}^{'} = (λ z) ρ_{n}$ . The propagated wavefield $ψ_{z} (r_{n}^{'})$ extends to a maximum radius of

R^{'} = (λ z) P = \frac{λ z}{2 R} N,

with sampling interval

〈 Δ r^{'} 〉 = \frac{R^{'}}{N} = \frac{λ z}{2 R} = \frac{λ z}{2 N 〈 Δ r 〉} .

Obviously, the output plane real space sampling interval

〈 Δ r^{'} 〉

might differ from the input plane interval

〈 Δ r 〉

. As an example, consider a zone plate with diameter

D

and outermost zone width

d r_{zp}

, where again the paraxial approximation gives

D d r_{zp} = λ f

. For propagation of a wavefield from the zone plate to the focus, or

z = f

, Eq. (37) becomes

〈 Δ r^{'} 〉 = \frac{D d r_{zp}}{2 R} = \frac{D d r_{zp}}{2 N 〈 Δ r 〉} .

If we want to keep the sampling interval constant, or

〈 Δ r^{'} 〉 = 〈 Δ r 〉

, we need to choose the number of sampling points

N

to be

N = \frac{D d r_{zp}}{2 {(〈 Δ r 〉)}^{2}} .

Equation (37) can be used to choose

N

for maintaining

〈 Δ r^{'} 〉 = 〈 Δ r 〉

at other propagation distances

z

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 $z_{0}$ , 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 $M$ on a fine sampling interval $〈 Δ r^{'} 〉$ 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).

Rapid calculation of paraxial wave propagation for cylindrically symmetric optics

Abstract

1. INTRODUCTION

2. QUASI-DISCRETE HANKEL TRANSFORM AND THE SAMPLING THEOREM

3. RAPID CALCULATION USING PARTIAL TRANSFORMS

4. COMPARISONS AGAINST ANALYTICAL CALCULATIONS

5. COMBINED PROPAGATION WITHIN AND BEYOND $z_{0}$

6. CONCLUSION

Funding

Acknowledgment

REFERENCES

Cited By

Figures (8)

Equations (39)

Journal of the Optical Society of America A

Abstract

1. INTRODUCTION

2. QUASI-DISCRETE HANKEL TRANSFORM AND THE SAMPLING THEOREM

3. RAPID CALCULATION USING PARTIAL TRANSFORMS

4. COMPARISONS AGAINST ANALYTICAL CALCULATIONS

5. COMBINED PROPAGATION WITHIN AND BEYOND z0

6. CONCLUSION

Funding

Acknowledgment

REFERENCES

Cited By

Figures (8)

Equations (39)

Journal of the Optical Society of America A

5. COMBINED PROPAGATION WITHIN AND BEYOND $z_{0}$