Ray mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation

TL;DR

This paper presents a numerical method for designing double freeform surfaces for collimated beam shaping that surpasses the paraxial approximation, enabling complex irradiance and phase control in laser beam applications.

Contribution

It introduces a combined approach using optimal mass transport and optimization to solve the PDE for non-paraxial beam shaping with double freeform surfaces.

Findings

Effective design of complex beam shapes demonstrated

Method surpasses paraxial approximation limitations

Applications include single-lens and two-lens systems

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 an optimization. The presented approach enables the efficient calculation of compact double freeform surfaces for complex target irradiances. This is demonstrated by applying it to the design of a single-lens and a two-lens system.

Figures

Equations

$$det(\nabla \mathbf{u}(x,y))I_T(\mathbf{u}(x,y))=I_S(x,y).$$

$$d(I_S,I_T{)}^2=in{f}_{\mathbf{u}\in M}\int |\mathbf{u}(\mathbf{x})-\mathbf{x}{|}^2{I}_S(\mathbf{x})d\mathbf{x},$$

$${\partial }_y{u}_x(x,y)-{\partial }_x{u}_y(x,y)=0,$$

$$\displaystyle\mathbf{s}_{1}=\begin{pmatrix}0\\ 0\\ z_{\textrm{I}}(x,y)\end{pmatrix},\ \ \ \mathbf{s}_{2}=\begin{pmatrix}u_{x}(x,y)-x\\ u_{y}(x,y)-y\\ z_{\textrm{II}}(u_{x},u_{y})-z_{\textrm{I}}(x,y)\end{pmatrix},\ \$$

$${\mathbf{s}}_3=\left(\begin{array}{c}0\\ 0\\ z_T-z_{\text{II}}(u_x,u_y)\end{array}\right),{\mathbf{s}}_4=\left(\begin{array}{c}{u}_x(x,y)-x\\ u_y(x,y)-y\\ z_T\end{array}\right).$$

$${\mathbf{n}}_{\text{I}}\cdot (\nabla \times {\mathbf{n}}_{\text{I}})=0$$

$${\mathbf{n}}_{\text{I}}=n_1{\hat{\mathbf{s}}}_1-n_2{\hat{\mathbf{s}}}_2,$$

$${\mathbf{s}}_2(\nabla \times {\mathbf{s}}_1)=n_1\frac{\{{\mathbf{s}}_2\times [({\mathbf{s}}_2\nabla ){\mathbf{s}}_2]{\}}_3}{\mathbf{n}\cdot {\mathbf{s}}_2}-{\mathbf{s}}_2(\nabla \times {\mathbf{s}}_3)+{\mathbf{s}}_2(\nabla \times {\mathbf{s}}_4)$$

$$\mathbf{v}\nabla z_{\text{I}}(x,y)=n_1\cdot \frac{\mathbf{v}\cdot \left[\left({\mathbf{v}}^{\perp }\cdot \nabla \right){\mathbf{v}}^{\perp }\right]}{\mathbf{n}\cdot {\mathbf{s}}_2}+\mathbf{v}\nabla z_{\text{II}}(u_x,u_y)-(z_{\text{II}}(u_x,u_y)-z_{\text{I}}(x,y))\nabla \mathbf{v}.$$

$$\nabla (z-z_{\text{I}}(x,y))\stackrel{!}{=}\frac{{\mathbf{n}}_{\text{I}}(x,y)}{({\mathbf{n}}_{\text{I}}(x,y){)}_z}\iff \left(\begin{array}{c}{\partial }_x{z}_{\text{I}}(x,y)\\ {\partial }_y{z}_{\text{I}}(x,y)\end{array}\right)\stackrel{!}{=}\left(\begin{array}{c}\frac{n_2\cdot (u_x-x)}{|{\mathbf{s}}_2|\cdot ({\mathbf{n}}_{\text{I}}{)}_z}\\ \frac{n_2\cdot (u_y-y)}{|{\mathbf{s}}_2|\cdot ({\mathbf{n}}_{\text{I}}{)}_z}\end{array}\right)\propto {\mathbf{v}}^{\perp },$$

$$\mathbf{v}\nabla z_{\text{I}}(x,y)=0,$$

$$n_1\cdot \frac{\mathbf{v}\cdot \left[\left({\mathbf{v}}^{\perp }\cdot \nabla \right){\mathbf{v}}^{\perp }\right]}{\mathbf{n}\cdot {\mathbf{s}}_2}+n_2\frac{g(u_x,u_y)}{({\mathbf{n}}_{\text{II}}{)}_z\cdot |{\mathbf{s}}_2|}-(z_{\text{II}}(u_x,u_y)-z_{\text{I}}(x,y))\nabla \mathbf{v}=0,$$

