## Abstract

Numerous applications require the simultaneous redistribution of the irradiance and phase of a laser beam. The beam shape is thereby determined by the respective application. An elegant way to control the irradiance and phase at the same time is from double freeform surfaces. In this work, the numerical design of continuous double freeform surfaces from ray-mapping methods for collimated beam shaping with arbitrary irradiances is considered. These methods consist of the calculation of a proper ray mapping between the source and the target irradiance and the subsequent construction of the freeform surfaces. By combining the law of refraction, the constant optical path length, and the surface continuity condition, a partial differential equation (PDE) for the ray mapping is derived. It is shown that the PDE can be fulfilled in a small-angle approximation by a mapping derived from optimal mass transport with a quadratic cost function. To overcome the restriction to the paraxial regime, we use this mapping as an initial iterate for the simultaneous solution of the Jacobian equation and the ray mapping PDE by a root-finding algorithm. The presented approach enables the efficient calculation of double freeform lenses with small distances between the freeform surfaces for complex target irradiances. This is demonstrated by applying it to the design of a single-lens and a two-lens system.

© 2017 Optical Society of America

## 1. INTRODUCTION

In recent years, the manufacturing of freeform surfaces has become increasingly feasible. These freeform surfaces offer an elegant way of simultaneous irradiance and phase control. Therefore, the development of numerical algorithms for the calculation of continuous freeform surfaces for control of the irradiance and/or the phase of a beam is of great interest.

In this work, the problem of designing continuous double freeform surfaces in a geometrical optics approximation is considered, in which two collimated beams of arbitrary irradiance are mapped onto each other. The simultaneous control of the phase thereby requires the second freeform surface. Hence, in contrast to the single freeform surface design for irradiance control, additional difficulties arise, since the surfaces can in general not be calculated independently. Several methods for phase and irradiance control with double freeform surfaces have been proposed in the literature.

One of the first design methods for the control of two sets of wavefronts by coupled freeform surfaces is the simultaneous multiple surface (SMS) method, which was developed by Benitez and Miñano [1]. The surfaces are thereby constructed from generalized Cartesian ovals and by applying constant optical path length (OPL) conditions [2]. The design method can be utilized for numerous applications in imaging and nonimaging optics [2–4].

In Ref. [5], the author reports to have solved the problem of laser beam shaping with double freeform surfaces by solving a Monge–Ampère-type equation. Details about the algorithm are not given.

Zhang *et al.* [6] and Chang *et al.* [7] solve the design problem by describing it in the form of a Monge–Ampère-type partial differential equation (PDE), discretizing the equation by finite differences and then solving the resulting nonlinear equation system by the Newton method. The design method can be applied to a variety of wavefront shapes [7]. Limitations and the performance of the method are unfortunately not discussed.

An alternative approach to construct freeform surfaces for irradiance and phase control is from ray-mapping methods [8–19]. These methods are based on the separation of the design process into two separate steps: the calculation of an *integrable* ray mapping between the source and the target irradiance, and the subsequent construction of the *continuous* freeform surfaces from the mapping. The integrability thereby ensures the continuity of the freeform surfaces and the mapping of the input irradiance onto the ouptut irradiance. Since the integrability depends on the physical properties of the optical system, it is in general a nontrivial task to find such a mapping.

As shown in several publications, there is a strong relation between the inverse problem of nonimaging optics and optimal mass transport (OMT) [20–25]. The cost function, which has to be applied to a certain optical configuration, is thereby problem-specific. For example, the mapping of two collimated beams with arbitrary irradiance onto each other with double freeform *mirrors* is described by a quadratic cost function [21] and can be solved by corresponding numerical schemes [19]. The same problem statement with double freeform *lenses*, which is considered in this work, is described by a different cost function, which depends on the OPL between the freeform surfaces as it was shown by Rubinstein and Wolansky [23].

The investigations presented here are inspired by several publications [9–14], in which the authors applied the quadratic OMT cost function to calculate a ray mapping to deal with the lens design problem. With this ray mapping, designs have been demonstrated of both single freeform surfaces for irradiance control for collimated input beams and point sources [9,10,14], and that of double freeform surfaces for irradiance and phase control [11–13]. As demonstrated for illumination control with single freeform surfaces in Ref. [17] and for collimated beam shaping with double freeform surfaces in Ref. [18], the design problems are thereby restricted to a paraxial approximation. Design methods for irradiance and phase control with continuous double freeform lenses using the quadratic cost function ray mapping like [11,12,18] are therefore inherently restricted to the paraxial regime. This means that either large enough distances between the freeform surfaces have to be considered or the surface continuity has to be given up, which leads to more difficult manufacturing processes and possibly to stronger diffraction effects. The main purpose of this work is therefore to present a detailed and efficient numerical algorithm to overcome these restrictions to the paraxial regime. Additionally, it is argued theoretically that the convergence of the presented algorithm can be ensured, which is of key importance for freeform design methods dealing with highly nonlinear PDEs.

