Unified mathematical framework for a class of fundamental freeform optical systems

: We present a unified mathematical framework for sixteen fundamental optical systems. The systems have a parallel or point source and a parallel, point, near-field or far-field target. These choices give eight configurations if we use reflectors only and take the minimum number of freeform surfaces required. Similarly, we get eight lens systems if we only use lens surfaces. The mathematical model for each system is based on Hamilton’s characteristic functions and conservation of luminous flux. Some configurations lead to standard or generalized Monge-Ampère equations. The remaining systems are described by so-called generated Jacobian equations.


Introduction
In this paper we present a unified mathematical framework for sixteen systems. These basic configurations have a parallel or point source and a parallel, point, near-field or far-field target. These choices give eight systems, if we use reflectors only and take the minimum number of freeform surfaces required. The systems are shown in Fig. 3. Similarly, we get eight basic lens systems if we only use lens surfaces, see Fig. 4. Our goal is to find the shape and location of the surfaces that convert a given source distribution into a desired target distribution. All systems can be modeled by a nonlinear partial differential equation (PDE). Many papers in freeform illumination optics cover only one particular optical system. Here we derive models for all sixteen using our unified framework.
A second approach is to solve the Monge-Kantorovich transport problem that corresponds to the design problem using optimization strategies, see Angenent et al. [37], Benamou et al. [38], and Haber et al. [39]. Doskolovich et al. reduce the transport problem to a linear assignment problem [40][41][42]. Oliker et al. have introduced the methods of supporting paraboloids and ellipsoids, now known as the supporting quadric method (SQM) [43][44][45][46][47][48][49][50][51]. SQM combines geometrical techniques with methods from calculus of variations. Instead of dealing directly with the complicated nonlinear PDEs, it produces weak solutions that are the envelope of a set of e.g. paraboloids of revolution. By construction there exists a supporting paraboloid at every point of the computed optical surface. The surfaces found are not necessarily differentiable. The method provides both proofs of existence of solutions as well as a way to compute them numerically.
Third are ray-mapping techniques. In these methods a ray mapping is computed which is subsequently used to compute the surface using the law of reflection or Snell's law. This method is used by Bösel and Gross [52,53] and Feng et al. [54][55][56].
We have organized this paper as follows. In Section 2 we formulate the freeform design problem for a general optical system, introduce notation and recap the main properties of Hamilton's characteristic functions. Next we present three mathematical models. In Sections 3 and 4 we discuss the eight reflector and eight lens systems, respectively. We sketch our least-squares algorithm to numerically solve the mathematical models in Section 5. We draw conclusions in Section 6.

