1. Introduction

The problem of collimated beam shaping consists in the design of a refractive optical element transforming a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront. In the general case, an optical element that performs this transformation has two working surfaces. As a rule, the problem under study is addressed within geometrical optics and has a wide range of applications, such as illumination, laser beam processing and optical lithography. The problem is easy to solve only in the axisymmetric case. In this case, the problem is reduced to a single spatial dimension, enabling the optical surfaces to be calculated by solving ordinary differential equations [1–6]. In the general two-dimensional case, the design of an optical element is essentially more challenging.

In [7–9], the problem of collimated beam shaping was reduced to finding a solution to an elliptic nonlinear partial differential equation (PDE). The PDE is solved using finite-difference methods, in which the derivatives are replaced by their finite-difference approximations. Thus, the solution of the PDE is reduced to the solution of a system of nonlinear equations with respect to the values of the sought-for function, which describes the first optical surface and is defined on a certain grid. The resulting system of equations is solved using Newton’s method. The second optical surface is then reconstructed analytically from the first surface by applying the requirement of a constant optical path length. General shortcomings of this approach include high computational complexity of solving a system of nonlinear equations as well as the choice of the initial approximation to the solution of the system and the formulation of the boundary conditions.

The considered problem can also be addressed by using another, mathematically less rigorous approach based on the so-called ray-mapping methods [10–12]. In this case, the ray mapping, i.e. the mapping between the coordinates of the rays on the wavefronts of the incident and the output beams is found by solving a standard Monge–Ampère (MA) equation. The MA equation is solved either using the above-mentioned finite-difference method based on reducing the task to the solution of a system of nonlinear equations, or using an explicit finite-difference method based on the replacement of the original equation by the non-stationary parabolic Monge–Ampère equation [11,12]. The optical surfaces are calculated on the basis of the computed ray mapping using numerical integration. In addition to having high computational complexity, the described ray-mapping methods can offer only an approximate solution because the standard MA equation is valid only in the paraxial approximation [9,13].

In the basic theoretical works [14,15], the problem of collimated beam shaping was formulated as a Monge–Kantorovich mass transportation problem, and the corresponding cost functions were derived. It is important to note that these cost functions are non-quadratic, meaning that in this case the ray mapping is not described by a standard MA equation. In a discrete version, the mass transportation problem can be formulated as a linear assignment problem (LAP) [13,16]. Such an approach is utilized in a recent work by the present authors [13] for solving an inverse problem of calculating the eikonal of a light field.

In this work, we for the first time apply this LAP-based approach to a problem of collimated beam shaping. Due to the fact that the element under study has two working optical surfaces, the application of the LAP-based method for the solution of this problem has a number of specific features and differs from [13]. As examples, we design refractive optical elements transforming a plane circular beam into plane beams with uniform irradiance distribution in a rectangular and a triangular domains in a nonparaxial case. These examples demonstrate high efficiency of the proposed LAP-based approach: the computation time for square grids of 143 × 143 points on a standard PC is of about 3 minutes. In contrast, the calculation of such a mapping using finite-difference methods for the solution of a PDE of Monge–Ampère type requires solving a system of 143² = 20449 nonlinear equations and is by several orders of magnitude slower [17]. Besides, the computation of the ray mapping by solving a mass transportation problem makes it possible to design elements with piecewise-smooth continuous optical surfaces, as distinct from smooth optical surfaces resulting from the solution of an elliptic nonlinear PDE [7–12]. It widens the range of the problems that can be solved, in particular, enabling the generation of uniform irradiance distributions in complex domains with non-smooth boundaries. As a demonstration of such a capability, we design an element with a piecewise-smooth optical surface transforming a circular beam into a cross-shaped beam.

2. Formulation of the problem

In this section, we formulate the problem within the geometrical optics framework using the associated notions of wavefronts, rays, and irradiances. The geometrical optics approximation implies that the characteristic dimensions of the optical surfaces and generated light distributions are much larger than the light wavelength.

