## Abstract

We present a fast Hankel transform (FHTn) method for direct numerical evaluation of electromagnetic (EM) field propagation through an axially symmetric system. Comparing with the vector-based plane-wave spectrum (VPWS) method, we present an alternative approach to implement the fast Hankel transform which does not require an additional coordinate transformation for Fourier transform. The proposed FHTn method is an efficient approach for numerical evaluation of an arbitrary integer order of the Hankel transform (HT). As an example to demonstrate the effectiveness of the proposed method, we apply the FHTn technique to the analysis of cylindrical EM field propagation through a diffractive microlens.

© 2002 Optical Society of America

## 1. Introduction

It is known that the plane-wave spectrum (PWS) method or alternatively known as an angular spectrum method, is an efficient technique for studying the propagation of electromagnetic wave from one plane to another by determining the Fourier transform (FT) of the incident field and a multiplicative transfer function over a propagation distance [1]. In numerical analysis of beam propagation through an axially symmetric system, the Fourier transform is replaced by the Hankel transform (HT) with separable variables in polar coordinates.

In 1999, Shi and Prather [2] reported a vector-based PWS (VPWS) method for the study of the propagation of a cylindrical electromagnetic field that involves reduplicate numerical evaluations of arbitrary integer orders of HT other than the only 0th-orders in many cases. In 1992, Magni et al [3] proposed an efficient fast Hankel transform algorithm (FHATHA) to implement numerical evaluations of the 0th-order HT. However, the FHATHA method has limited applications in the VPWS method, so it has to require an additional step of coordinate transformation when studying the propagation of the cylindrical wave in the VPWS method. A fast Hankel transform of n-th order has been reported by Ferrari et al [4] in 1999. However, when the order of the transform increases, the sampling frequency of the transform function must be changed in Ferrari’s algorithm. A numerical integration is involved in the algorithm. Since the analytical solution of the integration in Ferrari’s algorithm can not be given directly, numerical evaluation approach should be introduced and thus the computational speed is limited. The FHTn algorithm of this paper does not require the computations of derivatives and numerical integrals, and there is no sampling frequency problem.

In this paper, we present a new approach for the implementation of a fast Hankel transform for direct evaluation of an arbitrary integer order of HT using the VPWS method. The proposed FHTn method is an efficient technique for vector analysis of circular symmetric diffractive optical elements (DOEs) by studying the propagation of an EM field through an arbitrary distance, which can be determined by the finite-difference time-domain (FDTD) method [5] in the near field region.

As an example, the proposed FHTn method is employed here to perform a numerical experiment to demonstrate its effectiveness in the study of the propagation of an EM field emerging from a diffractive microlens to its focal plane.

## 2. Theory

It is well known that Fourier transform of a circular symmetric function with separable variables can be written as γth-order HT as follows

where *J*_{γ}
is the γth-order Bessel
function and *C* is a constant. We use
*H*_{γ}
{*f*(*r*)}
to represent the γth-order HT of a function
*f*(*r*), so Eq. (1) can be re-written as

The FHTn method begins with the Bessel function property:

