1. Introduction

Generation of light beams with tailored irradiance distributions is of great interest for various applications including the problems of transforming and focusing laser beams, the formation of complex images in the augmented reality systems and holographic displays, the problems of structural lighting, optical trapping, microscopy, etc. [1–5]. One of the promising ways for this is the use of diffractive optical elements (DOEs), which in recent years again attracted considerable attention due to the advances in the design and fabrication methods. The problem of DOE design considered in the present work consists in calculating a phase diffractive microrelief providing the generation of a required irradiance distribution in a certain given plane [6]. Since the DOE microrelief height is proportional to the phase function of the light field generated by the DOE, the problem of calculating such elements is considered as the problem of calculating the DOE phase function, which is often referred to as the phase retrieval problem from two known intensity distributions. For solving this problem in the Fresnel and Fraunhofer approximations, various iterative Fourier transform algorithms (IFTAs) have been proposed, including the “classic” Gerchberg–Saxton (GS) algorithm and a wide range of its modifications [6–11].

Despite the large number of promising applications of DOEs [1–11], the existing iterative methods for DOE design (IFTAs) have several important disadvantages. In particular, the DOE microrelief calculated using an iterative algorithm has, as a rule, complex and irregular shape, often resembling white noise. Such a relief is difficult to fabricate, and the fabrication errors lead to uncontrolled scattering of light on the microrelief and to the formation of speckle structures, which significantly impair the quality of the generated distribution. In addition, despite the long history (the GS algorithm was proposed in 1972), the existing IFTAs still have significant drawbacks associated with slow convergence and stagnation of the iterative process at local minima.

In order to tackle these problems, “hybrid” methods based on the combination of IFTAs and geometrical-optics methods were proposed in recent years [2,12,13]. This approach is promising due to the fact that a phase function calculated in the geometrical optics approximation from the condition of focusing to a given region is smooth or piecewise-smooth. The use of such a phase function as an initial approximation usually provides fast convergence of the IFTA. With such an initial approximation, as a rule, the problem of stagnation of the algorithm at local minima does not arise. It is important to note that in this case, the phase function obtained using an IFTA retains a quasi-regular form, which significantly simplifies the fabrication of the DOE. Moreover, there exist modifications of the iterative algorithms, which allow one to control the “smoothness” of the obtained solution [2].

The key component of hybrid methods [2,12,13] is the method for calculating the geometrical-optics phase function. The problem of calculating this function can be reduced to solving a nonlinear differential equation (NDE) of elliptic type (the Monge–Ampére equation) [12–14]. The calculation of the geometrical-optics phase function in [2,12,13] was based on the numerical solution of this NDE using finite-difference methods [14–16]. Note that these methods are essentially iterative as well, since they reduce the problem to solving a system of nonlinear equations with respect to the values of the phase function on a certain grid, which is then solved using an iterative method. The described approach to the calculation of the geometrical-optics phase function imposes a number of significant restrictions on the class of the irradiance distributions that can be generated. Indeed, the formulation of the problem of calculating a geometrical-optics phase function as the problem of solving an NDE assumes that the calculated phase function is smooth. A smooth phase function can only implement a continuous ray mapping and therefore does not enable the formation of an irradiance distribution defined on a disconnected region or on a region with non-smooth boundaries.

In this work, we propose a method for calculating geometrical-optics phase functions enabling generating required irradiance distributions, including the distributions defined on disconnected regions or on regions with complex non-smooth boundaries. In the method, the problem of calculating the phase function is formulated in a semi-discrete form as a problem of maximizing a concave function, for solving which a gradient method is used. The expressions for the gradient are obtained in an analytical form. The proposed method is simple to implement and is based on just several main formulas. In contrast to the geometrical-optics NDE-based methods in [2,12,13], the proposed method enables calculating continuous and piecewise-smooth phase functions and thus allows one to generate “discontinuous” irradiance distributions defined on disconnected regions. As an example, a geometrical-optics phase function is calculated, which provides the generation of a discontinuous irradiance distribution in the form of a hexagram on a zero background. The utilization of the obtained geometrical-optics solution as an initial approximation in the IFTAs demonstrates the possibility of calculating DOEs with a quasi-regular or piecewise-smooth microrelief structure in the framework of the scalar diffraction theory. The calculation results are confirmed by the results of the experimental studies including the fabrication of a DOE and the investigation of the generated irradiance distributions.

