An algorithm is reported for the design of a phase-only diffractive optical element (DOE) that reshapes a beam focused using a high numerical aperture (NA) lens. The vector diffraction integrals are used to relate the field distributions in the DOE plane and focal plane. The integrals are evaluated using the chirp-z transform and computed iteratively within the Method of Generalized Projections (MGP) to identify a solution that simultaneously satisfies the beam shaping and DOE constraints. The algorithm is applied to design a DOE that transforms a circularly apodized flat-top beam of wavelength λ to a square irradiance pattern when focused using a 1.4-NA objective. A DOE profile is identified that generates a 50λ×50λ square irradiance pattern having 7% uniformity error and 74.5% diffraction efficiency (fraction of focused power). The diffraction efficiency and uniformity decrease as the size of the focused profile is reduced toward the diffraction limited spot size. These observations can be understood as a manifestation of the uncertainty principle.
© 2008 Optical Society of America
Beam shaping diffractive optical elements (DOEs) are becoming an essential part of many optical systems. A DOE can be used to modify the amplitude, phase, and polarization of an incident beam so that it focuses into a targeted irradiance distribution at the image plane [1–3]. Beam shaping can enhance performance in optical lithography, laser-based materials processing, direct laser writing, surgical applications, and optical data storage . For many applications the optimum irradiance distribution consists of a flat-top profile having a defined geometry within the focal plane. Such irradiance patterns can be characterized in terms of the diffraction efficiency, κ, and uniformity error, δ. The diffraction efficiency quantifies the fraction of total optical power directed into the targeted region of interest and the uniformity error provides a measure of flatness in the irradiance distribution across that region.
Many excellent scalar techniques have been reported for designing beam shaping DOEs. These approaches are based on methods that include geometric mapping [5, 6], analytical solution , iterative processes [8–11], and genetic optimization . Although exceptional results have been achieved with these algorithms, they are all based on scalar diffraction theory and as such are only valid in the paraxial domain of diffractive optics . For systems with high numerical aperture (NA), depolarization effects are significant , so vectorial diffraction theory must be used in the DOE design process. This becomes particularly challenging because the overall beam shape is determined by the summed irradiance of the x-, y-, and z-polarized electric fields. Although the field components are orthogonal, they are not entirely independent because each is reshaped by a common DOE. As a result, the DOE must be designed to reshape the field components collectively so the irradiance of the total field matches the targeted beam shape. One report of vectorial beam shaping has appeared, but the method was not applied to high-NA systems . Given that high-NA systems are being increasingly employed in frontier technologies, further applications of beam shaping will be stymied unless accurate methods for vectorial beam shaping are developed.
In this work we report a vectorial beam shaping algorithm that can be used to design phase-only DOEs for use under high-NA conditions. The algorithm was developed by incorporating the vector diffraction integrals  into the Method of Generalized Projections (MGP) . The diffraction integrals are used to interrelate the DOE phase profile and the resulting vectorial electric field in the focal plane. The integrals are evaluated using the chirp-z transform  to improve computational speed and accuracy. Iterative projection of constraints in the pupil and focal planes progressively forces the simulation toward a DOE phase profile that generates the targeted beam shape. The new algorithm is applied to the problem of designing a phase-only DOE that transforms a circularly apodized flat-top input beam into a square flat-top irradiance distribution when focused using a 1.4-NA objective. In beam shaping, high diffraction efficiency and low uniformity error are known to be mutually exclusive characteristics that must be considered jointly in optimizing DOEs . In this work, we also investigate how κ and δ change as the size of the focused beam profile approaches the diffraction-limited spot size.
2.1. Vector diffraction integrals
The optical geometry of the focusing system is depicted in Fig. 1. A DOE and an aberration free aplanatic lens which fulfils the sine condition and has focal length f are positioned such that their optical axes are collinear with the z-axis of a Cartesian coordinate system whose origin is located at the Gaussian focus of the lens. The numerical aperture of the lens is NA=1.4 in all calculations. Monochromatic linearly polarized plane waves, with electric field vector parallel to the x-axis, propagate along the z-axis, passing through the DOE and entering the pupil of the lens. The light focuses into a medium of refractive index n=1.516. In the absence of the DOE, this situation is consistent with common applications of high-NA oilimmersion objective lenses.
At any point within the optical system the irradiance is I=(1/2)ncε0|E|2. The reshaped beam is the spatial map of the focused irradiance If(xf, yf) for all P. The speed of light and electric permittivity in vacuum are c and ε0, respectively. The wave number of light transmitted through the lens is kt=2π/λ=[k 2 x+k 2 y+k 2 z]1/2, with kx, ky, and kz being the plane wave components, and λ is the wavelength within the medium. The NA of the lens system sets k max=ktNA/n. The function T(kx, ky)exp[iΦ(kx, ky)] describes the transmission amplitude (T) and phase (Φ) of the DOE. The amplitude of the incident electric field Ein is assumed to be spatially constant, so this term was brought outside the integral. The focal plane is divided into a region of interest Ω and its complement Ωc. The region Ω wholly contains and bounds the targeted beam shape It.
It is helpful to cast the vector diffraction integrals into a form consisting of dimensionless variables by normalizing to the lateral extent of the input beam Iin and the targeted beam profile It . The Cartesian coordinates of the aperture plane (xa, ya) can be related to the x- and y-components of the wave vector by
where R is the radius of the entrance pupil, and u and v represent the normalized kx and ky components of the wave vector, respectively. Likewise, the focal plane coordinates (xf, yf) are normalized by D=mλ (see Fig. 1), giving the transformed focal plane coordinates (ξ, η):
The size of It can be scaled conveniently by m multiples of λ. All free parameters of the system can be combined into the single variable, β, given by :
Where q(u, ν)=[1-(NA/n)2(u 2+ν 2)]1/2, and ωx=k max Dξ=βξ and ωy=k max Dη=βη have been introduced so the diffraction integrals in Eq. (5) can be evaluated as a Fourier transform.
2.3. Chirp-z transform
In this work, the double integrals of Eq. (5) were evaluated using a two-dimensional (2D) chirp-z transform (CZT) with 512×512 sampling points in both the aperture and the focal planes, irrespective of the magnitude of R and D. Although 2D fast Fourier transform can be used, CZT is computationally faster and better suited for the present situation because it internalizes zero-padding and allows the spacing of sampling points in the aperture and focal planes to be set independently . This greatly reduces the number of sampling points required when R and D differ substantially in magnitude, as in the present case.
Evaluating Eq. (5) effectively propagates the field Ex(u, v) forward into Ex(ξ, η), Ey(ξ, η), and Ez(ξ, η). The fundamental operation of beam shaping involves applying constraints associated with It to the field in the focal plane and then calculating a new DOE phase profile that comes closest to generating the reshaped field. A new DOE transmission profile is obtained by “backward propagating” to the pupil plane through an inverse-CZT applied to the reshaped field.
Because the forward and inverse CZTs are applied to a finite region of the DOE and focal planes, Gibbs artifacts are generated . If these numerical errors are allowed to accumulate, they can degrade the uniformity of If or even cause the algorithm to diverge from a solution. Gibbs artifacts were suppressed by applying a Kaiser window to the amplitude of the focal field profiles immediately after they were computed using CZT .
2.4. Method of generalized projections
For a given input beam Iin we seek a DOE phase function Φ that generates a focal plane field distribution such that the sum of the x-, y-, and z-polarized irradiance Ix+Iy+Iz=If, matches the targeted irradiance It for all (ξ, η). An exact match is generally not possible because it is not known a priori that an arbitrary It is a solution to the wave equation . This is the case for the present example because the targeted square irradiance profile requires a discontinuous drop in the field at the interface of Ω and Ωc. The problem is further complicated because Φ(u, ν) affects each of Ex, Ey, and Ez, so the field components are not truly independent. Evaluating Eq. (5) and integrating Ix, Iy, and Iz over the entire focal plane shows that their fractional power content is 0.74, 0.01 and 0.25, respectively, and these values are independent of Φ(u, ν). So, high-NA beam shaping demands that Ix, Iy, and Iz are optimized collectively. This problem cannot be solved analytically, so iterative numerical techniques must be employed. The MGP is particularly well suited to the current problem because it can find solutions that closely satisfy sets of inconsistent and non-physical constraints . In the MGP the optical field is repeatedly propagated forward and backward between the DOE and focal planes, and constraints associated with both domains are applied in each iteration until a satisfactory solution is found.
2.5. Starting conditions
The goal is to design a phase-only DOE, so the initial transmission amplitude T0 is set to unity. The rate of convergence and quality of the solution can be greatly improved by initiating the vector diffraction algorithm using a well chosen starting DOE phase profile, Φ 0(u, ν). Geometrical optics based methods can be used to identify suitable starting DOEs . Geometrical transformations have been applied successfully to obtain Φ 0(u, ν) analytically when Iin and It are either separable or axially symmetric ; however, the beam shaping example studied here is neither separable nor axially symmetric. To overcome this problem, a procedure described by Aagedal et al.  was employed to obtain Φ 0(u, ν) as a combination of two separate DOEs that together achieve the required geometric beam transformation. The first element converts the axially symmetric Iin into a standard Gaussian beam, which is separable and axially symmetric. The second element converts the Gaussian beam into a square-shaped super-Gaussian beam. The resulting Φ 0(u, ν) does not adequately reshape Iin to It under high-NA, but it provides a good starting point for the vector diffraction algorithm.
2.6. Algorithm flow
The analytically calculated starting DOE, Φ0(u, ν), is substituted into Eq. (5). The diffraction integrals are then evaluated using the CZT giving |Ex|exp(iϕx), |Ey|exp(iϕy), and |Ez|exp(iϕz) in the focal plane, and the corresponding irradiance distribution is If. The diffraction efficiency and uniformity error are calculated for the solution as
Under these definitions, a “perfect” solution would have κ=1 and δ=0.
Second, the constraint of the targeted beam shape It is applied. For the present example, It is defined as
with the limits of Ω set to -0.25≤ξ,η≤0.25. Within the beam shaping area, total power is conserved and the irradiance is homogenized by setting the latter to its average value:
The x-polarized field within Ω is reshaped as
whereas the following are left unchanged: |Ex(ξ, η)| outside Ω; |Ey(ξ, η)| and |Ez(ξ, η)| across all Ω+Ωc; and the phases ϕ(ξ, η) of all field components in Ω+Ωc. γ is an adjustable scalar that augments |Ex(ξ, η)| in Ω relative to that in Ωc. This operation provides a means for slowly pulling energy from Ωc into Ω . In this work γ=1.03 was used for all iterations.
Third, an inverse CZT is applied to the reshaped Ex(ξ, η) to generate a complex DOE function T i+1(u, ν)exp[iΦ i+1(u, ν)]. Given that a phase-only DOE is required, we apply this constraint by resetting the transmission amplitude to T0 while retaining the phase. The new complex DOE transmission function becomes T0(u, ν)exp[iΦ i+1(u, ν)]. The electric field is then propagated forward again using the new DOE and the reshaped beam it generates is evaluated based on κ and δ. This process continues until the algorithm converges to a suitable solution or until a fixed number of iterations are completed.
3. Results and discussion
Equation (10) is intentionally configured so that the beam shaping constraint is only directly applied to the x-component of the field amplitude lying within Ω. This arrangement provides amplitude freedom outside the region of interest and phase freedom across the entire focal plane that help the algorithm converge to a solution [11, 24]. Additionally, it solves the problem of how to reshape three independent field components that are effectively coupled through a common DOE. The intended beam shape is applied repeatedly to the x-polarized field, as it contains the majority of the focused power, and only it is propagated backward to obtain the DOE phase function for the next iteration. The y- and z-polarized components of the focal field are reshaped indirectly when they are calculated by forward propagation through the new DOE. Repeated iterations effectively pull Ex(ξ, η), Ey(ξ, η), and Ez(ξ, η) toward distributions that collectively satisfy Ix+Iy+Iz→It.
We now discuss the results obtained when the vector diffraction algorithm was used to design a phase-only DOE that transforms a circularly apodized flat-top input beam into a focused square irradiance pattern of area 50λ×50λ (D=100λ). Figure 2 shows the normalized focal irradiance distributions of the three polarization components, Ix, Iy, and Iz, and the total focal irradiance distribution If generated by the DOE phase profile of Fig. 3. This DOE and the associated irradiance distributions were obtained after 600 iterations. The overall beam shape is square as intended with δ=7% and κ=74.5%, indicating that it has good uniformity and power confinement within the region of interest.
In contrast, the irradiance distributions of the constituent polarizations are non-uniform. I x most resembles the targeted profile, but appears doubly concave, as though squeezed along the x-axis. Although Ix is non-zero across the coordinate axes, Iy and Iz have node(s) at these positions where their field amplitudes drop to zero. Iy is most complex, appearing approximately four-fold symmetric with power concentrated in the corners of Ω. Iz exhibits two-fold symmetry with a single nodal plane lying along the y-axis. The regions of high irradiance in Iy and Iz fill in around the edges of the x-polarized profile making the total irradiance distribution If uniform and square. These profiles show that the vector diffraction algorithm successfully generates a DOE for which all polarization components of the field are reshaped concurrently to achieve a targeted irradiance distribution under high-NA focusing.
Figure 4 shows how the uniformity error and diffraction efficiency change during the calculation. The diffraction efficiency progressively increases because the parameter γ=1.03 causes power to transfer from Ωc into Ω with each iteration. On the other hand, the uniformity error drops rapidly and reaches an apparent plateau after circa 500 iterations. It is known that high uniformity and high diffraction efficiency are mutually exclusive characteristics in beam shaping . As a result, attempting to improve the diffraction efficiency beyond the level of 74% achieved at approximately 600 iterations caused the uniformity to erode. Obtaining solutions that are optimized in terms of both δ and κ could be achieved by extending the present vector diffraction algorithm through Tikhonov regularization theory . It is noteworthy that the diffraction efficiency and uniformity are very poor for the first iteration. This results because the starting DOE was designed using a geometrical transformation method, which does not account for the vector character of the field. It underscores the importance then of using vector diffraction theory to achieve accurate beam shaping under high-NA conditions, and it demonstrates the improvement that can be achieved in beam shaping using the present vectorial approach.
The problem of shaping a beam whose size approaches the diffraction limit was examined by repeating the calculations described above for D=50λ, 25λ, 10λ, and 5λ. The resulting focused irradiance distributions and the corresponding DOE phase profiles are shown in Fig. 5. Comparing the irradiance distributions reveals that the intended beam transformation can be achieved, even for targeted beam profiles having an edge-length of D/2=2.5λ (see Fig. 5(G)). The DOEs themselves have approximate four-fold symmetry with respect to rotation about the z-axis, as expected for a square target beam shape (consider also Fig. 3). The DOEs are also comprised of many concentric rings of steadily increasing phase reminiscent of a Fresnel lens. These concentric phase rings effectively negate some of the focusing power of the high-NA lens, so understandably their number and radial density decreases as the target beam size is reduced toward the diffraction limit.
As D decreases, the reshaped beam degrades in uniformity and sharpness at the boundary of Ω. The sharpness of the irradiance profiles was characterized empirically by fitting the normalized individual distributions to a super-Gaussian of the form
The order of the super-Gaussian, N, provides a measure of the sharpness of the beam profile at the boundary of the region of interest. Rather than applying Eq. (11) to all of Ω+Ωc, the fitting was restricted to a region within which If exceeds 0.5. This procedure yields fits that more accurately describe the steepness of the profiles at the interface between Ω and Ωc because it does not include irradiance fluctuations in Ωc that are necessarily part of any real solution to the wave equation. The beam shaping parameters κ, δ, and N obtained for each value of D are collated in Table 1. The data show that the beam uniformity, diffraction efficiency, and edge sharpness all deteriorate as the size of the reshaped beam is reduced toward the diffraction limit.
The results in Fig. 5 suggest that it becomes increasingly difficult to achieve a targeted irradiance distribution when the beam size becomes comparable to the diffraction limit, as has also been observed for scalar beam shaping . This phenomenon can be understood as a manifestation of the uncertainty principle . The irradiance distribution It that can be achieved is inherently limited by the finite spatial bandwidth of the input beam Iin and the limited range of wave vectors over which focusing occurs, as quantified by the NA. The limit in the achievable beam shape can be expressed as 
The particular usefulness of β now becomes apparent. Given that β is determined by all free parameters of the system (R, D, λ, and NA), it provides a single measure of the difficulty of the beam shaping problem. Values of β2ΔIinΔIt for the beam profiles of Fig. 5 are also included in Table 1. Because β2ΔIinΔIt depends quadratically on D, it drops rapidly within this series and is most comparable to unity at D=5λ, for which the reshaped beam quality is poorest. These data and Eq. (12) imply then that If will differ increasingly from It as the targeted beam size is decreased toward the diffraction limited spot size, with all other parameters kept fixed.
A vector diffraction algorithm was developed for designing phase-only DOEs that reshape beams focused under high-NA conditions. The algorithm accounts for depolarization effects that occur under high-NA focusing by relating the DOE complex transmittance function and the electric field in the focal plane using the vector diffraction integrals. The algorithm was applied in the design of a phase-only DOE that reshapes the focused irradiance distribution of a circularly apodized flat-top input beam into a uniform square profile when focused using a 1.4-NA objective lens. We observe that beam uniformity, diffraction efficiency, and edge sharpness all degrade as the size of the targeted flat-top beam is decreased. This suggests that beam shaping becomes increasingly difficult as the area of the targeted irradiance distribution approaches that of the diffraction limited spot size. There are many possibilities for extending the method reported here. A wider range of focused beam shapes could be considered by appropriately modifying the constraints used to define the targeted irradiance distribution. The search for solutions could be made more general by including other free parameters. For example, one could allow for freedom in the size of the homogenized area, so that the targeted size is optimized along with diffraction efficiency and uniformity. Three-dimensional beam shaping could be achieved by applying constraints in multiple planes. This work may also be useful for extending methods employed in the design of phase masks for high-resolution photo-lithography.
We wish to thank Marcel Leutenegger for helpful discussions pertaining to this work. We also thank the reviewers for their many helpful and insightful comments.
References and links
1. C. Dorrer and J. D. Zuegel, “Design and analysis of binary beam shapers using error diffusion,” J. Opt. Soc. Am. B 24, 1268–1275 (2007). [CrossRef]
3. L. A. Romero and F. M. Dickey, “Lossless laser beam shaping,” J. Opt. Soc. Am. A 13, 751–760 (1995). [CrossRef]
4. F. M. Dickey, S. C. Holswade, and D. L. Shealy, Laser Beam Shaping Applications (Taylor & Francis Group, Boca Raton, 2006).
5. O. Bryngdhal, “Geometrical transformations in optics,” J. Opt. Soc. Am. 64, 1092–1099 (1974). [CrossRef]
6. D. L. Shealy and S. H. Chao, “Geometric optics-based design of laser beam shapers,” Opt. Eng. 42, 3123–3138 (2003). [CrossRef]
7. H. Aagedal, M. Schmid, S. Egner, J. Müller-Quade, T. Beth, and F. Wyrowski, “Analytical beam shaping with application to laser-diode arrays,” J. Opt. Soc. Am. A 14, 1549–1553 (1997). [CrossRef]
8. M. Johansson and J. Bengtsson, “Robust design method for highly efficient beam-shaping diffractive optical elements using iterative-Fourier-transform algorithm with soft operations,” J. Mod. Opt. 47, 1385–1398 (2000). [CrossRef]
9. V. V. Kotlyar, P. G. Seraphimovich, and V. A. Soifer, “An iterative algorithm for designing diffractive optical elements with regularization,” Opt. Laser Eng. 29, 261–268 (1998). [CrossRef]
10. J. S. Liu and M. R. Taghizadeh, “Iterative algorithm for the design of diffractive phase elements for laser beam shaping,” Opt. Lett. 27, 1463–1465 (2002). [CrossRef]
11. F. Wyrowski, “Diffractive optical elements: iterative calculations of quantized, blazed phase structures,” J. Opt. Soc. Am. A 7, 961–969 (1990). [CrossRef]
12. G. Zhou, X. Yuan, P. Dowd, Y. L. Lam, and Y. C. Chan, “Design of diffractive phase elements for beam shaping: hybrid approach,” J. Opt. Soc. Am. A 18, 791–800 (2001). [CrossRef]
14. T. G. Jabbour and S. M. Kuebler, “Vector diffraction analysis of high numerical aperture focused beams modified by two- and three-zone annular multi-phase plates,” Opt. Express 14, 1033–1043 (2006). [CrossRef] [PubMed]
16. E. Wolf, “Electromagnetic diffraction in optical systems I. An integral representation of the image field,” Proc. Roy. Soc. A 253, 349–357 (1959). [CrossRef]
17. A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984). [CrossRef]
19. H. Kim, B. Yang, and B. Lee, “Iterative Fourier transform algorithm with regularization for the optimal design of diffractive optical elements,” J. Opt. Soc. Am. A 21, 2353–2365 (2004). [CrossRef]
20. R. Kant, “An analytical solution of vector diffraction for focusing optical systems with Seidal aberrations I. Spherical aberration, curvature of field, and distortion,” J. Mod. Opt. 40, 2293–2310 (1993). [CrossRef]
21. F. M. Dickey and D. L. Shealy, Laser Beam Shaping (Marcel Dekker, Inc., New York, 2000). [CrossRef]
23. R. Piestun and J. Shamir, “Synthesis of three-dimensional light fields and applications,” Proceedings of the IEEE 90, 222–244 (2002). [CrossRef]
24. F. Wyrowski, “Diffraction efficiency of analog and quantized digital amplitude holograms: analysis and manipulation,” J. Opt. Soc. Am. A 7, 383–393 (1990). [CrossRef]
25. J. S. Liu, A. J. Caley, and M. R. Taghizadeh, “Symmetrical iterative Fourier-transform algorithm using both phase and amplitude freedoms,” Opt. Commun. 267, 347–355 (2006). [CrossRef]