Consider an incident beam with a plane wavefront perpendicular to the z axis and with the irradiance distribution I₀(u), u ∈ G, where u = (u₁, u₂) are the Cartesian coordinates in the plane z = 0 [Fig. 1(a)]. The beam is being transformed by a refractive optical element with the refractive index n > 1, which is surrounded by a medium with the refractive index n₀ = 1. The element has two optical surfaces f and g. The problem to be solved consists in the design of the optical surfaces generating an output beam with a plane wavefront perpendicular to the z axis and a prescribed irradiance distribution I(x), x ∈ D, where x = (x₁, x₂) are the Cartesian coordinates in a certain plane z = f₀ located after the optical element. Note that the function I(x) should satisfy the following normalization condition:

 

Fig. 1 (a) Geometry of the problem of the design of an optical element. (b) Approximation of a circular domain G and a triangular domain D with equal-flux cells. Dots correspond to the centers of the cells.

Download Full Size | PDF

Let x = T(u) = (x₁(u₁, u₂), x₂(u₁, u₂)) denote a ray mapping that represents the Cartesian coordinates x of a ray in the output plane z = f₀ through the coordinates of the ray u in the plane z = 0 [Fig. 1(a)]. In the geometrical optics approximation, we can use the energy conservation law along the ray tubes. According to this law, the mapping T: G → D conserves the light flux. The light flux conservation law can be written in the following integral form [13,14]:

Note that for a ray to remain parallel to the z axis after the refraction by the optical surfaces described by the functions z = f(u) and z = g(x), the normal vectors to the surfaces calculated at the points u and x(u) = T(u) have to be parallel. In this case, the optical surfaces can be locally considered as a plane-parallel plate, which shifts a ray without changing its propagation direction. Using the vector form of the refraction law, the partial derivatives of the functions z = f(u) and z = g(x) can be easily derived in the form [14]

It has been shown [14] that a mapping x = T(u) that corresponds to the solution of the problem in question satisfies the integrability conditions in Eq. (4) and provides an extremum of the following functional:

As mentioned above, the optical surfaces in the vicinity of the points u₀ and x(u₀) locally correspond to a plane-parallel plate that shifts the ray going from the point u₀ to the point x(u₀) in the output plane. It can be easily shown that the cost function in Eq. (6) turns to infinity when, at given parameters δ and n it is impossible to shift the ray by the magnitude |s₀| = |u₀ − T(u₀)|. If the cost function turns to infinity on a set with nonzero measure, there may be no solution to the problem. Note that with an increase in the parameter γ the achievable value of the shift also increases. In particular, we can increase the achievable shift by increasing δ and, hence, the thickness of the optical element.

One can show that a mapping x = T(u) that minimizes the functional (5) corresponds to a Galilean configuration of rays. The maximum of the functional corresponds to a mapping with a Keplerian configuration of rays.

In the paraxial approximation (i.e. under the condition $|\mathbf{x}(\mathbf{u})-\mathbf{u}{|}^2\ll {\gamma }^2)C(\mathbf{s})\approx -\gamma +\frac{|\mathbf{s}{|}^2}{2\gamma }$, and the mapping x(u) = T(u) is computed by solving the mass transportation problem with a quadratic cost function. In this case, should there exist a differentiable mapping x(u), its computation is reduced to finding a solution to a standard Monge–Ampère equation [9–12].

It is important to note that the direct determination of the ray mapping from the solution of a mass transportation problem in Eqs. (5) and (6) makes it possible to operate with mappings that are differentiable not everywhere, but almost everywhere. In this case, the optical surfaces [functions z = f (u) and z = g(x)] are continuous and piecewise-smooth, as distinct from smooth surfaces obtained as the solution of elliptic PDEs [7–9]. This makes it possible to extend the range of problems that can be solved, in particular, to generate prescribed irradiance distributions on a zero-intensity background in domains with complex and non-smooth boundaries.

3. Reduction to the linear assignment problem