The paper is organized as follows. In Section 2, the statement of the problem of calculating a DOE in the geometrical optics approximation is presented. This problem is formulated as a problem of calculating an eikonal function (or a geometrical-optics phase function) of the light field defined in a certain plane, from the condition of generating a prescribed irradiance distribution in another given plane. In Section 3, the proposed gradient method for calculating the eikonal function is described. Section 4 is dedicated to the problem of calculating the DOE phase function in the Fresnel approximation of the scalar diffraction theory. In Section 5, the results of the calculation of the geometrical-optics phase function for generating an irradiance distribution in the form of a hexagram are presented, as well as the results of the optimization of this function using the IFTAs. Finally, in Section 6, the experimental results of the fabrication and investigation of a DOE generating a hexagram-shaped irradiance distribution are presented.

2. Problem statement in the geometrical-optics approximation

As it was noted above, the problem of designing a DOE is usually formulated as a problem of the calculation of a phase function of the light field defined in a certain input plane from the condition of generating a prescribed irradiance distribution in another given plane [6]. However, in the geometrical optics approximation, it is more convenient to work not with the phase function, but with the eikonal function, which equals the phase function divided by the wavenumber $k=2\pi / \lambda$, where $\lambda$ is the wavelength. This is due to the fact that the finiteness of the wavelength of the incident radiation is not taken into account in the geometrical optics approximation.

Let us consider the inverse problem of calculating the eikonal function in more detail. Let, in a certain region $G$ in the plane $z=0$, an irradiance distribution $I_0 (\mathbf {u}),\, \mathbf {u}\in G$ and an eikonal function $\Phi (\mathbf {u}),\, \mathbf {u}\in G$ be defined, where $\mathbf {u}=(u_1 ,u_2 )$ are the Cartesian coordinates in the plane $z=0$ (Fig. 1).

 

Fig. 1. Geometry of the problem of calculating the eikonal function.

Download Full Size | PDF

The eikonal function $\Phi (\mathbf {u})$ defines the direction of the ray $\vec {e}(\mathbf {u})=\left (\nabla \Phi (\mathbf {u}),\sqrt {1-(\nabla \Phi (\mathbf {u}))^{2} } \right )$ coming out of a certain point $\mathbf {u}$ of the plane $z=0$ in the positive direction of the $z$ axis (we assume that the rays propagate in a medium with a unit refractive index). Here, $\nabla \Phi (\mathbf {u})=\left (\partial \Phi / \partial u_1 , \partial \Phi / \partial u_2 \right )$ and the upper arrow denotes a three-dimensional vector. The rays define a mapping $\gamma :G\to D$ between the points of the region $G$ of the input plane $z=0$ and the points of a certain region $D$ located in the plane $z=f>0$. Let us denote by $\mathbf {x}=(x_1 ,x_2 )$ the Cartesian coordinates in the plane $z=f$. Then, the “ray” mapping $\gamma$ can be defined in the following way. The point $\mathbf {x}=\gamma (\mathbf {u})\in D$ is defined as the intersection point of the plane $z=f$ and the straight line passing through the point $\mathbf {u}=(u_1 ,u_2 )$ with the direction vector $\vec {e}(\mathbf {u})$. In the coordinate form, the mapping $\gamma :G\to D$ reads as

Let us formulate the inverse problem of calculating the eikonal function from the condition of generating a prescribed irradiance distribution. Let, in a certain region $G$ of the plane $z=0$ and a region $D$ of the plane $z=f$, irradiance distributions $I_0 (\mathbf {u}),\, \mathbf {u}\in G$ and $I(\mathbf {x}),\, \mathbf {x}\in D$ be defined. It is necessary to calculate an eikonal function $\Phi (\mathbf {u}),\, \mathbf {u}\in G$ defining a ray mapping $\gamma :G\to D$, which satisfies the light flux conservation law (2) with the given distributions $I_0 (\mathbf {u})$ and $I(\mathbf {x})$.

