## Abstract

The arrangement of binary subwavelength structures is a promising alternative to the conventional multiheight level technique to generate computer generated holograms (CGHs). However, the current heuristic design approach leads to a slight mismatch between the target signal and experimental data. To evaluate this deviation, a diffractive beam splitter design is investigated rigorously using a finite-difference time-domain (FDTD) method. Since the use of a rigorous Maxwell-equation solver like FDTD requires a massive computational effort, an alternative scalar approach, a fast Fourier transform beam propagation method (FFT-BPM), is investigated with a substantial higher computing speed, showing still a good agreement with the FDTD simulation and experimental data. Therefore, an implementation of this scalar approach into the CGH design process offers the possibility to significantly increase the accuracy.

© 2013 OSA

## 1. Introduction

Computer generated holograms (CGH) are widely applied in many areas of modern optics as beam shaping, high-precision interferometric testing of optical surfaces, image encryption or pattern generation. The efficiency of phase-only computer generated holograms can be enhanced by increasing the number of phase levels. Especially diffractive optical elements (DOE) with a high numerical aperture (e.g. $\ge 0.3$) in the visible range require a surface profile with a high aspect ratio and small feature sizes. Because of alignment errors and a limited resist resolution, the fabrication process often allows only the realization of a small number of phase steps [1–3]. Just for this case, the recently developed binary effective medium (EM) approach is applicable, which enables a high-resolution patterning due to an advanced binary technology. An effective medium can be described as a binary grating structure with a period small compared to the wavelength of light leading to only zero-order diffraction. Such a medium behaves analogously to a homogeneous slab. The artificial effective refractive index of the medium can be varied between the index of the incident medium (supposedly 1) and the refractive material by changing the filling factor of the grating. So far, the use of the EM approach for DOE structures bases on a certain periodical repetition of the subwavelength structures (SWS) [4–9]. In [10], we could demonstrate that such a periodical repetition is not stringently necessary. The phase delay at a certain point in the CGH phase pattern is realized by a single non-repeated SWS. Therefore, an arbitrary phase distribution can be generated by a lateral arrangement of various binary SWS. This approach enables the use of an efficient and fast binary technology chain which enables the generation of EM CGHs with high diffraction angles and a good efficiency due to a strong twin-image suppression [11–13].

Although the previously demonstrated EM CGHs show a good performance, there exists a discrepancy between the designed target signal and the measured intensity distribution [10,11]. Until now, this deviation could not be exactly quantified with rigorous simulation methods in consequence of the complexity of the former diffractive elements. In this paper, we investigate the impact of the used approximations and assumptions which are made within the heuristic EM CGH design process for a particular beam splitter design. The element and the corresponding design approach are presented in the second section of the paper. The simplicity of the element enables a rigorous handling, in our case the use of the finite-difference time-domain (FDTD) method. This allows a direct comparison of the near and far field distributions assumed in the EM CGH design with those given by the FDTD simulation as shown in the third section. Since the FDTD is very accurate but requires a massive computational effort, an alternative scalar approach, the fast Fourier transform beam propagation method (FFT-BPM), is introduced in the fourth section. Due to its fast Fourier transform (FFT) implementation, this method benefits from a low computational effort and a very short computing time. To evaluate appropriate sampling parameters for the FFT-BPM, the results are compared with those of the FDTD simulation as shown in the 5th section. In the 6th section, the used technology chain to fabricate the binary diffractive element is briefly introduced. A comparison of the EM design output with both simulation methods and experimental data is shown and discussed in the 7th section.

## 2. Heuristic design approach of an effective medium computer generated hologram

