## Abstract

We propose a version of the supporting quadric method for calculating a refractive optical element with two working surfaces for collimated beam shaping. Using optimal mass transportation theory and generalized Voronoi cells, we show that the proposed method can be regarded as a gradient method of maximizing a concave function, which is a discrete analogue of the Lagrange functional in the corresponding mass transportation problem. It is demonstrated that any maximum of this function provides a solution to the problem of collimated beam shaping. Therefore, the proposed method does not suffer from “trapping” at a local extremum, which is typical for gradient methods. We present design examples of refractive optical elements illustrating high performance of the method.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The problem of designing freeform optical elements generating light beams with prescribed irradiance distribution and wavefront has a wide range of applications including the design of lighting and backlight systems, laser beam processing and shaping, and optical lithography, among others. In the present work, we consider the problem of collimated beam shaping in the geometrical optics approximation. This problem consists in calculating a refracting (or reflecting) optical element, which transforms a given incident beam with a plane wavefront into an output beam with a prescribed irradiance distribution and a plane wavefront. Simultaneous control of both the irradiance distribution and the wavefront of the output beam requires the use of an optical element with two “working” optical surfaces.

In the geometrical optics approximation, the problem of collimated beam shaping can be reduced to finding a solution to a complex, fully nonlinear differential equation (NDE) of Monge–Ampère type. For the solution of this NDE, various finite-difference methods have been proposed [1–7]. The formulation of the problem of designing an optical element as an NDE assumes that the calculated surfaces of the optical element are smooth. This imposes a number of restrictions on the class of irradiance distributions that can be generated. In particular, an optical element with smooth surfaces cannot generate an irradiance distribution defined on a disconnected or a multiply connected region, or on a region with complex non-smooth boundary [8,9].

In the basic theoretical works [10–13], the problem of collimated beam shaping (PCBS) was formulated as a Monge–Kantorovich mass transportation problem (MTP) for the cases of a refractive optical element [10,12,13] and a system of two mirrors [11]. This MTP describes a problem of calculating an integrable ray mapping (i. e. a mapping relating the coordinates of the rays incident on an optical element to the coordinates of the refracted (or reflected) rays in the target plane), which provides the formation of an output beam with the required irradiance distribution. The calculation of an optical element in this MTP-based approach can be reduced to finding a solution to a linear programming problem [9] or a linear assignment problem [8]. In contrast to the NDE-based formulation, the MTP-based formulation of the PCBS enables calculating optical elements with continuous piecewise-smooth optical surfaces generating prescribed irradiance distributions defined on disconnected regions and on regions with complex non-smooth boundaries [8,9].

One of the widely used methods of designing refracting and reflecting optical surfaces is the supporting quadric method (SQM) [14–20]. It was initially proposed by V. Oliker and co-workers for solving the problem of calculating mirrors generating prescribed irradiance distributions [14–16]. The method was then generalized to the case of refracting optical surfaces [18–20]. The main advantages of the SQM are its versatility and simplicity. In the method, the prescribed irradiance distribution is approximated by a discrete distribution defined on a finite set of $N$ points. Then, the optical surface is represented as a set of $N$ segments of quadrics, each focusing the incident beam to one of these points. Depending on the problem being solved, paraboloids, ellipsoids, hyperboloids or more complex surfaces (e. g., Cartesian ovals [18]) are used as quadrics. In particular, in the problem of calculating a mirror generating a given discrete irradiance distribution in the far field, the mirror surface is represented as a set of segments of paraboloids. The calculation of the parameters of the quadrics is carried out iteratively. Note that in the above-mentioned problem of calculating a mirror, the convergence of this iterative method has been strictly proved [14].

In the present work, we consider a version of the SQM for designing a refractive optical element with two working surfaces for collimated beam shaping. The proposed method is simple to implement and is based on just a few main equations. We demonstrate that despite its simplicity, the method has firm theoretical foundations. Using the MTP-based formulation of the problem, we show that the proposed SQM is, in its essence, a gradient descent method of maximizing a concave function being a discrete analogue of the Lagrange functional in the MTP. It is important to note that any maximum (any critical point) of this function provides a solution to the PCBS. Therefore, the proposed version of the SQM does not suffer from the “false” convergence to a local maximum, which is typical for gradient methods. We illustrate the capabilities of the proposed method by designing optical elements confirming its high performance. The presented examples demonstrate the possibility of calculating optical elements with continuous piecewise-smooth surfaces generating prescribed irradiance distributions in doubly connected regions.

The paper is organized as follows. In Section 2, we formulate the problem of nonimaging optics being solved. In Section 3, we represent the surfaces of the designed element as envelopes of quadrics. Section 4 is dedicated to the description of the proposed version of the SQM. In Sections 5, we establish the connection between the SQM and the gradient method of maximizing a certain concave function associated with the MTP. In Section 6, we present two design examples illustrating the efficiency of the method. Section 7 concludes the paper.

## 2. Problem statement