For solving this problem, it is convenient to represent the eikonal function in a special form. Let us call the “lens eikonal” the eikonal of a converging spherical beam with the focus at the point $\mathbf {x}=(x_1 ,x_2 )$ of the plane $z=f$, written in the plane $z=0$. The lens eikonal function maps the whole region $G$ to the point $\mathbf {x}=(x_1 ,x_2 )$ and has the form [17]

Next, we assume the functions $\Phi (\mathbf {u})$ and $\Psi (\mathbf {x})$ to be continuous and piecewise-smooth. One can show [17] that in this case, the relations (6)–(8) remain valid with the difference that the mapping $\gamma :G\to D$ becomes multivalued on a set of measure zero. However, this change does not affect the light flux conservation law of Eq. (2).

3. Gradient method for calculating the eikonal function

In this section, we consider a gradient method for solving the problem of calculating the eikonal function. Let the required irradiance distribution $I(\mathbf {x}),\, \mathbf {x}\in D$ be approximated by a discrete distribution (a sum of Dirac delta functions) $I_\textrm {d} (\mathbf {x})=\sum _{i=1}^{N}I_i \delta (\mathbf {x}-\mathbf {x}_i ) ,\, i=\overline {1,N}$ defined on $N$ points $\mathbf {x}_i \in D$, which are the nodes of a certain grid covering the region $D$. In this case, the function $\Psi (\mathbf {x})$ is replaced by a set of $N$ values $\Psi _i ,\, i=\overline {1,N}$. It is these values that will be varied in the method. In this case, the eikonal function $\Phi (\mathbf {u})$ becomes a piecewise-smooth function consisting of fragments of lens eikonals of Eq. (3) providing the focusing to the points $\mathbf {x}_i ,\, i=\overline {1,N}$ of the plane $z=f$. According to Eq. (6), this function has the form

In the present work, for finding the maximum of the function $V(\Psi _1 ,\ldots , \Psi _N )$, we will use a gradient method. In this case, the values $\Psi _i, \, i=\overline {1,N}$ are calculated iteratively according to the rule

The gradient method of Eq. (15) has a simple physical interpretation. Indeed, from Eq. (10), it is easy to understand that with a decrease in $\Psi _i$ (at other $\Psi _j$ values, $j\ne i$, being fixed), the region $C(\mathbf {x}_i ; \Psi _1 ,\ldots ,\Psi _N )$ increases and, consequently, the light flux $I_{\textrm {est},i}$ described by Eq. (11) also increases. Therefore, if at the point $\mathbf {x}_i$, the required light flux $I_i >I_{\textrm {est},i}$, it is necessary to increase the $\Psi _i$ value, and to decrease it in the opposite case. In the gradient method of Eq. (15), the amount of increase or decrease in the $\Psi _i$ value is chosen proportionally to the difference $(I_i - I_{\textrm {est},i} )$ between the required and calculated values of the light flux.

The implementation of the gradient method requires calculating the $I_{\textrm {est},i}$ values in Eq. (11), which are expressed through integrals over the regions $C(\mathbf {x}_i ;\Psi _1 ,\ldots ,\Psi _{N} )$. For the calculation of these values, the ray tracing method can be used, which corresponds to the rectangle method for calculating the integrals of Eq. (11). In this case, in the region $G$ in the plane $z=0$, $N_\textrm {tr}$ rays are generated, going out of the nodes of a certain grid $\mathbf {u}_j, \, j=\overline {1,N_\textrm {tr} }$ introduced in the region $G$. These rays “carry” the light flux values proportional to the irradiance values at the grid nodes. According to Eqs. (7) and (9), for every ray going out of the point $\mathbf {u}_j$, the index of the lens eikonal $i(j)=\mathop {\arg \max }_i \left (\Psi _i -\rho (\mathbf {u}_j ,\mathbf {x}_i )\right )$ is found, which has the maximal value at the point $\mathbf {u}_j$. This ray is directed to the corresponding point $\mathbf {x}_i$. Then, at the points $\mathbf {x}_i \in D$ of the plane $z=f$, the light flux values are calculated as sums of the light flux values of the rays arriving to these points.