Therefore, we investigate the design by ray-mapping methods of double freeform surfaces that map two collimated beams with arbitrary irradiance onto each other *beyond* the paraxial approximation. To overcome the restriction to the paraxial regime, which is necessary for the construction of compact systems, the design problem will be modeled by two coupled PDEs. This involves on one hand the Jacobian equation, expressing the local energy conservation, and on the other hand a ray-mapping PDE, enforcing the surface continuity and the constant OPL. The PDEs will then be solved by a root-finding scheme with the OMT mapping from the quadratic cost function as the initial iterate, leading to a construction approach for the freeform surfaces.

To do so, the work is structured as follows. In Section 2, by using the law of refraction, the constant OPL condition and a surface continuity condition, a PDE for an integrable ray mapping, is derived. Together with the Jacobian equation it builds a system of PDEs for the determination of the mapping components. It is argued that the PDE system is fulfilled within the paraxial approximation by the quadratic cost function OMT map. In Section 3, a method for solving the PDE system for general distances between the freeform surfaces is presented. It is based on discretizing the PDEs with finite differences and solving the resulting system of nonlinear equations by a standard root-finding scheme with the quadratic cost function OMT map as the initial iterate. A summary of the design algorithm and a detailed discussion of the implemenation is presented in Section 4, followed by the application of the presented method to the design of a single-lens and a two-lens system in Section 5. Finally, in Section 6, a short discussion of the results is presented.

## 2. FREEFORM DESIGN IN PARAXIAL APPROXIMATION

#### A. Energy Conservation and Cost Functions

In Ref. [17], a design method was presented for the construction of a single freeform surface for a collimated input beam with irradiance ${I}_{S}(x,y)$ and an arbitrary illumination pattern ${I}_{T}(x,y)$ on a target plane. It was shown that in the paraxial approximation, the design process can be decoupled into two separate steps. In the first step, a ray mapping $\mathbf{u}(x,y)=({u}_{x}(x,y),{u}_{y}(x,y))$ is calculated from the theory of OMT, and in the second step the freeform surface is constructed from the mapping.

There are several basic physical principles that a ray mapping needs to fulfill. First, to map the source irradiance ${I}_{S}(x,y)$ onto the target irradiance ${I}_{T}(x,y)$, the ray mapping should be energy conserving. The local energy conservation is expressed through the Jacobian equation

Second, in the case of freeform illumination optics, the mapping should allow the calculation of*continuous*freeform surfaces. As shown in several publications, these so called

*integrable*ray mappings are related to problem-specific cost functions representing different optical settings, where one has to distinguish between point sources and/or collimated beams, mirrors and/or lenses, and so on [20–25]. The cost function defines a metric between the source distribution ${I}_{S}(x,y)$ and the target illumination pattern ${I}_{T}(x,y)$, and therefore represents an additional constraint to the underdetermined Eq. (1). In the case of a single freeform surface for the redistribution of collimated input beams, the quadratic cost function

As shown by Rubinstein and Wolansky, the cost function for collimated beam shaping with double freeform *lenses* takes a different form than Eq. (2) [23]. The authors propose to minimize the corresponding cost function by a steepest-descent algorithm to get the ray mapping [23], but unfortunately, a numerically stable implementation is a nontrivial problem.

Due to its applicability in the paraxial approximation (see below) and the availability of numerous published stable numerical schemes for its calculation, the quadratic cost function OMT mapping serves as an initial iterate for the root-finding scheme presented below. It will therefore build the basis of the design approach presented in Sections 3 and 4.

#### B. Ray-Mapping Condition

We follow the approach from Refs. [17,18] by expressing the basic geometry according to Fig. 1 in terms of the collimated input and output vector fields ${\mathbf{s}}_{1}$ and ${\mathbf{s}}_{3}$, the refracted vector field ${\mathbf{s}}_{2}$ and the ray-mapping vector field ${\mathbf{s}}_{4}$

Since the goal is to calculate at least continuous freeform surfaces, we have to apply the surface continuity condition