The geometry of the considered problem is presented in Fig. 1(a). A light beam, which has a plane wavefront parallel to the plane $z=0$, impinges on an optical element with two refracting surfaces [Fig. 1(a)]. The irradiance distribution of the incident beam is defined by $I(\mathbf {x})$, $\mathbf {x} = (x_1, x_2)\in G$, where $G$ is a certain domain in the plane $z=0$. The first surface $R_1$ of the optical element is defined as $z=u_1 (\mathbf {x})$, $\mathbf {x}\in G$. It is convenient to define the second surface $R_2$ with respect to a certain “output” plane $z=f>0$, in which the irradiance distribution of the transmitted beam will be defined. Let us denote by $\mathbf {y} = (y_1, y_2)$ the Cartesian coordinates in this plane and define the second surface through the function $u_2 (\mathbf {y})$ as $z=f-u_2 (\mathbf {y})$. We assume that the refractive index of the medium between the surfaces [at $u_1 (\mathbf {x})\le z\le f-u_2 (\mathbf {x})$] equals $n_1$, and the refractive index of the surrounding medium [at $z < u_1 (\mathbf {x})$ and at $f-u_2 (\mathbf {x}) < z$] equals $n_0$.

An incident ray from a point $\mathbf {x}\in G$ of the plane $z=0$ is successively refracted by the surfaces $R_1$ and $R_2$ and arrives at some point $\mathbf {y}$ of the “target” domain $D$ in the output plane $z=f$ [Fig. 1(a)]. In this sense, the surfaces $R_1$ and $R_2$ define a ray mapping $\gamma : G\to D$, which relates the coordinates of the refracted rays $\mathbf {y}=\gamma (\mathbf {x})$ to the coordinates of the incident rays. This mapping determines the irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ generated in the plane $z=f$ according to the light flux conservation law:

where $J_{\gamma } (\mathbf {x})$ is the Jacobian of the mapping $\gamma$. The light flux conservation law can also be written in the integral formThe problem of collimated beam shaping, which we consider in the present work, is formulated in the following way. Assuming that the irradiance distribution of the incident beam $I(\mathbf {x})$, $\mathbf {x}\in G$ and the required irradiance distribution $L(\mathbf {y})$, $\mathbf {y}\in D$ in the output plane $z=f$ [Fig. 1(a)] are given, it is necessary to find such functions $u_1(\mathbf {x})$ and $u_2 (\mathbf {y})$ so that the transmitted light beam would have the given irradiance distribution $L(\mathbf {y})$ and a plane wavefront parallel to the plane $z=f$.

## 3. Representation of the optical surfaces

Let us consider a ray from a point $\mathbf {x}$ of the plane $z=0$, which, after being refracted by the surfaces $R_1$ and $R_2$, arrives at a point $\mathbf {y}=\gamma (\mathbf {x})$ of the plane $z=f$. Since the wavefronts of both the incident and the output beams are plane, the optical path length $F$ between these points is constant and does not depend on $\mathbf {x}$. This condition can be written as

Let us show that the surfaces $R_1$ and $R_2$ of the optical element can be represented as envelope surfaces of a special kind. Indeed, for differentiable surfaces $R_1$ and $R_2$, the ray mapping $\mathbf {y}=\gamma (\mathbf {x})$ can be described in the following way. The ray coming from a point $\mathbf {x}\in G$, after being refracted at the point $\left (\mathbf {x}, u_1(\mathbf {x})\right )$ of the surface $R_1$, goes to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$. Therefore, the surface $R_1$ can be represented as the envelope of a family of surfaces, each focusing the incident plane beam to one of the points $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ of the surface $R_2$ [Fig. 1(b)] [9].

The equation of the refracting surface $z=u_1^{\mathbf {y}} (\mathbf {x})$ focusing the incident beam to a certain point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ can be obtained by requiring the optical path length from the points $\mathbf {x}$ of the plane $z=0$ to the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$ to be constant and equal to $F-n_0 u_2 (\mathbf {y})$. Note that Eqs. (4) and (5) describe exactly this condition, and therefore the function $u_1^{\mathbf {y}} (\mathbf {x})$ has the form

Substituting Eq. (5) into Eq. (6), one can easily obtain that the surface $z=u_1^{\mathbf {y}} (\mathbf {x})$ is an ellipsoid of revolution if $n_0 <n_1$, and a hyperboloid of revolution if $n_0 >n_1$ [21]. The axis of this ellipsoid or hyperboloid is parallel to the $z$ axis, and one of its foci coincides with the point $\left (\mathbf {y}, f-u_2 (\mathbf {y})\right )$.Considering $\mathbf {y} = (y_1, y_2)$ as parameters, we can treat the equation $z=u_1^{\mathbf {y}} (\mathbf {x})$ as a two-parameter family of surfaces. The first surface of the designed optical element, which is described by the function $u_1 (\mathbf {x})$, is the envelope of this family [9] and is defined by the following equations [22]:

Similarly, the second surface $R_2$ can be considered as the envelope of a family of surfaces