Above, we considered the gradient method for calculating the eikonal function $\Phi (\mathbf {u})$ [or the geometrical-optics phase function $\phi _\textrm {go} (\mathbf {u})=k\Phi (\mathbf {u})$] in a general nonparaxial case. In the problems of designing DOEs, it is usually assumed that the paraxial approximation condition is fulfilled. In this case, all the equations presented in this section retain their form, but instead of the distance $\rho (\mathbf {u},\mathbf {x})=\sqrt {f^{2} +(\mathbf {x}-\mathbf {u})^{2} }$, the following “paraxial” expression should be used:

4. DOE calculation in the Fresnel approximation

Let us now consider the problem of generating a light field with a prescribed irradiance distribution in the framework of the scalar diffraction theory. Let us denote by $w_0 (\mathbf {u})=A_0 (\mathbf {u})\exp \left (i\phi (\mathbf {u})\right )$ the complex amplitude of the light field in the plane $z=0$, where $\phi (\mathbf {u})$, $\mathbf {u}\in G$ is the phase function of the DOE, and $A_0 (\mathbf {u})=\sqrt {I_0 (\mathbf {u})}$, $\mathbf {u}\in G$ is the amplitude of the field at $z=0$. Then, in the Fresnel approximation, the complex amplitude of the field in the plane $z=f$ will be represented by the following integral:

In the Fresnel approximation, the problem of designing a DOE is formulated as the problem of calculating a phase function $\phi (\mathbf {u})$, which provides the generation of a light field with a prescribed irradiance $I(\mathbf {x})=\left |w(\mathbf {x})\right |^{2}$, $\mathbf {x}\in D$ in the plane $z=f$. Let us note that if the integral (17) is approximately calculated using the stationary phase method [20], the resulting irradiance $I(\mathbf {x})=\left |w(\mathbf {x})\right |^{2}$ will coincide with the irradiance distribution calculated in the geometrical optics approximation (in the paraxial case) [12]. The stationary point equation coincides with the paraxial equation of the ray. The stationary phase method has an $\textrm {O}(1 / k)$ error when applied to the calculation of the integral of Eq. (17). Accordingly, the geometrical-optics phase function $\phi _\textrm {go} (\mathbf {u})$ calculated in the approximation $k\to \infty$ $(\lambda \to 0)$ can be considered as an approximate solution from the point of view of the scalar diffraction theory.

For the rigorous solution of the problem, different IFTAs are used, including the Gerchberg–Saxton (GS) or error-reduction algorithm, as well as its various modifications [6–11]. In the GS algorithm, the phase function $\phi (\mathbf {u})$ is calculated iteratively, and each iteration includes the following four simple steps:

5. Design examples

In this section, the results of calculating a geometrical-optics phase function for generating a “discontinuous” irradiance distribution defined on a disconnected hexagram-shaped region [see Fig. 2(c)] are presented, as well as the results of the investigation of the “hybrid” approach to the DOE design, in which the geometrical-optics solution is used as the initial approximation in the iterative Fourier transform algorithms.

 

Fig. 2. (a), (b) Eikonal function (in mm) calculated from the condition of generating a uniform irradiance distribution in a hexagram-shaped region (c). Black lines in (b) show the derivative discontinuities of the eikonal function. (d) Normalized irradiance distribution generated in the target plane at the eikonal function shown in (a) and (b).

Download Full Size | PDF