To ensure a rigorous handling, we choose a simple diffractive beam splitter design with a low resolution of 20x20 pixels and three phase level configuration. The pixel size of 380 nm therefore leads to a 7.6 µm x 7.6 µm CGH unit cell. The beam splitter creates nine irregularly arranged spots with a maximum diffraction angle of 26.6°. Apart from constant pre-factors, the far field intensity distribution of an optical element is equal to the square of the Fourier transform (FT) of the scalar complex field right behind the element. A phase-only element does only tailors the incident wave front by introducing a local phase delay defined by the designed phase function. As an approximation, during the design procedure absorption effects are neglected, which means that the amplitude remains constant at one over the whole range. To design the phase function to a given target signal [Fig. 1(a) ], a standard iterative Fourier transform algorithm (IFTA) [14] is used. Commonly, the thin element approximation (TEA [15],) is applied to transfer the calculated phase function [Fig. 1(b)] into a physical structure, e.g. a surface profile. When the use of TEA is justified, the results can be obtained immediately with almost no computational effort. In this case, the far field intensity distribution of the real diffractive element only shows subtle variances to the IFTA target signal. The applicability of TEA is in fact limited to an element feature size, which should not be too small compared to the wavelength and a thickness which should be as low that diffraction within the structure is neglectable. A typical value for the feature size to wavelength ratio is greater than about 10 to ensure a valid approximation [16–18]. In [19] could be demonstrated that even structures with feature sizes in the range of the wavelength can still be described by TEA. In our case, TEA is assumed to transfer the resulting IFTA phase function into an effective refractive index distribution with effective indices between one and the index of the substrate material $n$. Considering TEA, this distribution has a constant element depth

where $N$ is the number of phase steps and $\lambda $ the used wavelength.For our setup with an operation wavelength of 532 nm and fused silica as the substrate material ($n$ = 1.46), the element depth is 770 nm. Then, an EM approach is applied to transfer the effective index distribution into a SWS pattern. To calculate the lateral extension of the SWS for different effective indices, a calibration curve is used which relates the effective index to the filling factor of a one or two dimensional subwavelength grating. As a first approximation, the effective index is almost directly proportional to the fraction of filled space in the unit cell. The exact shape of the calibration curve for different SWS geometries (e.g. square or circular hole or pillar structures) can be calculated numerically by a rigorous coupled wave analysis (RCWA) algorithm. Using these assumptions, the designed phase function obtained from the IFTA can be transferred into the corresponding arrangement of binary SWS [Fig. 1(c)]. Note that with this transfer into such a 2D arrangement, the lateral periodic boundary condition assumed for the calculation of the calibration curve is no longer satisfied. A detailed description of the heuristic design approach for EM CGHs can be found in [10]. For our design, we utilize a pillar geometry which leads to a subwavelength pillar distribution which is challenging but compatible with our fabrication process with a minimal pillar size of 256 nm and a minimal spacing between two adjacent structures of 62 nm.

## 3. Rigorous simulation of an effective medium computer generated hologram

To perform a rigorous simulation of the structure, we use the commercial RSoft Inc. FDTD solver named Fullwave. We choose periodic boundary conditions in $x$ and $y$ direction and an absorbing perfectly matched layer (PML) [20] that terminates the computational grid in $z$ direction. The incident field is given by a plane wave propagating in the positive $z$ direction with a constant amplitude of one. Both polarization states, TE (electric field vector in $x$ direction) and TM (electric field vector in $y$ direction) are investigated. The near field pattern is monitored 30 nm behind the structure (Fig. 2 ). In the following, the near and far field distributions of the FDTD simulations are compared to the IFTA target signal considering an ideal TEA element. Figure 3(a) shows the complex scalar near field distribution after propagation through an ideal TEA element. The amplitude of the incident plane wave remains unchanged at one and the phase is shifted exactly in the way given by the IFTA design. The FDTD simulation of the EM CGH pattern [Fig. 3(b)] reveals a markedly different behavior. The amplitude distribution shows clear waveguiding effects characterized by local field enhancements which are about three times larger than the initial amplitude of the input field. They are caused by the higher index of the pillars surrounded by a low index material (in this case air). The FDTD near field phase distribution has slight resemblance to the target phase function. The sharp edges and small structures of the surface profile are not visible as a phase retardation of the same shape in the optical field behind the element, what is to be expected because those are much smaller than the wavelength. The high spatial frequency components necessary to resolve the sharply outlined structural details are evanescent and so attenuated during propagation through the element. These evanescent components still carry a significant amount of energy which cannot be considered correctly while using TEA.

The FDTD far field intensity distribution in k-space

The cross correlation between TE and TM is ${C}_{TE-TM}\text{}$ = 0.999 and the summed intensity deviation over all 10 spots is only 0.01${I}_{0}$, whereas ${C}_{TEA-TE/TM}$ is 0.835 with an overall intensity deviation of 0.273${I}_{0}$. These results demonstrate the lack of accuracy considering EM CGHs as thin elements.

## 4. Fast Fourier transform beam propagation method