each focusing a plane beam, which propagates “backwards” from the plane $z=f$, to one of the points $\left (\mathbf {x}, u_1 (\mathbf {x})\right )$ of the first surface $R_1$. In this case, the function $u_2 (\mathbf {y})$ is the envelope of the family of surfaces of Eq. (8) with respect to the parameters $\mathbf {x} = \left (x_1, x_2 \right )\in G$ and is defined by the following equations:In what follows, we will consider a particular case of the envelope of Eqs. (7a) and (7b), when the condition of Eq. (7b) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-u_2 (\mathbf {y})$ with respect to $\mathbf {y}$, but to its minimum. Such a choice provides the so-called Galilean configuration of rays [10]. In this case, the envelope surface defined by Eqs. (7a) and (7b) takes the form

Equations (10) and (11) make it possible to formulate the problem of collimated beam shaping as a problem of calculating such a function $u_2 (\mathbf {y})$, so that the corresponding ray mapping (11) satisfies the light flux conservation law of Eq. (1) or Eq. (2). Let us note that not every mapping satisfying the light flux conservation law is integrable, i. e. can be implemented by a pair of optical surfaces. In order for the mapping $\gamma$ to be integrable, it is necessary for the function $u_2 (\mathbf {y})$ defining it to exist and to satisfy Eq. (7b). Therefore, below we will refer to Eq. (7b) as the integrability condition.

## 4. Supporting quadric method

The supporting quadric method (SQM) [14–20] is widely used for designing mirrors and refracting surfaces generating prescribed irradiance distributions. In the method, the required continuous irradiance distribution is approximated by a discrete distribution consisting of a finite set of points. The designed optical surface is assumed to be piecewise-smooth and consisting of segments of quadrics. The parameters of the quadrics are found iteratively from the condition of obtaining the required light flux values at the points of the discrete distribution.

In this section, we consider a version of the SQM for solving the problem of collimated beam shaping [Fig. 1(a)]. Strictly speaking, a rigorous discussion of the SQM requires the introduction of the so-called “weak solution” concept, which makes it possible to correctly define the solution in the cases, in which the ray mapping is discontinuous on a set of measure zero [9]. However, for the sake of simplicity, we will concentrate on the practical side of the method. Let the prescribed continuous irradiance distribution $L(\mathbf {y})\in D$ be approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta (\mathbf {y}-\mathbf {y}_i)$ (a sum of Dirac delta functions) concentrated at $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $u_2 (\mathbf {y})$ in Eqs. (10) and (11) is replaced by a set of $N$ numbers $u_{2, i}$, $i=1, \ldots , N$. It is these values that are varied in the method. In this case, refractive surface $R_1$ is a piecewise-smooth surface consisting of segments of refractive lenses of Eq. (6) (ellipsoids or hyperboloids) focusing the incident plane beam to the points $(\mathbf {y}_i, f-u_{2, i})$, $i=1, \ldots , N$. This surface is schematically shown in Fig. 1(b) with bold segments of the colored curves showing the ellipsoids focusing the incident light to the individual points of the second optical surface. According to Eq. (10), the function $u_1 (\mathbf {x})$ defining the surface $R_1$ has the form

In practical calculations, the integral in the right-hand side of Eq. (15) can be calculated numerically using the ray tracing method, which is schematically shown in Fig. 1(b). In this approach, the incident beam is approximated by a set of $N_\textrm {tr}$ rays going from the nodes of a certain grid $\mathbf {x}_j$, $j=1, \ldots , N_\textrm {tr}$ in the plane $z=0$. These rays “carry” the light flux values $I_j$, $j=1, \ldots , N_\textrm {tr}$, which are proportional to the values of the irradiance of the incident beam at the grid nodes. According to Eq. (13), for each ray coming from a certain point $\mathbf {x}_j$, we find an index $i(j)=\mathop {\arg \min }\nolimits _i \left [{\mathcal {K}}(\mathbf {x}_j, \mathbf {y}_i)-u_{2,i} \right ]$ corresponding to the lens, which the ray intersects first [Fig. 1(b)]. This ray is directed to the corresponding point $\left (\mathbf {y}_{i(j)}, f-u_{2,i(j)} \right )$. Then, at the points of the second optical surface, the light flux values are calculated, which are equal to the sums of the light flux values of the rays focused to these points. In Fig. 1(b), these values are shown by numbers near each of the points of the “discretized” second optical surface. It is easy to see that with an increase in the $u_{2, i}$ value (in the case of fixed $u_{2, j}, j\ne i$), the domain $C(\mathbf {y}_i)$ increases, so that the light flux $L_{\textrm {est}, i}$ increases as well. Thus, if at some point $u_{2, i}$, the required light flux is greater than the obtained one [$L_{\textrm {d}} (\mathbf {y}_i)=L_i >L_{\textrm {est}, i}$], it is necessary to increase the $u_{2i}$ value (and to decrease it in the opposite case). It is natural to choose the magnitude of this increase or decrease proportional to the difference between the required and calculated light flux values $\left (L_i -L_{\textrm {est}, i} \right )$ [20]. In this way, the values $u_{2, i}$ can be found iteratively. At each iteration, the generated distribution $L_{\textrm {est}, i}$, $i=1, \ldots , N$ is calculated, and the $u_{2, i}$ values are updated according to the rule