For obtaining the geometrical-optics phase function, the eikonal function was calculated, defined on a circular region $G$ in the plane $z=0$, from the condition of generating a uniform irradiance distribution in a hexagram-shaped region $D$ [Fig. 2(c)] in the plane $z=f$. For calculating the eikonal function, the proposed gradient method presented in Section 3 was used with uniform irradiance distributions $I_0 (\mathbf {u}) = I_0, \, \mathbf {u}\in G$ and $I(\mathbf {x})=I, \, \mathbf {x}\in D$. The calculation was carried out with the following parameters: radius of the region $G$ (DOE aperture) $R=5\textrm {~mm}$, distance to the focal plane $f=100\textrm {~mm}$, length of the sides of the internal and external triangles constituting the hexagram $d_1 = 10\textrm {~mm}$ and $d_2 = 15\textrm {~mm}$, respectively, and the thickness of the hexagram lines $w=0.45\textrm {~mm}$.

In accordance with the proposed method, the required hexagram-shaped irradiance distribution was approximated by a discrete distribution containing $N=16768$ points. The points were located at the nodes of a square grid with the step $\Delta = 0.05\textrm {~mm}$ covering the hexagram region. The eikonal function was calculated iteratively using the gradient method of Section 3. At each iteration, the light flux values at the points of the discrete distribution were calculated using the ray tracing technique, and then the values $\Psi _i, \, i=\overline {1,N}$ were corrected using Eq. (15). In the ray tracing method, the incident beam was approximated by a set of $N_\textrm {tr} = 4N$ rays defined on a square grid in the plane $z = 0$. In the calculation, about 20000 iterations were made. The total time required for the calculation of the eikonal function on a desktop computer (Core i9-7940X CPU) amounted to approximately 2.5 hours.

As it was shown in Section 3, the proposed method (15) corresponds to the gradient method for maximizing the concave function $V(\Psi _1, \ldots ,\Psi _N)$ of Eq. (13). To illustrate the convergence of the method, the upper panel of Fig. 3(a) shows a graph of the function $V_\textrm {opt} - V(\Psi _1, \ldots ,\Psi _N)$ vs. the iteration number, where $V_\textrm {opt}$ is the final value of this function obtained at the last iteration. From Fig. 3(a), it is evident that the difference $V_\textrm {opt} - V(\Psi _1, \ldots ,\Psi _N)$ monotonically decreases. Since the function $V(\Psi _1, \ldots ,\Psi _N)$ does not explicitly represent the error of the generated irradiance distribution, in the lower panel of Fig. 2(a) we showed the convergence plot for the normalized root-mean-square deviation (NRMSD) of the generated distribution from the required one

 

Fig. 3. (a) Convergence rate of the gradient method: the difference $V_\textrm {opt} - V(\Psi _1, \ldots ,\Psi _N)$ (upper panel) and the NRMSD of the generated distribution from the required one (lower panel) vs. the iteration number. (b) The NRMSD plots vs. the iteration number for the GS algorithm with geometrical-optics initial phase (blue line), for the GS algorithm with the smoothing $\sigma =7.4\textrm {~}\mu \textrm {m}$ (red line), and for the GS algorithm with quadratic initial phase (orange line).

Download Full Size | PDF

The eikonal function $\Phi (\mathbf {u})$ calculated using the gradient method is shown in Figs. 2(a) and 2(b). The corresponding generated irradiance distribution shown in Fig. 2(d) was calculated using the ray tracing technique and is quite close to the required distribution: the normalized root-mean-square deviation (NRMSD) of the generated distribution from the required one [Fig. 2(c)] is only 4.5%. Let us note that the obtained eikonal function is continuous and piecewise-smooth. The sharp bends (derivative discontinuities) of the eikonal function, which are clearly visible in Fig. 2(a) and are shown with black lines in Fig. 2(b), are caused by the discontinuous character of the generated distribution. It is important to note that in the case of such irradiance distributions, the methods for calculating the geometrical-optics phase functions used in [2,12,13] and based on the numerical solution of the Monge–Ampère equation are either inapplicable or numerically unstable. Let us also note that in the recent work [2], in which the prospects of using the geometrical-optics phase function as an initial approximation for IFTAs are discussed and the results of calculating a freeform hologram generating a discontinuous distribution corresponding to an image of a robot are presented, the calculation of the geometrical-optics phase function was, in fact, not carried out. As an initial approximation for the IFTA in [2], only a quadratic initial phase was used, providing the formation of a defocused light spot with the dimensions approximately coinciding with the size of the generated image of a robot.

