Comparison of iterative angular spectrum and optimal rotation angle methods in designing beam-fanners

Seyyed H. Kazemi; Mir M. Mirsalehi; Amir R. Attari

doi:10.1364/OE.17.014825

1. Introduction

Several methods have been used to design diffractive optical elements (DOEs). These methods can be generally divided to unidirectional and bidirectional methods. The unidirectional methods combine the analysis of DOE with optimization techniques. Examples of analysis methods are scalar phase only approximation and rigorous electromagnetic methods such as finite difference time domain (FDTD) [1] and boundary element method [2, 3]. The analysis methods often use additional methods to model the wave propagation from the output plane (the plane just behind the DOE) to the observation plane, e.g. Fresnel and Fraunhofer propagation approach or angular spectrum (AS) approach [4]. The optimization methods, such as micro-genetic [1] or simulated annealing [3], search the space of unknown DOE parameters using fitness function. Another unidirectional method is the optimal rotation angle (ORA) [5, 6]. In this method, the phase of each cell in the DOE is changed to maximize the intensities at a limited number of spots.

In bidirectional methods, analysis and inverse-analysis are used iteratively to improve the DOE structure. Gerchberg–Saxton method [7], iterative Fourier transform algorithm (IFTA) [8], the method proposed by Di [9], vectorial beam shaping by Jabbour [10] and iterative angular spectrum algorithm (IASA) [11] are various kinds of bidirectional methods. The Gerchberg–Saxton method and the IFTA method are used for far and near field cases. The method proposed by Di uses FDTD and thin element approximation (phase only approximation) as analysis and inverse-analysis methods, respectively and AS approach for wave propagation simulation. In the IASA, the optical field in the output plane is calculated using phase only approximation. This field is then propagated to the observation plane by AS approach.

In this paper, we compare the IAS and ORA methods in designing FADOEs as 2-D beamfanners. Since the original IAS method was presented for 1-D case, we have extended it for the 2-D case. The extended iterative angular spectrum (IAS) method is described in section 2 that is followed by a description of the ORA method. The transfer functions of these methods are compared in section 4. Both methods, as well as the FDTD method, are used to design and analyze 1-to-4 and 1-to-6 beamfanners in section 5 which is followed by conclusions.

2. Iterative angular spectrum method

The geometry of a 2-D finite aperture DOE is shown in Fig. 1 . The input plane that includes the aperture is located in front of the DOE and the observation plane is located at z=z _obs. It is assumed that a plane wave propagating along the z-axis interacts with the DOE at the input plane and it is expected to produce the desired pattern at the observation plane.

Fig. 1 Geometry of 2D finite aperture DOE and the two regions used in the IAS method.

Download Full Size | PDF

The space between the input and the observation plane is divided to two regions: DOE region and AS region. The DOE region extends from the input plane to the output plane (P₀) and includes DOE, while the AS region extends from the output plane located at z=0 to the observation plane (P₁). The thickness of DOE region is equal to the maximum etch depth (d _max). The refraction index of the dielectric material used is n ₂, while the surrounding medium is assumed to be air (n ₁=1).

The transmission function of the DOE is

t(x,y) = τ exp[j ϕ (x, y)] Π (x / L_{x}, y / L_{y}),

where τ=2n ₁/(n₁+n ₂) is Fresnel transmission coefficient, Π(x/L_x, y/L_y) is equal to 1 for any point inside a rectangular aperture of L_x×L_y, and equal to zero for other points. The phase function according to phase only approximation is φ(x,y)=k ₀Δnd(x,y), where Δn is the difference between the two refractive indices, k ₀=2π/λ, λ is the wavelength in free space and d(x,y) is the etch depth. We assume an incident plane wave, u _inc, that is propagating normal to the aperture. The wave at the output plane of the first region is u ₀=t(x,y)u _inc. The angular spectrum of u ₀, i.e., U ₀, is defined as [4]

U_{0} (f_{x}, f_{y}) = \iint_{P_{0}} u_{0} (x, y) \exp [- j 2 π (f_{x} x + f_{y} y)] d x d y .

Using the fast Fourier transform (FFT) technique, U ₀ can be computed. Assuming a square aperture, L_x=L_y=L, the distance between two adjacent points in spatial frequency domain is δf=δf_x=δf_y=1/L. To improve the resolution in the frequency domain, we can compute the FFT on an area larger than the DOE area. This is particularly useful if the desired intensity profile in the observation plane extends beyond the DOE area.

The field distribution in the observation plane, u ₁, can be obtained by

u_{1} (x, y) = \int \int_{- \infty}^{\infty} U_{0} \exp [j k_{2} z_{o b s} \sqrt{1 - {(λ_{2} f_{x})}^{2} - {(λ_{2} f_{y})}^{2}}] \exp [j 2 π (f_{x} x + f_{y} y)] d f_{x} d f_{y},