Problem formulation
We consider a general optical system such as the one sketched in Fig. 1. The goal of freeform design is to find the shape and location of one or more surfaces that convert a given source distribution into a desired target distribution. The surfaces can be either reflector or lens surfaces. The source distribution (emittance or intensity) is denoted by f (x), where x ∈ X are either spatial or directional coordinates and X is the source domain. Likewise, the target distribution (illuminance or intensity) is g(y), y ∈ Y, where y ∈ Y are spatial or directional coordinates and Y is the target domain. The optical map m : X → Y connects x with y, so y = m(x).
We use spatial coordinates for a parallel source, parallel target and near field target. We assume standard Cartesian coordinates in these cases. We use directional coordinates for a point source, point target and far field target. It is convenient to use so-called stereographic coordinates. These are a parametrization of the unit sphere: If we denote the direction vector of an emitted ray bŷ s = (s 1 , s 2 , s 3 ), then its stereographic projection onto the equatorial plane z = 0 from the south pole is given by x = (x 1 , x 2 ) = (s 1 , s 2 )/(1 + s 3 ), see Fig. 2. Note that the hat (ˆ) indicates thatŝ has length one. The inverse projection is given byŝ = (2x 1 , 2x 2 , 1 − |x| 2 )/(1 + |x| 2 ). In the same way, the direction vectort = (t 1 , t 2 , t 3 ) of a target ray has stereographic projection y = (y 1 , y 2 ). At the target, the projection can be either from the south or the north pole.
In Sections 2.3-2.5 we present models for the optical map. These models follow from the evaluation of one of the four Hamilton's characteristic functions [57,58]. This gives a geometrical description of the optical system. We briefly state the definitions and main properties of the characteristic functions. Consider a light ray that intersects the source plane z = z s at (two-dimensional) spatial coordinates q s and the target plane z = z t at spatial coordinates q t . The subscripts s and t refer to source and target. Hamilton's point characteristic V = V(q s , q t ) is defined as the optical path length between (q s , z s ) and (q t , z t ), so V = ∫ C n ds, with C the path connecting the two points and n the local refractive index. Let p s be the vector with the first two components of the momentum at the source plane and, in the same way, p t at the target plane. We have p s = n(s 1 , s 2 ) and p t = n(t 1 , t 2 ). The point characteristic V satisfies In our models, we use V if both x and y are Cartesian coordinates, so if x = q s and y = q t . The mixed characteristic W is defined as W(q s , p t ) = V(q s , q t ) − q t · p t . Note that W depends on two variables only: Using Eq. (1), we have ∂W/∂q t = ∂V/∂q t − p t = 0. It satisfies We will use W if x is Cartesian (x = q s ) and y is the stereographic projection oft. The mixed characteristic of the second kind W * is defined as W * (p s , q t ) = V(q s , q t ) + q s · p s . It depends on two variables only, because ∂W * /∂q s = 0. It satisfies We will use W * if x is the stereographic projection ofŝ and y is Cartesian (y = q t ). Finally, the angular characteristic T is defined as T(p s , p t ) = V(q s , q t ) + q s · p s − q t · p t . It depends on two variables only, because ∂T/∂q s = ∂T/∂q t = 0. It satisfies We will use T if x and y are the stereographic projections ofŝ andt, respectively.

Conservation of luminous flux
Let A be a part of the source domain, so A ⊂ X. Because luminous flux is conserved, we know that all light emitted from A should reach the image set m(A) in the target domain. This leads to the conservation law ∬ Here, J X (x) equals one if we use Cartesian coordinates in the source domain. If we use stereographic coordinates, it is the Jacobian of the stereographic projection of the direction vector s of the emitted rays, so In the same way J Y (y) equals one for Cartesian coordinates in the target domain and J Y (y) = 4/(1 + |y| 2 ) 2 for stereographic coordinates. Dm is the Jacobian matrix of m. We assume that det(Dm(x))>0 in Eq. (5). For far-field targets we use the far-field approximation.
Because Eq. (5) needs to hold for any subset A of X, we find the flux balance Equation (7) will lead to elliptic PDEs. For all models, we impose the condition m(X) = Y stating that all light from the source domain must be transferred to the target domain.

Model 1: standard Monge-Ampère equation
The first model we consider for the optical map will allow us to write m as the gradient of a function u 1 , so m = ∇u 1 . The function u 1 is related to the location of the optical surface. We will show that u 1 satisfies a second-order nonlinear PDE, the so-called Monge-Ampère equation.
A typical example of an optical system that can be described by this model is the parallel-tofar-field reflector sketched in Fig. 3(d). A detailed derivation of the model for this system can be found in [1,Sec. 3.2].
For the source we use Cartesian coordinates, so J X = 1. For the target we use stereographic coordinates, so J Y (y) = 4/(1 + |y| 2 ) 2 . We assume that the source domain is located at z = 0 and the target domain at z = −L. Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. With Eq. (2) it follows that Hamilton's mixed characteristic W is a function of p t (and hence y) only. We define the location of the reflector as 3 , which can be rewritten as The first term in Eq. (8) depends on y only and we call it u 2 (y). For the right-hand side of Eq. (8), we use p t /(1 − t 3 ) = y. This leads to A possible solution of Eq. (9) reads The second equation in Eq.
We call this PDE the standard Monge-Ampère (SMA) equation.