for a general surface normal vector field $\mathbf{n}(x,y)$ to both freeform surfaces ${z}_{\mathrm{I}}(x,y)$ and ${z}_{\mathrm{II}}(x,y)$ and their normal vector fields ${\mathbf{n}}_{\mathrm{I}}(x,y)$ and ${\mathbf{n}}_{\mathrm{II}}(x,y)$, respectively. Thereby the normal vector fields can be expressed with the law of refraction in terms of the normalized incoming and refracted vector fields ${\widehat{\mathbf{s}}}_{1}$ and ${\widehat{\mathbf{s}}}_{2}$ and the refractive indices ${n}_{1}$ and ${n}_{2}$. Hence Eq. (5) can be written asThe left-hand side (LHS) of Eq. (8) represents the dot product of the projected gradient of the first surface $\nabla {z}_{\mathrm{I}}(x,y)=({\partial}_{x}{z}_{\mathrm{I}}(x,y),{\partial}_{y}{z}_{\mathrm{I}}(x,y))$ and the direction perpendicular to the ray deflection $(\mathbf{u}(x,y)-\mathbf{Id})$. Therefore, a nonvanishing RHS of Eq. (8) contradicts the law of refraction, which states that the incoming, the refracted, and the normal vector have to lie in the same plane. This can be seen directly by using the relation

A further relation can be derived by applying the chain rule to the equation ${\nabla}_{\mathbf{u}}{z}_{\mathrm{II}}({u}_{x},{u}_{y})=({\partial}_{{u}_{x}}{z}_{\mathrm{II}}({u}_{x},{u}_{y}),{\partial}_{{u}_{y}}{z}_{\mathrm{II}}({u}_{x},{u}_{y}))=\frac{{\mathbf{n}}_{\mathrm{II}}({u}_{x},{u}_{y})}{{({\mathbf{n}}_{\mathrm{II}})}_{z}}$ [27]. This provides us with the gradient $\nabla {z}_{\mathrm{II}}(x,y)$, which is used to rewrite the second term of the RHS of Eq. (8).

Hence from the continuity condition [Eq. (5)] and the law of refraction [Eq. (6)] follow the system of equations

These Eqs. (10) and (11), together with Eq. (1), build a system of PDEs for the unknown mapping $\mathbf{u}(x,y)$ and the surfaces ${z}_{\mathrm{I}}(x,y)$ and ${z}_{\mathrm{II}}(x,y)$.

To decouple the design process into separate steps as described in the beginning of this section, one needs to find a ray mapping fulfilling the condition [Eq. (11)] that is nontrivial without any *a priori* knowledge about the freeform surfaces. For single and double freeforms, this can be done by considering the small-angle approximation $({z}_{\mathrm{II}}({u}_{x},{u}_{y})-{z}_{\mathrm{I}}(x,y))\gg |\mathbf{u}(x,y)-\mathbf{Id}|$ leading to a vanishing first and second term in Eq. (11), since in this limiting case we get [using Eq. (4)] ${\mathbf{n}}_{\mathrm{I}}\xb7{\mathbf{s}}_{2}\sim ({z}_{\mathrm{II}}({u}_{x},{u}_{y})-{z}_{\mathrm{I}}(x,y))$, ${({\mathbf{n}}_{\mathrm{II}})}_{z}\xb7|{\mathbf{s}}_{2}|\sim ({z}_{\mathrm{II}}({u}_{x},{u}_{y})-{z}_{\mathrm{I}}(x,y))$, and the numerators are only dependent on the map $\mathbf{u}(x,y)$. Additionally, using the mapping from the quadratic cost function defined by Eq. (3), which gives $\nabla \mathbf{v}=0$, the condition [Eq. (11)] is fulfilled. Hence, the surfaces can be calculated by using Eq. (10) with appropriate boundary conditions. The boundary conditions can be derived by considering the law of refraction [Eq. (9)] on the boundaries of ${I}_{S}(x,y)$ and ${I}_{T}(x,y)$ [17]. As discussed in Ref. [17], this leads to a path-independent integration of Eq. (9) to calculate the surface ${z}_{\mathrm{I}}(x,y)$.

An alternative way to derive the validity of the quadratic cost function in the paraxial approximation results is by utilizing an expansion of the Rubinstein–Wolansky cost function [23] for small angles.

The double freeform design process can further be simplified by also using the constant OPL condition

Since ${\mathbf{n}}_{\mathrm{I}}\xb7{\mathbf{s}}_{2}$ and ${({\mathbf{n}}_{\mathrm{II}})}_{z}\xb7|{\mathbf{s}}_{2}|$ depend on ${z}_{\mathrm{II}}({u}_{x},{u}_{y})-{z}_{\mathrm{I}}(x,y)$, Eq. (11) can be written as a PDE for the components of mapping function $\mathbf{u}(x,y)$ only. Therefore, Eqs. (1) and (11) build a system of PDEs for the functions ${u}_{x}(x,y)$ and ${u}_{y}(x,y)$. Both equations build the basis of the design process for double freeform surfaces described in Sections 3 and 4.