where k ₂=n ₂ k ₀ and λ₂=λ/n ₂ are the wave-number and wavelength in the AS region, respectively [4]. Using the inverse fast Fourier transform (IFFT), u ₁ can be calculated as u ₁=IFFT(HU ₀), where

H(f_{x}, f_{y}) = \exp [j k_{2} z_{o b s} \sqrt{1 - {(λ_{2} f_{x})}^{2} - {(λ_{2} f_{y})}^{2}}]

is the angular spectrum transfer function. The maximum spatial frequency is f _max=1/(2δ), where δ is the minimum feature size of the DOE. To consider all propagating waves, the condition f _max>1/λ₂, which is equivalent to δ <λ₂/2, should be satisfied.

If the intensity at the observation plane satisfies the desired specification, the design process is terminated and the obtained DOE is used. Otherwise, the amplitude (not phase) of the field at the observation plane, u ₁, is modified such that the difference between the obtained amplitude and the desired amplitude is reduced. The modified field distribution is then propagated backward to find the corresponding field distribution at the DOE output plane. Also, using the phase of the back going filed, the DOE structure is modified to satisfy the new field distribution. This modified DOE is quantized and used to repeat the above procedure.

3. Optimal rotation angle method

In this section, the ORA method is briefly described. For detailed information, the reader is referred to [5] and [6]. The geometry used in this method is shown in Fig. 2 . The output of the DOE is assumed to be on the left plane and it is divided to small cells. The output field at the center of the k-th cell is

u_{0, k} = A_{k}^{inc} \exp [j (ϕ_{k}^{inc} + ϕ_{k})],

where

A_{k}^{inc}

and

ϕ_{k}^{inc}

are the amplitude and phase of the incident field in the center of the k-th pixel, (x_k , y_k), respectively, and φ_k is the amount of phase modulation produced by the DOE. Considering the m-th spot in the observation plane, located at (x_m , y_m , z _obs), the field at this spot can be represented as

u_{1, m} = \sum_{k} u_{0, k} h_{k, m},

where h_k,m is a complex value quantity that relates the field at the k-th cell of the DOE to the field at the m-th spot of the observation plane. Using the Helmholtz-Kirchhoff integral, it is obtained as [5]

h_{k, m} = \frac{1}{4 π} [- j k_{2} - z_{o b s} (j k_{2} - 1 / r_{k m}^{c}) / r_{k m}^{c}] \frac{4 \exp (j k_{2} r_{k m}^{c})}{r_{k m}^{c}} \frac{\sin ({\tilde{k}}_{x} \frac{a}{2}) \sin ({\tilde{k}}_{y} \frac{b}{2})}{{\tilde{k}}_{x} {\tilde{k}}_{y}},

where

r_{k m}^{c}

=[(x_k-x_m)²+(y_k-y_m)²+z _obs ²]^½,

{\tilde{k}}_{x} = k_{2} (x_{k} - x_{m}) / r_{k m}^{c},

{\tilde{k}}_{y} = k_{2} (y_{k} - y_{m}) / r_{k m}^{c},

and a and b are the x and y dimensions of each DOE cell, respectively.

Fig. 2 Geometry of DOE and the spots used in the ORA method.

Download Full Size | PDF

According to Bengtesson's ORA method, the intensities of the observation spots will be maximized if the phase of each DOE cell is changed by Δφ_k and this procedure is repeated until the phases of the cells are converged. For the formula of Δφ_k, the reader is reffered to [5].

4. Transfer functions

As mentioned in the previous sections, ORA and IAS methods are scalar techniques that compute the field distribution at the observation plane in two different ways. We compare the transfer functions of these two methods as z _obs is changed. We use Eq. (4) as the transfer function of IAS method. For the ORA method, we start with the Helmholtz-Kirchhoff integral [4] and obtain the spatial impulse response of the propagation medium between DOE and observation plane as

h (x, y; z_{o b s}) = \frac{1}{4 π} [- j k_{2} - \frac{z_{o b s}}{r} (j k_{2} - \frac{1}{r})] \frac{\exp (j k_{2} r)}{r},

where r=(x ²+y ²+z _obs ²)^1/2. The above spatial impulse response is a continuous function of space, while h_k,m in Eq. (7) is a discrete quantity. Equation (7) can be obtained from Eq. (8) using sampling theorem [4]. Next, u ₁ can be represented by using the 2-D convolution as

u_{1} (x, y) = \iint_{P_{0}} u_{0} (η, ξ) h (x - η, y - ξ; z_{o b s}) d η d ξ .