Model 2: generalized Monge-Ampère equation
In the second model, the optical map m is no longer the gradient of a function u 1 . It is a (complicated) function of x and ∇u 1 , where u 1 is again related to the location of the optical surface. We will show that u 1 satisfies a generalized Monge-Ampère (GMA) equation.
A typical example of an optical system that can be described by this model is the point-to-farfield reflector sketched in Fig. 3(h). A detailed derivation of the model for this system can be found in [20].
For both the source and the target we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 and J Y (y) = 4/(1 + |y| 2 ) 2 . We assume that the source domain is located at z = 0 and the target domain at z = −L. The reflector surface has parametrization r(ŝ) = u(ŝ)ŝ. Because we have a point source, we find q s = 0. With Eq. (4) it follows that Hamilton's angular characteristic T is a function of p t (and hence y) only. Evaluating T gives T(p t ) = (1 −ŝ ·t)u(ŝ) − Lt 3 , which can be rewritten as log(T(p t ) + Lt 3 ) = log(1 −ŝ ·t) + log(u(ŝ)). Transforming to stereographic coordinates leads to where The function c is known as the cost function in optimal transport theory. Note that Eq. (13) is a generalization of Eq. (9). A possible solution of Eq. (13) reads The second equation in Eq. (14) implies In principle one can compute m from Eq. (15), but this is difficult and we proceed in a different manner. Substituting y = m(x) and once more differentiating with respect to x leads to the following equation for Dm Here D 2 u 1 and D xx c are Hessian matrices, D xy c is the matrix with mixed second-order derivatives. We have now found the matrix equation It follows that det(P) = det(C) det(Dm). Note that P is a 2 × 2-matrix, so det(−P) = det(P). Using the definition of P and Eq. (7) we obtain the flux balance which is a GMA equation.