Before we give any details, we want to discuss the condition [Eq. (11)] briefly for freeform mirrors.

### 1. Freeform Mirrors

For mirrors, the refractive indices in Eq. (6) have to be replaced by ${n}_{1}\equiv {n}_{2}\equiv -1$, and we get ${\mathbf{n}}_{\mathrm{I}}{\mathbf{s}}_{2}=-({\mathbf{n}}_{\mathrm{II}})|{\mathbf{s}}_{2}|$. Therefore, Eq. (11) reduces to

Hence, in contrast to the single-lens, single-mirror and double-freeform lens systems, the design problem can be solved by the quadratic cost function without any additional assumptions like the paraxial approximation.

## 3. FREEFORM LENS DESIGN BEYOND PARAXIAL APPROXIMATION

As mentioned in the previous section, Eqs. (1) and (11) are the basis of the design approach presented in the following. Since Eq. (11) is exactly fulfilled by the mapping with the quadratic cost functions for an infinite distance between the freeform surfaces or an infinite OPL, respectively, Eq. (11) represents a correction to Eq. (3) for finite OPLs. Hence, for finite distances between ${z}_{\mathrm{I}}(x,y)$ and ${z}_{\mathrm{II}}(x,y)$, we are searching for corrections $\mathrm{\Delta}\mathbf{u}(x,y)=(\mathrm{\Delta}{u}_{x}(x,y),\mathrm{\Delta}{u}_{y}(x,y))$ with $\mathrm{\Delta}\mathbf{u}(x,y)\stackrel{\mathrm{OPL}\to \infty}{\to}0$ to the ray mapping defined by the quadratic cost function, which we will denote by ${\mathbf{u}}^{\infty}(x,y)$ in the following. Hence, after writing Eq. (11) in terms of Eq. (3) plus a perturbation term, we want to solve the system of equations

The scalability of the mapping correction $\mathrm{\Delta}\mathbf{u}(x,y)$ with the parameter ${\mathrm{OPL}}_{\text{red}}$ thereby suggests solving Eqs. (15) and (16) by a root-finding method, since ${\mathrm{OPL}}_{\text{red}}$ can be reduced step by step with the solution $\mathrm{\Delta}\mathbf{u}(x,y)$ from the previous step as an initial iterate. This ensures the convergence of the root finding and can be done until the desired design goal or ${\mathrm{OPL}}_{\text{red}}$, respectively, is reached.

Additionally, we want to apply boundary conditions to Eqs. (15) and (16). To do so, we use standard boundary conditions for OMT problems by demanding that the boundary of the support of ${I}_{S}(x,y)$ is mapped onto the boundary of the support of ${I}_{T}(x,y)$. In the case of mapping two unit squares onto each other, such as in the example in Section 5, we therefore apply

We discretize Eqs. (15) and (16) using the standard central finite differences for the derivatives of $\mathrm{\Delta}{u}_{x}(x,y)\to {(\mathrm{\Delta}{u}_{x})}_{i;j}$ and $\mathrm{\Delta}{u}_{y}(x,y)\to {(\mathrm{\Delta}{u}_{y})}_{i;j}$ with $i=1,\dots ,N$; $j=1,\dots ,N$ at the inner points $i=2,\dots ,N-1$; $j=2,\dots ,N-1$ and second-order finite differences at the boundary. For the derivatives of, e.g., $\mathrm{\Delta}{u}_{x}(x,y)$, we get

## 4. DESIGN ALGORITHM

The design algorithm for the construction of double freeform surfaces presented in the previous section is summarized in the following. Since some of the steps offer freedom to choose between different applicable numerical techniques, we add some remarks below which are important for the examples in Section 5.

- 1. Calculate the OMT map with quadratic cost function ${\mathbf{u}}^{\infty}(x,y)$ between the distributions ${I}_{S}(x,y)$ and ${I}_{T}(x,y)$.
- 3. Choose the physical parameters (${n}_{1},{n}_{2},{\mathrm{OPL}}_{\text{red}}$) of the system and solve the system of nonlinear equations with a predefined tolerance and the initial iterate ${(\mathrm{\Delta}{u}_{x})}_{i;j}=0$ and ${(\mathrm{\Delta}{u}_{y})}_{i;j}=0$, with $i;j=1,\dots ,N$.