It is important to note that according to the light beam reversibility principle, the considered PCBS can be formulated in a “reverse” form. In this case, the propagation direction is reversed and the input and required irradiance distributions are interchanged. In such a reverse formulation, one can discretize the irradiance distribution $I(\mathbf {x})$ [instead of $L(\mathbf {y})$] and then use an equation similar to Eq. (16) for the calculation of the values $u_{1,i}$ representing the first surface.

## 5. Connection of the SQM with a gradient method

#### 5.1 Mass transportation problem

In Ref. [10], it was shown that in the considered collimated beam shaping problem, the calculation of a ray mapping satisfying the light flux conservation law (1) can be reduced to finding a solution to a Monge–Kantorovich mass transportation problem (MTP) with the cost function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ of Eq. (5). The strict formulation of this fact is given by the following theorem:

**Theorem 1.** Let the mapping $\gamma : G\to D$ provide a minimum to the functional

**Proof.** Let us consider the Lagrange functional for the considered variational problem [25]:

The variations of the functional ${\mathcal {H}}(T, \lambda )$ can be easily calculated. From Eq. (20), one can easily show that the condition $\delta _{\lambda } {\mathcal {H}}=0$ of zero variation with respect to $\lambda$ is equivalent to the light flux conservation law [26]. To calculate the variation $\delta _{T} {\mathcal {H}}$ with respect to $T$, let us split the integral in the right-hand side of Eq. (20) into two integrals and make the change of variables $\mathbf {y}=T(\mathbf {x})$ in the second integral:

The obtained equations are equivalent to the integrability condition (7b). According to the conditional extremum theorem [25], if the mapping $\gamma$ provides a conditional extremum to the problem of Eqs. (17) and (18), then there exists a function $\lambda (\mathbf {y})$, which satisfies Eqs. (22). This means that the mapping $\gamma$ is integrable, i. e. can be implemented by a system of two refracting surfaces described by the functions $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})\equiv \lambda (\mathbf {y})$, and the function $u_1 (\mathbf {x})$ is defined through $u_2 (\mathbf {y})$ using Eq. (10). This concludes the proof of the theorem.

From the presented proof, it is evident that the solution of the collimated beam shaping problem can be obtained from any critical point of the Lagrange functional (20), and not only from the point corresponding to the minimum of the functional (17). In particular, the mapping $\gamma$, which provides a maximum to the functional (17), also gives a solution of the problem. From the form of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$ [see Eq. (5)], it is clear that the conditions of minimization and maximization of the functional (17) correspond to the Galilean and Keplerian configurations of rays, respectively [8,9].

Let us note that the ray mapping appears under the integral in the Lagrange functional (20). This makes it possible to use the MTP approach for calculating refracting surfaces implementing not only smooth mappings $\gamma$, but also mappings, which are discontinuous on a set of measure zero. In the case of a discontinuous mapping $\gamma$, the optical surface $u_1 (\mathbf {x})$ calculated using Eq. (10) will be piecewise-smooth.

#### 5.2 Modification of the Lagrange functional

From the proof of Theorem 1, it follows that if a pair $(T, \lambda )$ is a critical point of the Lagrange functional ${\mathcal {H}} (T, \lambda )$, then the mapping $\mathbf {y}=T(\mathbf {x})$ satisfies the light flux conservation law and is defined by a system of refracting surfaces $u_1 (\mathbf {x})$ and $u_2 (\mathbf {y})$, where $u_2 (\mathbf {y})=\lambda (\mathbf {y})$, and $u_1 (\mathbf {x})$ is found from Eq. (10). At a critical point, the mapping $\mathbf {y}=T(\mathbf {x})$ and the function $\lambda (\mathbf {y})$ are related by the integrability condition (22). Therefore, the mapping can be expressed through the function $\lambda (\mathbf {y})$ using Eq. (22). Similarly to the discussion above, where we restricted our consideration to the envelopes defined by Eq. (10), in what follows, we consider the critical points of a special kind. Namely, we assume that the condition of the equality of the derivatives to zero in Eq. (22) corresponds not just to some critical point of the function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})-\lambda (\mathbf {y})$, but to a minimum point. In this case, the mapping $T$ can be written in a form similar to the expression for the ray mapping (11):

#### 5.3 Generalized Voronoi cells and gradient descent method

In this section, we establish the connection between the functional (24) and the supporting quadric method (16). For this, let us consider the form of the functional (24) when the continuous irradiance distribution $L(\mathbf {y})$ is approximated by a discrete distribution $L_{\textrm {d}} (\mathbf {y})=\sum _{i=1}^{N} L_i \delta \left (\mathbf {y}-\mathbf {y}_i \right )$ defined on $N$ points $\mathbf {y}_i$, $i=1, \ldots , N$. In this case, the function $\lambda (\mathbf {y})$ is replaced by a vector $\boldsymbol \lambda = (\lambda _1, \ldots , \lambda _N)$, where $\lambda _i =\lambda (\mathbf {y}_i)$, and the mapping $T^{\lambda }$ becomes a mapping of the domain $G$ to a discrete set of points $\mathbf {y}_i$, $i=1, \ldots , N$.