The FFT-BPM or split-step Fourier method was first introduced by M. D. Feit and J. A. Fleck [21]. Due to the use of a fast Fourier transform algorithm, it is a very powerful numerical technique to solve the scalar Helmholtz equation. While the method is applied to model various tasks of propagation problems, it is a standard tool to calculate the mode properties of optical fibers [22–24]. Although the validly region of this method is limited to paraxial propagation $({k}_{z}^{2}\gg {k}_{x}^{2}+{k}_{y}^{2})$ [25], we demonstrate that the FFT-BPM provides a good approximation for light wave propagation through an EM CGH with subwavelength feature sizes and thus also high spatial frequency components. This scalar approach offers a very high computing speed compared with more rigorous methods like FDTD or RCWA. A requirement for the applicability of the method is the use of polarization independent structures, which is quite common for many kinds of CGHs. To perform a FFT-BPM, the structure is divided into a certain number of sub-steps along the propagation direction $z$. The spatial refractive index distribution $n(x,y)=\overline{n}\pm \Delta n(x,y)$ is described in terms of an averaged refractive index $\overline{n}$ and a local index change $\Delta n(x,y)$. The propagation operator itself is spitted into a homogeneous medium propagation with $\overline{n}$ via a standard angular spectrum of plane waves (ASPW) approach [15] and a thin element index shift of $\pm (2\pi \cdot \Delta z\cdot (\overline{n}-1))/\lambda $ to centrally implement the phase delay between the higher and the lower index (Fig. 5 ). The propagation of a scalar field $U(x,y,z)$ about $\Delta z$ can thus be calculated with the following expression:

## 5. Fast Fourier transform beam propagation method versus finite-difference time-domain method

To verify the applicability of the FFT-BPM for EM CGHs, we compare the computed near and far field patterns with results of the FDTD simulation. The FFT-BPM calculates a spatial near field $U(x,y,z)$ at the output plane which is located 30 nm behind the structure. To calculate the far field intensity distribution

The intensity deviation to the FDTD summed over all 10 main spots is 0.039${I}_{0}$, whereas the major part of 0.016${I}_{0}$ is in the zeroth order. Compared to the TEA approach with ${C}_{TEA-FDTD}=$ 0.835 and an overall intensity deviation to the FDTD simulation of 0.273${I}_{0}$, this is a significant improvement of accuracy. In order to evaluate the increase in computing speed, the following example is given as a rough estimation: While for the presented beam splitter configuration the FFT-BPM only takes 2.8 seconds on a state of the art desktop computer (single x86 core @ 3.8 GHz), the FDTD needs about 30 minutes on a SMP workstation of the same generation (32 x86 cores @ 2 GHz). Note that numerical conditions are equal for both calculations. The step width is 10 nm in the direction of propagation and the spatial resolution is 1024x1024 sampling points. Regarding the computational effort, the FFT-BPM does not only save a factor of about 10^{4} floating point operations in this particular example, it also scales better with the number of sampling points concerning memory and computation speed. This enables the use of this method for more complicates problems where an FDTD simulation is not possible because of memory issues. Further an optimization of the designed phase function might be possible with this approach in the future.

## 6. Element fabrication and characterization

The fabrication effort for binary EM CGHs is much lower than for multiheight level elements fabricated with multistep lithography which requires several exposure and etching steps and thus suffers from alignment errors [1,2]. Our binary process has only two main steps, one lithographic and one etching step. The lithographic step is performed by a Vistec SB350OS electron-beam writer which enables a very fast writing speed with high precision due to its variable shaped-beam writing strategy. Therefore, a common positive-tone chemically amplified electron beam resist (FEP 171, Fuji Film) is used. This resist enables a writing speed of approximately one square millimeter per 25 seconds for the given sample. The developed 250 nm thick resist pattern is then transferred into a 50 nm thick chromium layer to enhance the selectivity of the following reactive ion beam etching process into the fused silica substrate. The whole technology chain is capable to produce the required high aspect ratio with a depth to groove width ratio of 12.5 [Fig. 9(a)
]. For a good experimental handling, the 7.6 µm x 7.6 µm unit cell is replicated across a sample area of 5 mm x 5 mm. The experiment is performed at normal incidence with a collimated TE polarized 532 nm diode laser [Fig. 9(b)]. The transmitted intensity of the main spots is detected using an integrating sphere. The measured normalized intensity of the $m$^{th} spot ${I}_{m}={P}_{m}/{P}_{0}$ is calculated from the measured optical power ${P}_{m}$ of the $m$^{th} spot and the overall transmitted optical power ${P}_{0}$.