There are numerous publications presenting numerical methods for the calculation of OMT maps for quadratic cost functions ${\mathbf{u}}^{\infty}(x,y)$. In the example section below, we choose the numerical scheme developed by Sulman *et al.* [28]. It provides an efficient calculation of ${\mathbf{u}}^{\infty}(x,y)$, but has some drawbacks, such as the restriction to square supports of the irradiance distributions ${I}_{S}(x,y)$ and ${I}_{T}(x,y)$. In *our* implementation of Sulman’s algorithm, we recognized instabilities if the distributions show large irradiance gradients. For such distributions in the example in Section 5 we therefore choose a background irradiance ${I}_{T}(x,y)>0$ to ensure the convergence. An OMT method that is less restrictive can be found, for example, in Ref. [29].

For solving the system of nonlinear equations in step 3 we have used the *fsolve*() function from MATLAB’s 2015b optimization toolbox with the trust-region-reflective solver. Providing *fsolve*() with the structure of the Jacobian matrix of the object function allows an efficient calculation of the solution of the nonlinear equation system even for a large number of variables. The scalability with the parameter ${\mathrm{OPL}}_{\text{red}}$ of the distance of the initial iterate to the solution of Eqs. (15) and (16) suggests that the root finding could be accelerated by using, e.g., the Newton algorithm. In practice, we never experienced the necessity of iteratively reducing ${\mathrm{OPL}}_{\text{red}}$ to ensure the convergence of the algorithm. It was always possible to directly set ${\mathrm{OPL}}_{\text{red}}$ to the final design goal.

Solving the equation system in the form of Eqs. (15) and (16), we recognized oscillations in the solution $\mathrm{\Delta}\mathbf{u}(x,y)$. These seem to arise due to the irradiance root finding by the Jacobian Eq. (15), which in contrast to Eq. (16), is already fulfilled to a high degree by the initial iterate $\mathrm{\Delta}\mathbf{u}(x,y)=0$. Hence, we replace the $\mathrm{det}(\nabla {\mathbf{u}}^{\infty}(x,y)){I}_{T}(\mathbf{u}(x,y))$ term in Eq. (15) by

To calculate the first freeform, we use Eq. (9) with Eq. (13), which can be integrated along an arbitrary path to give the surface sag. The basic assumption for this path-independent integration is the satisfaction of Eq. (16), as mentioned in Section 2. Since Eq. (16) cannot be perfectly fulfilled numerically, the path integration leads to an accumulation of errors and therefore to deviations from the predefined plane wavefront.

For the second surface, one could use Eq. (13) directly, which gives this surface on a scattered grid. Since we want to test the design algorithm with ray-tracing software and due to advantages for possible freeform manufacturing processes, we calculate the second surface on a regular grid. This can either be done by scattered data interpolation, which might lead to inaccuracies in the irradiance pattern, or by applying the ray-tracing formula

## 5. EXAMPLES

In the following, the presented design algorithm is applied to two example target distributions. The first example consists of redistributing a Gaussian input beam with a waist of $1/\sqrt{2}$*a.u*. onto the letters “IAP,” and in the second example, the same Gaussian will be mapped onto the test image “house” (see Fig. 2). In “IAP,” the difficulties arise mainly due to the steep gradients, whereas “house” shows numerous gray levels. For the resolution of both images, $250\text{\hspace{0.17em}}\mathrm{pixels}\times 250\text{\hspace{0.17em}}\mathrm{pixels}$ are chosen.

Since input and on unit squares are of the same size, the freeform surfaces have the shape of squares with the side length 1*a.u*. For the example “IAP,” we choose as the physical parameters of the system, the refractive indices ${n}_{1}=1.5$ of the lenses, ${n}_{2}=1$ of the surrounding medium and the desired reduced OPL ${\mathrm{OPL}}_{\text{red}}=-0.2$*a.u*. are used. Therefore, a two-lens system, such as in Fig. 1(b) is considered. A distance between both freeform surfaces of approximately 0.4*a.u*. can be estimated from Eq. (13), since, due to symmetry reasons, the corner points of the irradiances are mapped onto each other. The reference ray for calculating the OPD intersects ${z}_{\mathrm{I}}(-0.5,-0.5)$ and ${z}_{\mathrm{II}}(-0.5,-0.5)$.

For the example “house,” we calculate a single-lens system defined by the parameters ${n}_{1}=1$, ${n}_{2}=1.5$, and ${\mathrm{OPL}}_{\text{red}}=0.2$*a.u*., leading to an approximate distance of 0.6*a.u*. between the freeform surfaces.