**Definition.** Let us define the generalized weighted Voronoi cells as the following domains $C^{\boldsymbol \lambda } (\mathbf {y}_i)$: a point $\mathbf {x}\in G$ belongs to the domain $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ if and only if

In the case of a discrete distribution $L_{\textrm {d}} (\mathbf {y})$, the functional ${\mathcal {H}}(\lambda )$ becomes a function defined on ${\mathbb {R}}^{N}$:

In the Appendix, we consider the calculation of the partial derivatives of the function (26) for an arbitrary function ${\mathcal {K}}(\mathbf {x}, \mathbf {y})$. According to the results presented in the Appendix, these derivatives have the form

Note that since the function $\lambda (\mathbf {y})$ coincides with the function $u_2 (\mathbf {y})$ defining the surface $R_2$, then $u_{2,i}=\lambda _i$ and the generalized Voronoi cells $C^{\boldsymbol \lambda } (\mathbf {y}_i)$ coincide with the domains $C (\mathbf {y}_i)$ [Eq. (14)], on which the SQM quadrics are defined, which focus the incident beam to the points $\mathbf {y}_i$, $i=1, \ldots , N$. The values $L_{\textrm {est}, i} = \int _{C^{\boldsymbol \lambda } (\mathbf {y}_i)} I(\mathbf {x})\mathrm {d}{\mathbf {x}}$ are equal to the light flux focused to the points $\mathbf {y}_i$, $i=1, \ldots , N$. According to Eq. (27), at a critical point of the function $h(\boldsymbol \lambda )$, the following conditions hold:

It is important to note that the above theoretical analysis is different from the theoretical results of papers [10,12], which, in the discrete case, allow one to reduce the PCBS to a linear programming problem (LPP). In contrast, in the present work we reduce the PCBS to an unconstrained optimization problem, which consists in maximizing the concave nonlinear function $h(\boldsymbol \lambda )$ being a discrete analogue of the modified Lagrange functional (24). This approach, in comparison with the LPP-based formulations in Refs. [10,12], turns out to be much simpler in terms of the computational complexity. Indeed, the LPP representing the PCBS has large dimensions and, in most cases, cannot be solved using the commonly available LPP solvers (e.g. SOPLEX, MOSEK, or linprog MATLAB routine). For example, if both surfaces of the optical element are defined on grids with the dimensions of $100\times 100$ containing $N=10000$ points each (in this case, the incident beam and the generated collimated beam are also defined on grids with the same dimensions), the number of constraints in the mentioned LPP will be equal to $N^{2}=10000^{2}=10^{8}$ [10]. The LPP with such a large number of constrains is extremely complex for the existing LPP solvers. Thus, the “direct” application of the LPP-based technique of Refs. [10,12] is possible only for rather coarse meshes (with the sizes of about $50\times 50$), which is insufficient for many practical problems. The design examples presented in the next section demonstrate that the proposed SQM version of Eq. (16) [the gradient descent method of maximizing the function (26)] enables designing high-performance optical elements generating collimated beams with required irradiance distributions defined on the grids with sizes $100\times 100$ and $150\times 150$.

## 6. Design examples

#### 6.1 Generation of a “square-frame” beam with uniform irradiance

As a first example, let us consider the calculation of an optical element transforming a beam with a circular cross-section and uniform irradiance

For the calculation of the optical surfaces, the required irradiance distribution of the output beam (30) was approximated by a discrete distribution defined on a square grid with spacing $\Delta _y = 0.0173\textrm {~mm}$. In this case, the discrete distribution contains $N=100^{2}$ points. For the calculation of the second surface (of the values $u_{2,i}$, $i=1,\ldots ,N$), we used the SQM version described in Section 4. Under this approach, at each iteration the light flux at the points $\left (\mathbf {y}_i, f-u_{2,i} \right )$, $i=1,\ldots ,N$ is calculated using the ray tracing technique, and then the values $u_{2,i}$ are updated using Eq. (16). In the ray tracing method, the incident beam was approximated by a set of $N_\textrm {tr} =2N=20000$ rays defined on a square grid in the plane $z=0$. In the considered example, 25000 iterations were performed. The total calculation time on a desktop computer (Intel Core i9-7940X CPU) was of about 6 minutes.

After the calculations, the second surface defined by a discrete set of values $u_{2,i}$, $i=1,\ldots ,N$ was approximated using the least-squares approach by a system of two-dimensional B-splines [30]. Then, from the calculated second surface, the first surface was reconstructed using Eq. (10). The obtained optical surfaces are shown in Fig. 3(a). It is interesting to note that in the vicinity of the origin the first surface has sharp bends (derivative discontinuities) and resembles a tetrahedral pyramid. This is due to the fact that the central part of the first surface maps the central part of the incident beam (the neighborhood of the point $\mathbf {x}=0$) to the inner contour of the square-frame. This mapping is discontinuous. The possibility to design piecewise-smooth surfaces implementing discontinuous ray mappings is an important advantage of the SQM. In this case, the existing methods based on finding a numerical solution to a partial differential equation of elliptic type [1–7] are either inapplicable or numerically unstable.