Figure 4(a) shows the geometrical-optics phase function $\phi _\textrm {go} (\mathbf {u}) = \textrm {mod}_{2\pi M} \left [k\Phi (\mathbf {u})\right ]$ at the wavelength $\lambda =550\textrm {~nm}$ ($k = 2\pi /\lambda = 11.424\textrm {~}\mu \textrm {m}^{-1}$). For the sake of visual clarity, the phase function is shown modulo $8\pi$ (at $M=4$). Due to the continuity of the eikonal function, the boundaries of the zones of the phase function are continuous, but contain sharp bends caused by the sharp bends of the eikonal function. Figure 4(b) shows the irradiance distribution $I(\mathbf {x})=\left |w(\mathbf {x})\right |^{2}$ generated in the Fresnel approximation in the case of the geometrical-optics phase function. The irradiance distribution was calculated numerically using a fast Fourier transform routine. Figure 4(b) demonstrates a reasonably good quality of the generated distribution, albeit with some diffraction effects being present (see the inset). The NRMSD of the generated distribution [Fig. 4(b)] from the prescribed uniform hexagram-shaped distribution amounts to 17.5%.

 

Fig. 4. (a), (c), (e) Geometrical-optics phase function (a) and the phase functions obtained using the GS algorithm (c) and the GS algorithm with phase smoothing (e) (all shown modulo $8\pi$). (b), (d), (f) Normalized irradiance distributions generated in the target plane in the case of the phase functions of Figs. (a), (c), and (e), respectively.

Download Full Size | PDF

Figures 4(c) and 4(d) show the results of the optimization of the phase function using the Gerchberg–Saxton algorithm [with the phase being updated using Eq. (20)] and the irradiance distribution obtained after 30 iterations, respectively. According to the convergence plot in Fig. 3(b), such a number of iterations turned out to be sufficient for the convergence of the algorithm. The irradiance distribution in Fig. 4(d) is close to the required one [Fig. 2(c)], and its NRMSD from the prescribed uniform hexagram-shaped distribution is only 2.1%. At the same time, the resulting phase function [Fig. 4(c)], although retained a quasi-regular form, acquired a high-frequency component [compare the insets to Figs. 4(a) and 4(c)]. In particular, it can be seen in the inset to Fig. 4(c) that the boundaries of the zones have become “indented”, and characteristic grating-type modulations have appeared inside the zones. The appearance of these irregularities complicates the DOE fabrication process.

To obtain a piecewise-smooth phase function, we optimized the phase in Fig. 4(c) using the GS algorithm with phase smoothing [Eq. (19)] with $\sigma =7.4\textrm {~}\mu \textrm {m}$. The corresponding convergence plot is shown in Fig. 3(b) as a continuation of the convergence plot for the conventional GS algorithm. After 70 iterations with smoothing, we obtained the phase function and the corresponding irradiance distribution shown in Figs. 4(e) and 4(f), respectively. The NRMSD of the generated distribution shown in Fig. 4(f) from the required one amounts to 7.4%. This value exceeds the value obtained in the previous case. At the same time, the resulting phase function [see Fig. 4(e) and the inset to it] is piecewise-smooth. Let us note that by choosing the parameter $\sigma$, it is possible to achieve a required tradeoff between the NRMSD of the generated distribution and the smoothness of the obtained phase function. The used value $\sigma =7.4\textrm {~}\mu \textrm {m} \approx 13\lambda$ was chosen in accordance with the technological capabilities of the equipment available to the authors and used for the subsequent DOE fabrication. Besides, since the chosen smoothing is one order of magnitude larger than the design wavelength, one can expect that the scalar diffraction theory provides an accurate estimation of the DOE performance, whereas the effects of rigorous electromagnetic theory can be neglected.