We want to point out that it is possible with the same algorithm to design systems with caustic surfaces between the freeform surfaces for the presented examples, which can be done by mirroring the initial map ${\mathbf{u}}^{\infty}(x,y)$ at the point of origin. This transformation preserves the property [Eq. (3)] and therefore $\nabla \mathbf{v}=0$ in Eq. (11) so that the arguments given in the previous sections hold analogously for this map. Using the mirrored intial map for the root finding and the mirrored boundary conditions, Eq. (18) leads therefore to the additional solution. We want to remark that the arguments given in the previous sections, which led to the presented design algorithm, only hold for ${\mathbf{u}}^{\infty}(x,y)$ and its mirror map, since $\nabla \mathbf{v}=0$ for both of them. This does not rule out the existence of additional solutions to the freeform design problem, such as saddle-shaped solutions with vertically or horizontally crossing rays.

The initial map ${\mathbf{u}}^{\infty}(x,y)$ can also be scaled and shifted by constants, which corresponds to a scaling of the size of the target irradiance and shifting of its position, without changing Eq. (3). These transformations are of course limited by physical boundaries, e.g., in form a total internal reflection. Additionaly, single-lens systems and two-lens systems can be calculated. Hence, the design method offers the freedom to choose between different optical configurations without a recalculation of the initial map ${\mathbf{u}}^{\infty}(x,y)$, which is one of the useful properties of the presented method.

We want to emphasize that the transformations above (scaling and shifting of the target irradiance or mirroring ${\mathbf{u}}^{\infty}(x,y)$) do not change the presented design algorithm. In our experience, the transformations do not introduce any additional numerical complications and are limited by physical constraints. Therefore, we are concentrating in the following on the case of noncrossing rays and source and target irradiances of the same size. The calculation of the initial map ${\mathbf{u}}^{\infty}(x,y)$ with our implementation of Sulman’s algorithm took 188 sec for “IAP” and 575 sec for “house” on an Intel Core i3 at $2\times 2.4$ Ghz with 16GB RAM.

After fixing the physical properties, the mapping has to be optimized according to Sections 3 and 4 to obtain the solution of the equation systems (15) and (16). To do that, one could use a starting value smaller/larger than the design goal ${\mathrm{OPL}}_{\text{red}}=-0.2$*a.u*. or ${\mathrm{OPL}}_{\text{red}}=0.2$*a.u*., respectively, to ensure the convergence to the solution of Eqs. (15) and (16). However, in practice we did not experience any significant benefit from using smaller/larger starting values for ${\mathrm{OPL}}_{\text{red}}$ than the design goal. For the tolerance of MATLAB’s *fsolve*(), the value of ${10}^{-6}$ was used. For both examples, the root-finding processes with 125,000 design variables took about 67 sec for “IAP” and 71 sec for “house.”

In Fig. 3 (“IAP”) and Fig. 4 (“house”), Eqs. (15) and (16) are plotted with the initial map ${\mathbf{u}}^{\infty}(x,y)$ and the optimized map $\mathbf{u}(x,y)$, respectively, to compare the quality of both mappings. Whereas the solution of the Jacobian Eq. (15) and therefore the local energy-preserving property remains nearly the same, the solution of the mapping Eq. (16) improves drastically in both cases.

The construction of the surfaces is done by the integration of Eqs. (9) and (13), first from $(x,y)=(-0.5,-0.5)$ along the $x$ direction and then along the $y$ direction. According to the arguments given in the previous sections, the corresponding surface will fulfill Eq. (5) and therefore be continuous. The system layouts with the freeform surfaces for both examples together with a few rays can be seen in Fig. 5.

To evaluate the improvement by the optimized map, the freeform surfaces are calculated from the initial mapping ${\mathbf{u}}^{\infty}(x,y)$ and the optimized map $\mathbf{u}(x,y)$ and imported into ray-tracing software for the calculation and the comparison of the illumination patterns and the wavefronts. The results of the ray tracing with $200\xb7{10}^{6}$ rays can be seen in Fig. 6 (“IAP”) and Fig. 7 (“house”). For quantitative comparisons of the results produced by ${\mathbf{u}}^{\infty}(x,y)$ and $\mathbf{u}(x,y)$, we calculated the rms, the peak-to-valley value $pv$ of the absolute difference $\mathrm{\Delta}{I}_{T}(x,y)\u2254({I}_{T}(x,y)-{I}_{T,RT}(x,y))$ between predefined distribution ${I}_{T}(x,y)$ from Fig. 2 and the ray-tracing illumination patterns ${I}_{T,RT}(x,y)$, the correlation coefficient ${\mathrm{corr}}_{{I}_{T}}$ between ${I}_{T}(x,y)$ and ${I}_{T,RT}(x,y)$, and the rms values of the optical path difference. The corresponding values are listed in Tables 1 and 2. The energy efficiency for the example “IAP” was 99.9% and for “house” 99.7%.