Figures 3(b) and 3(c) show the calculated normalized irradiance distributions generated by the designed optical element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. Note that these distributions were calculated on a denser grid with the spacing $\Delta _{y} /4=0.0043\textrm {~mm}$ by tracing $10^{7}$ rays. The presented simulation results demonstrate the formation of a square-frame beam with the required dimensions and almost uniform irradiance. The normalized root-mean-square deviations of the obtained irradiance distributions from the uniform one amount to 3.3% and 3.8% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. Let us note that the dimensions of the square-frame beam at different distances from the optical element remain the same, therefore, the wavefront of the generated output beam is plane. Indeed, the calculated root-mean-square (RMS) deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $9\textrm {~nm}$ at $z=10\textrm {~mm}$.

#### 6.2 Generation of an annular beam with non-uniform irradiance

As the second example, we designed an optical element transforming the incident beam (29) with the radius $R=1\textrm {~mm}$ into an annular beam with a non-uniform irradiance distribution

The calculated surfaces of the optical element are shown in Fig. 4(a). As in the first example, the first surface has sharp bends near the origin of coordinates, which is due to the discontinuity of the ray mapping, which maps the neighborhood of the point $\mathbf {x}=0$ to the inner boundary of the ring. Figures 4(b) and 4(c) show the simulated normalized irradiance distributions generated by the designed element in the output plane $z=f=10\textrm {~mm}$ and in the plane $z=30\textrm {~mm}$. As in the previous example, the distributions were calculated on a grid with $0.0043\textrm {~mm}$ spacing by tracing $10^{7}$ rays. The presented results demonstrate the formation of an annular beam with the required dimensions. The normalized root-mean-square deviations of the calculated irradiance distributions from the required distribution defined by Eq. (31) amount to 6.2% and 6.7% at $z=10\textrm {~mm}$ and $z=30\textrm {~mm}$, respectively. The wavefront of the generated beam is almost plane and the calculated root-mean-square deviation of the optical path length from the constant design value $F=11\textrm {~mm}$ amounts to $13\textrm {~nm}$ at $z=10\textrm {~mm}$.

#### 6.3 Discussion of the convergence criteria

The SQM convergence is usually controlled by using criteria that represent the difference between the calculated and prescribed discrete distributions (e. g., the RMS deviation) [15,20]. However, such a choice is not optimal in the problem under consideration. Indeed, according to the results presented in Section 5.3, the SQM of Eq. (16) is a gradient method of maximizing a concave function $h(\boldsymbol \lambda )$ [Eq. (26)], which is a discrete counterpart of the Lagrange functional. To illustrate the concavity of this function, we chose two random vectors, $\boldsymbol \lambda _1$ and $\boldsymbol \lambda _2$, and plotted the function along a line segment connecting them: $h(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, $t\in [0, 1]$ [see the upper panel in Fig. 5(a)]. For comparison, the lower panel of Fig. 5(a) shows the RMS deviation between the vectors consisting of the prescribed, $\mathbf {L}$, and the calculated, $\mathbf {L}_\textrm {est}(t \boldsymbol \lambda _1 + (1-t)\boldsymbol \lambda _2)$, light flux values. The presence of several local minima and noisy behavior of the RMS deviation make it a poorer criterion for controlling the convergence of the iterative process than $h(\boldsymbol \lambda )$ or a certain monotonic function of it.

Let us illustrate this fact using as an example the optical element generating a square-frame beam (Fig. 3). The upper panel of Fig. 5(b) shows the plot of the quantity $h_\textrm {opt} - h(\boldsymbol \lambda )$, where $h(\boldsymbol \lambda )$ is calculated using Eq. (26) at each SQM iteration at $\boldsymbol \lambda =(u_{2,1}, \ldots ,u_{2,N})$, and $h_\textrm {opt}$ is the final value of this function obtained at the last iteration. From Fig. 5(b), it is evident that the difference $h_\textrm {opt} - h(\boldsymbol \lambda )$ monotonically decreases. An abrupt decrease in this value at iteration 5001 is caused by a decrease in the step $\alpha$ in Eq. (16): the first 5000 iterations were performed with the step $\alpha =1$, and then the step was decreased to $\alpha = 0.1$. A similar convergence plot for the RMS criterion is shown in the lower panel of Fig. 5(b). In contrast to $h(\boldsymbol \lambda )$, the RMS deviation graph exhibits a counterintuitive increase starting from the 5001st iteration, and is much more noisy. Therefore, the RMS criterion is not suitable for controlling the convergence, and it is the $h(\boldsymbol \lambda )$ function that should be used for this purpose.

## 7. Conclusion

We proposed a version of the supporting quadric method for designing refractive optical elements for collimated beam shaping. We demonstrated that this method is, in its essence, the gradient method of maximizing a concave function, which is a discrete analogue of the Lagrange functional in a corresponding mass transportation problem. We established that any maximum of this function provides a solution to the problem of collimated beam shaping. Therefore, the proposed method does not suffer from the convergence to a local extremum, which is typical for the “conventional” gradient methods. Another important advantage of the method is simplicity, since its implementation is based on only two simple equations: Eqs. (13) and (16). It is important to note that instead of the simple gradient descent method of Eq. (16), one can use more sophisticated convex optimization algorithms (e. g. Nesterov’s accelerated gradient method) to obtain a higher convergence rate.

As examples, we designed refractive optical elements, which transform a circular incident beam into a square-frame beam with uniform illuminance and to an annular beam with non-uniform illuminance. The presented examples demonstrate the possibility of designing optical elements with continuous piecewise-smooth surfaces and show high performance of the proposed method.

The authors believe that the connection between the SQM and the gradient method of maximizing a function related to a certain MTP can also be obtained for other problems of nonimaging optics, which can be formulated as an MTP. In particular, similar versions of the SQM can be formulated for the problems of calculating mirrors and refracting surfaces generating prescribed irradiance distributions in the far field. The differences will only appear when defining the generalized weighted Voronoi cells, the form of which is determined by the cost function, depending on the problem of nonimaging optics being solved.

## Appendix

Let us consider the calculation of the derivatives of the function $h(\boldsymbol \lambda )$ in Eq. (26). For the sake of convenience, let us introduce the function

Let $\boldsymbol \delta _m =(0,\ldots ,0,\delta ,0,\ldots ,0)$ be a vector of dimension $N$ having a single ($m$-th) nonzero component. Let us consider the difference

Let us divide the last integral in (34) into three integrals over the regions (35):

Since $C^{\boldsymbol \lambda + \boldsymbol \delta } (\mathbf {y}_m ) \subset C^{\boldsymbol \lambda } (\mathbf {y}_m)$, both minima defining the values $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m )$ and $M(\mathbf {x}, \boldsymbol \lambda )$ [see Eq. (32)] are reached at the point $\mathbf {y}_m$:

