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

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 ₀ = [2exp(α)+exp(2α)]/{[1+exp(α)]² [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,…,2N-1 according to Eq. (7), therefore the values of g_γ (ρ_m ) from N to 2N-1 are insignificant in mathematics. The parameters r ₀ and α can be chosen arbitrarily in theory, however, for better numerical results, we take the same values of r ₀ =[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

 

Fig. 1. (a)First-order Hankel transform g(y) of f(r)=r (0≤r≤1) for a constant C = 20π ; (b) Difference (error) between the analytical transform and the results obtained with the numerical first-order fast hankel transformations.

Download Full Size | PDF

 

Fig.2. (a)Second-order Hankel transform g(y) of _{f(r)=r⁴ (0≤r≤1)} for a constant C = 20π ; (b) Difference (error) between the analytical transform and the results obtained with the numerical second-order fast hankel transformations.

Download Full Size | PDF

whose analytical solutions can be easily determined. In this numerical experiment, we implemented the calculations on the Hankel transforms H ₁{r} and H ₂{r ⁴} 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 ₁(ρ)and g ₂(ρ). 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 ₀ can be determined with those in the near field z = 0 as follows

where j = √-1 , $k_z={\left[k_0^2-\left(k\prime k_0\right)\right]}^{\frac{1}{2}},$, k ₀ = 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 ₀Δ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.

 

Fig. 3 Electric field magnitude at the emergent plane of the lens.

Download Full Size | PDF

 

Fig. 4 Electric field magnitude in the focal plane of the lens.

Download Full Size | PDF

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.