Freeform Optics: Optimal Transport, Minkowski Method, and Monge–Ampère-Type Equations

1. Description and Background

A freeform optical surface, simply stated, refers to a surface whose shape lacks rotational symmetry. The use of such surfaces allows generation of complex, compact, and highly efficient imaging systems. Ever since lenses without symmetries were used in World War I in periscopes, the engineering and design of freeform optical surfaces have gone through remarkable evolution, with applications in a wide range of areas, including medical devices, clean energy technology, military surveillance equipments, mobile displays, remote sensing, and several other areas of imaging and nonimaging optics that can benefit from distributing light from a source to a target in a controlled fashion using spatially and energy-efficient systems. See 17 and the references therein.

Mathematically, the design of freeform optical surfaces is an inverse problem related to optimal transportation theory and leads to a class of nonlinear partial differential equations (PDE) called generated Jacobian equations for which the Monge–Ampère equation is a prototype. The modeling of the problem is based on the systematic application of the laws of reflection and refraction in geometric optics along with energy conservation principles.

For completeness we recap the laws of geometric optics. Consider a reflecting or refracting surface $\mathcal{S}$ and suppose a unit vector in the direction of $x$ is incident on the surface $\mathcal{S}$ at a point where the normal to $\mathcal{S}$ is given by the unit vector $\nu$. Reflection and refraction happen in tandem. (Fig. 1). However, for simplicity of the models we treat them separately. If $r$ is a unit vector in the direction of the reflected ray, the law of reflection (angle of incidence is equal to angle of reflection) in vector form says:

$$\begin{equation} r = x - 2\langle x, \nu \rangle \nu . {\tag{1.1}} \end{equation}$$

If medium $I$ and medium $II$ are two homogeneous isoptropic media with respective refractive indices $n_1$ and $n_2$ and $m$ is the direction of the refracted ray into medium $II$, the law of refraction (Snell’s law), which states $\dfrac{\sin \theta _i}{\sin \theta _t} = \dfrac{n_2}{n_1}$ where $\theta _i$ is the angle between $x$ and $\nu$ (angle of incidence) and $\theta _t$ is the angle between $m$ and $\nu$ (angle of refraction) in vector form, states that:

$$\begin{equation} m = \dfrac{1}{\kappa }(x - \lambda \nu ) {\tag{1.2}} \end{equation}$$

where $\kappa = n_2/n_1$ is the relative refractive index, $\lambda = \lambda (\kappa , x, \nu )$ and $\nu$ points into medium $II$.

Figure 1.

Laws of reflection and refraction.

When medium I is optically denser than medium II or equivalently $\kappa <1$, the refracted ray bends away from the normal. As a result, there is a critical value $\theta _c$ of the angle of incidence for which there is no value of $\theta _t$ unless $\sin \theta _i \leq \sin \theta _c = \kappa$. Thus for $\kappa < 1$ we impose the condition that

$$\begin{equation} x\cdot m \geq \kappa . {\tag{1.3}} \end{equation}$$

Similar physical restriction applies for $\kappa > 1$.

By using equation 1.1 or 1.2 we can define a map $T_{\mathcal{S}}: \Omega \to \Omega ^*$ that sends points in $\Omega$ (incident rays) to points in $\Omega ^*$ (reflected rays if we use 1.1 or refracted rays if we use 1.2) according to laws of geometric optics.

Generally, in an optical surface design problem the data consists of two bounded domains $\Omega$, $\Omega ^*$ contained in $\mathbb{R}^n$ and two probability measures $\rho _1, \rho _2$ defined on $\Omega$ and $\Omega ^*$ respectively, satisfying the energy conservation statement

$$\begin{equation} \rho _1(\Omega ) = \rho _2(\Omega ^*). \tag{1.4} \end{equation}$$

$\Omega$ and $\Omega ^*$ correspond to the directions of incident light and the directions of light after reflection or refraction respectively. Likewise $\rho _1$ and $\rho _2$ are the intensity of light from the source $\Omega$ and a desired illumination intensity on the target $\Omega ^*$, respectively. We regard $(\Omega , \rho _1)$ as a source pair and $(\Omega ^*, \rho _2)$ as a target pair.

The objective is then to find an optical surface (lens/mirror) $\mathcal{S}$ such that

$$\begin{equation} \rho _1(T_{\mathcal{S}}^{-1} (F)) =\rho _2(F) {\tag{1.5}} \end{equation}$$

for all $F \subset \Omega ^*$.

In practice, $\mathcal{S}$ is a reflecting mirror or refracting lens and for $m \in \Omega ^*$, $T_{\mathcal{S}}^{-1}(m)$ represents the set of all directions from the source that will, after reflection or refraction off of the lens $\mathcal{S}$, propagate in the direction of $m$. We use the notation $\mathcal{T}_{\mathcal{S}}$ for $T_{\mathcal{S}}^{-1}$. The map $\mathcal{T}_{\mathcal{S}}$ is called the ray-tracing map of $\mathcal{S}$.