Model 3: generated Jacobian equation
In the third model, we cannot separate the variables as we did in the left-hand side of Eqs. (9) and (13). Instead we will find a more complicated relation between u 1 (x) and u 2 (y). We will show that u 1 satisfies a generated Jacobian equation (GJE). A typical example of an optical system that can be described by this model is the parallel-tonear-field reflector sketched in Fig. 3(c). A detailed derivation of the model for this system can be found in [1,Sec. 3.3].
For both the source and target we use Cartesian coordinates, so J X = J Y = 1. We assume that the source domain is located at z = 0 and the target domain at z = −L. Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. With Eq. (1) it follows that Hamilton's point characteristic V is a function of q t (and hence y) only. We define the location of the reflector as Solving this equation for u 1 leads to the following generalization of Eq. (13): where u 2 (y) = V(y) and G(x, y, v) = 1 2 v − |x − y| 2 /(2v). The function G is the so-called generating function. Because ∂G/∂v>0, the function G(x, y, ·) has an inverse for fixed x and y. We denote this inverse by H(x, y, ·), so u 2 (y) = H(x, y, u 1 (x)). A possible solution of Eq. (19) reads The second equation in Eq. (20) implies In principle one can compute m from Eq. (21), but this is difficult and we proceed in a different manner. We introduce˜︁ H(x, y) := H(x, y, u 1 (x)), substitute y = m(x) and differentiate once more with respect to x to find the following equation for Dm We have now found the matrix equation H(x, m(x)).
It follows that det(P) = det(C) det(Dm). Using the definition of P and Eq. (7) we obtain det(D xx˜︁ H(x, m(x) which is a GJE.

Basic reflector systems
In this section we discuss the eight basic reflector systems sketched in Fig. 3. For all of the systems we assume that the source is located at z = 0.

Parallel-to-parallel reflector
We consider the parallel source to parallel target reflector system sketched in Fig. 3(a). This system can be described by Model 1 from Section 2.3. A detailed derivation of the model for this system can be found in [32]. For both the source and target we use Cartesian coordinates, so J X = J Y = 1. We assume that the target is located at z = L. The first reflector has equation z = u 1 (x).

Parallel-to-point reflector
In this section we consider the parallel source to point target reflector system sketched in Fig. 3(b). This system can be described by Model 2 from Section 2.4. Note that the system is the same as the point source to parallel target system, that we will consider in Section 3.5, when we interchange source and target coordinates. A detailed derivation of the model for the latter system can be found in [59]. For the source we use Cartesian coordinates, so J X = 1. For the target we use stereographic coordinates, so J Y (y) = 4/(1 + |y| 2 ) 2 . We assume that the target is located at z = L. Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. The target is a point, so q t = 0. With Eq. (2) it follows that Hamilton's mixed characteristic W is constant. If we introduce the reduced optical path length β = W − L, we get Eq. (13) with c(x, y) = log Repeating the steps described in Section 2.4 leads to the GMA equation

Parallel-to-near-field reflector
We consider the parallel source to near field reflector system sketched in Fig. 3(c). This system can be described by Model 3 from Section 2.5, and in fact it is the system discussed in that section.

Parallel-to-far-field reflector
We consider the parallel source to far field reflector system sketched in Fig. 3(d). This system can be described by Model 1 from Section 2.3, and in fact it is the system discussed in that section.

Point-to-parallel reflector
In this section we consider the point source to parallel target reflector system sketched in Fig. 3(e). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [59]. Note that the system is the same as the parallel source to point target system, described in Section 3.2, when we interchange source and target coordinates. For the source we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 . For the target we use Cartesian coordinates, so J Y = 1. We assume that the target is located at z = L. Because we have a point source, we find q s = 0. Target rays are parallel and perpendicular to z = L, so p t = 0. With Eq. (3) it follows that Hamilton's mixed characteristic W * is constant. If we introduce the reduced optical path length β = W * − L, we get Eq. (13) with Repeating the steps in Section 2.4 leads to the GMA equation

Point-to-point reflector
In this section we consider the point source to point target reflector system sketched in Fig. 3(f). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [33,Sec. 3.4]. For both the source and the target we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 and J Y (y) = 4/(1 + |y| 2 ) 2 . We assume that the target is located at z = L. Because we have both a point source and a point target, we find q s = q t = 0. With Eq. (4) it follows that Hamilton's angular characteristic T is constant. We get Eq. (13) with Repeating the steps in Section 2.4 leads to the GMA equation

Point-to-near-field reflector
In this section we consider the point source to near field reflector system sketched in Fig. 3(g). This system can be described by Model 3 from Section 2.5.
For the source we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 . For the target we use Cartesian coordinates, so J Y = 1. We assume that the target is located at z = L. Because we have a point source, we find q s = 0. With Eq. (3) it follows that Hamilton's mixed characteristic W * is a function of q t (and hence y) only. Evaluating W * we obtain Eq. (19) with Repeating the steps in Section 2.5 we find the GJE det (D xx˜︁ H(x, m(x)

Point-to-far-field reflector
In this section we consider the point source to far field reflector system sketched in Fig. 3(h). This system can be described by Model 2 from Section 2.4, and in fact it is the system discussed in that section.

Basic lens systems
In this section we discuss the eight basic lens systems sketched in Fig. 4. For all of the systems we assume that the source is located at z = 0, the target at z = L and that the lens has refractive index n>1.

Parallel-to-parallel lens
We consider the parallel source to parallel target lens system sketched in Fig. 4(a). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [31]. For both the source and target we use Cartesian coordinates, so J X = J Y = 1. Because both emitted and target rays are parallel and perpendicular to z = 0 and z = L, we have p s = p t = 0. With Eq. (1) it follows that Hamilton's point characteristic V is constant. If we introduce the reduced optical path length β = V − L, we get Eq. (13) with Repeating the steps in Section 2.4 leads to the GMA equation

Parallel-to-point lens
In this section we consider the parallel source to point target lens system sketched in Fig. 4(b). This system can be described by Model 3 from Section 2.5. Note that the system is the same as the point source to parallel target system, that we will consider in Section 4.5, when we interchange source and target coordinates. For the source we use Cartesian coordinates, so J X = 1. For the target we use stereographic coordinates, so J Y (y) = 4/(1 + |y| 2 ) 2 . Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. The target is a point, so q t = 0. With Eq. (2) it follows that Hamilton's mixed characteristic W is constant. If we introduce the reduced optical path length β = W − L, we get Eq. (19) with Repeating the steps in Section 2.5 we find the GJE

Parallel-to-near-field lens
We consider the parallel source to near field lens system sketched in Fig. 4(c). Note that the first lens surface is flat and perpendicular to the emitted rays, so it does not change their direction. This system can be described by Model 3 from Section 2.5. For both the source and target we use Cartesian coordinates, so J X = J Y = 1. Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. With Eq. (1) it follows that Hamilton's point characteristic V is a function of q t (and hence y) only. Evaluating V we obtain Eq. (19) with Repeating the steps in Section 2.5 we find the GJE det (D xx˜︁ H(x, m(x)

Parallel-to-far-field lens
We consider the parallel source to far field lens system sketched in Fig. 4(d). Note that the first lens surface is flat and perpendicular to the emitted rays, so it does not change their direction. This system can be described by Model 1 from Section 2.3. A detailed derivation of the model for this system can be found in [33,Sec. 3.1.2]. For the source we use Cartesian coordinates, so J X = 1. For the target we use stereographic coordinates, so J Y (y) = 4/(1 + |y| 2 ) 2 . Because emitted rays are parallel and perpendicular to z = 0, we have p s = 0. With Eq. (2) it follows that Hamilton's mixed characteristic W is a function of p t (and hence y) only. Evaluating W we get Eq. (9) where u 1 (x) defines the location of the freeform lens surface via z = u 1 (x). Repeating the steps in Section 2.3 we find the SMA Eq. (12).

Point-to-parallel lens
In this section we consider the point source to parallel target lens system sketched in Fig. 4(e). This system can be described by Model 3 from Section 2.5. Note that the system is the same as the parallel source to point target system, described in Section 4.2, when we interchange source and target coordinates.
For the source we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 . For the target we use Cartesian coordinates, so J Y = 1. Because we have a point source, we find q s = 0. Target rays are parallel and perpendicular to z = L, so p t = 0. With Eq. (3) it follows that Hamilton's mixed characteristic W * is constant. If we introduce the reduced optical path length β = W * − L, we get Eq. (19) with Repeating the steps in Section 2.5 we find the GJE det(D xx˜︁ H(x, m(x))) = det(C) 4f (x) g(m(x))(1 + |x| 2 ) 2 . (42)

Point-to-point lens
In this section we consider the point source to point target lens system sketched in Fig. 4(f). This system can be described by Model 3 from Section 2.5. For both the source and the target we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 and J Y (y) = 4/(1 + |y| 2 ) 2 . Because we have both a point source and a point target, we find q s = q t = 0. With Eq. (4) it follows that Hamilton's angular characteristic T is constant. We get Eq. (19) with and S := Repeating the steps in Section 2.5 we find the GJE det (D xx˜︁ H(x, m(x)

Point-to-near-field lens
In this section we consider the point source to near field lens system sketched in Fig. 4(g). This system can be described by Model 3 from Section 2.5.
For the source we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 . For the target we use Cartesian coordinates, so J Y = 1. The first lens surface is spherical and does not change the direction of the rays. The second lens surface is freeform. Because we have a point source, we find q s = 0. With Eq. (3) it follows that Hamilton's mixed characteristic W * is a function of q t (and hence y) only. Evaluating W * we obtain Eq. (19) with with S := √︁ (nv − p s · y) 2 − (n 2 − 1)(v 2 − |y| 2 ). Repeating the steps described in Section 2.5 we find the GJE det(D xx˜︁ H(x, m(x)

Point-to-far-field lens
In this section we consider the point source to far field lens system sketched in Fig. 4(h). This system can be described by Model 2 from Section 2.4. A detailed derivation of the model for this system can be found in [21]. For both the source and the target we use stereographic coordinates, so J X (x) = 4/(1 + |x| 2 ) 2 and J Y (y) = 4/(1 + |y| 2 ) 2 . The first lens surface is spherical and does not change the direction of the rays. The second lens surface is freeform and has parametrization r(ŝ) = u(ŝ)ŝ. Because we have a point source, we find q s = 0. With Eq. (4) it follows that Hamilton's angular characteristic T is a function of p t (and hence y) only. Evaluating T we get Eq. (13) with Repeating the steps in Section 2.4 leads to the GMA Eq. (18).

Numerical method for finding the optical surface(s)
Our numerical method for finding the optical surface(s) is an iterative least-squares solver. It is based on the PDE of SMA, GMA or GJE type together with the condition m(X) = Y. Our implementation is based on the following breakdown in substeps. We compute P, b and m such that Here, b maps the boundary of the source domain to the boundary of the target domain. Equations (48) are solved in an iterative procedure: We let n = 0 and choose an initial guess m 0 . Next, let J I (m, P) = 1 2 ∬ X ∥CDm − P∥ 2 dx and compute P n+1 such that it minimizes J I (m n , P). This leads to a constrained minimization problem-the contraint is Eq. (48b). To enforce the boundary condition Eq. (48c), we define J B (m, b) = 1 2 ∫ ∂X |m − b| 2 ds and compute b n+1 such that it minimizes J B (m n , b). We finally introduce J(m, P, b) = α J I (m, P) + (1 − α) J B (m, b) where 0<α<1 and compute m n+1 such that it minimizes J(m, P n+1 , b n+1 ). This leads to an elliptic PDE for m that we solve with a finite volume method.
For Models 1 and 2, we repeat the steps above until convergence. Upon convergence of the iterative procedure, the location of the surface is calculated from the optical map using Eq. (15). For Model 3, we need to do an additional step during the iteration to compute u 1 . To do so, we use Eq. (21) and let I(u 1 , m) = 1 2 ∬ X |∇ x H(x, m, u 1 ) + H v (x, m, u 1 )∇u 1 | 2 dx, where H v is the partial derivative of H with respect to its third argument. We compute u n+1 1 such that it minimizes I(u 1 , m n+1 ).

Conclusions
In this paper we presented a unified mathematical framework for sixteen basic optical systems, eight reflector and eight lens systems. Each system has a parallel or point source and a parallel, point, near-field or far-field target. The mathematical model for each system is based on Hamilton's characteristic functions and conservation of luminous flux.
Three systems lead to the standard Monge-Ampère (SMA) equation (Model 1). This holds true for the parallel-to-parallel reflector, the parallel-to-far-field reflector and the parallel-to-far-field lens system.
More complicated optical systems cannot be described by Model 1 and need a generalization. Six of them give a generalized Monge-Ampère (GMA) equation (Model 2). Systems in this category are the parallel-to-point reflector, the point-to-parallel reflector, the point-to-point reflector, the point-to-far-field reflector, the parallel-to-parallel lens and the point-to-far-field lens.
The remaining seven systems cannot be described by either Model 1 or 2. They need a yet more general formulation using so-called generating functions. This in turn leads to generated Jacobian equations (GJEs, Model 3). All systems with near-field targets belong to this category.
Note that the three models are proper generalizations in the sense that all Model 1 systems can also be described by Model 2, and Model 2 systems can be formulated using generating functions (Model 3) too.
Disclosures. The authors declare no conflicts of interest. Data availability. No data were generated or analyzed in the presented research.