These tables and the images in Figs. 6 and 7 demonstrate a significant improvement of quality of the illumination patterns and the wavefronts for both examples. Fractions of the rms value of $\mathrm{\Delta}{I}_{T}$ thereby result from ray-tracing noise, which can be estimated by averaging over several ray tracings with different seeds to lower the ${\mathrm{rms}}_{\mathrm{\Delta}{I}_{T}}$ values by approximately $0.1\xb7{10}^{-6}$.

The freeform calculation with the initial maps ${\mathbf{u}}^{\infty}(x,y)$ shows strong deviations from the predefined specifications, whereas the optimized maps $\mathbf{u}(x,y)$ show high-quality illumination patterns and a wavefront uniformity far beyond the diffraction limit. Remaining deviations from the ideal wavefront (besides fundamental numerical limitations) are from Eq. (16), which is not fulfilled exactly, and therefore leads to an error accumulation along the integration path of Eqs. (9) and (13). The precision of the illumination pattern, on the other hand, is mainly limited by the precision of ${\mathbf{u}}^{\infty}(x,y)$. The main deviations from the predefined distribution ${I}_{T}(x,y)$ result from steep gradients, which can be seen especially for the example “IAP” in Fig. 6(d). This is in agreement with Fig. 3(b) and observations that were made in Ref. [17].

To demonstrate the applicability of the design algorithm to output beam sizes different from the input beam, we give the results for another example. Therefore, we choose the Gaussian input beam from the previous examples and the target irradiance distribution “house.” In contrast to the previous examples, we change the area of the target distribution to a square with a side length of 2*a.u*. and the reduced OPL to ${\mathrm{OPL}}_{\text{red}}=-0.8$*a.u*., which corresponds to ${z}_{\mathrm{II}}(-1,-1)-{z}_{\mathrm{I}}(-0.5,-0.5)=1.8598$*a.u*., according to Eq. (13). Since we have already calculated ${\mathbf{u}}^{\infty}(x,y)$ in the previous examples, we can simply scale the initial map by a factor of 2 and subsequently calculate $\mathbf{u}(x,y)$ from it. The calculation time of the root finding remains unchanged. The corresponding system layout can be seen in Fig. 8 and the results from the ray tracing in Fig. 9 and Table 3. Also for this example, the designed double freeform lens gives a high-quality target irradiance pattern [Fig. 9(b)] and a uniform wavefront [Fig. 9(f)], which shows the applicability of the presented algorithm to design examples with input and output beams of different size.

## 6. CONCLUSION

A design method for the calculation of continuous double freeform lenses for collimated beam shaping with complex irradiance patterns beyond the paraxial approximation was presented. The method is based on the ray-mapping condition [Eq. (16)], which was derived from the law of refraction and the surface continuity condition in Section 2 and builds together with the Jacobian Eq. (15) a system of nonlinear PDEs for the unknown ray mapping $\mathbf{u}(x,y)$.

Due to the satisfaction of Eqs. (15) and (16) for infinite lens distances by the mapping from OMT with the quadratic cost function [Eq. (2)], this mapping serves as an ideal initial iterate for a root-finding approach for solving the system of Eqs. (15) and (16). As was shown by approximating Eqs. (15) and (16) by finite differences and using a standard root-finding scheme from MATLAB’s optimization toolbox, one can ensure a fast convergence to the solution of Eqs. (15) and (16). This was demonstrated by applying the presented method to several design examples with complex target distributions and validating the results by ray tracing. The double freeform surfaces showed thereby high accuracy for the irradiance patterns and the wavefront, which was assessed by the calculation of the corresponding rms values of the normalized differences.

Further improvements can be made by using OMT methods for more complex boundaries of the source and target distributions, which requires the replacement of Sulman’s method (with e.g., [29]) and the generalization of Eq. (18) to more complex boundary shapes. The scalability of the distance of the initial map ${\mathbf{u}}^{\infty}(x,y)$ by ${\mathrm{OPL}}_{\text{red}}$ to the solution of Eqs. (15) and (16) suggests the application of, e.g., the Newton algorithm for a faster root finding.

In our future research, we want to generalize the presented approach to double freeform surfaces for wavefronts different from the plane case, such as, e.g., spherical wavefronts.

## Funding

Bundesministerium für Bildung und Forschung (BMBF) (O3WKCK1D, 031PT609X).

## Acknowledgment

The authors want to thank Ralf Hambach and Mateusz Oleszko for valuable discussions, Johannes Stock for comments on the manuscript, and David Musick for a grammar and spelling check.