The physical meaning of equation 1.5 is that the prescribed intensity of reflected or refracted light on the set $F \subset \Omega ^*$, which is $\rho _2(F)$, is equal to the intensity of light that comes from the source, which is $\rho _1(\mathcal{T}_{\mathcal{S}} (F))$. Depending on whether $\mathcal{S}$ is a reflector or a refractor, the design problem is called a reflector/refractor problem.

The models used to obtain numerical and analytical solutions for reflector/refractor problems are developed by using one of the following three approaches. One of the approaches is to use variational method where the problem is cast as an optimal mass transportation problem from $\rho _1$ to $\rho _2$ for an appropriate cost function and use Kantorovich duality, 20. Another method is to derive the second order nonlinear PDE of Monge–Ampère-type satisfied by the surface defining function from the energy conservation relation 1.5 between the light energy emitted by the source and the light energy received by the target and analyze the PDE, 11. A third approach is to exploit mainly the inherent geometric features of the problem and use a classical approach of Minkowski and Aleksandrov which was used in solving the Minkowski problem 16 of finding a closed convex surface whose Gaussian curvature is a given function of the exterior unit normal.

In the discussion below, we will focus on refractor problems and briefly describe the aforementioned approaches by using prototype refractor problems. We will state the problems explicitly and show how the available methods are used to solve the corresponding design problems. We will also mention problems that are worth studying in the future. We point out to the reader that unless mentioned otherwise, we consider $\kappa <1$.

2. Optimal Transport: The Far Field Refractor Problem

Several freeform lens design problems have been successfully modeled by using optimal transport framework on the sphere. See 41020. The models have also been adopted by the optics community. Among these problems is the far field refractor problem.

2.1. Statement of the far field refractor problem with point source

In this problem, we are given two domains $\Omega , \Omega ^*$ contained in the unit sphere $S^{n-1}$ of $\mathbb{R}^n$ (in practice $n=3$), $|\partial \Omega | =0$; a punctual source of light located at the origin $\mathcal{O}$, surrounded by media I, from which a monochromatic ray of light is issued in each direction $x \in \Omega$ with intensity density function given by $g(x)$ for $g \in L^1(\Omega ), g > 0$ a.e. on $\Omega$ and a prescribed intensity density distribution $f \in L^1(\Omega ^*)$ which satisfies

$$\begin{equation} \int _{\Omega } g(x) dx = \int _{\Omega ^*} f(x) dx . {\tag{2.1}} \end{equation}$$

The objective is to find a refracting surface parametrized as $\mathcal{S} = \{\rho (x)x : x \in \Omega \}, \rho >0$, between medium $I$ and medium $II$ such that all rays emitted from the point $\mathcal{O}$ with directions $x \in \Omega$ are transported via refraction by the surface into media $II$ with directions in $\Omega ^*$ in such a way that the prescribed illumination pattern $f \in L^1(\Omega ^*)$ is achieved on $\Omega ^*$.

In particular, if we let $d\rho _1 = g dx$ and $d\rho _2 = f dx$ where $dx$ is the surface measure on $S^{n-1}$, this means $\rho _1$ and $\rho _2$ satisfy the mass balance condition $\rho _1 (\Omega ) = \rho _2 (\Omega ^*)$ and a weak solution to the far field refractor problem is defined as an interface $\mathcal{S}$ so that the associated map $T_{\mathcal{S}}$ satisfies

$$\begin{equation} T_{\mathcal{S}}(\Omega ) \subset \Omega ^* \qquad \text{and} \qquad \rho _2 = (T_{\mathcal{S}})_\#\rho _1. {\tag{2.2}} \end{equation}$$

We recall for $F \subset \Omega ^*$ that the measure $(T_{\mathcal{S}})_\#\rho _1 (F) \coloneq \rho _1(T_{\mathcal{S}}^{-1}(F))$ is the push forward measure. Because of 1.3, we assume $x \cdot m \geq \kappa$ for $x \in \Omega$ and $m \in \Omega ^*$.

By using Snell’s law 1.2 it can be shown that the semi-ellipsoid of revolution (ellipse in 2D case) with one focus at $\mathcal{O}$, given in polar representation as (10)

$$\begin{equation*} E(m,b) = \left\{\dfrac{b}{1-\kappa m \cdot x}x: \ x \in S^{n-1},\ m \cdot x \geq \kappa \right\} \end{equation*}$$

for some $b>0$ and $m \in \Omega ^*$, has a uniform refracting property. This means that if $E(m,b)$ is an interface between media $I$ and $II$, all rays issued from $\mathcal{O}$ in a direction $x$ with $x\cdot m\geq \kappa$ will be refracted in the direction $m$ after hitting $E(m,b)$. Motivated by this, a solution for the far field refractor problem is sought among surfaces $\mathcal{S} = \{\rho (x)x : x \in \Omega \}$ which are supported from outside by semi-ellipsoids at each point. That is, for each point $\rho (x_o)x_o \in \mathcal{S}$ there exists a semi-ellipsoid $E(m,b)$ with $m \in \overline{\Omega ^*}$ such that $\rho (x) \leq \dfrac{b}{1-\kappa m \cdot x}$ with equality at $x_o$. We now discuss how optimal transport technique is used to find such a solution. First, a brief review of optimal transport.