$$n_1|{\mathbf{s}}_1|+n_2|{\mathbf{s}}_2|+n_1|{\mathbf{s}}_3|\stackrel{!}{=}const\equiv OPL.$$

$$z_{\text{II}}(u_x,u_y)-z_{\text{I}}(x,y)=-\frac{n}{n^2-1}OP{L}_{red}\mp \frac{1}{n^2-1}\sqrt{OP{L}_{red}^2+(n^2-1)|\mathbf{u}(x,y)-\mathbf{x}{|}^2}$$

$$-\frac{(v_x^2+v_y^2)}{{\mathbf{n}}_{\text{I}}\cdot {\mathbf{s}}_2}\nabla \mathbf{v}-(z_{\text{II}}(u_x,u_y)-z_{\text{I}}(x,y))\nabla \mathbf{v}\stackrel{!}{=}0,$$

$$det(\nabla \mathbf{u}(x,y))I_T(\mathbf{u}(x,y))=I_S(x,y),$$

$${\partial }_y{u}_x-{\partial }_x{u}_y-\frac{n^2-1}{OP{L}_{red}^2+(n^2-1)|\mathbf{u}(x,y)-Id|}$$

$$\cdot [(u_x-x{)}^2{\partial }_y{u}_x-(u_y-y{)}^2{\partial }_x{u}_y+(u_x-x)(u_y-y)({\partial }_y{u}_y-{\partial }_x{u}_x)]=0,$$

$$\mathbf{u}(x,y)={\mathbf{u}}^{\infty }(x,y)+\Delta \mathbf{u}(x,y)$$

$$\Delta u_x(-0.5,y)=\Delta u_x(0.5,y)\stackrel{!}{=}0,y\in [-0.5,0.5],$$

$$\Delta u_y(x,-0.5)=\Delta u_y(x,0.5)\stackrel{!}{=}0,x\in [-0.5,0.5],$$

$${\partial }_x(\Delta u_x)\to \frac{1}{2\Delta x}[(\Delta u_x{)}_{i;j+1}-(\Delta u_x{)}_{i;j-1}],{\partial }_y(\Delta u_x)\to \frac{1}{2\Delta y}[(\Delta u_x{)}_{i+1;j}-(\Delta u_x{)}_{i-1;j}]$$

$${\partial }_x(\Delta u_x)\to -\frac{1}{2\Delta x}[3(\Delta u_x{)}_{i;j+2}-4(\Delta u_x{)}_{i;j+1}+(\Delta u_x{)}_{i;j}],$$

$${\partial }_y(\Delta u_x)\to -\frac{1}{2\Delta y}[3(\Delta u_x{)}_{i+2;j}-4(\Delta u_x{)}_{i+1;j}+(\Delta u_x{)}_{i;j}]$$

$$det(\nabla {\mathbf{u}}^{\infty }(x,y))I_T(\mathbf{u}(x,y))\stackrel{(\text{ref}eq:1)}{\to }\frac{I_S(x,y)}{I_T({\mathbf{u}}^{\infty }(x,y))}{I}_T(\mathbf{u}(x,y)).$$

$${\hat{\mathbf{s}}}_2=\frac{n_1}{n_2}{\hat{\mathbf{s}}}_1+\left\{-\frac{n_1}{n_2}\cdot \hat{\mathbf{n}}\cdot {\hat{\mathbf{s}}}_1+\sqrt{1-{\left(\frac{n_1}{n_2}\right)}^2[1-(\hat{\mathbf{n}}\cdot {\hat{\mathbf{s}}}_1{)}^2]}\right\}\hat{\mathbf{n}}$$

$$|{\mathbf{s}}_2|=OP{L}_{red}+\frac{n_1}{n_2}[z_{\text{II}}(x_m,y_m)-z_{\text{I}}(x_s,y_s)].$$

$${\mathbf{s}}_2(x_s,y_s,z_{II}(x_m,y_m))\stackrel{!}{=}(x_m-x_s,y_m-y_s,z_{II}(x_m,y_m)-z_I(x_s,y_s))$$

Full text

***\linkablechristoph.boesel@uni-jena.de **

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. Several methods for phase and irradiance control with double freeform surfaces have been proposed in literature.

One of the first design methods for the mapping of two 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, 3].