For the second integral in (36), the minima defining the values $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m)$ and $M(\mathbf {x},\boldsymbol \lambda )$ are reached at different points $\mathbf {y}_i$. Nevertheless, for every point $\mathbf {x} \in C^{\boldsymbol \lambda } (\mathbf {y}_m ) \backslash C^{\boldsymbol \lambda + \boldsymbol \delta _m} (\mathbf {y}_m )$, the weighted distance ${\mathcal {K}}(\mathbf {x}, \mathbf {y}_i)-\lambda _i$ to an arbitrary point $\mathbf {y}_i$ either did not change, or changed by $|\delta |$. Accordingly, the minimum of these distances $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m )$ changed by no more than $|\delta |$. Taking into account the sign of $\delta$, we can write

The third integral equals zero, because for any point $\mathbf {x}$ from the integration domain both minima $M(\mathbf {x},\boldsymbol \lambda + \boldsymbol \delta _m)$ and $M(\mathbf {x},\boldsymbol \lambda )$ are reached at the same $\mathbf {y}_i ,\, \, i\ne m$ and, consequently, their difference is zero.

Finally, using Eqs. (36), (37) and (39), one can easily obtain the derivatives of the function $h(\boldsymbol \lambda )$ [Eq. (26)] in the form of Eq. (27).

## Funding

Ministry of Science and Higher Education of the Russian Federation (State assignment to FSRC "Crystallography and Photonics" RAS); Russian Foundation for Basic Research (18-07-00982, 18-29-03067).

## Acknowledgments

The development of the design method and the investigation of the designed elements was supported by Russian Foundation for Basic Research, and the implementation of the software for the optical element design was supported by Ministry of Science and Higher Education of the Russian Federation.

## Disclosures

The authors declare no conflicts of interest.

## References

**1. **Y. Zhang, R. Wu, P. Liu, Z. Zheng, H. Li, and X. Liu, “Double freeform surfaces design for laser beam shaping with Monge–Ampère equation method,” Opt. Commun. **331**, 297–305 (2014). [CrossRef]

**2. **N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, “Computation of double freeform optical surfaces using a Monge–Ampère solver: Application to beam shaping,” Opt. Commun. **439**, 251–259 (2019). [CrossRef]

**3. **S. Chang, R. Wu, A. Li, and Z. Zheng, “Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge–Ampère type equation,” J. Opt. **18**(12), 125602 (2016). [CrossRef]

**4. **C. Bösel, N. G. Worku, and H. Gross, “Ray-mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation,” Appl. Opt. **56**(13), 3679–3688 (2017). [CrossRef]

**5. **X. Mao, J. Li, F. Wang, R. Gao, X. Li, and Y. Xie, “Fast design method of smooth freeform lens with an arbitrary aperture for collimated beam shaping,” Appl. Opt. **58**(10), 2512–2521 (2019). [CrossRef]