Thus, the successive optimization of the geometrical-optics phase function in the iterative GS algorithm and in the GS algorithm with smoothing ($\sigma =7.4\textrm {~}\mu \textrm {m}$) enabled obtaining a piecewise-smooth phase function providing the generation of the required discontinuous hexagram-shaped irradiance distribution with an error not exceeding 8%.

In order to demonstrate the advantages of using an exact geometrical-optics solution as an initial approximation in the IFTA, we also iteratively calculated the phase function with the following initial approximation:

Figure 5(a) shows the phase function obtained using the GS algorithm in the case of the quadratic initial phase of Eq. (21) (at $r=9\textrm {~mm}$). This phase function also provides the generation of the required irradiance distribution [Fig. 5(b)], but is highly irregular. Note that if a random or constant phase is used in the GS algorithm as an initial approximation, the resulting phase functions become even more irregular and resemble realizations of an uncorrelated random process like the white noise. It is also important to note that at the initial phase of Eq. (21) or in the case of a random initial approximation, the GS algorithm requires a significantly larger number of iterations [see the convergence plots in Fig. 3(b)]. The NRMSD of the obtained irradiance distribution [Fig. 5(b)] amounts to 9.3%, which exceeds the deviation obtained in the GS algorithm with phase smoothing (19) in the case of the geometrical-optics initial phase [Fig. 4(f)]. It is also important to note that a further optimization of the irregular phase function shown in Fig. 5(a) using the GS algorithm with phase smoothing [Eq. (19)] with the previously used smoothing value $\sigma =7.4\textrm {~}\mu \textrm {m}$ does not lead to satisfactory results. In this case, the NRMSD of the obtained irradiance distribution (not shown in the paper) exceeds 30%.

 

Fig. 5. (a) Phase function modulo $8\pi$ obtained using the GS algorithm with the quadratic initial phase of Eq. (21). (b) Normalized irradiance distribution generated in the target plane.

Download Full Size | PDF

In the opinion of the authors, the proposed gradient method for calculating a geometrical-optics phase function expands the capabilities of hybrid methods for the DOE design, especially, in the problem of generating discontinuous irradiance distributions. In this approach, the phase distributions obtained using the IFTA turn out to have a regular structure and be piecewise-smooth. This simplifies the DOE fabrication process and reduces the stray light and speckle effects associated with the light scattering on an irregular microrelief.

6. Experimental results

In order to verify the correctness of the presented design, we fabricated a DOE with the phase function of Fig. 4(e). Let us remind that this phase function was obtained from the geometrical-optics phase function using the GS algorithm followed by the GS algorithm with phase smoothing. The height of the DOE microrelief was calculated from the phase function in the thin optical element approximation as

The DOE was fabricated using the direct laser writing technique at a circular laser writing station CLWS-2014 [21]. The diffractive relief was patterned into a $4\textrm {~}\mu \textrm {m}$ layer of positive photoresist FP-3535, which was spin-coated on a quartz substrate. Figure 6(a) shows a photograph of the fabricated DOE, Fig. 6(b) shows an image of the central part of the DOE microrelief obtained using an optical microscope, and Fig. 6(c) demonstrates the central part of the DOE microrelief measured using the white light interferometer NewView 7300. The results of interferometric measurements of various fragments of the DOE microrelief showed that the maximum height differs from the specified value $h_\textrm {max}=3.62\textrm {~}\mu \textrm {m}$ by 100–150 nm.

 

Fig. 6. Photograph of the fabricated DOE (a), image of the central part of the DOE microrelief obtained in an optical microscope (b), and the central part of the DOE microrelief measured in the white light interferometer NewView 7300 (c).

Download Full Size | PDF

