1. Introduction

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 [4]. 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 [7], iterative processes [8–11], and genetic optimization [12]. 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 [13]. For systems with high numerical aperture (NA), depolarization effects are significant [14], 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 [15]. 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 [16] into the Method of Generalized Projections (MGP) [17]. 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 [18] 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 [19]. In this work, we also investigate how κ and δ change as the size of the focused beam profile approaches the diffraction-limited spot size.

2. Theory

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.

 

Fig. 1. Optical setup of the beam shaping problem. The aperture represents the input pupil of the objective lens. The focal plane is divided into two regions. Ω represents the region of interest that contains and bounds the targeted beam shape. Its complement Ω^c represents the remainder of the focal plane.

Download Full Size | PDF

The electric field at an arbitrary point P(x_f, y_f) in the focal plane (z_f=0) can be calculated using the vector diffraction integrals [16, 20] given as

At any point within the optical system the irradiance is I=(1/2)ncε₀|E|². The reshaped beam is the spatial map of the focused irradiance I_f(x_f, y_f) for all P. The speed of light and electric permittivity in vacuum are c and ε₀, respectively. The wave number of light transmitted through the lens is k_t=2π/λ=[k ² _x+k ² _y+k ² _z]^1/2, with k_x, k_y, and k_z being the plane wave components, and λ is the wavelength within the medium. The NA of the lens system sets k _max=k_tNA/n. The function T(k_x, k_y)exp[iΦ(k_x, k_y)] describes the transmission amplitude (T) and phase (Φ) of the DOE. The amplitude of the incident electric field E_in 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 I_t.

2.2. Normalization

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 I_in and the targeted beam profile I_t [21]. The Cartesian coordinates of the aperture plane (x_a, y_a) 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 k_x and k_y components of the wave vector, respectively. Likewise, the focal plane coordinates (x_f, y_f) are normalized by D=mλ (see Fig. 1), giving the transformed focal plane coordinates (ξ, η):

The size of I_t can be scaled conveniently by m multiples of λ. All free parameters of the system can be combined into the single variable, β, given by [21]:

Substituting Eqs. (3) and (4) into Eq. (1) and algebraically manipulating gives the following normalized vectorial diffraction integrals for the field at the focal plane:

Where q(u, ν)=[1-(NA/n)²(u ²+ν ²)]^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 [22]. 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 E_x(u, v) forward into E_x(ξ, η), E_y(ξ, η), and E_z(ξ, η). The fundamental operation of beam shaping involves applying constraints associated with I_t 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 [18]. If these numerical errors are allowed to accumulate, they can degrade the uniformity of I_f 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 [18].

2.4. Method of generalized projections

For a given input beam I_in 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 I_x+I_y+I_z=I_f, matches the targeted irradiance I_t for all (ξ, η). An exact match is generally not possible because it is not known a priori that an arbitrary I_t is a solution to the wave equation [11]. 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 E_x, E_y, and E_z, so the field components are not truly independent. Evaluating Eq. (5) and integrating I_x, I_y, and I_z 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 I_x, I_y, and I_z 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 [23]. 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 T₀ 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, Φ ₀(u, ν). Geometrical optics based methods can be used to identify suitable starting DOEs [5]. Geometrical transformations have been applied successfully to obtain Φ ₀(u, ν) analytically when I_in and I_t are either separable or axially symmetric [7]; however, the beam shaping example studied here is neither separable nor axially symmetric. To overcome this problem, a procedure described by Aagedal et al. [7] was employed to obtain Φ ₀(u, ν) as a combination of two separate DOEs that together achieve the required geometric beam transformation. The first element converts the axially symmetric I_in 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 Φ ₀(u, ν) does not adequately reshape I_in to I_t under high-NA, but it provides a good starting point for the vector diffraction algorithm.

2.6. Algorithm flow

The analytically calculated starting DOE, Φ₀(u, ν), is substituted into Eq. (5). The diffraction integrals are then evaluated using the CZT giving |E_x|exp(iϕ_x), |E_y|exp(iϕ_y), and |E_z|exp(iϕ_z) in the focal plane, and the corresponding irradiance distribution is I_f. 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 I_t is applied. For the present example, I_t 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: |E_x(ξ, η)| outside Ω; |E_y(ξ, η)| and |E_z(ξ, η)| across all Ω+Ω^c; and the phases ϕ(ξ, η) of all field components in Ω+Ω^c. γ is an adjustable scalar that augments |E_x(ξ, η)| in Ω relative to that in Ω^c. This operation provides a means for slowly pulling energy from Ω^c into Ω [8]. In this work γ=1.03 was used for all iterations.

Third, an inverse CZT is applied to the reshaped E_x(ξ, η) 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 T₀ while retaining the phase. The new complex DOE transmission function becomes T₀(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 E_x(ξ, η), E_y(ξ, η), and E_z(ξ, η) toward distributions that collectively satisfy I_x+I_y+I_z→I_t.

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, I_x, I_y, and I_z, and the total focal irradiance distribution I_f 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.

 

Fig. 2. (A) Calculated irradiance distribution resulting when a circularly apodized flat-top input beam of radius R is passed through the phase-only DOE shown in Fig. 3 and focused using a 1.4-NA objective. The DOE was designed to reshape the beam into a flat-top square irradiance pattern of area 50λ×50λ (D=100λ). (B)–(D) Irradiances of the constituent x-, y-, and z-polarized components of the total field. Each profile is normalized to the peak of I_f.

Download Full Size | PDF

 

Fig. 3. DOE phase profile that generates the focal irradiance distributions shown in Fig. 2. The phase is plotted in units of radians.

Download Full Size | PDF

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 I_x is non-zero across the coordinate axes, I_y and I_z have node(s) at these positions where their field amplitudes drop to zero. I_y is most complex, appearing approximately four-fold symmetric with power concentrated in the corners of Ω. I_z exhibits two-fold symmetry with a single nodal plane lying along the y-axis. The regions of high irradiance in I_y and I_z fill in around the edges of the x-polarized profile making the total irradiance distribution I_f 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 [19]. 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 [19]. 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.

 

Fig. 4. (A) Evolution of the beam shaping diffraction efficiency and (B) uniformity error versus iteration number for D=100λ.

Download Full Size | PDF

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.

 

Fig. 5. (Left) Normalized focused irradiance distributions and (right) corresponding DOE phase profiles obtained when the vector diffraction algorithm was used to reshape the input beam to a focused square flattop irradiance profile for which (A, B) D=50λ, (C, D) D=25λ, (E, F) D=10λ, (G, H) and D=5λ. The DOE phase is plotted in units of radians.

Download Full Size | PDF

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 I_f 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.

Table 1. Performance parameters for the reshaped beams shown in Figs. 2 and 5 as function of the focused beam parameter D.

View Table

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 [25]. This phenomenon can be understood as a manifestation of the uncertainty principle [21]. The irradiance distribution I_t that can be achieved is inherently limited by the finite spatial bandwidth of the input beam I_in 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 [21]

where

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 β²ΔI_inΔI_t for the beam profiles of Fig. 5 are also included in Table 1. Because β²ΔI_inΔI_t 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 I_f will differ increasingly from I_t as the targeted beam size is decreased toward the diffraction limited spot size, with all other parameters kept fixed.

4. Conclusion

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.