**6. **C. Bösel and H. Gross, “Double freeform illumination design for prescribed wavefronts and irradiances,” J. Opt. Soc. Am. A **35**(2), 236–243 (2018). [CrossRef]

**7. **Z. Feng, B. D. Froese, C.-Y. Huang, D. Ma, and R. Liang, “Creating unconventional geometric beams with large depth of field using double freeform-surface optics,” Appl. Opt. **54**(20), 6277–6281 (2015). [CrossRef]

**8. **L. L. Doskolovich, D. A. Bykov, E. S. Andreev, E. A. Bezus, and V. Oliker, “Designing double freeform surfaces for collimated beam shaping with optimal mass transportation and linear assignment problems,” Opt. Express **26**(19), 24602–24613 (2018). [CrossRef]

**9. **V. Oliker, L. L. Doskolovich, and D. A. Bykov, “Beam shaping with a plano-freeform lens pair,” Opt. Express **26**(15), 19406–19419 (2018). [CrossRef]

**10. **J. Rubinstein and G. Wolansky, “Intensity control with a free-form lens,” J. Opt. Soc. Am. A **24**(2), 463–469 (2007). [CrossRef]

**11. **T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge–Kantorovich mass transfer problem and Fermat’s principle,” J. Math. Mech. **53**(5), 1255–1278 (2004). [CrossRef]

**12. **V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Arch. Ration. Mech. Anal. **201**(3), 1013–1045 (2011). [CrossRef]

**13. **J. Rubinstein, G. Wolansky, and Y. Weinberg, “Ray mappings and the weighted least action principle,” J. Math. Ind. **8**(1), 6 (2018). [CrossRef]

**14. **L. A. Caffarelli, S. A. Kochengin, and V. I. Oliker, “On the numerical solution of the problem of reflector design with given far-field scattering data,” Ser. Contemp. Appl. Math. **226**, 13–32 (1999). [CrossRef]

**15. **S. Kochengin and V. Oliker, “Computational algorithms for constructing reflectors,” Comput. Visualization Sci. **6**(1), 15–21 (2003). [CrossRef]

**16. **V. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in * Trends in Nonlinear Analysis*, M. Kirkilionis, S. Krömker, R. Rannacher, and F. Tomi, eds. (Springer, 2003).

**17. **F. R. Fournier, W. J. Cassarly, and J. P. Rolland, “Fast freeform reflector generation using source-target maps,” Opt. Express **18**(5), 5295–5304 (2010). [CrossRef]

**18. **D. Michaelis, P. Schreiber, and A. Bräuer, “Cartesian oval representation of freeform optics in illumination systems,” Opt. Lett. **36**(6), 918–920 (2011). [CrossRef]

**19. **V. Oliker, “Controlling light with freeform multifocal lens designed with supporting quadric method (SQM),” Opt. Express **25**(4), A58–A72 (2017). [CrossRef]

**20. **L. L. Doskolovich, M. A. Moiseev, E. A. Bezus, and V. Oliker, “On the use of the supporting quadric method in the problem of the light field eikonal calculation,” Opt. Express **23**(15), 19605–19617 (2015). [CrossRef]

**21. **L. L. Doskolovich, A. Y. Dmitriev, E. A. Bezus, and M. A. Moiseev, “Analytical design of freeform optical elements generating an arbitrary-shape curve,” Appl. Opt. **52**(12), 2521–2526 (2013). [CrossRef]

**22. **L. P. Eisenhart, * A Treatise on the Differential Geometry of Curves and Surfaces* (Schwarz, 2008).

**23. **D. A. Bykov, L. L. Doskolovich, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam,” Appl. Opt. **58**(33), 9131–9140 (2019). [CrossRef]

**24. **L. L. Doskolovich, A. A. Mingazov, D. A. Bykov, and E. A. Bezus, “Formulation of the inverse problem of calculating the optical surface for an illuminating beam with a plane wavefront as the Monge–Kantorovich problem,” Comput. Opt. **43**(5), 705–713 (2019). [CrossRef]

**25. **L. C. Evans, “Partial differential equations and Monge–Kantorovich mass transfer,” in * Current Developments in Mathematics*, R. Bott, A. Jaffe, D. Jerison, G. Lusztig, I. Singer, and S.-T. Yau, eds. (International, 1999).

**26. **L. L. Doskolovich, D. A. Bykov, A. A. Mingazov, and E. A. Bezus, “Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions,” Opt. Express **27**(9), 13083–13097 (2019). [CrossRef]

**27. **Q. Mérigot, “A multiscale approach to optimal transport,” Comput. Graph. Forum. **30**(5), 1583–1592 (2011). [CrossRef]

**28. **F. Aurenhammer, F. Hoffmann, and B. Aronov, “Minkowski-type theorems and least-squares clustering,” Algorithmica **20**(1), 61–76 (1998). [CrossRef]

**29. **G. Wolansky, “On semi-discrete Monge–Kantorovich and generalized partitions,” J. Optim. Theory App. **165**(2), 359–384 (2015). [CrossRef]

**30. **C. de Boor, * A Practical Guide to Splines* (Springer-Verlag, 2001).