We use the 2-D fast Fourier transform (FFT) to calculate the ORA transfer function in the spatial frequency domain. The phase of both transfer functions are the same, but their amplitudes are different.

The differences in the amplitude can be divided in two types. First, for z_obs smaller than 10λ the frequencies near edge are amplified. Second, for z _obs greater than 10λ, the spatial frequency bandwidth of the ORA transfer function becomes smaller (See Fig. 3 ). The first difference produces high spatial frequency components in the DOE structures designed by ORA method, which can destructively affect the field distribution. Figure 4 shows the effect of these high spatial frequency components in two typical DOE profiles. The second difference eliminates high spatial frequencies that are important specially in cases that the spots are located at a large distance from each other compared with the DOE dimensions.

Fig. 3 The middle cross sections of the transfer functions for IAS (^...) and ORA (^___) methods at z_obs=10λ, 30λ and 100λ

Download Full Size | PDF

Fig. 4 The middle cross sections of typical DOE profiles designed by (a) IAS and (b) ORA methods for z_obs=10λ case.

Download Full Size | PDF

Therefore, the designs obtained by the IAS method for the near field cases are superior to the corresponding designs obtained by the ORA method. For a specific pattern, as the distance between the observation plane and the DOE is increased, the results obtained by the two methods become more similar.

5. Design and analysis

We investigate the effect of high spatial frequency components in the design of 1-to-4 and 1-to-6 beamfanners by IAS and ORA methods. The DOE is assumed to be a square of dimension D=10λ₀ made of silicon with a refractive index of n ₂=3.4 and λ₀ is set to 5μm. The spatial quantization steps along the three directions (x, y and z), i.e., a, b, and δ, and the etch depth quantization step are chosen equal to λ₀/25=.0.2μm. As mentioned in section 2, to improve the frequency resolution, we put additional zeroes around the above matrix to increase its dimension by a factor of 10. As a result, the spatial frequency resolution will be δ_f=1/(100λ₀) and therefore, the transmission function is represented by a 2500×2500 complex-valued matrix.

In designing 1-to-4 and 1-to-6 beamfanners by the ORA method, only 4 and 6 spots of equal intensities are defined at the observation plane, while in the IAS method, the desired intensity profile at the observation plane is used. We used a Gaussian shape as the desired intensity profile at the observation plane. The Gaussian functions were defined as exp{-[(x-x_m)²+(y-y_m)²]/(2w ₀ ²)]}, where (x_m,y_m) corresponds to the center of the m-th spot. In our simulations, we have used w ₀= 0.14λ. This width is very small compared to the distance between the spot points (3λ to 8λ). Therefore, the two distribution profiles are approximately the same in numerical calculations.

To analyze the designed DOEs, we used a 3-D FDTD algorithm with unsplit step PML absorbing boundary condition. The Analyzed region is divided to the total field and the scattered filed regions, as shown in Fig. 5 . We have used the technique introduced by Sullivan for implementing the PML [12, 13].

Fig. 5 Cross section of the analyzed region.

Download Full Size | PDF

The mesh size for the FDTD simulation was chosen equal for all three directions, i.e., δ_x=δ_y=δ_z=δ. Also, δt=δ/(2c), where c is the velocity of light in free space, was chosen as the time step. The depths of PML and scattered field were set to 7δ. The incident field was assumed to be a TEM wave with a sharp Gaussian time varying amplitude, E _xinc=exp{-(t-t _d)²/(2w _t ²)}, where t _d and w _t are the time delay and the spread parameter of the Gaussian pulse, respectively. Using the FDTD method, the distributions of the electric and magnetic fields were obtained for 200 time steps.

After the fields were found at P ₀ plane by the FDTD method, the angular spectrum technique was used to obtain the fields at the observation plane, i.e., P ₁. The distribution of light intensity, I, at P ₁ is then calculated. The diffraction efficiency defined as

DE = \frac{\iint_{W} I(x,y) d x d y}{\iint_{P_{1}} I(x,y) d x d y}

is then computed, where W is the detectors window which consists of four (six) circular areas of radius equal to λ₀, centered at the four (six) spots positions.

The diffraction efficiency (DE) obtained by FDTD method are provided in Table 1 for 1-to-4 and 1-to-6 beamfanners, where S₄ and S₆ are separation distance of adjacent spots in the 1-to-4 and 1-to-6 beamfanners, respectively. It can be seen, that the IAS method provides better designs than the ORA.

Table 1. The diffraction efficiency (DE) obtained by FDTD method for designed 1-to-4 and 1-to-6 beamfanners.

View Table

6. Conclusions