Zhang et al. [4] and Shengqian et al. [5] solve the design problem by describing it in the form of a Monge-Ampère type 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 [5].

An alternative approach to construct freeform surfaces for irradiance and phase control is from ray mapping methods[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. These methods are based on the seperation 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 surfacesfrom 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)[6, 7, 8, 9, 10, 11]. 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 [7] and can be solved by corresponding numerical schemes[23]. 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[9].

The investigations presented here are inspired by several publications [13, 14, 15, 16, 17, 18], 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 [13, 14, 18], and that of double freeform surfaces for irradiance and phase control[15, 16, 17]. As demonstrated for illumination control with single freeform surfaces in Ref. [21] and for collimated beam shaping with double freeform surfaces in Ref. [22] the design problems are thereby restricted to a paraxial approximation.

Here we first investigate the design by ray-mapping methods of double freeform surfaces which 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 PDE’s. 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 an optimization scheme with the OMT mapping from the quadratic cost function as the intial 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 lens lens distances is presented. It is based on discretizing the PDEs with finite differences and solving the resulting system of nonlinear equations by a standard optimization 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

2.1 Energy conservation and cost functions

In Ref. [21] 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 raymapping $\mathbf{u}(x,y)=(u_{x}(x,y),u_{y}(x,y))$ is calculated from the theory of optimal mass transport, 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. Firstly, 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 Eq.

[TABLE]

Secondly, in 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[6, 7, 8, 9, 10, 11]. 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

[TABLE]

which is valid in the paraxial approximation, was studied by the authors[21]. A key property there was the vanishing curl

[TABLE]

characterizing the quadratic cost function in Eq. (2) [25].

As shown by Rubinstein and Wolansky, the cost function for collimated beam shaping with double freeform lenses takes a different form than Eq. (2)[9]. The authors propose to minimize the corresponding cost function by a steepest descent algorithm to get the ray mapping[9], but unfortunatly 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 optimization scheme presented below. It will therefore build the basis of the design approach presented in sections 3 and 4.

2.2 Ray mapping condition

We follow the approach from Ref. [21, 22] 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}$:

[TABLE]

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

[TABLE]

to both freeform surfaces $z_{\textrm{I}}(x,y)$ and $z_{\textrm{II}}(x,y)$ and their normal vector fields $\mathbf{n}_{\textrm{I}}(x,y)$ and $\mathbf{n}_{\textrm{II}}(x,y)$, respectively. Thereby the normal vector fields can be expressed with the law of refraction

[TABLE]

in terms of the normalized incoming and refracted vector fields $\hat{\mathbf{s}}_{1}$ and $\hat{\mathbf{s}}_{2}$ and the refractive indices $n_{1}$ and $n_{2}$. Hence Eq (5) can be written as

[TABLE]

and therefore, by plugging in the vector fields given in Eq. (4) and defining $\mathbf{v}:=(\mathbf{u}(x,y)-\mathbf{Id})^{\perp}=(-(u_{y}(x,y)-y),u_{x}(x,y)-x)$, we get

[TABLE]

Comparing this equation with the case of a single freeform surface[21], we see that the $\mathbf{v}\nabla z_{\textrm{II}}(u_{x},u_{y})$-term on the right hand side (RHS) arises due to the second freeform surface $z_{\textrm{II}}(x,y)$ and is therefore connected to the recollimation of the refracted vector field.