The discrete version of the mass transportation problem in Eqs. (5) and (6) can be formulated as a linear assignment problem (LAP) [16]. Indeed, let us assume that the domains G and D are divided into N equal-flux cells (or approximated by N cells), with the equality ${\displaystyle {\int }_{g_i}{I}_0(\mathbf{u})\mathrm{d}\mathbf{u}}={\displaystyle {\int }_{d_j}I(\mathbf{x})\mathrm{d}\mathbf{x}=e}$ being valid for any pair of cells g ⊂ G, d_j ⊂ D. For the sake of illustration, Fig. 1(b) shows approximations of a circular domain G and a triangular domain D by grids of N = 196 square cells. While constructing these approximations, the grid steps were chosen so that N points fall within the domains G and D. This approach based on the choice of the rectangular grid steps is suitable for generating uniform irradiance distributions in the domains of interest. In the case of nonuniform distributions I₀(u) and I(x), an equal-flux set of N cells in the domain G (or D) can be obtained using one of the approaches utilized in the so-called source-target mapping methods based on the construction of an equal-flux grid in the source and target domains [18–20]. In particular, an equal-flux set of N = M₁ M₂ cells can be constructed in the following way. The domain G can be divided by parallel lines u₁ = u_1,i, i = 1, …, M₁ into M₁ equal-flux sub-domains. Afterwards, each sub-domain G_i is further divided into M₂ sub-domains by parallel lines u₂ = u_2,j, j = 1, …, M₂.

For the grids (approximations) introduced in the domains G and D, all energy-conserving mappings G → D can be described by permutations (i₁, …, i_N) of N numbers, which determine to which cells d_i of the domain D the cells g₁, …, g_N are mapped. In terms of the LAP, the mapping of the cell g_i to the cell d_j can be interpreted as the assignment of the j-th task to the i-th agent. According to Eq. (6), the cost of the tasks is described by the matrix

4. Reconstruction of the optical surface

Due to the discrete character of the LAP in Eq. (8), the resulting discrete mapping ${\mathbf{x}}_i=\widehat{T}({\mathbf{u}}_i)$, i = 1, …, N describes the sought-for continuous mapping x(u) = T(u) with an error. Besides, if the domain D has a non-smooth boundary, the mapping T(u) may have discontinuities. In this case, the reconstruction of the optical surfaces by direct numerical integration of Eqs. (3) may result in a significant error. Because of this, in this section we propose a method for reconstructing optical surfaces from a computed discrete mapping. Within the proposed method, the second optical surface [function g(x)] is assumed to be smooth, while the first surface [function f(u)] is assumed to be continuous and piecewise-smooth. This case corresponds to the most challenging example of an optical element described in Subsection 5.3 and designed to generate a cross-shaped beam.

4.1. Reconstruction of the upper optical surface

Let us first consider the calculation of the function g(x), which describes the second (upper) optical surface of the element. Because the function g(x) is assumed to be smooth, it can be approximated by two-dimensional B-splines [23].