## 7. Comparison of different simulation methods and experiment

In Fig. 10 a direct comparison of the main spots for the different simulation results and the experimental data is illustrated. Table 1 additionally shows the total intensity deviation in relation to the FDTD simulation. The measurement verifies the FDTD simulation very well. The intensity deviation summed over all 10 main spots is only about 0.044${I}_{0}$. A reason for this mismatch could be the slight differences between theoretical and fabricated structural parameters. As visible in Fig. 9(a), the technological process leads to a slight rounding of the edges and a notable line edge roughness. Further, shadowing effects in the thin grooves partly lead to an insufficient etching depth in this region. Comparing the different simulation methods, it is obvious that the simple heuristic TEA design based on an effective medium approach leads to a substantial error whose magnitude always depends on the particular design. Already the near fields in Fig. 3 show that the pillars do not really behave like a homogeneous effective medium pixel. The used RCWA calibration curve bases on periodic boundary conditions which are not valid for a pillar embedded in an EM CGH superstructure. The internal diffraction occurring within the element cannot be considered with TEA. Thus, the heuristic design approach only leads to a weak approximation of the far field intensity distribution for our beam splitter design. Especially, the massive loss in the zeroth order cannot be described.

The FDTD is a powerful and accurate simulation tool, but it requires a massive computational effort. This barely allows the computation of even larger EM CGH structures as the here demonstrated example. The scalar FFT-BPM, however, allows a fast computation with an acceptable error deviation for polarization independent elements. The calculated FFT-BPM near field distribution of the presented beam splitter show a good accordance to the FDTD simulation and also the far field pattern is in good agreement with the FDTD and experimental data. The total intensity deviation between measurement and FFT-BPM summed over the main 10 spots is only about 0.059${I}_{0}$, and is thus only slightly larger than the deviation between measurement and FDTD. Consequently, since the error induced by the fabrication is in the same order of magnitude, the FFT-BPM provides a sufficient solution to evaluate the performance of EM CGHs.

## 8. Conclusion

As demonstrated in this paper, the current heuristic TEA design based on an effective medium approach only leads to a weak approximation of the generated far field distribution of the presented beam splitter design. The results of the FFT-BPM, however, show a good agreement to the FDTD simulation and experimental data. Therefore, while the integration into CGH design algorithms and optimization with a rigorous Maxwell equation solver via FDTD or RCWA is not possible due to the high computational effort, it becomes within reach with the FFT-BPM. Therefore, the common IFTA-based design approach has to be optimized by an appropriate FFT-BPM step, which will be the next step within the presented framework.

## Acknowledgments

The authors are grateful to W. Gräf, C. Fuchs, H. Schmidt and T. Käsebier (Institute of Applied Physics, Friedrich-Schiller-University Jena) who enable the challenging fabrication process. Special thanks to H.-C. Eckstein (Fraunhofer Institute for Applied Optics and Precision Engineering, IOF, Jena) for the helpful discussions and for providing his software SToNE Lab where I could implement my own algorithms. The authors thank D. Lehr and S. Kroker (Institute of Applied Physics, Friedrich-Schiller-University Jena) for their careful proofreading of the manuscript. Further, the financial support of the German Federal Ministry of Education and Research (BMBF, grant number 03IS2101A) and the Thuringian Ministry for Education, Science and Culture (TMBWK, grant number PE116-1) is gratefully acknowledged.

## References and links

**1. **B. Goebel, L. L. Wang, and T. Tschudi, “Multilayer technology for diffractive optical elements,” Appl. Opt. **35**(22), 4490–4493 (1996). [CrossRef]

**2. **J. M. Miller, M. R. Taghizadeh, J. Turunen, and N. Ross, “Multilevel-grating array generators: fabrication error analysis and experiments,” Appl. Opt. **32**(14), 2519–2525 (1993). [CrossRef]

**3. **M. Banasch, L.-C. Wittig, and E.-B. Kley, “Fabrication tolerances of binary and multilevel Computer Generated Holograms (CGHs) with submicron Pixel Size,” *MOC´04–10th Microoptics Conference, Germany* (2004).

**4. **E. Noponen and J. Turunen, “Binary high-frequency-carrier diffractive optical elements: electromagnetic theory,” J. Opt. Soc. Am. A **11**(3), 1097–1109 (1994). [CrossRef]