The left hand side (LHS) of Eq. (8) represents the dot product of the projected gradient of the first surface $\nabla z_{\textrm{I}}(x,y)$ =$(\partial_{x}z_{\textrm{I}}(x,y),\partial_{y}z_{\textrm{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

[TABLE]

leaving us with the condition that the RHS of Eq. (8) has to be equal to zero.

A further relation can be derived by applying the chain rule to the Eq. $\nabla_{\mathbf{u}}z_{\textrm{II}}(u_{x},u_{y})=(\partial_{u_{x}}z_{\textrm{II}}(u_{x},u_{y}),\partial_{u_{x}}z_{\textrm{II}}(u_{x},u_{y}))\stackrel{{\scriptstyle!}}{{=}}\frac{\mathbf{n_{\textrm{II}}}(u_{x},u_{y})}{(\mathbf{n_{\textrm{II}}})_{z}}$[24]. This provides us with the gradient $\nabla z_{\textrm{II}}(x,y)$, which is used to rewrite the second term of the RHS of Eq. (8).

Hence from the continuity condition (5) and the law of refraction (6) follow the system of Eqs.

[TABLE]

with $g(u_{x},u_{y})=-v_{x}^{2}\partial_{x}u_{y}+v_{y}^{2}\partial_{y}u_{x}+v_{x}v_{y}(\partial_{x}u_{x}-\partial_{y}u_{y})$ and similar Eqs. by considering Eq. (5) and Eq. (6) for the second freeform surface $z_{\textrm{II}}(x,y)$.

These Eq. (10) together with Eq. (1) build a system of PDEs for the unknown mapping $\mathbf{u}(x,y)$ and the surfaces $z_{\textrm{I}}(x,y)$ and $z_{\textrm{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 (10b) which 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_{\textrm{II}}(u_{x},u_{y})-z_{\textrm{I}}(x,y))\gg|\mathbf{u}(x,y)-\mathbf{Id}|$ leading to a vanishing first and second term in (10b). Addionally, using the mapping from the quadratic cost function defined by Eq. (3), the condition (10b) is fulfilled and the surfaces can be calculated by using Eq. (10a) with appropriate boundary conditions. The boundary conditions can be derived by considering the law of refraction (9) on the boundaries of $I_{S}(x,y)$ and $I_{T}(x,y)$ [21]. As discussed in Ref. [21] , this leads to a path independent integration of Eq. (9) to calculate the surface $z_{\textrm{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 [9] for small angles.

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

[TABLE]

By plugging in Eq. (4), the Eq. (11) can be solved for

[TABLE]

with $n:=n_{1}/n_{2}$ and $OPL_{red}:=(OPL-n_{1}z_{T})/n_{2}$ between the first and second surface. The sign in Eq. (12) depends on whether we have a single-lens ($OPL_{red}>0$; $n<1$: negative sign) or a two-lens system ($OPL_{red}<0$; $n>1$: positive sign). According to Eq. (12), for single-lens systems the mapping values are restricted by the relation $|\mathbf{u}(x,y)-\mathbf{x}|^{2}<OPL_{red}^{2}/|n^{2}-1|$.

Since $\mathbf{n}\cdot\mathbf{s}_{2}$ and $(\mathbf{n_{\textrm{II}}})_{z}\cdot|\mathbf{s}_{2}|$ depend on $z_{\textrm{II}}(u_{x},u_{y})-z_{\textrm{I}}(x,y)$, Eq. (10b) can be written as a PDE for the components of mapping function $\mathbf{u}(x,y)$ only. Therefore, the Eqs. (1) and (10b) build a system of PDEs for the functions $u_{x}(x,y)$ and $u_{y}(x,y)$. Both Eqs. 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 (10b) briefly for freeform mirrors.

2.2.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}_{\textrm{I}}\mathbf{s}_{2}=-(\mathbf{n}_{\textrm{II}})|\mathbf{s}_{2}|$. Therefore, Eq. (10b) reduces to

[TABLE]

Hence the two-refractor problem with collimated beams can be solved if the mapping fulfills $\nabla\mathbf{v}=\partial_{y}u_{x}-\partial_{x}u_{y}=0$, which is the case for the quadratic cost function defined by Eq. (1) and Eq. (3). This was proven in a mathematically rigorous way by Glimm and Oliker [7].

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, the Eqs. (1) and (10b) are the basis of the design approach presented in the following. Since Eq. (10b) is exactly fulfilled by the mapping with the quadratic cost functions for an infinite distance between the freeform surfaces or an infinite OPL, respectively, the Eq. (10b) represents a correction to Eq. (3) for finite OPLs. Hence, for finite distances between $z_{\textrm{I}}(x,y)$ and $z_{\textrm{II}}(x,y)$ we are searching for corrections $\Delta\mathbf{u}(x,y)=(\Delta u_{x}(x,y),\Delta u_{y}(x,y))$ with $\Delta\mathbf{u}(x,y)\stackrel{{\scriptstyle OPL\rightarrow\infty}}{{\rightarrow}}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. (10b) in terms of Eq. (3) plus a perturbation term, we want to solve the system of equations

[TABLE]

with

[TABLE]

and the given function $\mathbf{u}^{\infty}(x,y)$ for the correction $\Delta\mathbf{u}(x,y)$. The scalability of the mapping correction $\Delta\mathbf{u}(x,y)$ with the parameter $OPL_{red}$ thereby suggests to solve Eq. (14) by an optimization method, since $OPL_{red}$ can be reduced step by step with the solution $\Delta\mathbf{u}(x,y)$ from the previous step as an initial iterate. This ensures the convergence of the optimization and can be done until the desired design goal or $OPL_{red}$, respectively, is reached.

Additionally, we want to apply boundary conditions to Eq. (14). 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, like in the example section 5, we therefore apply

[TABLE]

implying that the edges of the boundary of $I_{S}(x,y)$ are mapped onto the opposing edges of the boundary of $I_{T}(x,y)$.

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

[TABLE]

for the inner points and

[TABLE]

on the boundary. A system of $2\cdot N^{2}$ nonlinear equations for the unknows $(\Delta u_{x})_{i;j}$ and $(\Delta u_{y})_{i;j}$ is left, which can be solved by standard numerical methods. In the next section, we give an overview of the design algorithm and a detailed discussion of the implementation.

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.

Calculate the optimal mass transport map with quadratic cost function $\mathbf{u}^{\infty}(x,y)$ between the distributions $I_{S}(x,y)$ and $I_{T}(x,y)$. 2. 2.

Discretize the system of Eq. (14) with e.g. finite differences (17) and (18) and apply boundary conditions. 3. 3.

Choose the physical paramters ($n_{1},n_{2},OPL_{red}$) of the system and solve the system of nonlinear Eqs. with a predefined tolerance and the inital iterate $(\Delta u_{x})_{i;j}=0$ and $(\Delta u_{y})_{i;j}=0$, with $i;j=1,...,N$. 4. 4.

Calculate the freeform surfaces from Eqs. (9) and (12).

There are numerous publications presenting numerical methods for the calculation of optimal mass transport 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.[26]. It provides an efficient calculation of $\mathbf{u}^{\infty}(x,y)$, but has some drawbacks like 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 section 5 we therefore choose a background irradiance $I_{T}(x,y)>0$ to ensure the convergence.

For solving the system of nonlinear Eqs. in step $3$ we have used the fsolve() function from MATLABs 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 Eq. system even for a large number of variables. The scalability with the parameter $OPL_{red}$ of the distance of the initial iterate to the solution of Eq. (14) suggests that the optimization could be accelerated by using e.g. the Newton algorithm.

Solving the Eq. system in the form (14), we recognized oscillations in the solution $\Delta\mathbf{u}(\mathbf{x})$. These seem to arise due to the irradiance optimization by the Jacobian Eq. (14a), which in contrast to Eq. (14b), is already fulfilled to a high degree by the initial iterate $\Delta\mathbf{u}(\mathbf{x})=0$. Hence, we replace the $det(\nabla\mathbf{u}^{\infty}(x,y))I_{T}(\mathbf{u}(x,y))$-term in (14a) by

[TABLE]

This replacement contains the wrong assumption that the Jacobian Eq. is exactly fulfilled by $\mathbf{u}^{\infty}(\mathbf{x})$, but in our experience this leads to only minor changes in the local energy preserving property of the mapping by the optimization. The optimized map will therefore solve the Jacobian Eq. to (nearly) the same degree as $\mathbf{u}^{\infty}(\mathbf{x})$, but also Eq. (14b) to a high degree. Alternatively, instead of using (19), one could also redefine $I_{T}(x,y)$ with (1) and $\mathbf{u}^{\infty}(\mathbf{x})$ to optimize solely for Eq. (14b).

To calculate the first freeform we use Eq. (9) with Eq. (12), 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. (14b) as mentioned in section 2. Since Eq. (14b) 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. (12) 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 by applying the ray tracing formula

[TABLE]

and Eq. (11) in the form

[TABLE]

Hereby $(x_{s},y_{s})$ are the unknown initial coordinates of the incoming vector $\mathbf{\hat{s}}_{1}(x_{s},y_{s})$ and $(x_{m},y_{m})$ the predefined target coordinates. By interpolating the surface $z_{I}(x,y)$, which we will do bilinearly, the Eq. (20) gives the normalized refracted vector field $\mathbf{\hat{s}}_{2}(x,y)$. To calculate $z_{II}(x,y)$ on the desired grid point $(x_{m},y_{m})$ we use the refracted vector field $\mathbf{s}_{2}(x,y,z_{II}):=|\mathbf{s}_{2}|\cdot\mathbf{\hat{s}}_{2}(x,y)$ and solve the system

[TABLE]

at every grid point $(x_{m},y_{m})$ for the unknown values $x_{s}$, $y_{s}$ and $z_{II}(x_{m},y_{m})$.

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 grey levels. For the resolution of both images $250\times 250$ pixels are chosen.

Since input and output irradiances are defined on unit squares of the same size, the freeform surfaces have the shape of squares with the side length $1a.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 optical path length $OPL_{red}=-0.2a.u.$ are used. Therefore, a two-lens system like in Fig. 1b) is considered. An distance between both freeform surfaces of approximately $0.4a.u.$ can be estimated from Eq. (12) since, due to symmetry reasons, the corner points of the irradiances are mapped onto each other.

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

We want to point out that it would be possible to design a system with crossing rays for the presented examples. This can be done by mirroring $\mathbf{u}^{\infty}(x,y)$ at the point of origin, which preserves the property (3), and using the mirror map for the optimization and according boundary conditions. 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 the Eq. (3). 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)$.

The calculation of the initial map $\mathbf{u}^{\infty}(x,y)$ with our implementation of Sulmans algorithm took $188sec$ for “IAP” and $575sec$ for “house” on an Intel Core i3 at $2\times 2.4$Ghz with $16$GB RAM.

After fixing the physical properties, the mapping has to be optimized according to sections 3 and 4 to obtain the solution of the Eq. system (14). To do that, one could use a starting value smaller/larger than the design goal $OPL_{red}=-0.2a.u.$ or $OPL_{red}=0.2a.u.$, respectively, to ensure the convergence to the solution of (14). However, in practice we did not experience any significant benefit from using smaller/larger starting values for $OPL_{red}$ than the design goal. For the tolerance of MATLABs fsolve(), the value of $10^{-6}$ was used. For both examples the optimization processes with $125,000$ design variables took about $67sec$ for “IAP” and $71sec$ for “house”.

In Fig. 3 (“IAP”) and Fig. 4 (“house”) the Eqs. (14) are plotted with the initial map $\mathbf{u}^{\infty}(x,y)$ and the optimized map $\mathbf{u}(x,y)$, respectively. Whereas the solutions of the Jacobian Eq. (14a), according to the calculated rms values, remains nearly the same, the solution of the mapping Eq. (14b) improves drastically in both cases.

The construction of the surfaces is done by the integration of Eqs. (9) and (12), first from $(x,y)=(-0.5,-0.5)$ along the x-direction and then along the y-direction. 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\cdot 10^{6}$ rays can be seen in Fig. 6 (“IAP”) and Fig. 7 (“house”). The rms values of the normalized difference $\Delta I_{T}(x,y):=(I_{T}(x,y)-I_{T,RT}(x,y))/I_{T}(x,y)$ between predefined distribution $I_{T}(x,y)$ from Fig. 2 and the ray tracing illumination patterns $I_{T,RT}(x,y)$ as well as the rms values of the optical path difference show a significant improvement of quality of the illumination patterns and the wavefronts for both examples.

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 beyond the diffraction limit. Remaining deviations from the ideal wavefront result (besides fundamental numerical limitiations) are from Eq. (14b), which is not fulfilled exactly and therefore leads to an error accumulation along the integration path of Eqs. (9) and (12). 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)$ are resulting from steep gradients, which can be seen especially for the example “IAP” in Fig. 6d). This is in agreement with Fig. 3b) and observations which were made in Ref. [21] .