By substituting Eq. (3) into Eq. (1) with an assumption that the integration is implemented
numerically, we can simplify the integral and reduce it to the evaluation of the
integral over a small sampling interval. Consequently, in each discrete interval
*ξ*_{i}
< *r*
<
*ξ*
_{i+1}, the value
of *f*(*r*)/*r*^{γ}
can be determined by
*f*(*r*_{n}
)/${r}_{n}^{\mathrm{\gamma}}$
at a sampling point *r*
_{n}. Thus the integral over each
interval is given by

$$\phantom{\rule{12em}{0ex}}=\frac{f\left({r}_{n}\right)}{C\xb7{r}_{n}^{\gamma}\rho}\xb7\left[{\xi}_{n+1}^{\gamma +1}\xb7{J}_{\gamma +1}\left(\mathit{C\rho}\xb7{\xi}_{n+1}\right)-{\xi}_{n}^{\gamma +1}\xb7{J}_{\gamma +1}\left(\mathit{C\rho}\xb7{\xi}_{n}\right)\right]$$

Following the approach of Ref. [3], we can rewrite the HT integral of Eq. (1) as a discrete cross-correlation by summing the discrete integrals. The discrete cross-correlation can be numerically evaluated by

where *F* {} represents the Fourier transform (FT) and *F*
^{-1}{}the inverse FT. In Eq. (5), the functions
*j*
_{(γ+1)n}
and *φ*_{n}
are defined as

for *n* = 0, 1,…,2N-1,and

where *f*(*r*_{N}
) ≡ 0,
*F*_{n}
=
*f*(*r*_{n}
)/${r}_{n}^{\gamma}$ -
*f*(*r*
_{n+l})/${r}_{n+1}^{\gamma}$, *k*
_{0} =
[2exp(*α*)+exp(2*α*)]/{[1+exp(*α*)]^{2}
[1-exp(-2*α*)]} , and *N* is the total
number of sampling points. It is seen that the FHTn can be reduced to the FHATHA for
the case of *γ*=0 in Eq. (7). In numerical calculation we only need to reserve the values
of *g*_{γ}
(*ρ*_{m}
) from *m* = 0 to
*N*-1 in Eq. (5) because we pattern the value of
*φ*
_{n} to be zero for *n* =
*N*,
*N*+1,…,2*N*-1 according to Eq. (7), therefore the values of
*g*_{γ}
(*ρ*_{m}
) from *N* to
2*N*-1 are insignificant in mathematics. The parameters
*r*
_{0} and *α* can be chosen
arbitrarily in theory, however, for better numerical results, we take the same
values of *r*
_{0}
=[1+exp(*α*)]exp(-*αN*)/2
and *α* =
-ln[1-exp(-*α*)](*N*-1) as used in Ref [3].

## 3. Numerical experiments

To test the numerical accuracy of the proposed algorithm, we consider two special functions

whose analytical solutions can be easily determined. In this numerical experiment, we
implemented the calculations on the Hankel transforms
*H*
_{1}{*r*} and
*H*
_{2}{*r*
^{4}} with a constant
*C* = 20*π*. The numerical calculations
were completed within 0.01 second on a 450MHz PC with 128MByte RAM for 1024 sampling
points. Figures 1(a) and 2(a) show the plots of both numerical results *g*_{1}(*ρ*)and *g*_{2}(*ρ*). The errors between the numerical
results and the analytical solutions are plotted in Figures 1(b) and 2(b). The good agreement between the two results show that
the FHTn method is highly accurate, and it can thus be used as a direct Hankel
transform approach for studying the propagation through an axially symmetric
system.

## 4. Application for VPWS

We now apply the FHTn method to demonstrate its effectiveness in studying the EM
field distribution on the focal plane of a diffractive lens. The lens is normally
illuminated with an incident plane wave. For axially symmetric structures, the EM
components can be represented in terms of a cylindrical harmonic expansion in the
azimuthal direction. From Ref. [2], the expansion coefficients of the EM components in the
focal plane *z* = *z*
_{0} can be determined
with those in the near field *z* = *0* as follows

$$\phantom{\rule{11em}{0ex}}\times \phantom{\rule{.2em}{0ex}}\sum _{i=0}^{M}i\{\left[{E}_{r}^{m}(i\xb7\Delta r,0)+{E}_{\varphi}^{m}(i\xb7\Delta r,0)\right]$$

$$\phantom{\rule{12em}{0ex}}\times \phantom{\rule{.2em}{0ex}}{k}_{0}^{2}{H}_{m+1}\left\{{J}_{m+1}\left(k\prime \xb7l{k}_{0}\Delta \phantom{\rule{.2em}{0ex}}\rho \right)\mathrm{exp}\left(-j{k}_{z}{z}_{0}\right)\right\}$$

$$\phantom{\rule{8em}{0ex}}\phantom{\rule{.2em}{0ex}}+\left[{E}_{r}^{m}(i\xb7\Delta r,0)-{E}_{\varphi}^{m}(i\xb7\Delta r,0)\right]$$

$$\phantom{\rule{12em}{0ex}}\times \phantom{\rule{.2em}{0ex}}{k}_{0}^{2}{H}_{m-1}\left\{{J}_{m-1}\left(k\prime \xb7l{k}_{0}\Delta \phantom{\rule{.2em}{0ex}}\rho \right)\mathrm{exp}\left(-j{k}_{z}{z}_{0}\right)\right\}\}$$

$$\phantom{\rule{11em}{0ex}}\times \phantom{\rule{.2em}{0ex}}\sum _{i=0}^{M}i\{\left[{E}_{\varphi}^{m}(i\xb7\Delta r,0)+{E}_{r}^{m}(i\xb7\Delta r,0)\right]$$

$$\phantom{\rule{12em}{0ex}}\times \phantom{\rule{.2em}{0ex}}{k}_{0}^{2}{H}_{m+1}\left\{{J}_{m+1}\left(k\prime \xb7l{k}_{0}\Delta \phantom{\rule{.2em}{0ex}}\rho \right)\mathrm{exp}\left(-j{k}_{z}{z}_{0}\right)\right\}$$

$$\phantom{\rule{8em}{0ex}}+\phantom{\rule{.2em}{0ex}}\left[{E}_{\varphi}^{m}(i\xb7\Delta r,0)-{E}_{r}^{m}(i\xb7\Delta r,0)\right]$$

$$\phantom{\rule{12em}{0ex}}\times \phantom{\rule{.2em}{0ex}}{k}_{0}^{2}{H}_{m-1}\left\{{J}_{m-1}\left(k\prime \xb7l{k}_{0}\Delta \phantom{\rule{.2em}{0ex}}\rho \right)\mathrm{exp}\left(-j{k}_{z}{z}_{0}\right)\right\}\}$$

$$\phantom{\rule{2em}{0ex}}\times \phantom{\rule{.2em}{0ex}}\sum _{i=0}^{M}i\xb7\{{E}_{z}^{m}\left(i\xb7\Delta r,0\right)$$

$$\phantom{\rule{9em}{0ex}}\times \phantom{\rule{.5em}{0ex}}{k}_{0}^{2}{H}_{m}\left\{{J}_{m}\left(k\prime \xb7l{k}_{0}\Delta \phantom{\rule{.2em}{0ex}}\rho \right)\text{}\phantom{\rule{.2em}{0ex}}\times \phantom{\rule{.2em}{0ex}}\mathrm{exp}(j{k}_{z}\xb7{z}_{0})\right\}$$

where *j* = √-1 , ${k}_{z}={\left[{k}_{0}^{2}-\left(k\prime {k}_{0}\right)\right]}^{\frac{1}{2}},$, *k*
_{0} =
2*π*/*λ*, 0 ≤
*k*’ ≤ 1, Δ*r*
and Δ*ρ* are the sampling lengths in the
original and propagated planes, respectively, and λ is the wavelength of
illumination light. The index *i* and *l* represent
the values of the *i*-th sampling point in the near field and the
*l*-th in the focal plane, respectively. *M* is
the total number of sampling points in the near field, and the constant
*C* in the process of HT is equal to
*Nk*
_{0}Δ*r*.