Assume that the domain D, where g(x) is defined, is bounded by a rectangle $\overline{D}$ with the sides w₁ and w₂. Let x₁ = w₁· [m/(N₁ − 1) − 1/2, m = 0, …, N₁ − 1 be a set of N₁ equidistant knots on the interval [−w₁/2, w₁/2] of the x₁ axis. Using these knots, let us introduce basis splines (B-splines) of order q, denoting them by B_m(x₁), m = 1, …., N₁ + q − 2 [23]. In a similar way, let us introduce B-splines P_n(x₂) for the second coordinate x₂ ∈ [−w₂/2, w₂/2]. The function g(x) = g(x₁, x₂) can be represented using the introduced basis functions as

The coefficients p_m,n are determined from the condition that the values of the derivatives of the spline g(x) at the points ${\mathbf{x}}_i=\widehat{T}({\mathbf{u}}_i)$, i = 1, …, N constitute a least-squares approximation of the derivatives $\frac{\partial \widehat{g}}{\partial x_k}({\mathbf{x}}_i)$. The latter are derived from Eq. (3) on the basis of a discrete mapping $\widehat{T}$, i.e. by assuming $\mathbf{u}({\mathbf{x}}_i)={\widehat{T}}^{-1}({\mathbf{x}}_i)$ in Eq. (3). Thus, the coefficients p_m,n can be obtained by solving the following minimization problem:

Note that for a nonrectangular domain $D\subset \overline{D}$, the solution of the problem (10) may be numerically unstable. In this case, one should apply a regularized least-squares method, in which the term $\lambda {\displaystyle {\sum }_{m,n}{p}_{m,n}^2}$, where λ is a small number, is added to the minimized function S(p).

4.2. Reconstruction of the lower optical surface

Next, let us consider the calculation of the function f(u), which describes the first (lower) optical surface. As we mentioned earlier, the function f(u) is assumed to be piecewise-smooth, which does not allow us to approximate it by splines. We propose to reconstruct the function f(u) through the function g(x) using the representation of the optical surface from the supporting quadric method (SQM) [26–29]. The SQM is intended for designing a piecewise-smooth (refracting or reflecting) optical surface to focus the incident beam into a given array of points. In this case, the optical surface is composed of continuously stitched lens segments.

To apply the SQM to our problem, we assume that the second optical surface is defined on a sufficiently fine grid and is described by a set of points X_i = (x_i, g(x_i)), i = 1, …, N_g. Using the spline representation of the optical surface g(x) [Eq. (10)], the array of points X_i can be calculated on a uniform grid with an arbitrarily small step. In the vicinity of each point u ∈ G, the function f(u) corresponds to a lens segment that focuses the incident plane beam into one of the points X_i on the upper optical surface.

To derive the lens equation, let us use the condition of a constant optical path length from the points on the incident wavefront to the focal point [30,31]:

It is easy to show that the lens surface is described by an ellipsoid of rotation with the rotation axis parallel to the z -axis, eccentricity ε = 1/n, focal parameter p = δ/n, and one of the foci located at the point X_i [30, 31].

According to the SQM, the surface that focuses into the set of points X_i consists of a set of ellipsoid segments defined by Eq. (12). At each point u ∈ G, the surface is calculated from the following expression [26–29]:

Thus, we propose to calculate the function f(u) using Eqs. (12)–(13) through the second optical surface defined on a set of points X_i = (x_i, g(x_i)), i = 1, …, N_g. Note that Eq. (13) is also convenient for calculating the mapping T(u) at arbitrary values of u. Indeed, according to Eq. (13), we obtain

5. Design examples

5.1. Transformation of a circular beam to a rectangular beam

As a first example, we consider the design of an optical element transforming a beam with circular cross-section and uniform irradiance $I_0(\mathbf{u})=1/(\pi R^2)$, $\mathbf{u}\in G=\{(u_1,u_2)|u_1^2+u_2^2\le R^2\}$ to a beam with rectangular cross-section and uniform irradiance I(x) = 1/(w₁w₂), x ∈ D = {(x₁, x₂)| |x₁| ≤ w₁/2, |x₂| ≤ w₂/2}. The wavefronts of the incident and output beams are plane and perpendicular to the z axis. The element was designed for the following parameters: the refractive index of the element material n = 1.5, the incident beam radius R = 1 mm, the sizes of the output rectangular beam w₁ = 5 mm and w₂ = 2.5 mm, the output plane z = f₀ = 10 mm, and the parameter δ = 2.5 mm. For the considered parameters, the optical path length from the plane z = 0 to the plane z = f₀ is L = f₀ + δ = 12.5 mm. In this example, the irradiances I₀(u) and I(x) are centrally symmetric. Therefore, the central ray originating from the point u = (0, 0) will propagate through the optical element without refraction. This allows us to calculate the element thickness at u = (0, 0) as h₀ = δ/(n − 1) = 5 mm.

The mapping x(u) = T(u) was found by solving the LAP of Eq. (8). For the calculation of the mapping, the rectangle D was represented as a set of N = 143² = 20449 rectangular cells d_i with the size 0.0282 × 0.0141 mm². The circular domain G was also approximated by 20449 square cells g_i with the size 0.0124 × 0.0124 mm². For the solution of the LAP (8), the auction algorithm [22] implemented in [32] was used. For the considered example, the time of solving the LAP was less than 3 minutes. For comparison, let us note that the calculation of this mapping using finite-difference methods for solving PDE of Monge–Ampère type [7–9] requires the solution of a system of 20449 nonlinear equations and is by several orders of magnitude slower [17].

The optical surfaces of the element were reconstructed from the computed discrete mapping using the technique described above in Section 4. The function g(x) was represented as a two-dimensional B-spline at q = 4 and N₁ = N₂ = 10. The coefficients of the spline were found using the least squares method from the condition of minimizing expression (10). Then, the second optical surface defined as spline in Eq. (9) was recalculated on a uniform 400 × 400 rectangular grid (N_g = 400²). Then, using Eqs. (12) and (13), the first optical surface was calculated on a uniform 143 × 143 square grid. The resulting optical surfaces are shown in Fig. 2(a). For the sake of illustration, Figs. 3(a) and 3(b) depict the mapping T(u) calculated using Eq. (14). Figure 3(a) shows a square grid in the domain G, whereas Fig. 3(b) shows its image by the mapping T(u). It should be noted that in the considered example, the angles of ray deflection by the optical surfaces (which reach maximum at the corners of the rectangle D) are as large as 25°. Hence, the optical element is designed in a non-paraxial case.

 

Fig. 2 (a) Surfaces of an optical element transforming a circular beam into a beam with rectangular cross-section. (b, c) Normalized irradiance distributions generated by the element in the planes z = 10 mm and z = 30 mm calculated using TracePro. The irradiance cross-sections along the coordinate axes are shown at the top and at the right of the distributions.

Download Full Size | PDF

 

Fig. 3 The mappings T: G → D corresponding to the transformation of a circular beam (a) to a rectangular beam (b), a triangular beam (c), and a cross-shaped beam (d). Blue lines mark a straight line in the G domain and its images in the generated domains D.

Download Full Size | PDF

In order to verify the proposed design method, the performance of the optical element was numerically simulated using the commercial ray-tracing software TracePro [33]. To do this, the optical surfaces shown in Fig. 2(a) were exported to a computer-aided design software Rhinoceros [34], in which a 3D-model of the element was created [Fig. 2(a)]. Then this model was exported to TracePro for simulation. Figures 2(b) and 2(c) show the simulated normalized irradiance distributions generated by the optical element in the planes z = 10 mm and z = 30 mm. The simulation results demonstrate the formation of a rectangle with the required size and near-uniform irradiance. The normalized root-mean-square deviations (NRMSD) of the numerically simulated irradiance distributions from a constant value amount to 5.3% and 6.1% at z = 10 mm and z = 30 mm, respectively. The fact that the size of the rectangle remains unchanged at different distances from the element shows that the wavefront of the beam generated by the element is nearly plane. Actually, the root-mean-square deviation of the optical path length from a constant value at z = 10 mm does not exceed 1 nm.

5.2. Transformation of a circular beam to a triangular beam

As a second example, let us consider the design of an optical element transforming a uniform circular incident beam to a uniform triangular beam. The equilateral triangle (the domain D) has a 3-mm side, with the rest geometric and simulation parameters (R, f₀, δ, n, N, q, N₁, N₂, N_g) remaining the same. The designed optical surfaces and the mapping T(u) are depicted in Fig. 4(a) and Figs. 3(a) and 3(c), respectively. Figures 4(b) and 4(c) show the simulated normalized irradiance distributions generated by the optical element in the planes z = 10 mm and z = 30 mm obtained using TracePro. The simulation results show that a triangle-shaped irradiance domain with the required size is generated. The NRMSD of the calculated irradiance distributions from a uniform value amount to 5.6% and 6.4% at z = 10 mm and z = 30 mm, respectively. As in the previous example, the root-mean-square deviation of the optical path length from a constant value at z = 10 mm does not exceed 1 nm.

 

Fig. 4 (a) Surfaces of an optical element transforming a circular beam to a triangular one. (b, c) Normalized irradiance distributions generated by the element in the planes z = 10 mm and z = 30 mm calculated using TracePro. The irradiance cross-sections along the coordinate axes are shown at the top and at the right of the distributions.

Download Full Size | PDF

5.3. Transformation of a circular beam to a cross-shaped beam

The final and the most challenging example involves the design of an optical element transforming a uniform incident beam with circular cross-section into a beam with cross-shaped uniform irradiance distribution. In this case, the domain D corresponds to a cross obtained by a superposition of two mutually perpendicular rectangles with the sizes 2.8 × 1 mm² and 1 × 2.8 mm². The rest geometric and calculation parameters remain the same. The resulting optical surfaces are shown in Fig. 5(a). Note that the lower surface of the optical element has ‘breaks’ (sharp bends) along the bisector in each quadrant of the coordinate system. Figure 6 depicts a magnified fragment of the lower optical surface with a break, which is highlighted by a red rectangle in Fig. 5(a). These breaks are caused by the discontinuities in the mapping x(u)= T(u). Indeed, the rays that strike the lower optical surface on the opposite sides of the bisector are mapped to different parts of the cross. This is confirmed by the form of the mapping T(u) shown in Fig. 3(d). For instance, the image of a straight blue line shown in Fig. 3(a) has discontinuities near the inner corners of the cross and is composed of several curvilinear segments. It is worth noting that since the inverse mapping T⁻¹: D → G is continuous, the obtained upper optical surface is smooth.

 

Fig. 5 (a) Surfaces of an optical element transforming a circular beam to a cross-shaped beam. A magnified fragment of the lower optical surface shown in Fig. 6 is highlighted with a rectangle. (b, c) Normalized irradiance distributions generated by the element in the planes z = 10 mm and z = 30 mm calculated using TracePro. The irradiance cross-sections along the coordinate axes are shown at the top and at the right of the distributions.

Download Full Size | PDF

 

Fig. 6 A magnified fragment of the lower surface [f(u)] of the designed optical element [see Fig. 5(a)]. The break is marked with a red line.

Download Full Size | PDF

The simulation results obtained using TracePro and shown in Figs. 5(b) and 5(c) demonstrate that a high-quality cross-shaped irradiance distribution is generated. The NRMSD of the calculated irradiance distributions from a constant value are 7.1% and 7.9% at z = 10 mm and z = 30 mm, respectively. The root-mean-square deviation of the optical path length from a constant value is 1.1 nm at z = 10 mm. The authors believe that the somewhat larger NRMSD of the generated irradiance distributions from a constant value is due to the approximation error of the piecewise-smooth lower optical surface in the used CAD software Rhinoceros. This software adopts non-uniform rational basis splines (NURBS) for representing a surface from a point cloud. As a result, the sharp bends of the lower optical surface are somewhat smoothed.

The third example demonstrates one of the advantages of the proposed method, which is the capability of operating with discontinuous mappings. In this case, the known methods based on the numerical solution of an elliptic PDE [7–12] are either not applicable or numerically unstable.

6. Conclusion

We proposed a method for designing refractive optical elements for collimated beam shaping. The computation of a ray mapping is reduced to a mass transportation problem with a non-quadratic cost function. In the proposed method, the problem of computing a ray mapping is formulated as a linear assignment problem, which constitutes a discrete version of the corresponding mass transportation problem. A technique for reconstructing an optical surface from the computed ray mapping has also been proposed. As distinct from the methods based on the numerical solution of an elliptic PDE [7–12], the method proposed here is suitable for the design of optical elements with continuous piecewise-smooth optical surfaces.

We have designed refractive optical elements transforming a beam with circular cross-section into variously shaped (rectangular, triangular, and cross-shaped) beams with plane wavefronts. The proposed method is computationally efficient. It requires less than 3 minutes on a standard computer to compute a ray mapping on a 143 × 143 grid. High performance of the designed optical elements is confirmed by the numerical simulation results.