6 Conclusion

A design method for the calculation of compact continuous double freeform surfaces for collimated beam shaping with complex irradiance patterns was presented. The method is based on the ray mapping condition (14b), which was derived from the law of refraction and the surface continuity condition in section 2 and builds together with the Jacobian Eq. (14a) a system of nonlinear PDEs for the unknown ray mapping $\mathbf{u}(x,y)$.

Due to the satisfaction of Eq. (14) for infinite lens distances by the mapping from OMT with the quadratic cost function (2), this mapping serves as an ideal initial iterate for an optimization approach for solving the system of Eq. (14). As it was shown by approximating the Eqs. (14) by finite differences and using a standard optimization scheme from MATLAB’s optimization toolbox one can ensure a fast convergence to the solution of the Eqs. (14). This was demonstrated by applying the presented method to two 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. [27]) and the generalization of Eq. (16) to more complex boundary shapes. The scalability of the distance of the initial map $\mathbf{u}^{\infty}(x,y)$ by $OPL_{red}$ to the solution of the Eqs. (14) suggests the application of e.g. the Newton algorithm for a faster optimization.

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

Acknowledgements.

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. We also want to acknowledge the Federal Ministry of Education and Research Germany for financial support through the project fo+ (FKZ: O3WKCK1D) and KoSimO (FKZ:031PT609X).