In this example, since the incident wave is a normal illumination, only the
*m* = 1 mode is needed. Thus, only zero-, first-, and second-order
HTs are used in the calculation. For oblique incidence, higher-order HT will be
involved in the VPWS calculation, and in this case, the FHTn method can still be
used because it is a generic formula for an arbitrary integer order HT.

For convenience of comparison, we used the same eight-level diffractive lens that was
studied in Ref. [5]. The lens parameters include relative permittivity
∊_{r} = 2.25; wavelength λ = 1.0 μm;
diameter = 102.47 μm; *f*-number = 0.78; and 15 zones.
Firstly, we employed our axially symmetric finite-difference time-domain (FDTD)
program to obtain the near electric field and then used Eqs. (8) to (10) to calculate the electric field propagated onto the focal
plane. Figure 3 shows the near electric field magnitude at the
emergent plane of the lens as determined by the axially symmetric FDTD method. Figure 4 shows the numerical result of the electric field
distribution at the focal plane of the lens, which is the same as that determined
with the EM method presented in Ref. [5]. For this example, the VPWS method with FHTn took 41s and
only 1M byte to propagate, compared with 913s and 897M byte for the conventional EM
method. It should be noted that the proposed FHTn method is an efficient alternative
to implement the VPWS analysis since no coordinate transformations are required in
the calculation.

## 5. Summary

We have presented a fast arbitrary integer order FHTn method for direct numerical evaluation of the field propagation through an axially symmetric system. It is an efficient and fast technique for evaluating arbitrary integer order HT. Comparing with Prather’s VPWS method, the proposed method does not require a coordinate transformation for an axially symmetric problem.

## References and Links

**1. **J. W. Goodman, *Introduction to Fourier Optics* (McGraw-Hill, New York,1968), Chap. 2, pp. 11–13.

**2. **Shouyuan Shi and Dennis W. Prather, “Vector-based plane-wave spectrum method for the propagation of cylindrical electromagnetic fields,”Opt. Lett. **24**, 1445 (1999). [CrossRef]

**3. **Vittorio Magni, Giulio Cerullo, and Sandro De Silvestri, “High-accuracy fast Hankel transform for optical beam propagation,”J. Opt. Soc. Am. A. **9**, 2031–2033 (1992). [CrossRef]

**4. **José A. Ferrari, Daniel Perciante, and Alfredo Dubra, “Fast Hankel transform of n-th order,”J. Opt. Soc. Am. A. **16**, 2581–2582 (1999). [CrossRef]

**5. **Dennis W. Prather and Shouyuan Shi, “Formulation and application of the finite-difference time-domain method for the analysis of axially symmetric diffractive optical elements,”J. Opt. Soc. Am. A. **16**, 1131–1142 (1999). [CrossRef]