**5. **J. Mait, D. Prather, and M. Mirotznik, “Design of binary subwavelength diffractive lenses by use of zeroth-order effective-medium theory,” J. Opt. Soc. Am. A **16**(5), 1157–1167 (1999). [CrossRef]

**6. **P. Lalanne, S. Astilean, P. Chavel, E. Cambril, and H. Launois, “Design and fabrication of blazed binary diffractive elements with sampling periods smaller than the structural cutoff,” J. Opt. Soc. Am. A **16**(5), 1143–1156 (1999). [CrossRef]

**7. **C. Ribot, P. Lalanne, M. S. Lee, B. Loiseaux, and J. P. Huignard, “Analysis of blazed diffractive optical elements formed with artificial dielectrics,” J. Opt. Soc. Am. A **24**(12), 3819–3826 (2007). [CrossRef]

**8. **H. J. Hyvärinen, P. Karvinen, and J. Turunen, “Polarization insensitive resonance-domain blazed binary gratings,” Opt. Express **18**(13), 13444–13450 (2010). [CrossRef]

**9. **W. Yu, K. Takahara, T. Konishi, T. Yotsuya, and Y. Ichioka, “Fabrication of multilevel phase computer-generated hologram elements based on effective medium theory,” Appl. Opt. **39**(20), 3531–3536 (2000). [CrossRef]

**10. **W. Freese, T. Kämpfe, E.-B. Kley, and A. Tünnermann, “Design of binary subwavelength multi-phase level computer generated holograms,” Opt. Lett. **35**(5), 676–678 (2010). [CrossRef]

**11. **W. Freese, T. Kämpfe, W. Rockstroh, E.-B. Kley, and A. Tünnermann, “Optimized electron beam writing strategy for fabricating computer-generated holograms based on an effective medium approach,” Opt. Express **19**(9), 8684–8692 (2011). [CrossRef]

**12. **W. Freese, T. Kämpfe, E.-B. Kley, and A. Tünnermann, “Multi-phase-level diffractive elements realized by binary effective medium patterns,” Proc. SPIE **7591**(75910Z), 75910Z-1–75910Z-7 (2010). [CrossRef]

**13. **W. Freese, T. Kämpfe, E.-B. Kley, and A. Tünnermann, “Design and fabrication of a highly off-axis binary multi-phase level computer-generated hologram based on an effective medium approach,” Proc. SPIE **7927**(792710), 792710, 792710-7 (2011). [CrossRef]

**14. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttg.) **35**, 227–246 (1972).

**15. **J. W. Goodman, *Introduction to Fourier Optics, 3rd ed.* (Roberts & Company, 2005).

**16. **U. Levy, E. Marom, and D. Mendlovic, “Thin element approximation for the analysis of blazed gratings: simplified model and validity limits,” Opt. Commun. **229**(1-6), 11–21 (2004). [CrossRef]

**17. **A. Vasara, M. R. Taghizadeh, J. Turunen, J. Westerholm, E. Noponen, H. Ichikawa, J. M. Miller, T. Jaakkola, and S. Kuisma, “Binary surface-relief gratings for array illumination in digital optics,” Appl. Opt. **31**(17), 3320–3336 (1992). [CrossRef]

**18. **D. Pommet, M. Moharam, and E. Grann, “Limits of scalar diffraction theory for diffractive phase elements,” J. Opt. Soc. Am. A **11**(6), 1827–1834 (1994). [CrossRef]

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

**20. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method* (Boston: Artech House, 2000).

**21. **M. D. Feit and J. A. Fleck Jr., “Light propagation in graded-index optical fibers,” Appl. Opt. **17**(24), 3990–3998 (1978). [CrossRef]

**22. **S. Kumar, T. Srinivas, and A. Selvarjan, “Beam propagation method and its application to integrated optic structures and optical fibers,” J. Phys. **34**, 347–358 (1989).

**23. **M. D. Feit and J. A. Fleck Jr., “Calculation of dispersion in graded-index multimode fibers by a propagating-beam method,” Appl. Opt. **18**(16), 2843–2851 (1979). [CrossRef]

**24. **B. Hermansson, D. Yevick, and J. Saijonmaa, “Propagating-beam-method analysis of two-dimensional microlenses and three-dimensional taper structures,” J. Opt. Soc. Am. A **1**(6), 663–671 (1984). [CrossRef]

**25. **M. Fertig and K.-H. Brenner, “Vector wave propagation method,” J. Opt. Soc. Am. A **27**(4), 709–717 (2010). [CrossRef]