We compared the transfer functions of the ORA and IAS methods in the near field region. This comparison was made in spatial frequency domain. We showed that, there are two weak points in the ORA transfer functions including 1) stronger high spatial frequency components and 2) smaller frequency bandwidths. By choosing the distance between the observation plane and the FADOE plane (z _obs) as an important variable, we found, that if z _obs is less than or equal to 10λ, the high spatial frequency components appear stronger in the ORA method. This results in increasing the number of discontinuities in the DOE phase profile. Vice versa, by increasing z _obs, the ORA frequency bandwidth becomes smaller. This puts a limit in the separating distance between the spots.

To investigate these concepts, we designed several 1-to-4 and 1-to-6 beamfanners using both methods and compared their diffraction efficiencies. In all designs, the IAS method outperforms the ORA method in the near field region.

Acknowledgments

The authors are thankful to Iran Telecommunication Research Center (ITRC) for supporting this research.

References and links

1. J. Jiang and G. P. Nordin, “A rigorous unidirectional method for designing finite aperture diffractive optical elements,” Opt. Express 7(6), 237–242 (2000). [CrossRef] [PubMed]

2. D. W. Prather, M. S. Mirotznik, and J. N. Mait, “Boundary integral methods applied to the analysis of diffractive optical elements,” J. Opt. Soc. Am. A 14(1), 34–43 (1997). [CrossRef]

3. D. W. Prather, J. N. Mait, M. S. Mirotznik, and J. P. Collins, “Vector-based synthesis of finite aperiodic subwavelength diffractive optical elements,” J. Opt. Soc. Am. A 15(6), 1599–1607 (1998). [CrossRef]

4. J. W. Goodman, Introduction to Fourier Optics (Roberts & Company, Englewood, 2005).

5. J. Bengtsson, “Design of fan-out kinoforms in the entire scalar diffraction regime with an optimal-rotation-angle method,” Appl. Opt. 36(32), 8435–8444 (1997). [CrossRef]

6. J. Stigwall and J. Bengtsson, “Design of array of diffractive optical elements with inter-element coherent fan-outs,” Opt. Express 12(23), 5675–5683 (2004). [CrossRef] [PubMed]

7. R. W. Gerchberg and W. O. Saxton, “Practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttg.) 35, 237–250 (1972).

8. F. Wyrowski and O. Bryngdahl, “Iterative Fourier-transform algorithm applied to computer holography,” J. Opt. Soc. Am. A 5(7), 1058–1065 (1988). [CrossRef]

9. F. Di, Y. Yingbai, J. Guofan, T. Qiaofeng, and H. Liu, “Rigorous electromagnetic design of finite-aperture diffractive optical elements by use of an iterative optimization algorithm,” J. Opt. Soc. Am. A 20(9), 1739–1746 (2003). [CrossRef]

10. T. G. Jabbour and S. M. Kuebler, “Vectorial beam shaping,” Opt. Express 16(10), 7203–7213 (2008). [CrossRef] [PubMed]

11. S. D. Mellin and G. P. Nordin, “Limits of scalar diffraction theory and an iterative angular spectrum algorithm for finite aperture diffractive optical element design,” Opt. Express 8, 706–722 (2001). [CrossRef]

12. D. M. Sullivan, “An unsplit step 3-D PML for use with the FDTD method,” IEEE Microwave and Guided Wave Letters 7(7), 184–186 (1997). [CrossRef]

13. D. M. Sullivan, Electromagnetic Simulation Using the FDTD Method, IEEE Press series on RF and microwave technology, 2000.

1-to-4 Beamfanner					1-to-6 Beamfanner
	DE (%) @ S₄=8λ		DE (%) @ S₄=6λ		DE (%) @ S₆=4λ		DE (%) @ S₆=3λ
z _obs/λ	ORA	IAS	ORA	IAS	ORA	IAS	ORA	IAS
5	40.4	47.3	42.8	54.4	44.4	62.1	37.7	58.7
10	36.1	62.8	36.2	68.2	36.5	69.0	38.3	60.0
20	25.9	62.7	30.2	65.6	47.7	52.3	47.5	58.3
30	24.2	53.6	28.5	60.6	46.1	60.0	48.3	64.1
50	28.1	42.3	25.7	50.1	49.6	62.7	46.4	59.7
80	16.1	28.4	24.7	48.5	34.0	48.4	44.6	59.0
100	19.3	37.1	15.8	28.8	32.2	40.4	32.2	41.5
120	13.5	26.6	8.4	15.7	27.0	31.3	26.2	34.5

Comparison of iterative angular spectrum and optimal rotation angle methods in designing beam-fanners

Abstract

1. Introduction

2. Iterative angular spectrum method

3. Optimal rotation angle method

4. Transfer functions

5. Design and analysis

6. Conclusions

Acknowledgments

References and links

Cited By

Figures (5)

Tables (1)

Equations (10)

Optics Express