## REFERENCES

**1. **P. Benítez and J. C. Miñano, “Simultaneous multiple surface optical design method in three dimensions,” Opt. Eng. **43**, 1489–1502 (2004). [CrossRef]

**2. **J. Chaves, *Introduction to Nonimaging Optics*, 2nd ed. (CRC Press, 2015), pp. 321–406.

**3. **J. C. Miñano, P. Benítez, and A. Santamaría, “Free-form optics for illumination,” Opt. Rev. **16**, 99–102 (2009). [CrossRef]

**4. **J. Chaves, J. C. Miñano, and P. Benítez, “Afocal video-pixel lens for tricolor LEDs,” Proc. SPIE **5942**, 594203 (2005). [CrossRef]

**5. **H. Ries, “Laser beam shaping by double tailoring,” Proc. SPIE **5876**, 587607 (2005). [CrossRef]

**6. **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. Comm. **331**, 297–305 (2014). [CrossRef]

**7. **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**, 125602 (2016). [CrossRef]

**8. **V. I. Oliker, J. Rubinstein, and G. Wolansky, “Ray mapping and illumination control,” J. Photon. Energy **3**, 035599 (2013).

**9. **A. Bäuerle, A. Bruneton, P. Loosen, and R. Wester, “Algorithm for irradiance tailoring using multiple freeform optical surfaces,” Opt. Express **20**, 14477–14485 (2012). [CrossRef]

**10. **A. Bruneton, A. Bäuerle, P. Loosen, and R. Wester, “High resolution irradiance tailoring using multiple freeform surfaces,” Opt. Express **21**, 10563–10571 (2013). [CrossRef]

**11. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Beam shaping system design using double freeform optical surfaces,” Opt. Express **21**, 14728–14735 (2013). [CrossRef]

**12. **Z. Feng, L. Huang, G. Jin, and M. Gong, “Designing double freeform optical surfaces for controlling both irradiance and wavefront,” Opt. Express **21**, 28693–28701 (2013). [CrossRef]

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

**14. **Z. Feng, B. D. Froese, and R. Liang, “Freeform illumination optics construction following an optimal transport map,” Appl. Opt. **55**, 4301–4306 (2016). [CrossRef]

**15. **L. L. Doskolovich, E. S. Andreev, S. I. Kharitonov, and N. L. Kazansky, “Reconstruction of an optical surface from a given source-target map,” J. Opt. Soc. Am. A **33**, 1504–1508 (2016). [CrossRef]

**16. **L. L. Doskolovich, E. A. Bezus, M. A. Moiseev, D. A. Bykov, and N. L. Kazanskiy, “Analytical source-target mapping method for the design of freeform mirrors generating prescribed 2D intensity distributions,” Opt. Express **24**, 10962–10971 (2016). [CrossRef]

**17. **C. Bösel and H. Gross, “Ray mapping approach for the efficient design of continuous freeform surfaces,” Opt. Express **24**, 14271–14282 (2016). [CrossRef]

**18. **C. Bösel and H. Gross, “Ray mapping approach in double freeform surface design for collimated beam shaping,” Proc. SPIE **9950**, 995004 (2016). [CrossRef]

**19. **N. K. Yadav, J. H. M. ten Thije Boonkkamp, and W. L. IJzerman, A Least-Squares Method for the Design of Two-Reflector Optical Systems, CASA-report 1619 (2016).

**20. **T. Glimm and V. I. Oliker, “Optical design of single reflector systems and the Monge–Kantorovich mass transfer problem,” J. Math. Sci. **117**, 4096–4108 (2003). [CrossRef]

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

**22. **X.-J. Wang, “On the design of a reflector antenna II,” Calc. Var. **20**, 329–341 (2004).

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

**24. **T. Glimm, “A rigorous analysis using optimal transport theory for a two-reflector design problem with a point source,” Inverse Probl. **26**, 045001 (2010). [CrossRef]

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

**26. **S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, “Optimal mass transport for registration and warping,” Int. J. Comp. Vis. **60**, 225–240 (2004). [CrossRef]

**27. **J. Rubinstein and G. Wolansky, “Reconstruction of optical surfaces from ray data,” Opt. Rev. **8**, 281–283 (2001). [CrossRef]

**28. **M. M. Sulman, J. F. Williams, and R. D. Russel, “An efficient approach for the numerical solution of Monge–Ampère equation,” Appl. Numer. Math. **61**, 298–307 (2011). [CrossRef]

**29. **J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Numerical solution of the optimal transportation problem using the Monge–Ampère equation,” J. Comp. Physiol. **260**, 107–126 (2014). [CrossRef]