2.2. Review of optimal transport problems

Let $X$ and $Y$ be compact subsets of $\mathbb{R}^n$ and $c: X \times Y \to [0, \infty ]$. We refer to $c$ as the cost function. Denote by $\mathcal{P}(X)$ and $\mathcal{P}(Y)$ the set of probability measures on $X$ and $Y$ respectively. Given $\mu \in \mathcal{P}(X)$ and $\nu \in \mathcal{P}(Y)$ a measurable map $T: X \to Y$ satisfying $T_{\#} \mu = \nu$ is called a transport map from $\mu$ to $\nu$. Monge’s formulation of the optimal transport problem is to minimize

$$\begin{equation} T \to \int _{X} c(x, T(x)) d\mu (x) {\tag{2.3}} \end{equation}$$

among all transport maps from $\mu$ to $\nu$. A minimizer of 2.3 is called an optimal transport map. A generalization due to Kantorovich 18 of Monge’s problem is stated as

$$\begin{equation} \inf _{\gamma \in \Gamma (\mu , \nu )} \int _{X \times Y} c(x, y) d \gamma (x, y) {\tag{2.4}} \end{equation}$$

among all transport plans $\gamma \in \Gamma (\mu , \nu )$ where

$$\begin{equation*} \Gamma (\mu , \nu ) : = \{\gamma \in \mathcal{P}(X \times Y) : (\pi _1)_{\#} \gamma = \mu , (\pi _2)_{\#} \gamma = \nu \}. \end{equation*}$$

Here $\pi _1: X \times Y \to X$ and $\pi _2: X \times Y \to Y$ are projection maps. If $c$ is lower semicontinuous and bounded from below, problem 2.4 has a minimizer in $\Gamma (\mu , \nu )$. A minimizer of 2.4 is called an optimal plan.

The $c$-transform of a given function $\phi : X \to [-\infty , \infty )$ is defined as $\phi ^c : Y \to \mathbb{R}$ where

$$\begin{equation*} \phi ^c(y) = \inf _{x \in X} c(x,y) - \phi (x). \end{equation*}$$

Also $\zeta ^{\bar{c}}: X \to \mathbb{R}$, the $\bar{c}$-transform of $\zeta :Y \to [-\infty , \infty )$, is defined by

$$\begin{equation*} \zeta ^{\bar{c}}(x) = \inf _{y \in Y} c(x,y) - \zeta (y). \end{equation*}$$

A function $\phi \in C(X)$ is called $c$-concave if $\phi = \zeta ^{\bar{c}}$ for some $\zeta$. We also define for any $x \in X$ the $c$-superdifferential $\phi$ at $x$ as

$$\begin{equation*} \partial _c \phi (x) = \{y \in Y: \phi (z) \leq c(z,y) - c(x,y) + \phi (x), \forall z \in X\}. \end{equation*}$$

Problem 2.4 is a linear optimization problem with convex constraints and thus, it admits a dual problem of maximizing

$$\begin{equation} \left\{ \int _X \phi (x) \, d\mu + \int _Y \psi (y) \, d\nu : \phi (x) + \psi (y) \leq c(x,y) \right\} {\tag{2.5}} \end{equation}$$

over $\phi \in L^1(X)$ and $\psi \in L^1(Y)$ 18.

One of the basic results in optimal transport theory is that if $c$ is continuous and $X$ and $Y$ are compact, then the dual problem 2.5 has a maximizer of the form $(\phi , \phi ^c)$ where the map $\phi$, which is referred to as Kantorovich potential, is a $c$-concave function. Furthermore, if $\partial _c \phi (x)$ is single valued for $\mu -a.e.\ x$, then $T: = \partial _c \phi$ induces the optimal plan $\gamma$ for 2.4. That is, $\gamma = (Id, T)_\# \mu$ minimizes 2.4. Consequently, $\partial _c \phi (x)$ is an optimal transport map from $\mu$ to $\nu$.

2.3. The far field refractor problem as optimal transport

To transform the far field refractor problem 2.2 to an optimal transport problem, we introduce the logarithmic cost function $c: \Omega \times \Omega ^* \to [0, \infty ]$ given by

$$\begin{equation} c(x,m) = -\log (1 - \kappa x \cdot m). {\tag{2.6}} \end{equation}$$

$c(x,m)$ is continuous and bounded from below since $x\cdot m \geq \kappa$. Thus, if we consider the dual problem 2.5 corresponding to $c(x,m)$ given as

$$\begin{equation} \sup _{u \in C(\bar{\Omega }), v \in C(\bar{\Omega }^*)} \left\{ \int _{\Omega } u(x) d\rho _1 + \int _{\Omega ^*} v(m) d\rho _2 : u(x) + v(m) \leq c(x,m)\right\}, {\tag{2.7}} \end{equation}$$

there exists 10 a $c$-concave $\phi$ such that $(\phi , \phi ^c) \in C(\bar{\Omega }) \times C(\bar{\Omega }^*)$ is a maximizer of 2.7. Let $\phi \in C(\Omega )$ be $c$-concave for the cost function $c(x,m) = -\log (1 - \kappa x \cdot m)$. If we set

$$\begin{equation} \rho (x) = e^{\phi (x)}, {\tag{2.8}} \end{equation}$$

then for each $x_o \in \Omega$ there exists $b>0$ and $m \in \Omega ^*$ such that

$$\begin{equation} e^{\phi (x)} \leq \dfrac{b}{1-\kappa m \cdot x} \,\, \forall x \in \Omega \,\, \text{ and } \,\, e^{\phi (x_o)} = \dfrac{b}{1-\kappa m \cdot x_o}. {\tag{2.9}} \end{equation}$$

That is, at each point $\rho (x_o) x_o$ the surface $\mathcal{S} = \{ \rho (x) x: x \in \Omega \}$ is supported from outside by a semi-ellipsoid $E(m,b)$ for some $b>0$ and $m \in \Omega ^*$.

If $\mathcal{S}$ is smooth, then $\mathcal{S}$ and $E(m,b)$ have the same normals at $\rho (x_o)x_o$. By the refraction property of $E(m,b)$ a ray emanating from $\mathcal{O}$ in a direction $x_o$ will be refracted off $\mathcal{S}$ in the direction $m$ by the surface $\mathcal{S} = \{x e^{\phi (x)}: x \in \Omega \}$. Therefore, $T_{\mathcal{S}}(x_o) = m$. Thus, $\mathcal{S}$ satisfies, $T_{\mathcal{S}}(\Omega ) \subset \Omega ^*$. Also, from 2.9 we get

$$\begin{equation*} \phi (x) \leq -\log (1 - \kappa m \cdot x) + \log (1 - \kappa m \cdot x_o) + \phi (x_o) \qquad \text{for all} \qquad x \in \Omega \end{equation*}$$

and hence

$$\begin{equation*} T_{\mathcal{S}}(x_o) = m \in \partial _c \phi (x_o). \end{equation*}$$

Furthermore if $c$ is given by 2.6, for any $c$-concave $\phi$ the $c$-superdifferential mapping $\partial _c \phi$ is single valued for $\rho _1-a. e.\ x$. Thus if $(\phi , \phi ^c)$ is a maximizer of 2.7, $\partial _c \phi$ is a measure preserving map from $\rho _1$ to $\rho _2$. Therefore $\partial _c \phi = T_{\mathcal{S}} \rho _1$-a.e., and more importantly $\rho _2 = (T_{\mathcal{S}})_\#\rho _1$ proving that $\mathcal{S} = \{x e^{\phi (x)}: x \in \Omega \}$ will be a solution of the far field refractor problem. From this we deduce that solving the far field refractor problem corresponds to finding a maximizer of the dual problem to the optimal transport problem 2.7. It is further shown in 10 that if $\mathcal{S} = \{x \rho (x): x \in \Omega \}$ and $\tilde{\mathcal{S}} = \{x \tilde{\rho }(x): x \in \Omega \}$ are two solutions of the far field refractor problem, then $\rho = C \tilde{\rho }$ for some constant C.

The optimal transport approach can also be used in the case of anisotropic materials where optical properties vary according to the direction of propagation of light 6. In this type of media modeling the design problem is more complicated due to issues such as birefringence, a situation where the incident rays may be refracted into two rays. It is also used to develop numerical methods to approximate solutions 4 to freeform design problems. Recently, a more efficient method called entropic regularization is applied to the the reflector problem 2 and it is likely that this approach can be applied to refractor problems.

Finally, it is worth noting that the far field problem has not been studied in its full generality using optimal transport method. In practice, if a ray of light hits an interface between two media which have different optical properties, the ray is partly reflected back and partly transmitted (Fig. 1). It is an open problem whether or not this energy imbalance between the source and target can be modeled as a variational problem by using the optimal transport approach.

3. Minkowski Approach: The Near Field Refractor Problem

The Minkowski problem involves finding a closed convex surface whose Gaussian curvature at a point with exterior normal $\xi$ is given by a continuous positive function $k(\xi )$ 16. To solve this, the approach was first to define a set function on $S^{n-1}$ by $\sigma (M) = \int _M \dfrac{d\omega (\xi )}{k(\xi )}$ where $d\omega$ is a surface element on the unit sphere, and consider the general problem of finding a convex hypersurface for which $\sigma$ is a surface area function. That is, to find a convex hypersurface $F$ such that for a given subset $M'$ of $F$, the area of $M'$ is equal to $\sigma (M)$ where $M' = \{x \in F: \text{the normal to $F$ at $x$ is in } M\}$. Then partition the sphere into smaller domains $g_k$ with area $\sigma _k$ and obtain an approximation $\sigma ' = \sum _{k=1}^\ell \sigma _k \delta _{\xi _k} ; \xi _k \in g_k$ to $\sigma$. Replacing $\sigma$ by $\sigma '$ gives a discrete problem of existence of convex polyhedron with faces of areas $\sigma _k$ and outward normal $\xi _k$. Let $P_\ell$ be a solution of the discrete problem. As $\ell \to \infty$ (or the diameter of $g_k$ to $0$), it is proved that the sequence $P_\ell$ of surfaces converges to a surface which solves the Minkowski problem. This classical approach of Minkowski is used to obtain iterative methods to solve not only various geometric optics problems related to both reflection 3 and refraction, but also to develop numerical methods to solve PDEs 15. The approach is also generalized to approximate solutions to semi-discrete optimal mass transport problems and generated Jacobian equations. See 1 and the reference therein. It is also one of the methods used by the optics community to design freeform optical surfaces 14. Below we exhibit how this approach is used in the near field refractor problem.

3.1. Statement of the near field refractor problem with point source

For this problem we are given $\Omega \subset S^{n-1}$, $\Omega ^* \subset \Sigma$, a hyperplane in $\mathbb{R}^n$ with $|\partial \Omega | =0$ and the origin $\mathcal{O} \notin \Omega ^*$. For each $x \in \Omega$, a light ray is issued in the direction of $x$ from a point source of light located at $\mathcal{O}$. The illumination intensity of the source is given by a density function $g \in L^1(\Omega )$, $g > 0$ a.e. A prescribed irradiance distribution on the target set $\Omega ^*$ is given by a nonnegative density function $f \in L^1(\Omega ^*)$. $f$ and $g$ satisfy the mass balance 2.1. We assume that $\Omega$ is surrounded by medium I and $\Omega ^*$ is surrounded by medium II. Additional constraints are imposed on $\Omega$ and $\Omega ^*$ to make the problem physically feasible 9.

The near field refractor problem involves finding an interface $\mathcal{S}$ between media $I$ and $II$, parametrized radially as $\mathcal{S} = \{\rho (x) x : x \in \overline{\Omega }\}$, that redirects each ray with direction $x\in \Omega$ by refraction into $\Omega ^*$ so that a prescribed irradiance distribution $f$ is obtained on $\Omega ^*$. More precisely we will require

$$\begin{equation} \int _{\mathcal{T}_{\mathcal{S}} (F)} g(x) \, d\sigma (x) = \int _F f(y) \, d\sigma (y) {\tag{3.1}} \end{equation}$$

holds for all Borel sets $F\subset \Omega ^*$.

Similar to the semi-ellipsoids in the far field problem, there are surfaces that have uniform refracting property in the near field case. The Descartes ovoid refracts all rays to a single point (Fig. 2). Let $P \in \Omega ^*$ and $\kappa |P| < b < |P|$. The Descartes ovoid with one focus at the origin $\mathcal{O}$ is given as

$$\begin{equation*} O_b(P) = \{h(x,P,b)x : x \in S^{n-1}, x\cdot P \geq b \} \end{equation*}$$

where

$$\begin{equation} h(x,P,b) = \dfrac{(b-\kappa ^2 x \cdot P)- \sqrt {(b-\kappa ^2 x \cdot P)^2 - (1-\kappa ^2) (b^2 - \kappa ^2 |P|^2)}}{1 - \kappa ^2}. {\tag{3.2}} \end{equation}$$

If the relative refractive index of the material inside the oval $\{h(x,P,b)x : x \in S^{n-1} \}$ to that of outside is $\kappa$, then by using 1.2 it can be shown that any ray emanating from the origin $\mathcal{O}$ and having direction $x\in S^{n-1}$ with $x\cdot P\geq b$ will pass through the point $P$ after refracting off of the ovoid $O_b(P)$.

Figure 2.

Refraction property of oval; $\kappa = 0.7$, $P = (10,0), b=8$.

The Descartes ovoids will be used as the building blocks of the near field refractors and we look for a solution for the near field refractor problem among surfaces which are supported from above by an ovoid at each point.

It is known that the near field refractor problem can’t be cast as an optimal mass transport problem. However, Minkowski’s approach can be used to solve the problem.

3.2. The near field refractor problem using Minkowski’s approach

Once again let $d\rho _1 = g d\sigma (x)$ and $d\rho _2 = f d\sigma (y)$. The road map is as follows. First partition $\Omega ^*$ and approximate the measure $\rho _2$ on the target by a sequence of discrete measures $\mu _{\ell }$ concentrated at finite points in $\Omega ^*$ and that converge to $\rho _2$ as $\ell \to \infty$. Then solve the problem when the target distribution is $\mu _{\ell }$. Finally let $\ell \to \infty$ to obtain a solution to the problem.

Now, partitioning $\Omega ^* = \cup _{j=1}^{\ell } \Omega ^*_j$, one can define a discrete measure on $\Omega ^*$ by $\mu _\ell = \sum _{j=1}^{\ell } f_j \delta _{P_j}$ where $P_j \in \Omega ^*_j$, $\rho _2(\Omega _j^*) = f_j$. Clearly the mass balance $\mu _{\ell }(\Omega ^*) = \Sigma _{j=1}^{\ell } f_j = \rho _2(\Omega ^*)$ is satisfied. For the discrete measure $\mu _{\ell }$ on the target (see 8), given $\epsilon > 0$, and $b_1>0$ satisfying $\kappa |P_1| < b_1 \leq \kappa |P_1|+r_0 \dfrac{(1-\kappa )^2}{1+\kappa }$, ($r_o$ that depends on the geometry of the problem), one can find through an iterative scheme a vector $\mathbf{b}_{\epsilon } = (b_1, b_2^{\epsilon }, \cdots , b_\ell ^{\epsilon }) \in \prod _{i=1}^\ell \left(\kappa |P_i|,|P_i|\right)$ such that the multifaceted refractor $\mathcal{S}(\mathbf{b}_{\epsilon })$ defined by

$$\begin{equation} \mathcal{S}(\mathbf{b}_{\epsilon }) = \{\rho (x)x : x \in \overline{\Omega }\text{ and } \rho (x) = \min _{1 \leq i \leq \ell } h(x, P_i, b_i^{\epsilon })\} {\tag{3.3}} \end{equation}$$

satisfies

$$\begin{equation} \left|\int _{\mathcal{T}_{\S (\mathbf{b}_{\epsilon })}(P_i)} g(x) \, dx - f_i\right| \leq \epsilon ,\qquad 1\leq i\leq \ell . {\tag{3.4}} \end{equation}$$

This means that the surface $\mathcal{S}(\mathbf{b}_{\epsilon })$ solves the discrete near field refractor (Fig. 3) problem with an error of $\epsilon$ for the intensity distribution. For a given error $\epsilon$, the convergence in finite number of steps of this scheme is proved by obtaining appropriate Lipschitz estimates for the refractor measure $\int _{\mathcal{T}_{\mathcal{S}(\mathbf{b})}(P_i)} g(x) \, dx$ as a function of $\mathbf{b}$.

Figure 3.

Discrete near field refractor problem.

It can be shown that as $\epsilon \to 0$, $\mathbf{b}^{\epsilon } \to \mathbf{b}_\ell : = (b_1, b_2, \cdots , b_\ell )$ and that the multifaceted refractor $\mathcal{S}(\mathbf{b}_{\ell })$ defined by

$$\begin{equation} \mathcal{S}_\ell = \{\rho _{\ell }(x)x : x \in \overline{\Omega }\text{ and } \rho _{\ell }(x) = \min _{1 \leq i \leq \ell } h(x, P_i, b_i)\} {\tag{3.5}} \end{equation}$$

is an exact solution to the near field refractor problem if $\rho _2$ is given by the discrete measure $\mu _\ell$. It should be remarked here that the iterative process used to obtain $\mathbf{b}_{\epsilon }$ can also be implemented numerically. It is used to obtain an approximate solution for the near refractor problem when the irradiance distribution on the target is given by a discrete measure 8. The convergence of this process in more generality is discussed in 8. Notice that the freeform refracting surface which solves the near field refractor problem when $\rho _2$ is discrete measure is built from segments of Descartes ovoids.

Now that for each $\ell \in \mathbb{N}$, a multifaceted near field refractor $\mathcal{S}_\ell$ parametrized by $\rho _\ell$ is obtained corresponding to the discrete measure $\mu _\ell$ on the target, the next step is to let $\ell \to \infty$ and study the convergence. As discussed in 9, $\mu _\ell$ converges weakly to the measure $d\rho _2 = f d\sigma (y)$ and by using compactness argument, up to a subsequence $\rho _\ell$ converges to $\rho$. It is then shown that the surface $\mathcal{S} = \{\rho (x) x : x \in \Omega \}$ is the desired solution to the near field refractor problem from $d\rho _1 = g d\sigma (x)$ to $d\rho _2 = f d\sigma (y)$.

The Minkowski approach, also known in literature as the supporting quadric method, has been used to treat several other beam shaping freeform design problems. The fundamental principle is similar for all the problems. One first determines the surfaces with uniform reflection/refraction property, like the semi-ellipsoids or the Descartes ovoids, along with the length parameters and then uses those surfaces as building blocks to obtain the required freeform lens. More general far field refractor problems such as the one accounting for loss of energy due to internal reflection 7 and also in anisotropic media 6 have been studied by using this technique. Whether this approach can be extended to propagation in anisotropic media for the near field refractor problem and for the parallel refractor problem discussed in the next section remains open.

4. PDE Methods: Regularity of Solutions

The existence of a solution to a refractor design problem can be shown in most cases by exploiting the physical and geometric features of the specific design problem and using the Minkowski approach. Equally important is the regularity properties of these solutions. The analysis of regularity of surfaces that are solutions to problems in geometric optics is not only of deep mathematical interest, but also has physical significance since non-smoothness of solutions has a physical interpretation of aberrations and diffraction of light. It should be noted that some regularity properties can be easily obtained from the geometry of the problems. For instance for $\kappa < 1$ the solutions for the far field refractor problem are Lipschitz continuous by the fact that they are supported by semi-ellipsoids at each point. Further analytical properties of solutions can be obtained by exploiting the relations between solutions of refractor problems and solutions of nonlinear PDEs.

For some design problems, optimal transport theory could be used to obtain the related PDE. Indeed, for appropriate assumptions on the cost function $c(x,y): X\times Y \to \mathbb{R}$, it is known that the optimal transport map $T$ minimizing Monge’s problem solves a nonlinear PDE of Monge–Ampère-type. Since, as seen in Section 2, the far field refractor problem has an optimal transport formulation a similar nonlinear PDE can be deduced. In particular, it is known that the potential $\phi (x) = \log \rho (x)$ in 2.8 satisfies the corresponding equation, 10

$$\begin{equation*} \mathrm{det}\left( \nabla _{xx}^2 \phi (x) + \nabla ^2_{xx} c(x, T(x))\right) = \left|\mathrm{det}\nabla ^2_{xy} c(x,T(x))\right| \dfrac{g(x)}{f(T(x))} \end{equation*}$$

where for a function $h$ on $S^{n-1}$, $\nabla _x h(x)$ is the tangential gradient of $h$ at $x \in S^{n-1}$. Ma–Trudinger–Wang in 13 identified a key condition which depends on the fourth order derivatives of the cost $c$, to obtain some regularity results. Using this condition in 10 it was discussed that one can’t expect a $C^1$ regularity of the solution surface for the far field refractor problem when $\kappa < 1$. On the other hand for $\kappa >1$ the solutions are smooth provided the density functions $f$ and $g$ are smooth functions which are bounded and away from zero.

Not all refractor problems can be formulated as optimal transport problems. In that case a corresponding PDE can still be obtained by using energy conservation and by computing the Jacobi determinant of the ray tracing map explicitly. For instance, consider the parallel near field refractor problem. In this problem we have a source $\Omega$ which lies in medium I, and a target $\Omega ^*$ which lies in medium II and both are bounded domains in $\mathbb{R}^n$. For brevity we assume that $\Omega \subset \{x\in \mathbb{R}^{n+1} : x_{n+1} = 0\}$ and $\Omega ^* \subset \{x\in \mathbb{R}^{n+1} : x_{n+1} = k\}$ for some $k>0$. From each $x \in \Omega$, a ray of light is issued in a direction parallel to $e_{n+1}$—the unit vector in the direction of $x_{n+1}$ axis in $\mathbb{R}^{n+1}$. The intensity of this light is given by a nonnegative density function $g(x)$. The prescribed intensity at $y \in \Omega ^*$ is given by $f(y)$. Also $g$ and $f$ satisfy 2.1. A generalized solution to the problem is a refracting surface $\mathcal{S} = \{u(x) : x \in \Omega \}$ with $u:\Omega \to \mathbb{R}$ such that for the map $T_{\mathcal{S}} : \Omega \to \Omega ^*$ we have $T_{\mathcal{S}}(\Omega ) \subset \Omega ^*$ and that the energy conservation

$$\begin{equation*} \int _F g(x) \, dx = \int _{T_{\mathcal{S}}(F)} f(x) \, dx \end{equation*}$$

is satisfied for any measurable subset $F$ of $\Omega$.

If $dS_{\Omega }$ and $dS_{\Omega ^*}$ represent area elements on $\Omega$ and $\Omega ^*$ respectively, and assuming the ray tracing map is smooth, the energy conservation condition along with change of variables requires that

$$\begin{equation*} \det DT_{\mathcal{S}} = \dfrac{dS_{\Omega ^*}}{dS_{\Omega }} = \dfrac{g(x)}{f(T_{\mathcal{S}}(x))} \end{equation*}$$

where $DT_{\mathcal{S}}$ is the Jacobian matrix of $T_{\mathcal{S}}$. $DT_{\mathcal{S}}$ can then be explicitly expressed as in 11 in terms of $u$ and its gradient. In particular if $u \in C^2(\Omega )$ is a solution of the parallel near field refractor problem, then $u$ satisfies the Monge–Ampère-type PDE

$$\begin{equation*} |\det (D^2u + \mathcal{A}(\kappa , x, u, Du)) |= \left| \eta (\kappa , x, u, Du) \dfrac{g(x)}{f(T_{\mathcal{S}}(x))}\right|. \end{equation*}$$

The $C^2$ smoothness of weak solution when $f$ and $g$ are bounded and away from zero under some geometric conditions is proved in 11 by using standard technique from PDE theory.

A more comprehensive approach can be carried out by using what are called generated Jacobian equations.

Given $\Omega ,\Omega ^*$ domains in $\mathbb{R}^n$, a generating function $G:\Omega \times \Omega ^*\times \mathbb{R}^+\to \mathbb{R}^+$, with variables $G(x,y,v)$, satisfies the following conditions:

(a)

$G_v<0$ for all $(x,y)\in \Omega \times \Omega ^*$;

(b)

the map $(y,v)\to (G(x,y,v),G_x(x,y,v))$ is invertible for each $x\in \Omega$;

(c)

for each $y\in \Omega ^*$, $G(x,y,v)\to +\infty$ as $v\to 0^+$, uniformly in $x\in \Omega$;

(d)

$G$ is $C^1$ in $v$, is $C^2$ in $x,y$, and for each $\alpha >0$, $\sup _{\Omega \times \Omega ^*\times [0,\alpha ]}|G_x(x,y,v)|<\infty$.

By using the first condition, the dual generating function $H$ is defined by $G(x,y, H(x,y,v)) = v$ for all $(x, y, v)$ when the expression is well defined.

A generated Jacobian equation is a PDE of the form

$$\begin{equation} \det DY(x, u, Du) = \psi (x, u, Du) {\tag{4.1}} \end{equation}$$

where the vector field $Y(x,v,p)$ is generated from $G$ by

$$\begin{equation*} D_x G(x, Y, Z) = p \end{equation*}$$

and $Z(x,v,p) = H(x,Y(x,v,p),v)$. For brevity consider the case where

$$\begin{equation} \psi (x, u, Du) = \dfrac{g(x)}{f(Y(x, u, Du))} {\tag{4.2}} \end{equation}$$

for $g \in L^1(\Omega )$ and $f \in L^1(\Omega ^*)$.

For example, in an optimal transportation problem with cost $c(x,y)$ which satisfies compatible conditions, setting $G(x,y,v) = c(x,y) - v$ we get $G_x(x,y,v) = c_x(x,y)$. The vector field $Y(x,v,p)$ is generated by $c_x(x, Y(x,p)) = p$ and the equation corresponding to 4.1 will be a generalized Monge–Ampère equation given as $\det [D^2 u - A(\cdot , u, Du )] = B(\cdot , u, Du)$ where $A(x,z,p) = c_{xx}(x, Y(x,p))$ and $B(x, z, p) = \det c_{xy}(x, Y(x,p)) \psi (x,z,p)$. See 19.

In particular if $\mu = g(x) \, dx$ and $\nu = f(y) \, dy$ and the cost function is the quadratic cost $c(x,y) = \frac{1}{2} |x - y|^2$ ,we obtain the vector field $Y(x, p) = x-p$ from which 4.1 will be

$$\begin{equation} \det D(x-Du(x)) = \psi (x, u(x), Du(x)). \tag{4.3} \end{equation}$$

If $T$ is any transport map from $\mu$ to $\nu$, we know from energy relations that

$$\begin{equation} \det DT = \dfrac{g(x)}{f(T(x))}. {\tag{4.4}} \end{equation}$$

In particular if $T$ is the optimal map and $\phi$ is the potential, $T_{\phi }(x) = x - D\phi (x)$. Using this relation we observe 4.4 is of type 4.1 with $Y = x-p$. If we further express $T_{\phi } = Du$ we get the standard Monge–Ampère equation

$$\begin{equation*} \det D^2 u(x) = \dfrac{g(x)}{f(Du(x))}. \end{equation*}$$

A function $\phi :\Omega \to \mathbb{R}$ is $G$-convex if for each $x_o\in \Omega$ there exists $y_o\in \Omega ^*$ and $v_o\in \mathbb{R}^+$ such that $\phi (x)\geq G(x,y_o,v_o)$ for all $x\in \Omega$ with equality at $x=x_o$; i.e., $G(\cdot ,y_o,v_o)$ supports $\phi$ at $x_o$. For a $G$-convex function $\phi$, the $G$-subdifferential of $\phi$ at $x_0$ is defined by

$$\begin{equation*} \partial _G\phi (x_0)=\{y_0\in \Omega ^*: \text{there exists } v_0\in \mathbb{R}^+ \quad \text{such that $G(\cdot ,y_0,v_0)$ supports $\phi $ at $x_0$}\}. \end{equation*}$$

If $\phi$ is differentiable, then $\partial _G\phi (x_0)$ has exactly one element. Following 19 a $G$-convex function $\phi$ is called a generalized solution of 4.1 with $\psi$ given by 4.2 if

$$\begin{equation*} \int _E g(x) \, dx = \int _{\partial _G\phi (E)} f(x) \, dx \end{equation*}$$

for all $E \subset \Omega$.

For the near field refractor problem with point source, the corresponding generating functions will be

$$\begin{equation*} G(x, P, b) = \dfrac{1}{h(x, P, b)} \end{equation*}$$

with $h(x, P, b)$ defined as in 3.2. If $\mathcal{S} = \{\rho (x) x : x \in \overline{\Omega }\}$ is a solution to the near field refractor problem and $\phi (x) = \dfrac{1}{\rho (x)}$ then the $G$-subdifferential of $\phi$ and the inverse of the ray tracing mapping, $T_{\mathcal{S}}$ are related by

$$\begin{equation*} \partial _G\phi (x_o) = T_{\mathcal{S}} (x_o) \end{equation*}$$

and hence we have

$$\begin{equation*} \int _E g(x) \, dx = \int _{T_{\mathcal{S}}(E)} f(x) \, dx = \int _{\partial _G\phi (E)} f(x) \, dx. \end{equation*}$$

Thus, the near field refractor problem can be studied via the theory of generated Jacobian equations. For the particular case where $\kappa >1$ this is done in 12 where they proved the $C^{2, \alpha }$ regularity of solutions. The solutions to the reflector problem with collimated beam, the near field reflector problem with point source and other related problems in geometric optics can also be expressed as generated Jacobian equations 19. More discussion on the relationship between design of optical surfaces and generated Jacobian equations is in 5. However, for design problems involving both reflection and refraction in anisotropic media, the regularity theory using PDE methods or otherwise is an untouched territory.