The fabricated DOE was investigated in an optical experiment. The optical setup of the experiment is shown in Fig. 7 and operates in the following way. A tunable laser (a) generates a beam with the required wavelength, which is focused by a micro-objective (b) on a pinhole (c). Then, a lens (d) forms a collimated beam, which impinges on the investigated DOE (e). Due to a relatively large size of the irradiance distribution generated by the DOE (the length of the side of the outer triangles constituting the hexagram amounts to $15\textrm {~mm}$), digital camera Canon 1100d was used as the sensor (f). The size of its matrix $22.2\times 14.7\textrm {~mm}^{2}$ was sufficient for imaging the generated distribution without additional scaling. The used tunable laser NT-242 was sequentially tuned to different wavelengths in the vicinity of the design wavelength $\lambda =550\textrm {~nm}$. The wavelength was tuned in order to compensate for the systematic error in the DOE microrelief height, which, as mentioned above, amounted to $100$–$150\textrm {~nm}$. According to visual assessment, the best irradiance distribution was obtained at the wavelength $532\textrm {~nm}$. This distribution is shown in Fig. 8(a) and is in a reasonable agreement with the numerical simulation results shown in Fig. 4(f). The darkening in the upper part of the experimentally obtained distribution is due to an alignment error between the DOE and the camera lens, as a result of which the upper edge of the beam formed by the DOE was partially cut-off. A certain irregularity of the irradiance distribution on the hexagram lines and a slightly nonzero irradiance in the inner part of the hexagram are due to the microrelief fabrication errors, in particular, the “smoothing” of the vertical boundaries of the microrelief zones. According to the estimates made by the authors, instead of the vertical “walls”, the used direct laser writing technique generates inclined sections with an angle up to 30–40 degrees to the vertical.

 

Fig. 7. Optical setup used for measuring the irradiance distribution generated by the fabricated DOE: (a) tunable laser NT-242, (b) 20$\times$ micro-objective, (c) $15\textrm {~}\mu \textrm {m}$ pinhole, (d) collimating lens, (e) DOE, (f) digital camera Canon 1100d.

Download Full Size | PDF

 

Fig. 8. Irradiance distributions generated by the fabricated DOE at the wavelengths 532 nm (a), 669 nm (b), and 457 nm (c).

Download Full Size | PDF

Let us remind that the fabricated DOE implements the calculated phase function modulo $8\pi$ [$M=4$ in Eq. (22), see Fig. 4(e)] and has microrelief height, which is 4 times larger than that of a “usual” DOE implementing a phase function modulo $2\pi$. Such DOEs, similarly to the harmonic diffractive lenses [22–26], enable working with the radiation of several (the so-called harmonic) wavelengths related by the following expression:

Thus, the presented experimental results confirm the correctness of the calculations carried out in Section 5, and additionally demonstrate the possibility of using the fabricated DOE at several operating wavelengths.

7. Conclusion

In the present work, a geometrical-optics method was proposed for calculating the eikonal function (or the geometrical-optics phase function), which provides the generation of a required irradiance distribution that can be defined on a disconnected region or on a region with complex non-smooth boundaries. In the method, the problem of calculating the eikonal function was formulated in a semi-discrete form as a problem of maximizing a concave function. For solving the maximization problem, a gradient method was used and analytical expressions for the gradient were obtained. The proposed method is simple to implement and is based on just several main formulas. Let us also note that instead of the gradient descent method of Eq. (15), one can utilize more efficient accelerated gradient methods that enable faster convergence for convex optimization problems [27].

Using the proposed method, a geometrical-optics phase function was calculated providing the generation of a “discontinuous” irradiance distribution in the form of a hexagram. It was shown that the utilization of the obtained geometrical-optics solution as an initial approximation in an IFTA makes it possible to design DOEs with a quasi-regular microrelief structure. In particular, the use of the geometrical-optics phase function in the GS algorithm with smoothing ($\sigma =7.4\textrm {~}\mu \textrm {m}$) enabled obtaining a piecewise-smooth phase function providing the generation of the required hexagram-shaped distribution with a normalized RMS error not exceeding 8%. The resulting phase function is significantly simpler (from the fabrication point of view) than the phase functions obtained in the GS algorithm with a quadratic or random initial phase.

The calculation results were confirmed by the results of experimental investigations, which included the fabrication of a DOE generating the hexagram-shaped irradiance distribution, and the measurement of the generated distributions.